Momentum-dependent two-loop QCD corrections to the neutral Higgs-boson masses in the MSSM

Results are presented for the momentum dependent two-loop contributions of O(alpha_t alpha_s) to the masses and mixing effects in the Higgs sector of the MSSM. They are obtained in the Feynman-diagrammatic approach using a mixed on-shell/DRbar renormalization that can directly be matched onto the higher-order corrections included in the code FeynHiggs. The new two-loop diagrams are evaluated with the program SecDec. The combination of the new momentum dependent two-loop contribution with the existing one- and two-loop corrections in the on-shell/DRbar scheme leads to an improved prediction of the light MSSM Higgs boson mass and a correspondingly reduced theoretical uncertainty. We find that the corresponding shifts in the lightest Higgs-boson mass M_h are below 1 GeV in all scenarios considered, but can extend up to the level of the current experimental uncertainty. The results are included in the code FeynHiggs.


Introduction
The ATLAS and CMS experiments at CERN have recently discovered a new boson with a mass around 125.6 GeV [1,2]. Within the present experimental uncertainties this new boson behaves like the Higgs boson of the Standard Model (SM) [3]. However, the newly discovered particle can also be interpreted as the Higgs boson of extended models. The Higgs sector of the Minimal Supersymmetric Standard Model (MSSM) [4] with two scalar doublets accommodates five physical Higgs bosons. In lowest order these are the light and heavy CP-even h and H, the CP-odd A, and the charged Higgs bosons H ± . The measured mass value, having already reached the level of a precision observable with an experimental accuracy of about 500 MeV, plays an important role in this context. In the MSSM the mass of the light CP-even Higgs boson, M h , can directly be predicted from the other parameters of the model. The accuracy of this prediction should at least match the one of the experimental result.
The Higgs sector of the MSSM can be expressed at lowest order in terms of the gauge couplings, the mass of the CP-odd Higgs boson, M A , and tan β ≡ v 2 /v 1 , the ratio of the two vacuum expectation values. All other masses and mixing angles can therefore be predicted. Higher-order contributions can give large corrections to the tree-level relations [5,6]. An upper bound for the mass of the lightest MSSM Higgs boson of M h < ∼ 135 GeV was obtained [7], and the remaining theoretical uncertainty in the calculation of M h , from unknown higher-order corrections, was estimated to be up to 3 GeV, depending on the parameter region. Recent improvements have lead to a somewhat smaller estimate of up to ∼ 2 GeV [8,9] (see below).
Experimental searches for the neutral MSSM Higgs bosons have been performed at LEP [10,11], placing important restrictions on the parameter space. At Run II of the Tevatron the search was continued, but is now superseeded by the LHC Higgs searches. Besides the discovery of a SM Higgs-like boson the LHC searches place stringent bounds, in particular in the regions of small M A and large tan β [12]. At a future linear collider (ILC) a precise determination of the Higgs boson properties (either of the light Higgs boson at ∼ 125.6 GeV or heavier MSSM Higgs bosons within the kinematic reach) will be possible [13]. In particular a mass measurement of the light Higgs boson with an accuracy below ∼ 0.05 GeV is anticipated [14]. The interplay of the LHC and the ILC in the neutral MSSM Higgs sector has been discussed in Refs. [15,16].
For the MSSM 1 the status of higher-order corrections to the masses and mixing angles in the neutral Higgs sector is quite advanced. The complete one-loop result within the MSSM is known [21][22][23][24]. The by far dominant one-loop contribution is the O(α t ) term due to top and stop loops (α t ≡ h 2 t /(4π), h t being the top-quark Yukawa coupling). The computation of the two-loop corrections has meanwhile reached a stage where all the presumably dominant contributions are available [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]. In particular, the O(α t α s ) contributions to the selfenergies -evaluated in the Feynman-diagrammatic (FD) as well as in the effective potential (EP) method -as well as the O(α 2 t ), O(α b α s ), O(α t α b ) and O(α 2 b ) contributions -evaluated in the EP approach -are known for vanishing external momenta. An evaluation of the momentum dependence at the two-loop level in a pure DR calculation was presented in Ref. [40]. A (nearly) full two-loop EP calculation, including even the leading three-loop corrections, has also been published [41]. However, within the EP method all contributions are evaluated at zero external momentum, in contrast to the FD method, which in principle allows non-vanishing external momentum. Further, the calculation presented in Ref. [41] is not publicly available as a computer code for Higgs-mass calculations. Subsequently, another leading three-loop calculation of O(α t α 2 s ), depending on the various SUSY mass hierarchies, has been performed [42], resulting in the code H3m (which adds the three-loop corrections to the FeynHiggs result). Most recently, a combination of the full one-loop result, supplemented with leading and subleading two-loop corrections evaluated in the Feynmandiagrammatic/effective potential method and a resummation of the leading and subleading logarithmic corrections from the scalar-top sector has been published [8] in the latest version of the code FeynHiggs [7,8,17,26,43]. While previous to this combination the remaining theoretical uncertainty on the lightest CP-even Higgs boson mass had been estimated to be about 3 GeV [6,7], the combined result was roughly estimated to yield an uncertainty of about 2 GeV [8,9]; however, more detailed analyses will be necessary to yield a more solid result.
In the present paper we calculate the two-loop O(α t α s ) corrections to the Higgs boson masses in a mixed on-shell/DR scheme. Compared to previously known results [25,26,31] we evaluate here corrections that are proportional to the external momentum of the relevant Higgs boson self-energies. These corrections can directly be added to the corrections included in FeynHiggs. An overview of the relevant sectors and the calculation is given in Sect. 2, whereas in Sect. 3 we discuss the size and relevance of the new two-loop corrections. Our conclusions are given in Sect. 4.

The Higgs-boson sector of the MSSM
The MSSM requires two scalar doublets, which are conventionally written in terms of their components as follows, The Higgs boson sector can be described with the help of two independent parameters (besides the SM gauge couplings), conventionally chosen as tan β = v 2 /v 1 , the ratio of the two vacuum expectation values, and M 2 A , the mass of the CP-odd Higgs boson A. The bilinear part of the Higgs potential leads to the tree-level mass matrix for the neutral CPeven Higgs boson, in the (φ 1 , φ 2 ) basis and being expressed in terms of the parameters M Z , M A and the angle β. Diagonalization yields the tree-level masses m h,tree , m H,tree .
The higher-order corrected CP-even Higgs boson masses in the MSSM are obtained from the corresponding propagators dressed by their self-energies. The inverse propagator matrix in the (φ 1 , φ 2 ) basis is given by where theΣ(p 2 ) denote the renormalized Higgs-boson self-energies, p being the external momentum. The renormalized self-energies can be expressed through the unrenormalized self-energies, Σ(p 2 ), and counterterms involving renormalization constants δm 2 and δZ from parameter and field renormalization. With the self-energies expanded up to two-loop order, Σ =Σ (1) +Σ (2) , one has for the CP-even part at the i-loop level (i = 1, 2), The counterterms are determined by appropriate renormalization conditions and are given in the Appendix.
The renormalized self-energies in the (φ 1 , φ 2 ) basis can be rotated into the physical (h, H) basis where the tree-level propagator matrix is diagonal, via with the matrix which diagonalizes the tree-level mass matrix (2). The CP-even Higgs boson masses are determined by the poles of the (h, H)-propagator matrix. This is equivalent to solving the equation yielding the loop-corrected pole masses, M h and M H . Here we use the implementation in the code FeynHiggs [7,8,17,26,43], supplemented by the new momentum dependent O(α t α s ) corrections, as described in 2.4. Our calculation is performed in the Feynman-diagrammatic (FD) approach. To arrive at expressions for the unrenormalized self-energies and tadpoles at O(α t α s ), the evaluation of genuine two-loop diagrams and one-loop graphs with counterterm insertions is required. Example diagrams for the neutral Higgs-boson self-energies are shown in Fig. 1, and for the tadpoles in Fig. 2. For the counterterm insertions, described in subsection 2.2, one-loop diagrams with external top quarks/squarks have to be evaluated as well, as displayed in Fig. 3. The complete set of contributing Feynman diagrams has been generated with the program FeynArts [49] (using the model file including counterterms from Ref. [50]), tensor reduction and the evaluation of traces was done with support from the programs FormCalc [51] and TwoCalc [52], yielding algebraic expressions in terms of the scalar one-loop functions A 0 , B 0 [53], the massive vacuum two-loop functions [54], and two-loop integrals which depend on the external momentum. These integrals have been evaluated with the program SecDec [63,64], see subsection 2.3.

The scalar-top sector of the MSSM
The bilinear part of the top-squark Lagrangian, contains the stop mass matrix Mt, given by with Q t and T 3 t denote the charge and isospin of the top quark, A t is the trilinear coupling between the Higgs bosons and the scalar tops, and µ is the Higgsino mass parameter. Below we use M SUSY := Mt L = Mt R for our numerical evaluation. However, the analytical calculation has been performed for arbitrary Mt L and Mt R . Mt can be diagonalized with the help of a unitary transformation matrix Ut, parametrized by a mixing angle θt, to provide the eigenvalues m 2 t 1 and m 2 t 2 as the squares of the two on-shell top-squark masses. For the evaluation of the O(α t α s ) two-loop contributions to the self-energies and tadpoles of the Higgs sector, renormalization of the top/stop sector at O(α s ) is required, giving rise to the counterterms for one-loop subrenormalization (see Figs. 1,2). We follow the renormalization at the one-loop level given in Refs. [28,[44][45][46], where details can be found. In the context of this paper, we only want to emphasize that on-shell (OS) renormalization is performed for the top-quark mass as well as for the scalar-top masses. This is different from the approach pursued, for example, in Ref. [40], where a DR renormalization has been employed. Using the OS scheme allows us to consistently combine our new correction terms with the hitherto available self-energies included in FeynHiggs.
Finally, at O(α t α s ), gluinos appear as virtual particles only at the two-loop level (hence, no renormalization for the gluinos is needed). The corresponding soft-breaking gluino mass parameter M 3 determines the gluino mass, Mg = M 3 .

The program SecDec
The calculation of the momentum-dependent two-loop corrections to the Higgs-boson masses at order O(α t α s ) involves two-loop two-point functions with up to four different masses, in addition to the mass scale given by the external momentum p 2 . For two-loop diagrams of propagator type, analytical results in four space-time dimensions are known only sparsely if different masses are occurring in the loops [54][55][56][57][58][59][60][61][62]. The integrals which are lacking analytical results can be classified into four different topologies, shown in Figure 4. We have calculated these integrals numerically using the program SecDec [63,64], where up to four different masses in 34 different mass configurations needed to be considered, with differences in the kinematic invariants of several orders of magnitude. Figure 4: Topologies which have been calculated numerically using SecDec.
The program SecDec is a publicly available tool [65] to calculate multi-loop integrals numerically. Dimensionally regulated poles are factorized by sector decomposition as described in Refs. [66,67], while kinematic thresholds are handled by a deformation of the integration contour into the complex plane, as described e.g. in Ref. [64]. The numerical integration is done using the Cuba library [68].
The program has also been extended to be able to calculate tensor integrals of any rank [69], and to process efficiently the evaluation of large ranges of kinematic points using the "multinumerics" feature of the program, which is of particular importance for the calculation presented here. This feature allows to produce input files for large sets of kinematic points automatically, and to process the evaluation of these points in parallel if several cores or a cluster are available, without repeating the algebraic part of the sector decomposition, which can be done once and for all. The evaluation of a single phase space point for the most complicated topology, to reach a relative accuracy of at least 10 −5 , ranges between 0.01 and 100 seconds on an Intel core i7 processor, where the larger timings are for points very close to a kinematic threshold.

Evaluation and implementation in the program FeynHiggs
The resulting new contributions to the neutral CP-even Higgs-boson self-energies, containing all momentum-dependent and additional constant terms, are assigned to the differences Note the tilde (not hat) onΣ (2) (0) which signifies that not only the self-energies are evaluated at zero external momentum but also the corresponding counterterms, following Refs. [25,26]. A finite shift ∆Σ (2) (0) therefore remains in the limit p 2 → 0 due to δM (2) , but at p 2 = 0 inΣ (2) ; for details see Eqs. (22) and (24) in the Appendix.
The numerical evaluation to derive the physical masses for h, H as the poles (real parts) of the dressed propagators proceeds on the basis of Eq. (7) in an iterative way.
• In a first step, the squared masses M 2 h,0 , M 2 H,0 are determined by solving Eq. (7) excluding the new terms ∆Σ (2) ab (p 2 ) from the self-energies.
• In a second step, the shifts ∆Σ • In the third step, Eq. (7) is solved again, now including the constant shifts c h(H) ab in the self-energies, to deliver the refined masses M h (with c h ab ) and M H (with c H ab ). This procedure can be repeated for improving the accuracy; numerically it turns out that going beyond the first iteration yields only marginal changes.
The corrections of Eq. (11) are incorporated in FeynHiggs by the following recipe, which is more general and in principle applicable also to the case of the complex MSSM with CP violation.
1. Determine Higgs masses M h i ,0 without the momentum-dependent terms of Eq. (11); the index i = 1, . . . , 4 enumerates the masses of h, H, A, H ± in the real MSSM. This is done by invoking the FeynHiggs mass-finder.
2. Compute the shifts c h k ab = ∆Σ

Run
FeynHiggs' mass-finder again including the c h k ab as constant shifts in the selfenergies to determine the refined Higgs masses M h and M H .
This procedure could conceivably be iterated until full self-consistency is reached; yet the resulting mass improvements turn out to be too small to justify extra CPU time.
On the technical side we added an interface for an external program to FeynHiggs which exports relevant model parameters to the external program's environment, currently: where the Ut ,1i denote the elements of the stop mixing matrix, α s (m t ) the running strong coupling at the scale m t , and G F the Fermi constant. The renormalization scale is defined within FeynHiggs as µ R = m t · FHscalefactor. Invocation of the external program is switched on by providing its path in the environment variable FHEXTSE. The program is executed from inside a temporary directory which is afterwards removed. The output (stdout) is scanned for lines of the form 'se@m c r c i ' which specify the correction c r + ic i [with c r = Re(c h k ab ), c i = Im(c h k ab )] to self-energy se in the computation of mass m, where m is one of Mh0, MHH, MA0, MHp, and se is one of h0h0, HHHH, A0A0, HmHp, h0HH, h0A0, HHA0, G0G0, h0G0, HHG0, A0G0, GmGp, HmGp, F1F1, F2F2, F1F2. The latter three, if given, substitute where α is the tree-level 2 × 2 neutral-Higgs mixing angle in Eq. (6). Self-energies not given are assumed zero.

Numerical results
We show results for the subtracted two-loop self-energies ∆Σ (2) ab (p 2 ) given in Eq. (11), as well as for the mass shifts i.e. the difference in the physical Higgs-boson masses evaluated including and excluding the newly obtained momentum-dependent two-loop corrections. This quantity, in particular ∆M h for the light CP-even Higgs boson, can directly be compared with the current experimental uncertainty as well as with the anticipated future ILC accuracy of [14], The leading to stop mass values of With the introduction of the momentum dependence, thresholds occur in the self-energy diagrams when the external momentum p = p 2 , in the time-like region, is such that a cut of the diagram would correspond to on-shell production of the massive particles of the cut propagators. The resulting imaginary parts enter in the search for the complex poles of the inverse propagator matrix of the Higgs bosons. Therefore it is interesting to study the behaviour of the real and imaginary parts of the self-energies. In Fig. 5 we show the momentum dependent parts of the renormalized two-loop self-energies in the physical basis, Eq. (11) for two different values of tan β, tan β = 5 and tan β = 20, at a fixed A-boson mass M A = 250 GeV. The data points are not connected by a line in order to show that each numerical point is obtained from a calculation of the 34 analytically unknown integrals with the program SecDec. The inlays in Fig. 5 magnify the region p 2 ≤ (125 GeV) 2 , where one can observe that for p 2 → 0, the subtracted self-energies are not exactly zero. As mentioned in subsection 2.4, this is due to the fact that the on-shell renormalization condition for the A-boson self-energy is defined differently with regard to the calculation without momentum dependence. The resulting constant contributions are additionally suppressed by factors sin 2 β, sin β cos β and cos 2 β appearing in the counter terms δm φ 1 φ 2 and δm 2(2) φ 2 , respectively, according to Eqs. (24) in the Appendix. The imaginary part is independent of the A-boson mass, as this mass parameter solely appears in the counterterms of DR renormalized quantities and the δM 2(2) A counterterm, where only the real part contributes. Therefore, the imaginary parts displayed in Fig. 5 do not contain additional constant terms. As to be expected, the imaginary parts are zero below the tt production threshold at p = 2 m t , which results from the fact that the top mass is the smallest mass appearing in the loops. Beyond this threshold, the imaginary parts are growing substantially with increasing p 2 . From these observations, the mass shifts in the region below the first threshold at p = 2 m t are expected not to be large.
Similar results, now including a variation of M A are shown in Fig. 6. In the upper plot for ∆Σ hh and in the middle plot for ∆Σ hH the solid lines depict M A ∼ 100 GeV, while the dashed lines are for M A ∼ 900 GeV. In these plots the light shading covers the range for tan β = 5, while the dark shading for tan β = 20. In the lower plot for ∆Σ HH we show results for M A ∼ 100, 250, 600, 900 GeV as solid, dotted, dot-dashed, dashed lines, respectively (and shading has been omitted). For ∆Σ hh at low p values only a small variation with M A can be observed. For p and M A large, the contributions to the self-energy are bigger. In ∆Σ hH larger effects are observed at smaller M A for both small and large p values. For ∆Σ HH , on the other hand, at low p values, large effects can be observed for large M A due to the aforementioned counterterm contribution ∼ δM     The dependence of the light CP-even Higgs-boson mass on Mg is analyzed in Fig. 9 for tan β = 5, 20 and M A = 250 GeV. Here we show as dashed lines the results for M h,0 (i.e. without the newly obtained momentum dependent two-loop corrections) and as solid lines the results for M h (i.e. including the new corrections). While a maximum of the Higgs-boson mass can be observed around Mg ∼ 800 GeV, in agreement with the original definition of the m max h scenario [72], a downward shift by more than 4 GeV is found for Mg ∼ 5 TeV. Such a strong effect is due to a (squared) logarithmic dependence of the O(α t α s ) corrections evaluated at p 2 = 0, as given in Eq. (73) of Ref. [26]. In Fig. 9 it can be seen that the size of the momentum dependent two-loop corrections similarly grows with Mg, reaching ∼ 400 MeV, as was shown above in Fig. 8. Consequently, the logarithmic dependence of the light CP-even Higgs-boson mass on the gluino mass that was found analytically for the O(α t α s ) corrections at p 2 = 0, is now also found numerically for the momentum dependent two-loop corrections.

Scenario 2: light stops
Scenario 2 is oriented at the "light-stop scenario" of Ref. [70] 2 . We use the following parameters: leading to stop mass values of Scenario 2 is analyzed with the same set of plots shown for scenario 1 in the previous subsection. The effects of the new momentum dependent two-loop contributions on the renormalized Higgs-boson self-energies, ∆Σ ab (p 2 ), are shown in Fig. 10. As before, we show the results separately for the real and imaginary parts of the self-energies. An additional threshold beyond the top-mass threshold appears at p = 2 mt 1 , where the discontinuity stems from the derivative of the imaginary part of the B 0 function(s). Analogously to scenario 1, the largest contributions in the region below 200 GeV arise in the real part of ∆Σ hh amounting to about 15 GeV 2 at p = 125 GeV, where the dependence on the value of tan β is rather weak.
The dependence of ∆Σ ab (p 2 ) on M A is shown in Fig. 11, using the same line styles as in Fig. 6. The curves show the same qualitative behaviour as in Fig. 10, exhibiting again the new threshold at p = 2 mt 1 . In general, outside the threshold region the effects in scenario 2 are slightly smaller than in scenario 1.
We now turn to the effects on the physical neutral CP-even Higgs boson masses. In Fig. 12 we show the results for ∆M h (upper plot) and ∆M H (lower plot) as a function of M A (with the same line styles as in Fig. 7). As can be expected from the previous figures, the effects on M h and M H are in general slightly smaller in scenario 2 than in scenario 1, where ∆M h still reaches the anticipated ILC accuracy, see Eq. (14). For ∆M H around the threshold p = 2 mt 1 the largest shift of ∼ −1 GeV can be observed. However, this shift is still below the anticipated mass resolution at the LHC [71].

Conclusions
We have presented results for the leading momentum-dependent O(α t α s ) contributions to the masses of neutral CP-even Higgs-bosons in the MSSM. They are obtained by calculating the corresponding contributions to the dressed Higgs-boson propagators obtained in the Feynman-diagrammatic approach using a mixed on-shell/DR renormalization scheme. In the Higgs sector a two-loop renormalization has to be carried out for the mass of the neutral Higgs bosons and the tadpole contributions. Furthermore, renormalization of the top/stop sector at O(α s ) is needed entering at the two-loop level via one-loop subrenormalization. The diagrams were generated with FeynArts and reduced to a set of basic integrals with the help of FormCalc and TwoCalc. The two-loop integrals which are analytically unknown have been calculated numerically with the program SecDec.
We have analyzed numerically the effect of the new momentum-dependent two-loop corrections on the predictions for the CP-even Higgs boson masses. This is particularly important for the interpretation of the scalar boson discovered at the LHC as the light CP-even Higgs state of the MSSM. While currently a precision below the level of ∼ 500 MeV is reached, a reduction by about an order of magnitude can be expected at the future e + e − International Linear Collider (ILC).
In our numerical analysis we found that the effects on the light CP-even Higgs boson mass, M h , depend strongly on the value of the gluino mass, Mg. For values of Mg ∼ 1.5 TeV corrections to M h of about −50 MeV are found, at the level of the anticipated future ILC accuracy. For very large gluino masses, Mg 4 TeV, on the other hand, substantially larger corrections are found, at the level of the current experimental accuracy. Consequently, this type of momentum dependent two-loop corrections should be taken into account in precision analyses interpreting the discovered Higgs boson in the MSSM.
For the heavy CP-even Higgs boson mass, M H , the effects are mostly below current and future anticipated accuracies. Only close to thresholds, e.g. around p = 2 mt 1 , larger corrections to M H around ∼ 1 GeV are found.
The new results of O(α t α s ) have been implemented into the program FeynHiggs. A detailed description of our calculation will be presented in a forthcoming publication [73].

Appendix: Renormalization and counterterms
Renormalization and calculation of the renormalized self-energies is performed in the (φ 1 , φ 2 ) basis, which has the advantage that the mixing angle α does not appear and expressions are in general simpler.
Field renormalization is perfomed by assigning one renormalization constant for each doublet, which can be expanded to one-and two-loop order according to The field renormalization constants appearing in (4) are then given by The mass counterterms δm 2(i) ab in (4) are derived from the Higgs potential, including the tadpoles, by the following parameter renormalization, 2 , tan β → tan β 1 + δ tan β (1) + δ tan β (2) .
The parameters T 1 and T 2 are the terms linear in φ 1 and φ 2 in the Higgs potential. The renormalization of the Z mass M Z does not contribute to the O(α s α t ) corrections we are pursuing here; it is listed, however, for completeness.
The basic renormalization constants for parameters and fields have to be fixed by renormalization conditions according to a renormalization scheme. Here we choose the on-shell scheme for the parameters and the DR scheme for field renormalization and give the expressions for the two-loop part.
The tadpole coefficients are chosen to vanish at all orders; hence their two-loop counterterms follow from T (2) 1,2 + δT where T 1 , T in terms of the A-boson unrenormalized self-energy Σ AA . The appearance of a non-zero momentum in the self-energy goes beyond the O(α t α s ) corrections evaluated in Refs. [25,26,31].
For the renormalization constants δZ H 1 , δZ H 2 and δ tan β several choices are possible, see the discussion in [47]. As shown there, the most convenient choice is a DR renormalization of δ tan β, δZ H 1 and δZ H 2 , which reads at the two-loop level δZ (2) δ tan β (2) = δ tan β (2)DR = 1 2 δZ (2) The term in Eq. (23c) is in general not the proper expression beyond one-loop order even in the DR scheme. For our approximation, however, with only the top Yukawa coupling at the two-loop level, it is the correct DR form [48].
The two-loop mass counterterms in the self-energies (4) are now expressed in terms of the parameter renormalization constants determined above as follows, Note that the Z-mass counterterm is kept for completeness; it does not contribute in the approximation of O(α s α t ) considered here.