*[arabic]figure *[arabic]table [toc]toc.1toc.1ContentsContentstoc.1

DESY–17–066

CERN–TH–2017–081

Complete two-loop QCD contributions to the
lightest Higgs-boson mass in the MSSM
with complex parameters

Sophia Borowka1***email<EMAIL_ADDRESS>Sebastian Paßehr2,3†††email:
<EMAIL_ADDRESS>and Georg Weiglein3‡‡‡email<EMAIL_ADDRESS>
1 Theoretical Physics Department, CERN, Geneva, Switzerland

2Sorbonne Université, CNRS, Laboratoire de Physique Théorique et Hautes
Énergies (LPTHE),
4 Place Jussieu, F–75252 Paris CEDEX 05, France

3Deutsches Elektronensynchrotron DESY,
Notkestraße 85, D–22607 Hamburg, Germany

###### Abstract

Higher-order corrections to the MSSM Higgs-boson masses are desirable for
accurate predictions currently testable at the LHC. By comparing the
prediction with the measured value of the discovered Higgs signal, viable
parameter regions can be inferred. For an improved theory accuracy, we compute
all two-loop corrections involving the strong coupling for the Higgs-boson
mass spectrum of the MSSM with complex parameters. Apart from the dependence
on the strong coupling, these contributions depend on the weak coupling and
Yukawa couplings, leading to terms of
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ and
$\mathcal{O}{\left(\sqrt{\alpha_{q_{1}}}\sqrt{\alpha_{q_{2}}}\alpha_{s}\right)}$,
($q_{1,2}=t,b,c,s,u,d$). The full dependence on the external momentum and all
relevant mass scales is taken into account. The calculation is performed in
the Feynman-diagrammatic approach which is flexible in the choice of the
employed renormalization scheme. For the phenomenological results presented
here, a renormalization scheme consistent with higher-order corrections
included in the code FeynHiggs is adopted. For the evaluation of the results,
a total of 513 two-loop two-point integrals with up to five different mass
scales are computed fully numerically using the program SecDec. A comparison
with existing results in the limit of real parameters and/or vanishing
external momentum is carried out, and the impact on the lightest Higgs-boson
mass is discussed, including the dependence on complex phases. The new results
will be included in the public code FeynHiggs.

###### Contents

  1. 1 Introduction
  2. 2 The Higgs sector of the MSSM with complex parameters
    1. 2.1 Tree-level relations for masses and mixing
    2. 2.2 Masses and mixing beyond lowest order
  3. 3 Calculation of the renormalized two-loop self-energies
    1. 3.1 Two-loop counterterms
    2. 3.2 Sub-loop renormalization
    3. 3.3 Resummation of $\tan\beta$-enhanced terms
  4. 4 Numerical evaluation of the self-energies
  5. 5 Numerical results for the Higgs mass spectrum
    1. 5.1 Comparison with earlier results
    2. 5.2 Scenario 1: $m_{h}^{\text{mod}}$-like
    3. 5.3 Scenario 2: large $\tan\beta$
    4. 5.4 Scenario 3: low $M_{H}$
  6. 6 Conclusions
  7. References

##  1 Introduction

Since the discovery of a signal in the Higgs-boson searches at the LHC
Aad:2012tfa ; Chatrchyan:2012ufa  with a mass around $125$ GeV, it is a prime
goal to reveal the detailed nature of the new particle. While with the present
experimental and theoretical uncertainties the measured properties of the
detected particle are compatible with the expectations for the Higgs boson of
the Standard Model (SM) Sperka:2018 ; Khachatryan:2016vau , other
interpretations corresponding to very different underlying physics are also in
agreement with the data. A crucial question in this context is in particular
whether the observed particle is part of an extended Higgs sector that would
be associated with a more general theoretical framework beyond the SM.

Within the theoretically well motivated Minimal Supersymmetric extension of
the SM (MSSM), the observed particle can be interpreted as a light state
within a richer spectrum of scalar particles.111Within the MSSM it is usually
assumed that the observed particle is associated with the lightest neutral
Higgs boson of the model; see Ref. Bechtle:2016kui  for a recent update on the
viability of the interpretation in terms of the next-to-lightest neutral Higgs
boson of the MSSM. The Higgs-boson sector of the MSSM consists of two complex
scalar doublets leading to five physical Higgs bosons and three (would-be)
Goldstone bosons. At the tree-level, the physical states are given by the
neutral $CP$-even bosons $h,H$ and the $CP$-odd state $A$, together with the
charged $H^{\pm}$ bosons. The Higgs sector at lowest order can be parametrized
in terms of the $A$-boson mass $m_{A}$ and the ratio of the vacuum expectation
values of the scalar doublets, $\tan\beta=\left.v_{2}\middle/v_{1}\right.$.
The MSSM with complex parameters (cMSSM) is of particular interest since it
provides new sources of $CP$-violation in addition to the $CP$-violating phase
of the SM. Thereby the Higgs sector is $CP$-conserving at the tree level, but
potentially large loop contributions involving complex parameters from other
supersymmetric (SUSY) sectors can lead to an admixture of the $CP$-even states
$h,\,H,$ and the $CP$-odd $A$ resulting in the mass eigenstates
$h_{1},h_{2},h_{3}$ Pilaftsis:1998pe ; Demir:1999hj ; Pilaftsis:1999qt ;
Heinemeyer:2001qd ; Frank:2006yh . In this case $m_{A}$ is no longer a useful
input parameter; instead the mass of the charged Higgs boson $m_{H^{\pm}}$ is
used. Besides the input parameter $m_{A}$ or $m_{H^{\pm}}$ all other Higgs-
boson masses are predicted quantities in the MSSM. The Higgs-boson masses and
mixings in the neutral sector are strongly affected by loop contributions.
Especially for the experimentally measured Higgs boson at about $125$ GeV a
sufficiently high accuracy of the theoretical computation is essential for
drawing reliable conclusions on the viability of the investigated region of
parameter space.

A large amount of work has been invested into calculating higher-order
corrections to the mass spectrum within the MSSM with real parameters
Degrassi:2002fi ; Heinemeyer:1998np ; Heinemeyer:1998jw ; Heinemeyer:1999be ;
Haber:1990aw ; Ellis:1990nz ; Okada:1990vk ; Okada:1990gg ; Ellis:1991zd ;
Sasaki:1991qu ; Chankowski:1991md ; Brignole:1992uf ; Hempfling:1993qq ;
Casas:1994us ; Dabelstein:1994hb ; Carena:1995bx ; Carena:1995wu ;
Pierce:1996zz ; Haber:1996fp ; Heinemeyer:1998kz ; Zhang:1998bm ;
Espinosa:1999zm ; Carena:2000dp ; Espinosa:2000df ; Espinosa:2001mm ;
Degrassi:2001yf ; Martin:2001vx ; Brignole:2001jy ; Dedes:2002dy ;
Martin:2002iu ; Dedes:2003km ; Heinemeyer:2004xw ; Martin:2003qz ;
Martin:2003it ; Allanach:2004rh ; Heinemeyer:2004gx ; Martin:2005eg ;
Martin:2005qm ; Harlander:2008ju ; Kant:2010tf ; Brignole:2002bz ;
Harlander:2017kuc ; Borowka:2014wla ; Degrassi:2014pfa ; Borowka:2015ura ;
Draper:2013oza ; Vega:2015fna ; Lee:2015uza ; Hahn:2013ria ; Bahl:2016brp ;
Bahl:2017aev ; Athron:2016fuq ; Bagnaschi:2014rsa ; Bagnaschi:2017xid ;
Martin:2002wn ; Martin:2004kr ; Martin:2007pg  as well as the MSSM with
complex parameters Pilaftsis:1998pe ; Demir:1999hj ; Pilaftsis:1999qt ;
Choi:2000wz ; Ibrahim:2000qj ; Ibrahim:2002zk ; Carena:2000yi ;
Heinemeyer:2001qd ; Martin:2002wn ; Martin:2004kr ; Frank:2006yh ;
Martin:2007pg ; Heinemeyer:2007aq ; Hollik:2014wea ; Hollik:2014bua ;
Passehr:2017ufr ; Goodsell:2016udb . The largest loop contributions originate
from the Yukawa sector due to the size of the top-quark Yukawa coupling
$h_{t}$, where $\alpha_{t}=\left.h_{t}^{2}\middle/(4\pi)\right.$. At the two-
loop level QCD corrections enter. The dominant contribution at the two-loop
level is given by the $\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$ terms
which are known for the MSSM with complex parameters Heinemeyer:1998jw ;
Heinemeyer:1998np ; Heinemeyer:1999be ; Heinemeyer:2007aq . Also the
$\mathcal{O}{\left(\alpha_{t}^{2}+\alpha_{t}\alpha_{b}+\alpha_{b}^{2}\right)}$
corrections are known for the case of complex parameters Hollik:2014wea ;
Hollik:2014bua ; Passehr:2017ufr . Restricting to the case of real parameters,
the momentum-dependent $\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$
corrections Borowka:2014wla ; Degrassi:2014pfa ; Borowka:2015ura  and the
contributions for the case where all Yukawa couplings except the one of the
top quark are neglected Degrassi:2014pfa  are known. While the phases of the
complex parameters affect the predictions for the Higgs-boson masses,
production cross sections Liebler:2016ceh  and decays Heinemeyer:1998yj ;
Lee:2003nta ; Williams:2011bu , they also induce $CP$-violating effects that
are constrained by other experiments. These concern in particular the electric
dipole moments Baron:2013eja ; Pospelov:2005pr ; Afach:2015sja ; Baker:2006ts
; Baker:2007df ; Serebrov:2015idv ; Graner:2016ses ; Yamanaka:2017mef . For
the usual convention where the phase of the mass of the electroweakinos,
$\phi_{M_{2}}$, is set to zero without loss of generality, the phase of the
parameter $\mu$ is constrained to be very close to zero or $\pi$. The other
important phases of the gluino mass, $\phi_{M_{3}}$, and the trilinear soft-
breaking parameters of the stops, $\phi_{A_{t}}$, and sbottoms,
$\phi_{A_{b}}$, are much less constrained. In particular, the bounds on the
phases of the trilinear soft-breaking parameters are significantly weaker for
the third generation than for the second and first generation.

In this article the full two-loop QCD corrections to the Higgs-boson masses
are presented for the general case of the MSSM with complex parameters. They
contain all previously computed results for the MSSM with real or complex
parameters. The contributions are comprised of the
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ terms, involving the electroweak
gauge coupling $\alpha$, and the
$\mathcal{O}{\left(\sqrt{\alpha_{q_{1}}}\sqrt{\alpha_{q_{2}}}\alpha_{s}\right)}$
terms, involving the Yukawa couplings $\alpha_{q_{1}}$, $\alpha_{q_{2}}$,
where $q_{1,2}=t,b,c,s,u,d$. Terms with mixed up- and down-type Yukawa
couplings only appear in conjunction with $m_{H^{\pm}}$ as input parameter.
Mixed contributions of
$\mathcal{O}{\left(\sqrt{\alpha}\sqrt{\alpha_{q_{1}}}\alpha_{s}\right)}$
involving one gauge coupling and one Yukawa coupling do not appear in the
final result. The results obtained here for the MSSM can furthermore be used
as an approximation for higher-order contributions within the NMSSM, as
discussed in Refs. Drechsel:2016jdg ; Drechsel:2016htw ; Domingo:2017rhb . The
computation carried out below makes use of previously developed tools
Weiglein:1993hd ; Borowka:2014wla ; Hahn:2015gaa . The momentum-dependent two-
loop integrals appearing in the two-loop QCD corrections are evaluated with an
adapted version of SecDec 2 Carter:2010hi ; Borowka:2012yc ; Borowka:2013cma .
For the numerical analysis the new contributions are combined with the full
one-loop result Frank:2006yh  and the leading
$\mathcal{O}{\left(\alpha_{t}^{2}\right)}$ terms Hollik:2014wea ;
Hollik:2014bua  in the Feynman-diagrammatic approach for complex parameters,
available through the public program FeynHiggs Heinemeyer:1998np ;
Degrassi:2002fi ; Frank:2006yh ; Heinemeyer:1998yj ; Hahn:2010te . In deriving
the new contributions the renormalization scheme of Ref. Frank:2006yh  at the
one-loop level has been adopted and applied to the case of the
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ contributions. This ensures that
the obtained analytical results for the renormalized two-loop self-energies
can consistently be incorporated into FeynHiggs. In the results presented in
this paper no resummation of higher-order logarithmic contributions as
obtained in Refs. Hahn:2013ria ; Draper:2013oza ; Vega:2015fna ; Lee:2015uza ;
Bahl:2016brp ; Athron:2016fuq ; Bahl:2017aev  has been included. The
combination of resummed higher-order logarithmic contributions with the
results obtained in the present paper will be addressed in future work. In the
numerical analysis below, we show results for the masses of the three neutral
Higgs bosons of the MSSM with complex parameters and their phase dependence,
with a particular focus on those results which are phenomenologically most
relevant.

The paper is organized as follows: section 2 provides the theoretical
framework for the calculation and renormalization of the Feynman diagrams that
is used to arrive at expressions for the dressed propagators of the Higgs
sector up to the two-loop level. The calculation of the unrenormalized self-
energies and the construction of the two-loop counterterms are described in
section 3. Details on the numerical evaluation of the momentum-dependent two-
loop integrals are given in section 4, whereas the impact of the new
contributions on the Higgs-boson masses is discussed in section 5. The
conclusions are given in section 6.

##  2 The Higgs sector of the MSSM with complex parameters

###  2.1 Tree-level relations for masses and mixing

The two scalar $SU(2)$-doublets are conventionally expressed in terms of their
components as follows,

$\displaystyle\mathcal{H}_{1}$
$\displaystyle=\begin{pmatrix}v_{1}+\frac{1}{\sqrt{2}}(\phi_{1}-{i\mkern
1.0mu}\chi_{1})\\\ -\phi^{-}_{1}\end{pmatrix},$
$\displaystyle\quad\mathcal{H}_{2}$ $\displaystyle=\mathrm{e}^{{i\mkern
1.0mu}\xi}\begin{pmatrix}\phi^{+}_{2}\\\
v_{2}+\frac{1}{\sqrt{2}}(\phi_{2}+{i\mkern 1.0mu}\chi_{2})\end{pmatrix}\,,$
(1)

with the relative phase $\xi$. The Higgs potential can be written as a
polynomial in the field components,

$\displaystyle\begin{split}V_{H}&=-T_{\phi_{1}}\phi_{1}-T_{\phi_{2}}\phi_{2}-T_{\chi_{1}}\chi_{1}-T_{\chi_{2}}\chi_{2}\\\
&\quad+\frac{1}{2}\begin{pmatrix}\phi_{1},&\phi_{2},&\chi_{1},&\chi_{2}\end{pmatrix}\begin{pmatrix}\mathbf{M}_{\phi}&\mathbf{M}_{\phi\chi}\\\
\mathbf{M}_{\phi\chi}^{\dagger}&\mathbf{M}_{\chi}\end{pmatrix}\begin{pmatrix}\phi_{1},&\phi_{2},&\chi_{1},&\chi_{2}\end{pmatrix}^{\text{T}}+\begin{pmatrix}\phi^{-}_{1},&\phi^{-}_{2}\end{pmatrix}\mathbf{M}_{\phi^{\pm}}\begin{pmatrix}\phi^{+}_{1}\\\
\phi^{+}_{2}\end{pmatrix}+\dots\ ,\end{split}$ (2)

where terms of third and fourth power in the fields have been omitted, and the
relations $\phi^{-}_{1}=\left(\phi^{+}_{1}\right)^{\dagger}$ and
$\phi^{-}_{2}=\left(\phi^{+}_{2}\right)^{\dagger}$ have been used. Explicit
expressions for the tadpole coefficients $T$ and for the mass matrices
$\mathbf{M}$ can be found in Ref. Frank:2006yh . They are parametrized by the
phase $\xi$, the real SUSY-breaking quantities
$m_{1,2}^{2}=\tilde{m}_{1,2}^{2}+\lvert\mu\rvert^{2}$, and the complex SUSY-
breaking quantity $m_{12}^{2}$. With the help of a Peccei–Quinn transformation
Peccei:1977hh  the parameter $m_{12}^{2}$ can be redefined such that its phase
vanishes Dimopoulos:1995kn , leaving only the phase $\xi$ as a potential
source of $CP$-violation at tree level. The requirement of minimizing $V_{H}$
at the vacuum expectation values $v_{1}$ and $v_{2}$ is equivalent to the
requirement of vanishing tadpoles of the physical fields, which in turn
implies the condition $\xi=0$ at tree level. As a consequence, the Higgs
sector of the MSSM is $CP$-conserving at lowest order. This implies in Eq. (2)
that $\mathbf{M}_{\phi\chi}$ is equal to zero, and $\phi_{1,2}$ do not mix
with $\chi_{1,2}$ at tree-level.

The remaining $(2\times 2)$-matrices $\mathbf{M}_{\phi}$, $\mathbf{M}_{\chi}$,
$\mathbf{M}_{\phi^{\pm}}$ can be transformed into the mass eigenstate basis
with the help of orthogonal matrices $D(x)$, using the abbreviations
$s_{x}\equiv\sin{x}$, $c_{x}\equiv\cos{x}$,

$\displaystyle D{\left(x\right)}$
$\displaystyle=\begin{pmatrix}-s_{x}&c_{x}\\\ c_{x}&s_{x}\end{pmatrix},$
$\displaystyle\quad\begin{pmatrix}h\\\ H\end{pmatrix}$
$\displaystyle=D(\alpha)\begin{pmatrix}\phi_{1}\\\ \phi_{2}\end{pmatrix},$
$\displaystyle\quad\begin{pmatrix}A\\\ G\end{pmatrix}$
$\displaystyle=D(\beta_{n})\begin{pmatrix}\chi_{1}\\\ \chi_{2}\end{pmatrix},$
$\displaystyle\quad\begin{pmatrix}H^{\pm}\\\ G^{\pm}\end{pmatrix}$
$\displaystyle=D(\beta_{c})\begin{pmatrix}\phi^{\pm}_{1}\\\
\phi^{\pm}_{2}\end{pmatrix}.$ (3)

The Higgs potential in this basis can be expressed as follows,

$\displaystyle\begin{split}V_{H}&=-T_{h}\,h-T_{H}\,H-T_{A}\,A-T_{G}\,G\\\
&\quad+\frac{1}{2}\begin{pmatrix}h,&H,&A,&G\end{pmatrix}\mathbf{M}_{hHAG}\begin{pmatrix}h,&H,&A,&G\end{pmatrix}^{\text{T}}+\begin{pmatrix}H^{-},&G^{-}\end{pmatrix}\mathbf{M}_{H^{\pm}G^{\pm}}\begin{pmatrix}H^{+}\\\
G^{+}\end{pmatrix}+\dots\ \end{split}$ (4)

with the tadpole coefficients $T_{h,H,A,G}$ and the mass matrices

$\displaystyle\mathbf{M}_{hHAG}$
$\displaystyle=\begin{pmatrix}m^{2}_{h}&m^{2}_{hH}&m^{2}_{hA}&m^{2}_{hG}\\\
m^{2}_{hH}&m^{2}_{H}&m^{2}_{HA}&m^{2}_{HG}\\\
m^{2}_{hA}&m^{2}_{HA}&m^{2}_{A}&m^{2}_{AG}\\\
m^{2}_{hG}&m^{2}_{HG}&m^{2}_{AG}&m^{2}_{G}\end{pmatrix},\qquad\mathbf{M}_{H^{\pm}G\pm}\,=\begin{pmatrix}m^{2}_{H^{\pm}}&m^{2}_{H^{-}G^{+}}\\\
m^{2}_{G^{-}H^{+}}&m^{2}_{G^{\pm}}\end{pmatrix};$ (5)

explicit expressions for the entries are given in Ref. Frank:2006yh . The
tadpole terms in Eq. (4) are zero at the tree level, but they enter the
predictions for the Higgs-boson masses at higher orders. As mentioned above,
the ellipses denote terms of higher power in the fields which are not relevant
in our calculation.

After applying the minimization conditions to Eqs. (5), the mass matrices can
be brought into canonical form222We use a lower-case $m$ for the Higgs-boson
masses at the tree level.

$\displaystyle\mathbf{M}_{hHAG}^{(0)}$
$\displaystyle=\,\mathrm{diag}\left(m_{h}^{2},\,m_{H}^{2},\,m_{A}^{2},\,m_{G}^{2}\right),\quad\mathbf{M}_{H^{\pm}G^{\pm}}^{(0)}\,=\,\mathrm{diag}\left(m_{H^{\pm}}^{2},\,m_{G^{\pm}}^{2}\right),$
(6)

for $\beta=\beta_{n}=\beta_{c}$, with $\beta\in\left[0,\pi/2\right)$ given in
terms of the vacuum expectation values,

$\displaystyle\tan\beta$ $\displaystyle\equiv t_{\beta}=\frac{v_{2}}{v_{1}}\
,$ (7)

and for the second mixing angle $\alpha\in\left[-\pi/2,0\right)$ determined by

$\displaystyle\tan(2\alpha)$
$\displaystyle=\frac{m_{A}^{2}+m_{Z}^{2}}{m_{A}^{2}-m_{Z}^{2}}\,\tan(2\beta)\,.$
(8)

The Goldstone bosons are massless333The Goldstone bosons can acquire a non-
zero mass value by gauge fixing., $m_{G^{\pm}}=m_{G}=0$. The masses
$m_{H^{\pm}},m_{A},m_{h},m_{H}$ fullfil the relations

$\displaystyle m_{H^{\pm}}^{2}$ $\displaystyle=m_{A}^{2}+M_{W}^{2}\,,$ (9a)
$\displaystyle m_{h,\,H}^{2}$
$\displaystyle=\frac{1}{2}\left(m_{A}^{2}+M_{Z}^{2}\mp\sqrt{\left(m_{A}^{2}+M_{Z}^{2}\right)^{2}-4m_{A}^{2}M_{Z}^{2}\,c_{2\beta}^{2}}\right),$
(9b)

including the vector-boson masses $M_{W}$ and $M_{Z}$. Given the relation in
Eq. (9a), both $m_{A}$ and $m_{H^{\pm}}$ can be chosen as input parameter.

At lowest order, the irreducible two-point vertex functions of the neutral
Higgs sector

$\displaystyle\Gamma^{(0)}_{hHAG}(p^{2})$
$\displaystyle=\,i\,\Big{[}p^{2}\textbf{1}-\mathbf{M}^{(0)}_{hHAG}\Big{]}$
(10)

are diagonal, and the entries of the mass matrices in Eq. (6) provide the
poles of the diagonal lowest-order propagators

$\displaystyle\Delta^{(0)}_{hHAG}(p^{2})$
$\displaystyle=\,-\Big{[}\Gamma^{(0)}_{hHAG}(p^{2})\Big{]}^{-1}\,.$ (11)

###  2.2 Masses and mixing beyond lowest order

Going beyond leading order, the irreducible two-point functions are dressed by
adding the matrix $\mathbf{\hat{\Sigma}}_{hHAG}$ of the renormalized diagonal
and non-diagonal self-energies for the $h,H,A,G$ fields up to the considered
order,

$\displaystyle p^{2}\textbf{1}-\mathbf{M}^{(0)}_{hHAG}$
$\displaystyle\quad\to\quad
p^{2}\textbf{1}-\mathbf{M}^{(0)}_{hHAG}+\mathbf{\hat{\Sigma}}_{hHAG}(p^{2})\;\equiv\;p^{2}\textbf{1}-\mathbf{M}_{hHAG}(p^{2})\,,$
(12)

yielding the full renormalized two-point vertex function

$\displaystyle\hat{\Gamma}_{hHAG}(p^{2})=\,i\,\Big{[}p^{2}\textbf{1}-\mathbf{M}_{hHAG}\Big{]}\;.$
(13)

The latter generally contains a mixing of all fields with equal quantum
numbers. The dressed propagators are obtained by inverting the matrix
$\hat{\Gamma}_{hHAG}(p^{2})$.

Truncating the perturbative expansion at the two-loop level, the momentum-
dependent corrections to the neutral Higgs-boson mass matrices in Eq. (12) are
given by

$\displaystyle\mathbf{M}_{hHAG}^{(2)}(p^{2})$
$\displaystyle=\mathbf{M}_{hHAG}^{(0)}-\mathbf{\hat{\Sigma}}_{hHAG}^{(1)}(p^{2})-\mathbf{\hat{\Sigma}}_{hHAG}^{(2)}(p^{2})\
.$ (14)

For the MSSM with complex parameters, the one-loop self-energies are
completely known Frank:2006yh , and the leading two-loop
$\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$ and
$\mathcal{O}{\left(\alpha_{t}^{2}+\alpha_{t}\alpha_{b}+\alpha_{b}^{2}\right)}$
contributions have been obtained in the approximation of zero external
momentum Heinemeyer:2007aq ; Hollik:2014wea ; Hollik:2014bua ; Passehr:2017ufr
. In the case of the MSSM with real parameters also the momentum-dependent
corrections of $\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$ are known
Borowka:2014wla ; Degrassi:2014pfa ; Borowka:2015ura . The remaining QCD
contributions at the two-loop level are completed within this paper. These
contributions comprise terms of the
$\mathcal{O}{\left(\alpha_{x}\alpha_{s}\right)}$, where $\alpha_{x}$ is either
the gauge coupling $\alpha$ or the Yukawa coupling $\alpha_{q}$ with
$q=\\{u,d,s,c,b,t\\}$. We neglect CKM mixing for those contributions.

In order to obtain the physical Higgs-boson masses from the dressed
propagators at the considered order, it is sufficient to explicitly derive the
entries of the $(3\times 3)$-submatrix of Eq. (14) corresponding to the
$(h,H,A)$-components. A mixing of the neutral Higgs bosons with the Goldstone
boson, as well as Goldstone–$Z$ mixing, yields subleading two-loop
contributions to the Higgs-boson masses that are not of
$\mathcal{O}{\left(\alpha_{x}\alpha_{s}\right)}$.

The masses of the three neutral Higgs bosons are obtained from the real parts
of the complex poles of the $(h,H,A)$-propagator matrix. They are obtained as
the zeroes of the determinant of the renormalized two-loop two-point vertex
function,444We use an uppercase $M$ for the Higgs masses at higher order.

$\displaystyle\left.\operatorname{det}\hat{\Gamma}^{(2)}_{hHA}{\left(p^{2}\right)}\right|_{p^{2}\,=\,M_{j}^{2}\,-\,{i\mkern
1.0mu}\,M_{j}\,\Gamma_{j}}$ $\displaystyle=0,$
$\displaystyle\quad\hat{\Gamma}^{(2)}_{hHA}{\left(p^{2}\right)}$
$\displaystyle={i\mkern
1.0mu}\left[p^{2}{\textbf{1}}-\mathbf{M}_{hHA}^{(2)}{\left(p^{2}\right)}\right],$
$\displaystyle\quad j\in\\{h,H,A\\}\,,$ (15)

with $\mathbf{M}_{hHA}^{(2)}$ being the corresponding $(3\times 3)$-submatrix
of Eq. (14). The impact of the self-energies on the mixing and couplings of
the various Higgs bosons to other (MS)SM particles can be obtained with the
same formalism as described in Refs. Frank:2006yh ; Fuchs:2016swt .

##  3 Calculation of the renormalized two-loop self-energies

The renormalized two-loop self-energies can be written as

$\displaystyle\mathbf{\hat{\Sigma}}_{hHA}^{(2)}(p^{2})$
$\displaystyle=\mathbf{\Sigma}_{hHA}^{(2)}(p^{2})-\delta^{(2)}\mathbf{M}^{\mathbf{Z}}_{hHA}\,,$
(16)

where $\mathbf{\Sigma}_{hHA}^{(2)}$ denotes the unrenormalized self-energies
corresponding to the sum of genuine two-loop diagrams and one-loop diagrams
with counterterm insertions. The symbol
$\delta^{(2)}\mathbf{M}^{\mathbf{Z}}_{hHA}$ comprises all two-loop
counterterms resulting from parameter and field renormalization.

The contributing types of Feynman diagrams for the calculation of the full
two-loop QCD corrections entering Eq. (16) are depicted in Fig. 1. The
diagrams of the topologies $12$, $14$ and $15$ contribute only if all squarks
have the same flavor; couplings with different flavors vanish since the color
sum is equal to zero in that case. The diagrammatic calculation has been
performed with the help of FeynArts Kublbeck:1990xc ; Hahn:2000kx  in
generating the Feynman diagrams, and TwoCalc Weiglein:1993hd  and Reduze
vonManteuffel:2012np  for the two-loop trace evaluation and tensor reduction.
The one-loop renormalization constants have been computed with the help of
FormCalc Hahn:1998yk .

###  3.1 Two-loop counterterms

In order to obtain the renormalized self-energies in Eq. (16), counterterms
have to be introduced up to second order in the loop expansion, for the
tadpoles

$\displaystyle T_{i}$ $\displaystyle\rightarrow
T_{i}+\delta^{(1)}T_{i}+\delta^{(2)}T_{i}\,,\quad i=h,\,H,\,A\,,$ (17)

and for the mass matrices of Eq. (4)

$\displaystyle\mathbf{M}_{hHA}$
$\displaystyle\rightarrow\mathbf{M}_{hHA}^{(0)}+\delta^{(1)}\mathbf{M}_{hHA}+\delta^{(2)}\mathbf{M}_{hHA}\
,$ (18a) $\displaystyle\delta^{(k)}\mathbf{M}_{hHA}$
$\displaystyle=\begin{pmatrix}\delta^{(k)}m_{h}^{2}&\delta^{(k)}m_{hH}^{2}&\delta^{(k)}m_{hA}^{2}\\\
\delta^{(k)}m_{Hh}^{2}&\delta^{(k)}m_{H}^{2}&\delta^{(k)}m_{HA}^{2}\\\
\delta^{(k)}m_{Ah}^{2}&\delta^{(k)}m_{AH}^{2}&\delta^{(k)}m_{A}^{2}\end{pmatrix},$
(18b) $\displaystyle m^{2}_{H^{\pm}}$ $\displaystyle\rightarrow
m_{H^{\pm}}^{2}+\,\delta^{(1)}m^{2}_{H^{\pm}}+\,\delta^{(2)}m^{2}_{H^{\pm}}\,.$
(18c)

The two-loop counterterms of $\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ have
the same structure as the corresponding one-loop counterterms. They are listed
here for completeness and to fix our notation.

In order to ensure the correct form of the counterterms for the mass matrices,
the rotation angles $\beta_{n}$ and $\beta_{c}$ from Eqs. (3) have to be
distinguished from $\beta$ in Eq. (7). Whereas no renormalization is needed
for $\alpha$, $\beta_{n}$ and $\beta_{c}$, a counterterm associated with
$\beta$ of the form $\beta\to\beta+\delta\beta$ is required, in accordance
with the renormalization of $\tan\beta$,

$\displaystyle t_{\beta}$ $\displaystyle\rightarrow
t_{\beta}+\delta^{(1)}t_{\beta}+\delta^{(2)}t_{\beta}\ .$ (19)

In the resulting expressions for the counterterm matrices, the identification
$\beta_{c}=\beta_{n}=\beta$ can be made, see Ref. Frank:2006yh  for details of
the analogous treatment at the one-loop order (note that a different
convention for the counterterm of $t_{\beta}$ is used in Ref. Frank:2006yh ).
A complete list of the two-loop counterterms is given in the Appendix of Ref.
Hollik:2014bua .

Figure 1:  Types of two-loop self-energy diagrams for the neutral Higgs
bosons. One-loop counterterm insertions are denoted by a cross.
$\Phi_{i}=h,\,H,\,A$; $\;\tilde{q}\neq\tilde{o}$. Topologies $11$ and $13$
contain a one-point loop with a mass-less gluon and are therefore equal to
zero.

In addition to the parameter renormalization described previously, the field-
renormalization constants $\delta^{(1)}Z_{\mathcal{H}_{i}}$ and
$\delta^{(2)}Z_{\mathcal{H}_{i}}$ are introduced at the one-loop and two-loop
order (restricting the latter to the contributions of
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$) for each of the scalar doublets
of Eqs. (1) through the transformation

$\displaystyle\mathcal{H}_{i}\rightarrow\sqrt{Z_{\mathcal{H}_{i}}}\mathcal{H}_{i}$
$\displaystyle=\left[1+\frac{1}{2}\delta^{(1)}Z_{\mathcal{H}_{i}}+\frac{1}{2}\delta^{(2)}Z_{\mathcal{H}_{i}}\right]\mathcal{H}_{i}\,.$
(20)

The field-renormalization constants in the mass-eigenstate basis of Eqs. (3)
are obtained via

$\displaystyle\begin{pmatrix}h\\\ H\end{pmatrix}$ $\displaystyle\rightarrow
D(\alpha)\begin{pmatrix}\sqrt{Z_{\mathcal{H}_{1}}}&0\\\
0&\sqrt{Z_{\mathcal{H}_{2}}}\end{pmatrix}D(\alpha)^{-1}\begin{pmatrix}h\\\
H\end{pmatrix}$ $\displaystyle\equiv\mathbf{Z}_{hH}\begin{pmatrix}h\\\
H\end{pmatrix},$ (21a) $\displaystyle\begin{pmatrix}A\\\ G\end{pmatrix}$
$\displaystyle\rightarrow
D(\beta_{n})\begin{pmatrix}\sqrt{Z_{\mathcal{H}_{1}}}&0\\\
0&\sqrt{Z_{\mathcal{H}_{2}}}\end{pmatrix}D(\beta_{n})^{-1}\begin{pmatrix}A\\\
G\end{pmatrix}$ $\displaystyle\equiv\begin{pmatrix}Z_{AA}&Z_{AG}\\\
Z_{GA}&Z_{GG}\end{pmatrix}\begin{pmatrix}A\\\
G\end{pmatrix}\;\equiv\;\mathbf{Z}_{AG}\begin{pmatrix}A\\\ G\end{pmatrix},$
(21b) $\displaystyle\begin{pmatrix}H^{\pm}\\\ G^{\pm}\end{pmatrix}$
$\displaystyle\rightarrow
D(\beta_{c})\begin{pmatrix}\sqrt{Z_{\mathcal{H}_{1}}}&0\\\
0&\sqrt{Z_{\mathcal{H}_{2}}}\end{pmatrix}D(\beta_{c})^{-1}\begin{pmatrix}H^{\pm}\\\
G^{\pm}\end{pmatrix}$
$\displaystyle\equiv\mathbf{Z}_{H^{\pm}G^{\pm}}\begin{pmatrix}H^{\pm}\\\
G^{\pm}\end{pmatrix}.$ (21c)

The matrices $\mathbf{Z}_{ij}$ can be expanded as

$\displaystyle\mathbf{Z}_{ij}=\mathbf{1}+\delta^{(1)}\mathbf{Z}_{ij}+\delta^{(2)}\mathbf{Z}_{ij}\,.$
(22)

The required one-loop expressions for the entries in the
$\delta^{(1)}\mathbf{Z}_{ij}$-matrices are given in Ref. Frank:2006yh ; the
corresponding set of two-loop expressions is given in Ref. Hollik:2014bua .

The genuine two-loop counterterms $\delta^{(2)}\mathbf{M}^{\mathbf{Z}}_{hHA}$
of Eq. (16) can now be summarized as

$\displaystyle\begin{split}\delta^{(2)}\mathbf{M}^{\mathbf{Z}}_{hHA}&=\delta^{(2)}\mathbf{M}_{hHA}+\begin{pmatrix}\delta^{(2)}\mathbf{Z}_{hH}^{T}&\mathbf{0}\\\
\mathbf{0}&\delta^{(2)}Z_{AA}\end{pmatrix}\left(\mathbf{M}_{hHA}^{(0)}-p^{2}\textbf{1}\right)\\\
&\quad+\left(\mathbf{M}_{hHA}^{(0)}-p^{2}\textbf{1}\right)\begin{pmatrix}\delta^{(2)}\mathbf{Z}_{hH}&\mathbf{0}\\\
\mathbf{0}&\delta^{(2)}Z_{AA}\end{pmatrix}\ ,\end{split}$ (23)

where the required two-loop mass counterterms read

$\displaystyle\begin{split}\delta^{(2)}m_{h}^{2}&=c_{\alpha-\beta}^{2}\,\delta^{(2)}m_{A}^{2}+s_{\alpha+\beta}^{2}\,\delta^{(2)}m_{Z}^{2}+c_{\beta}^{2}\,\delta^{(2)}t_{\beta}\left(s_{2(\alpha-\beta)}\,m_{A}^{2}+s_{2(\alpha+\beta)}\,m_{Z}^{2}\right)\\\
&\quad+\frac{e\,s_{\alpha-\beta}}{2\,M_{W}\,s_{\text{w}}}\Biggl{[}\left(1+c_{\alpha-\beta}^{2}\right)\delta^{(2)}T_{h}+s_{\alpha-\beta}\,c_{\alpha-\beta}\,\delta^{(2)}T_{H}\Biggr{]}\
,\end{split}$ (24a)
$\displaystyle\begin{split}\delta^{(2)}m_{H}^{2}&=s_{\alpha-\beta}^{2}\,\delta^{(2)}m_{A}^{2}+c_{\alpha+\beta}^{2}\,\delta^{(2)}m_{Z}^{2}-c_{\beta}^{2}\,\delta^{(2)}t_{\beta}\left(s_{2(\alpha-\beta)}\,m_{A}^{2}+s_{2(\alpha+\beta)}\,m_{Z}^{2}\right)\\\
&\quad-\frac{e\,c_{\alpha-\beta}}{2\,M_{W}\,s_{\text{w}}}\Biggl{[}\left(1+s_{\alpha-\beta}^{2}\right)\delta^{(2)}T_{H}+c_{\alpha-\beta}\,s_{\alpha-\beta}\,\delta^{(2)}T_{h}\Biggr{]}\
,\end{split}$ (24b) $\displaystyle\delta^{(2)}m_{A}^{2}$
$\displaystyle=\delta^{(2)}m_{H^{\pm}}^{2}-\delta^{(2)}m_{W}^{2}\ ,$ (24c)
$\displaystyle\begin{split}\delta^{(2)}m_{hH}^{2}&=\frac{1}{2}\left(s_{2\left(\alpha-\beta\right)}\,\delta^{(2)}m_{A}^{2}-s_{2\left(\alpha+\beta\right)}\,\delta^{(2)}m_{Z}^{2}\right)-c_{\beta}^{2}\,\delta^{(2)}t_{\beta}\left(c_{2(\alpha-\beta)}\,m_{A}^{2}+c_{2(\alpha+\beta)}\,m_{Z}^{2}\right)\\\
&\quad+\frac{e}{2\,M_{W}\,s_{\text{w}}}\Biggl{[}s_{\alpha-\beta}^{3}\,\delta^{(2)}T_{H}-c_{\alpha-\beta}^{3}\,\delta^{(2)}T_{h}\Biggr{]}\
,\end{split}$ (24d) $\displaystyle\delta^{(2)}m_{hA}^{2}$
$\displaystyle=\frac{e}{2\,M_{W}\,s_{\text{w}}}\,s_{\alpha-\beta}\,\delta^{(2)}T_{A}\
,$ (24e) $\displaystyle\delta^{(2)}m_{HA}^{2}$
$\displaystyle=-\frac{e}{2\,M_{W}\,s_{\text{w}}}\,c_{\alpha-\beta}\,\delta^{(2)}T_{A}\
.$ (24f)

The entries of $\delta^{(2)}\mathbf{M}_{hHA}$ that are not listed here are
determined by symmetry. When replacing $\delta^{(2)}\to\delta$ they are
formally equal to the one-loop counterterms listed in Eqs. (53) of Ref.
Frank:2006yh  (up to the different convention for the counterterm of
$t_{\beta}$ used there).

The two-loop renormalization constants of Eqs. (17)–(24) are fixed by
extending the renormalization scheme of Ref. Frank:2006yh  from the one-loop
to the two-loop order:

  * •

The tadpole counterterms $\delta^{(2)}T_{i}$ are fixed by requiring that the
minimum of the Higgs potential is not shifted, which means that the tadpole
coefficients have to vanish at each order. At the two-loop level, the
condition reads

$\displaystyle T_{i}^{(2)}+\delta^{(2)}T_{i}$ $\displaystyle=0\ ,\quad
i=h,\,H,\,A\ ,$ (25)

where the $T_{i}^{(2)}$ denote the unrenormalized one-point functions at two-
loop order, see Fig. 4 for the contributing two-loop diagrams. The
aforementioned relation for the mixing angles $\beta_{n}=\beta_{c}=\beta$ is a
consequence of the tadpole conditions $T_{i}=0$ at lowest order.

Figure 2:  Types of two-loop tadpole diagrams contributing to $T^{(2)}_{i}$.
One-loop counterterm insertions are denoted by a cross. $\Phi_{i}=h,\,H,\,A$;
$\tilde{q}\neq\tilde{o}$. Topology $5$ contains a one-point loop with a mass-
less gluon and is therefore equal to zero.

Figure 3:  Types of two-loop self-energy diagrams for the charged Higgs
bosons. One-loop counterterm insertions are denoted by a cross. $\;q\neq o$,
$\;\tilde{q}\neq\tilde{o}$. Topologies $11$ and $13$ contain a one-point loop
with a mass-less gluon and are therefore equal to zero.

Figure 4:  Additional types of two-loop self-energy diagrams for the gauge
bosons besides the ones in analogy to Fig. 1 and 4.
$\;\tilde{q}\neq\tilde{o}$.

  * •

The charged Higgs-boson mass $m_{H^{\pm}}$ is the only independent mass
parameter of the Higgs sector and is used as an input quantity. Accordingly,
the corresponding mass counterterm is fixed by an independent renormalization
condition, chosen as on-shell, given by

$\displaystyle\Re\mathfrak{e}{\left[\hat{\Sigma}_{H^{\pm}}^{(2)}(m_{H^{\pm}}^{2})\right]}=0\text{
.}$ (26)

The renormalized charged-Higgs self-energy at the two-loop level can be
expressed in terms of the unrenormalized charged self-energy and its
respective counterterms

$\displaystyle\hat{\Sigma}_{H^{\pm}}^{(2)}(m_{H^{\pm}})=\Sigma_{H^{\pm}}^{(2)}(m_{H^{\pm}}^{2})-\delta^{(2)}m^{2}_{H^{\pm}}\,,$
(27)

leading to the mass counterterm

$\displaystyle\delta^{(2)}m_{H^{\pm}}^{2}$
$\displaystyle=\Re\mathfrak{e}{\left[\Sigma_{H^{\pm}}^{(2)}{\left(m_{H^{\pm}}^{2}\right)}\right]}$
(28)

when applying the on-shell condition. The contributing Feynman diagrams are
shown in Fig. 4. As we neglect flavor mixing $q$, $o$, $\tilde{q}$ and
$\tilde{o}$ always belong to the same generation. As a consequence, the
vertices with four squarks in topologies $12$, $14$ and $15$ are only non-zero
when all adjacent fields are of the same generation.

  * •

The field-renormalization constants of the Higgs mass eigenstates in Eq. (21)
are combinations of the basic doublet-field renormalization constants
$\delta^{(2)}Z_{\mathcal{H}_{1}}$ and $\delta^{(2)}Z_{\mathcal{H}_{2}}$, which
are fixed by the UV-divergent parts of the derivatives of the corresponding
self-energies,

$\displaystyle\delta^{(2)}Z_{\mathcal{H}_{1}}$
$\displaystyle=-\left[\frac{d\Sigma^{(2)}_{HH}(p^{2})}{dp^{2}}\right]_{\alpha=0}^{\text{div}},$
$\displaystyle\quad\delta^{(2)}Z_{\mathcal{H}_{2}}$
$\displaystyle=-\left[\frac{d\Sigma^{(2)}_{hh}(p^{2})}{dp^{2}}\right]_{\alpha=0}^{\text{div}}.$
(29)

  * •

Also $t_{\beta}$ is renormalized by a purely UV-divergent counterterm, which
was shown to be a convenient choice Freitas:2002um  (see also Refs.
Sperling:2013eva ; Sperling:2013xqa ). Alternative process-dependent
definitions for the renormalization of $t_{\beta}$ can be found in Ref.
Baro:2008bg . For the class of two-loop corrections of
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ the counterterm can be written as

$\displaystyle\delta^{(2)}t_{\beta}^{2}$
$\displaystyle=t_{\beta}^{2}\left(\delta^{(2)}Z_{\mathcal{H}_{2}}-\delta^{(2)}Z_{\mathcal{H}_{1}}\right).$
(30)

  * •

When neglecting momentum-dependent contributions and taking the gaugeless
limit, the purely UV-divergent two-loop counterterms
$\delta^{(2)}Z_{\mathcal{H}_{1}}$, $\delta^{(2)}Z_{\mathcal{H}_{2}}$ and
$\delta^{(2)}t_{\beta}$ cancel each other and are therefore not required for
renormalization, compare Ref. Passehr:2017ufr . If one of these two
limitations is dropped, $\delta^{(2)}Z_{\mathcal{H}_{1}}$,
$\delta^{(2)}Z_{\mathcal{H}_{2}}$ and $\delta^{(2)}t_{\beta}$ are necessary in
order to obtain a UV-finite result. In the corrections discussed in this
article these counterterms have to be taken into account as none of these
approximations is used.

It should also be noted that the chosen renormalization conditions for
$\delta^{(2)}Z_{\mathcal{H}_{2}}$ and $\delta^{(2)}t_{\beta}$ are not equal to
pure $\overline{\text{DR}}$ conditions, since the top-mass counterterm
$\delta^{(1)}m_{t}$ which enters in $\delta^{(2)}Z_{\mathcal{H}_{2}}$ is fixed
by an on-shell condition. The resulting differences between the two schemes
have been discussed in Degrassi:2014pfa ; Borowka:2015ura .

  * •

Renormalization of the $D$ terms in the Higgs–squark couplings which are
induced by the gauge coupling $g_{2}$, as well as the relation between the
charged and $CP$-odd Higgs masses require counterterms for the $Z$\- and
$W$-boson masses, $\delta^{(2)}M_{Z}^{2}$ and $\delta^{(2)}M_{W}^{2}$,
respectively. We treat $M_{W}$ and $M_{Z}$ as independent input parameters and
fix their renormalization constants by the on-shell conditions

$\displaystyle\Re\mathfrak{e}{\left[\hat{\Sigma}_{Z,W}^{(2)}(M_{Z,W}^{2})\right]}$
$\displaystyle=0\,,$ (31)

leading to

$\displaystyle\delta^{(2)}M_{Z,W}^{2}$
$\displaystyle=\Re\mathfrak{e}{\left[\Sigma_{Z,W}^{(2)}{\left(M_{Z,W}^{2}\right)}\right]}.$
(32)

Here $\Sigma_{Z,W}^{(2)}$ denote the transverse parts of the two-loop self-
energies of $Z$ and $W$, repectively.

Most of the Feynman diagrams contributing to the two-loop self-energies
$\Sigma_{Z,W}^{(2)}$ differ from the Higgs self-energies depicted in Figs. 1
and 4 only by the external fields. Pictorially, they can be otained by
replacing the neutral external Higgs fields by the $Z$, and the charged Higgs
field by the $W$ boson field. All additional topologies are depicted in Fig.
4.

###  3.2 Sub-loop renormalization

Apart from the genuine two-loop diagrams, the lowest-order QCD contributions
to the self-energies and tadpoles involve one-loop diagrams with insertions of
one-loop counterterms. This subrenormalization concerns masses and mixing of
the colored particles.

The required one-loop counterterms for subrenormalization arise from the quark
$q$ and scalar quark $\tilde{q}$ sectors. The squark mass matrices in the
$\big{(}\tilde{q}_{\text{L}},\,\tilde{q}_{\text{R}}\big{)}$ bases are given in
lowest order by

$\displaystyle\mathbf{M}_{\tilde{q}}$
$\displaystyle=\begin{pmatrix}m_{\tilde{q}_{\text{L}}}^{2}+m_{q}^{2}+M_{Z}^{2}\,c_{2\beta}(T_{q}^{3}-Q_{q}\,s^{2}_{\mathrm{w}})&m_{q}\left(A_{q}^{*}-\mu\,\kappa_{q}\right)\\\\[2.84544pt]
m_{q}\left(A_{q}-\mu^{*}\,\kappa_{q}\right)&m_{\tilde{q}_{\text{R}}}^{2}+m_{q}^{2}+M_{Z}^{2}\,c_{2\beta}\,Q_{q}\,s^{2}_{\mathrm{w}}\end{pmatrix},$
$\displaystyle\begin{split}\kappa_{t,c,u}&=\frac{1}{t_{\beta}}\,,\\\
\kappa_{b,s,d}&=t_{\beta}\,,\end{split}$ (33)

with $Q_{q}$ and $T^{3}_{q}$ denoting charge and isospin of
$q\in\\{u,c,t,d,s,b\\}$. For the sake of convenience we suppress repeating the
indices of the first and second generation in the following since
renormalization is analogous to the third generation. $SU(2)$-invariance
requires $m_{\tilde{t}_{\text{L}}}^{2}=m_{\tilde{b}_{\text{L}}}^{2}\equiv
m_{\tilde{Q}_{3}}^{2}$.

The squark mass eigenvalues can be obtained from unitary transformations,

$\displaystyle\mathbf{U}^{\tilde{q}}\mathbf{M}_{\tilde{q}}\mathbf{U}^{\tilde{q}\dagger}$
$\displaystyle=\mathrm{diag}{\left(m_{\tilde{q}_{1}}^{2},\,m_{\tilde{q}_{2}}^{2}\right)}.$
(34)

Since $A_{q}$ and $\mu$ are complex parameters, the unitary matrices
$\mathbf{U}^{\tilde{q}}$ can be described by the mixing angle
$\theta_{\tilde{q}}$ and an additional phase $\varphi_{\tilde{q}}$.

The independent parameters which enter the two-loop calculation through the
quark–squark sector are: the quark masses $m_{q}$, the soft SUSY-breaking
parameters $m_{\tilde{Q}_{i}}$ and $m_{\tilde{q}_{\text{R}}}$, and the complex
trilinear couplings $A_{q}=\lvert A_{q}\rvert\,\mathrm{e}^{i\phi_{A_{q}}}$.
These parameters have to be renormalized at the one-loop level,

$\displaystyle m_{q}$ $\displaystyle\rightarrow m_{q}+\delta^{(1)}m_{q}\,,$
$\displaystyle m_{\tilde{q}_{\text{L,R}}}^{2}$ $\displaystyle\rightarrow
m_{\tilde{q}_{\text{L,R}}}^{2}+\delta^{(1)}m_{\tilde{q}_{\text{L,R}}}^{2}\,,$
$\displaystyle A_{q}$ $\displaystyle\rightarrow A_{q}+\delta^{(1)}A_{q}\,,$
(35)

thus defining transformations
$\mathbf{M}_{\tilde{q}}\to\mathbf{M}_{\tilde{q}}+\delta^{(1)}\mathbf{M}_{\tilde{q}}$
for the mass matrices in Eq. (33). The other free parameter $\mu$, which is
related to the Higgsino sector, enters the self-energies as well. However, the
renormalization of $\mu$ does not receive one-loop corrections of
$\mathcal{O}{\left(\alpha_{s}\right)}$ and is therefore not part of the
contributions considered in this calculation.

The individual renormalization conditions for the colored sector are
formulated as follows:

  * •

Renormalization of the top quark mass is carried out in the on-shell scheme,
i. e.

$\displaystyle\delta^{(1)}m_{t}$
$\displaystyle=m_{t}\,\Re\mathfrak{e}{\left[\frac{1}{2}\left(\Sigma_{t}^{\text{L}(1)}{\left(m_{t}^{2}\right)}+\Sigma_{t}^{\text{R}(1)}{\left(m_{t}^{2}\right)}\right)+\Sigma_{t}^{\text{S}(1)}{\left(m_{t}^{2}\right)}\right]},$
(36)

where the quark self-energy is given in terms of its Lorentz decomposition

$\displaystyle\Sigma_{q}(p)$
$\displaystyle=\,\not{\\!p}\,\omega_{-}\,\Sigma_{q}^{\mathrm{L}}(p^{2})+\not{\\!p}\,\omega_{+}\,\Sigma_{q}^{\mathrm{R}}(p^{2})+m_{q}\,\Sigma_{q}^{\mathrm{S}}(p^{2})+m_{q}\gamma_{5}\,\Sigma_{q}^{\mathrm{PS}}(p^{2})$
(37)

with the left-, right-handed projectors
$\omega_{-,+}=\tfrac{1}{2}\left(1\mp\gamma_{5}\right)$.

The bottom mass is renormalized in the $\overline{\text{DR}}$ scheme (see
Refs. Dedes:2003km ; Heinemeyer:2004xw ; Heinemeyer:2010mm ) at the scale
$m_{t}^{\text{os}}$. The counterterm can be obtained by using the expression
in analogy to the counterterm for the top quark mass in Eq. (36) and
restricting to the UV-divergent contributions at the scale
$m_{t}^{\text{os}}$. The choice of a $\overline{\text{DR}}$ renormalization
for $m_{b}$ is convenient in order to incorporate a resummation of
$\tan\beta$-enhanced contributions to the relation between the bottom quark
mass and the bottom Yukawa coupling, see Sec. 3.3 below. The contributing
Feynman diagrams for the renormalization of $m_{t}$ and $m_{b}$ are depicted
in Fig. 5.

  * •

In order to fix the renormalization constants of the stop sector, we employ
the relation

$\displaystyle\delta^{(1)}\mathbf{M}_{\tilde{t}}$
$\displaystyle=\delta^{(1)}{\left(\mathbf{U}^{\tilde{t}\dagger}\,\mathrm{diag}{\left(m_{\tilde{q}_{1}}^{2},\,m_{\tilde{q}_{2}}^{2}\right)}\,\mathbf{U}^{\tilde{t}}\right)}=\mathbf{U}^{\tilde{t}\dagger}\begin{pmatrix}\delta^{(1)}m_{\tilde{t}_{1}}^{2}&\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}\\\
\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2\,*}&\delta^{(1)}m_{\tilde{t}_{2}}^{2}\end{pmatrix}\mathbf{U}^{\tilde{t}}\,.$
(38)

Thus we derive

$\displaystyle\delta^{(1)}m_{\tilde{t}_{\text{L}}}^{2}$
$\displaystyle=\sum\limits_{i=1}^{2}\lvert\mathbf{U}^{\tilde{t}}_{i1}\rvert^{2}\,\delta^{(1)}m_{\tilde{t}_{i}}^{2}+2\,\Re\mathfrak{e}{\left[\mathbf{U}^{\tilde{t}}_{21}\mathbf{U}^{\tilde{t}*}_{11}\,\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}\right]}-2\,m_{t}\,\delta^{(1)}m_{t}\,,$
(39a) $\displaystyle\delta^{(1)}m_{\tilde{t}_{\text{R}}}^{2}$
$\displaystyle=\sum\limits_{i=1}^{2}\lvert\mathbf{U}^{\tilde{t}}_{i2}\rvert^{2}\,\delta^{(1)}m_{\tilde{t}_{i}}^{2}+2\,\Re\mathfrak{e}{\left[\mathbf{U}^{\tilde{t}}_{22}\mathbf{U}^{\tilde{t}*}_{12}\,\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}\right]}-2\,m_{t}\,\delta^{(1)}m_{t}\,,$
(39b)
$\displaystyle\begin{split}\delta^{(1)}A_{t}&=\mathbf{U}^{\tilde{t}}_{11}\mathbf{U}^{\tilde{t}*}_{12}\frac{\delta^{(1)}m_{\tilde{t}_{1}}^{2}-\delta^{(1)}m_{\tilde{t}_{2}}^{2}}{m_{t}}+\mathbf{U}^{\tilde{t}}_{21}\mathbf{U}^{\tilde{t}*}_{12}\frac{\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}}{m_{t}}+\mathbf{U}^{\tilde{t}}_{22}\mathbf{U}^{\tilde{t}*}_{11}\frac{\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2\,*}}{m_{t}}\\\
&\quad-\left(A_{t}-\frac{\mu^{*}}{t_{\beta}}\right)\frac{\delta^{(1)}m_{t}}{m_{t}}\,.\end{split}$
(39c)

The counterterm $\delta^{(1)}A_{t}$ given in Eq. (39c) provides the
renormalization of the complex parameter $A_{t}$. It should be noted that the
counterterm in fact only contributes to the absolute value of $A_{t}$, while
the phase of $A_{t}$ remains unrenormalized.

The counterterms $\delta^{(1)}m_{\tilde{t}_{1}}^{2}$ and
$\delta^{(1)}m_{\tilde{t}_{2}}^{2}$ are fixed by on-shell conditions for the
top-squarks,

$\displaystyle\delta^{(1)}m_{\tilde{t}_{i}}^{2}$
$\displaystyle=\Re\mathfrak{e}{\left[\,\Sigma_{\tilde{t}_{ii}}^{(1)}{\left(m_{\tilde{t}_{i}}^{2}\right)}\right]},\quad
i=1,2\,,$ (40)

involving the diagonal $\tilde{t}_{1,2}$ self-energies, see Fig. 5. The
remaining counterterm $\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}$ is
fixed by the renormalization condition (see Ref. Heinemeyer:2007aq )

$\displaystyle\delta^{(1)}m_{\tilde{t}_{1}\tilde{t}_{2}}^{2}$
$\displaystyle=\frac{1}{2}\,\Re\mathfrak{e}{\left[\Sigma_{\tilde{t}_{12}}^{(1)}{\left(m_{\tilde{t}_{1}}^{2}\right)}+\Sigma_{\tilde{t}_{12}}^{(1)}{\left(m_{\tilde{t}_{2}}^{2}\right)}\right]}\
,$ (41)

which involves the non-diagonal squark self-energy shown in Fig. 5 with
incoming $\tilde{t}_{2}$ and outgoing $\tilde{t}_{1}$.

Figure 5:  Types of Feynman diagrams for the renormalization of the
quark–squark sector. $\;\tilde{q}\neq\tilde{o}$. Topology $5$ contains a one-
point loop with a mass-less gluon and is therefore equal to zero.

  * •

Between the gauge and mass eigenstates of the bottom squarks we employ an
analogous relation to Eq. (38). We derive

$\displaystyle\delta^{(1)}m_{\tilde{b}_{\text{L}}}^{2}$
$\displaystyle=\sum\limits_{i=1}^{2}\lvert\mathbf{U}^{\tilde{b}}_{i1}\rvert^{2}\,\delta^{(1)}m_{\tilde{b}_{i}}^{2}+2\,\Re\mathfrak{e}{\left[\mathbf{U}^{\tilde{b}}_{21}\mathbf{U}^{\tilde{b}*}_{11}\,\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}\right]}-2\,m_{b}\,\delta^{(1)}m_{b}\,,$
(42a) $\displaystyle\delta^{(1)}m_{\tilde{b}_{\text{R}}}^{2}$
$\displaystyle=\sum\limits_{i=1}^{2}\lvert\mathbf{U}^{\tilde{b}}_{i2}\rvert^{2}\,\delta^{(1)}m_{\tilde{b}_{i}}^{2}+2\,\Re\mathfrak{e}{\left[\mathbf{U}^{\tilde{b}}_{22}\mathbf{U}^{\tilde{b}*}_{12}\,\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}\right]}-2\,m_{b}\,\delta^{(1)}m_{b}\,,$
(42b)
$\displaystyle\begin{split}\delta^{(1)}A_{b}&=\Bigg{[}\mathbf{U}^{\tilde{b}}_{11}\mathbf{U}^{\tilde{b}*}_{12}\frac{\delta^{(1)}m_{\tilde{b}_{1}}^{2}-\delta^{(1)}m_{\tilde{b}_{2}}^{2}}{m_{b}}+\mathbf{U}^{\tilde{b}}_{21}\mathbf{U}^{\tilde{b}*}_{12}\frac{\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}}{m_{b}}+\mathbf{U}^{\tilde{b}}_{22}\mathbf{U}^{\tilde{b}*}_{11}\frac{\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2\,*}}{m_{b}}\\\
&\quad\hphantom{\Bigg{[}}-\left(A_{b}-\mu^{*}\,t_{\beta}\right)\frac{\delta^{(1)}m_{b}}{m_{b}}\Bigg{]}_{\overline{\text{DR}}}\,.\end{split}$
(42c)

As indicated by the subscript, we choose to renormalize $A_{b}$ in the
$\overline{\text{DR}}$ scheme, which has been shown to be convenient for
reasons of numerical stability Dedes:2003km ; Heinemeyer:2004xw ;
Heinemeyer:2010mm . The scale of $A_{b}$ is chosen to be $m_{t}^{\text{os}}$.

As a consequence of $SU(2)$ invariance the counterterm
$\delta^{(1)}m_{\tilde{b}_{\text{L}}}^{2}$ is not independent, but a derived
quantity which is fixed by the renormalization of the top–stop sector in Eq.
(39a), since

$\displaystyle\delta^{(1)}m_{\tilde{b}_{\text{L}}}^{2}=\delta^{(1)}m_{\tilde{Q}_{3}}^{2}=\delta^{(1)}m_{\tilde{t}_{\text{L}}}^{2}\,.$
(43)

Inserting Eq. (42a) and solving for $\delta^{(1)}m_{\tilde{b}_{1}}^{2}$ yields

$\displaystyle\delta^{(1)}m_{\tilde{b}_{1}}^{2}$
$\displaystyle=\frac{1}{\lvert\mathbf{U}^{\tilde{b}}_{11}\rvert^{2}}\left(\delta^{(1)}m_{\tilde{t}_{\text{L}}}^{2}-\lvert\mathbf{U}^{\tilde{b}}_{12}\rvert^{2}\,\delta^{(1)}m_{\tilde{b}_{2}}^{2}-2\,\Re\mathfrak{e}{\left[\mathbf{U}^{\tilde{b}}_{21}\mathbf{U}^{\tilde{b}*}_{11}\,\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}\right]}-2\,m_{b}\,\delta^{(1)}m_{b}\right).$
(44)

The other two counterterms $\delta^{(1)}m_{\tilde{b}_{2}}^{2}$ and
$\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}$ are fixed analogously as for
the stops:

$\displaystyle\delta^{(1)}m_{\tilde{b}_{2}}^{2}$
$\displaystyle=\Re\mathfrak{e}{\left[\,\Sigma_{\tilde{b}_{22}}^{(1)}{\left(m_{\tilde{b}_{2}}^{2}\right)}\right]}\,,$
(45a) $\displaystyle\delta^{(1)}m_{\tilde{b}_{1}\tilde{b}_{2}}^{2}$
$\displaystyle=\frac{1}{2}\,\Re\mathfrak{e}{\left[\Sigma_{\tilde{b}_{12}}^{(1)}{\left(m_{\tilde{b}_{1}}^{2}\right)}+\Sigma_{\tilde{b}_{12}}^{(1)}{\left(m_{\tilde{b}_{2}}^{2}\right)}\right]}\,.$
(45b)

Therefore in our scheme only $m_{\tilde{b}_{2}}$ is renormalized on-shell,
while the counterterm $\delta^{(1)}m_{\tilde{b}_{1}}^{2}$ is a derived
quantity according to Eq. (44).

###  3.3 Resummation of $\tan\beta$-enhanced terms

The Yukawa coupling of the bottom quark $h_{b}$ receives radiative corrections
proportional to $\tan\beta$. Those $\tan\beta$-enhanced contributions can be
resummed as described in Refs. Banks:1987iu ; Hall:1993gn ; Hempfling:1993kv ;
Carena:1994bv ; Carena:1999py ; Eberl:1999he ; Williams:2011bu . The resummed
contributions $\Delta_{b}$ are UV finite and generally yield complex numerical
results. For the numerical evaluation in Sec. 5, we use the version for
$\Delta_{b}$ at the one-loop order which is implemented in FeynHiggs and
outlined in the following. The largest $\tan\beta$-enhanced contributions can
be absorbed by using an effective bottom-quark mass, which is related to the
$\overline{\text{DR}}$-renormalized bottom quark mass in the MSSM as follows,

$\displaystyle
m_{b}^{\overline{\text{DR}},\text{MSSM}}{\left(m_{t}^{\text{os}}\right)}\simeq
m_{b,\text{eff}}$
$\displaystyle=\frac{m_{b}^{\overline{\text{DR}},\text{SM}}{\left(m_{t}^{\text{os}}\right)}}{\lvert
1+\Delta_{b}\rvert}\left(1-\delta_{b}\right),$ (46)

where $m_{b}^{\overline{\text{DR}},\text{SM}}(m_{t}^{\text{os}})$ is the
bottom mass in the $\overline{\text{DR}}$ renormalization scheme in the
Standard Model evaluated at the on-shell top mass. The $\tan\beta$-enhanced
contributions are captured in $\Delta_{b}$ and properly resummed by including
them in the denominator. The remaining parts of the scalar part of the
$\overline{\text{DR}}$-renormalized bottom self-energy
$\hat{\Sigma}_{b}^{\text{S}}$ which are not enhanced by $\tan\beta$ are
contained in $\delta_{b}$ such that

$\displaystyle\hat{\Sigma}_{b}^{\text{S}}{\left(0\right)}$
$\displaystyle=-\Delta_{b}-\delta_{b}\,.$ (47)

The expression $\Delta_{b}$ at the one-loop order contains contributions from
gluinos, charginos and neutralinos (ordered in decreasing numerical
importance) and reads

$\displaystyle\begin{split}\Delta_{b}&=\frac{2\,\alpha_{s}{\left(Q\right)}}{3\pi}\frac{M_{3}^{*}}{m_{b}}\sum\limits_{i=1}^{2}\mathbf{U}^{\tilde{b}}_{i1}\mathbf{U}^{\tilde{b}*}_{i2}\,B_{0}{\left(0,\lvert
M_{3}\rvert^{2},m_{\tilde{b}_{i}}^{2}\right)}\\\
&\quad+\frac{\alpha{\left(Q\right)}}{4\pi}\sum\limits_{g=1}^{3}\sum\limits_{i,j=1}^{2}\frac{m_{\tilde{\chi}^{\pm}_{i}}}{m_{b}}c_{\text{L}}c_{\text{R}}\lvert\mathbf{C}_{g3}\rvert^{2}\,B_{0}{\left(0,m_{\tilde{\chi}^{\pm}_{i}}^{2},m_{\tilde{u}^{g}_{j}}^{2}\right)}\\\
&\quad-\frac{\alpha{\left(Q\right)}}{8\pi}\sum\limits_{i=1}^{4}\sum\limits_{j=1}^{2}\frac{m_{\tilde{\chi}^{0}_{i}}}{m_{b}}n_{\text{L}}n_{\text{R}}B_{0}{\left(0,m_{\tilde{\chi}^{0}_{i}}^{2},m_{\tilde{b}_{j}}^{2}\right)}\
.\end{split}$ (48)

The couplings $\alpha_{s}$ and $\alpha$ are running parameters and are
evaluated at the scale
$Q=\sqrt{\vphantom{Q}m_{\tilde{b}_{1}}m_{\tilde{b}_{2}}}$. The symbol
$\mathbf{C}$ depicts the CKM matrix, and $u^{g},\,\tilde{u}^{g}$ are the $g$th
generation up-type quarks and squarks, whereas
$B_{0}{\left(0,m_{1},m_{2}\right)}$ and $B_{1}{\left(0,m_{1},m_{2}\right)}$
are one-loop functions.555These one-loop functions are given by
$\displaystyle B_{1}{\left(0,m_{1},m_{2}\right)}$
$\displaystyle=-\frac{1}{2}B_{0}{\left(0,m_{1},m_{2}\right)}+\frac{m_{2}-m_{1}}{2}D\\!B_{0}{\left(0,m_{1},m_{2}\right)}\
,$ $\displaystyle B_{0}{\left(0,m_{1},m_{2}\right)}$
$\displaystyle=\frac{A_{0}{\left(m_{1}\right)}-A_{0}{\left(m_{2}\right)}}{m_{1}-m_{2}}\
,$ $\displaystyle D\\!B_{0}{\left(0,m_{1},m_{2}\right)}$
$\displaystyle=\frac{m_{1}^{2}-m_{2}^{2}+2\,m_{1}\,m_{2}\log{\tfrac{m_{2}}{m_{1}}}}{2\,(m_{1}-m_{2})^{3}}\
,$ $\displaystyle A_{0}{\left(m\right)}$
$\displaystyle=m\left(1-\log{\tfrac{m}{\mu_{r}}}\right)\ .$ As mentioned
above, we otherwise neglect CKM mixing in the two-loop contributions that we
evaluate. The renormalization scale $\mu_{r}$ from the loop integrals drops
out. The coefficients $c_{\text{L,R}}$ and $n_{\text{L,R}}$ are given by

$\displaystyle c_{\text{L}}$
$\displaystyle=\frac{\mathbf{V}^{\tilde{\chi}*}_{i1}\,\mathbf{U}^{\tilde{u}^{g}}_{j1}}{s_{\text{w}}}-\frac{m_{u^{g}}\,\mathbf{V}^{\tilde{\chi}*}_{i2}\,\mathbf{U}^{\tilde{u}^{g}}_{j2}}{\sqrt{2}\,M_{W}\,s_{\beta}\,s_{\text{w}}}\,,$
$\displaystyle c_{\text{R}}$
$\displaystyle=\frac{m_{b}\,\mathbf{U}^{\tilde{\chi}*}_{i2}\,\mathbf{U}^{\tilde{u}^{g}*}_{j1}}{\sqrt{2}\,M_{W}\,c_{\beta}\,s_{\text{w}}}\,,$
(49a) $\displaystyle n_{\text{L}}$
$\displaystyle=\left(\frac{\mathbf{N}^{\tilde{\chi}*}_{i1}}{3\,c_{\text{w}}}-\frac{\mathbf{N}^{\tilde{\chi}*}_{i2}}{s_{\text{w}}}\right)\mathbf{U}^{\tilde{b}}_{j1}+\frac{m_{b}\,\mathbf{N}^{\tilde{\chi}*}_{i3}\,\mathbf{U}^{\tilde{b}}_{j2}}{M_{W}\,c_{\beta}\,s_{\text{w}}}\,,$
$\displaystyle n_{\text{R}}$
$\displaystyle=\frac{2\,\mathbf{N}^{\tilde{\chi}*}_{i1}\,\mathbf{U}^{\tilde{b}*}_{j2}}{3\,c_{\text{w}}}+\frac{m_{b}\,\mathbf{N}^{\tilde{\chi}*}_{i3}\,\mathbf{U}^{\tilde{b}*}_{j1}}{M_{W}\,c_{\beta}\,s_{\text{w}}}\,.$
(49b)

In order to obtain a full conversion of the bottom mass between the on-shell
scheme and the $\overline{\text{DR}}$ scheme in Eq. (46), those parts of the
bottom self-energy which are not enhanced by $\tan\beta$ are included in
$\delta_{b}$ and incorporated in the numerator of Eq. (46).

At the one-loop order they read

$\displaystyle\begin{split}\delta
b&=\frac{\alpha_{s}{\left(Q\right)}}{3\pi}\sum\limits_{i=1}^{2}B_{1}{\left(0,\lvert
M_{3}\rvert^{2},m_{\tilde{b}_{i}}^{2}\right)}+\frac{\alpha{\left(Q\right)}}{8\pi}\sum\limits_{g=1}^{3}\sum\limits_{i,j=1}^{2}\left[\left|c_{\text{L}}\right|^{2}+\left|c_{\text{R}}\right|^{2}\right]\lvert\mathbf{C}_{g3}\rvert^{2}\,B_{1}{\left(0,m_{\tilde{\chi}^{\pm}_{i}}^{2},m_{\tilde{u}^{g}_{j}}^{2}\right)}\\\\[-2.15277pt]
&\quad+\frac{\alpha{\left(Q\right)}}{16\pi}\sum\limits_{i=1}^{4}\sum\limits_{j=1}^{2}\left[\left|n_{\text{L}}\right|^{2}+\left|n_{\text{R}}\right|^{2}\right]B_{1}{\left(0,m_{\tilde{\chi}^{0}_{i}}^{2},m_{\tilde{b}_{j}}^{2}\right)}\
.\end{split}$ (50)

The parameters entering in $\Delta_{b}$ and $\delta_{b}$ are computed in the
limit of large $\tan\beta$. The chargino and neutralino masses and mixing
matrices are then obtained as

$\displaystyle\lim_{t_{\beta}\rightarrow\infty}\text{diag}{\left(m_{\tilde{\chi}^{\pm}_{1}},m_{\tilde{\chi}^{\pm}_{2}}\right)}$
$\displaystyle=\mathbf{U}^{\tilde{\chi}*}\mathbf{X}\mathbf{V}^{\tilde{\chi}\dagger}\,,$
$\displaystyle\lim_{t_{\beta}\rightarrow\infty}\text{diag}{\left(m_{\tilde{\chi}^{0}_{1}},m_{\tilde{\chi}^{0}_{2}},m_{\tilde{\chi}^{0}_{3}},m_{\tilde{\chi}^{0}_{4}}\right)}$
$\displaystyle=\mathbf{N}^{\tilde{\chi}*}\mathbf{Y}\mathbf{N}^{\tilde{\chi}\dagger}\,,$
(51)

where we use

$\displaystyle\mathbf{X}$
$\displaystyle=\lim_{t_{\beta}\rightarrow\infty}\mathbf{M}_{\chi^{\pm}}=\begin{pmatrix}M_{2}&\sqrt{2}M_{W}s_{\beta}\\\
0&\mu\end{pmatrix},$ $\displaystyle\mathbf{Y}$
$\displaystyle=\lim_{t_{\beta}\rightarrow\infty}\mathbf{M}_{\chi^{0}}=\begin{pmatrix}M_{1}&0&0&M_{Z}s_{\beta}s_{\text{w}}\\\
0&M_{2}&0&-M_{Z}s_{\beta}c_{\text{w}}\\\ 0&0&0&-\mu\\\
M_{Z}s_{\beta}s_{\text{w}}&-M_{Z}s_{\beta}c_{\text{w}}&-\mu&0\end{pmatrix}.$
(52)

Thereby the matrices $\mathbf{U}^{\tilde{\chi}}$ and
$\mathbf{V}^{\tilde{\chi}}$ yield a singular value decomposition for
$\mathbf{X}$, and the matrix $\mathbf{N}^{\tilde{\chi}}$ yields Takagi’s
factorization Takagi:1927  on $\mathbf{Y}$.

The sbottom masses in this limit are computed from the matrix given in Eq.
(33) at $A_{b}=0$. Since the bottom mass itself also enters that matrix, the
final solution for $m_{b,\text{eff}}$ is found iteratively.

By using Eq. (46) for the bottom mass in the one-loop contributions to the
Higgs masses, the leading higher-order corrections to the Higgs masses from
the bottom–sbottom sector are generated. The contributions of the
bottom–sbottom sector to the two-loop self-energies presented in this article
add further subleading shifts. It should be noted that the expression given in
Eq. (46), which employs the $\overline{\text{DR}}$ scheme in the MSSM, is
chosen such that no double counting of the terms contained in
$m_{b,\text{eff}}$ occurs at the two-loop level.

##  4 Numerical evaluation of the self-energies

The renormalized two-loop self-energies are expressed in terms of two-loop
two-point multi-scale integrals with non-zero external momenta. With the help
of TwoCalc Weiglein:1993hd  and Reduze vonManteuffel:2012np  all integrals can
be reduced to the four irreducible scalar two-loop topologies depicted in Fig.
6, and products of analytically well-known one-loop one- and two-point
functions.

The scalar two-loop integrals are defined as

$T_{i_{1}i_{2}\dots
i_{n}}(p^{2},m_{i_{1}}^{2},m_{i_{2}}^{2},\dots,m_{i_{5}}^{2})=\left(2\pi\mu_{r}\right)^{2(4-D)}\\\
\times\iint\frac{\text{d}^{D}q_{1}}{i\pi^{2}}\,\frac{\text{d}^{D}q_{2}}{i\pi^{2}}\frac{1}{(k_{i_{1}}^{2}-m_{i_{1}}^{2}+i\delta)(k_{i_{2}}^{2}-m_{i_{2}}^{2}+i\delta)\cdots(k_{i_{n}}^{2}-m_{i_{n}}^{2}+{i\mkern
1.0mu}\delta)}\text{ ,}$ (53)

where $p$ is the external momentum, $q_{i}$ are the loop momenta, $m_{i}$ the
masses of the propagators, $\mu_{r}$ is the renormalization scale and
$D=4-2\varepsilon$ the dimension. The ${i\mkern 1.0mu}\delta$ results from the
solutions of the field equations in terms of causal Green functions, while the
indices ${i_{1},i_{2},\dots i_{n}}$ label which $k_{i}$ and $m_{i}$ appear in
the propagators of the integral. The five different $k_{i}$ read

$\displaystyle k_{1}=q_{1},\hskip 10.0ptk_{2}=q_{1}+p,\hskip
10.0ptk_{3}=q_{2}-q_{1},\hskip 10.0ptk_{4}=q_{2},\hskip
10.0ptk_{5}=q_{2}+p\text{.}$ (54)

The irreducible two-loop integrals of Fig. 6 may depend on up to five
different internal mass scales taken from the following set,

$\displaystyle
m_{t},\,m_{b},\,m_{\tilde{t}_{1}},\,m_{\tilde{t}_{2}},\,m_{\tilde{b}_{1}},\,m_{\tilde{b}_{2}},\,m_{\tilde{g}}=|M_{3}|,$
(55)

in addition to a non-zero external momentum, taking the values of
$p^{2}=M_{h_{1}}^{2},\,M_{h_{2}}^{2},\,M_{h_{3}}^{2}$ when entering the
unrenormalized self-energies, or
$p^{2}=m_{H^{\pm}}^{2},\,m_{W}^{2},\,m_{Z}^{2}$ when entering the self-
energies through two-loop renormalization constants. Recently, a lot of
progress has been made towards describing and evaluating integrals of this
class analytically Bloch:2013tra ; Bloch:2014qca ; Adams:2013kgc ;
Remiddi:2013joa ; Adams:2015gva ; Adams:2015ydq ; Bloch:2016izu ;
Remiddi:2016gno ; Adams:2016xah ; Broedel:2017siw . However, to the best of
our knowledge, an implementation of the analytical results for all topologies
in Fig. 6 is not publicly available. We have therefore calculated these
integrals numerically using the program SecDec Carter:2010hi ; Borowka:2012yc
; Borowka:2013cma .

Figure 6: Irreducible two-loop topologies resulting from tensor reduction,
calculated numerically with the program SecDec. Some of the internal lines may
also be massless.

For the evaluation, the resulting new contributions to the neutral Higgs-boson
self-energies have been added to FeynHiggs via its interface to external
programs, see section 2.4 of Ref. Borowka:2014wla  for details. We have
extended the existing interface to the program SecDec in FeynHiggs to deal
with the 177 mass configurations of which 88 are computed at four different
kinematic points, 72 at two and the rest at one kinematic point. The
parameters entering the integrals are evaluated by FeynHiggs and passed on to
SecDec. It should be noted that the heavy growth of mass configurations with
respect to non-electroweak corrections is due to an increase in the number of
mass scales involved in the renormalized self-energies.

We constructed two independent integration setups to allow for consistency
checks of the numerical result. The two-point one-loop topologies entering the
self-energies up to $\mathcal{O}(\varepsilon)$ are known analytically. The
bulk of their implementation was previously tested in Ref. Borowka:2014wla
and compared with the authors of Ref. Degrassi:2014pfa . Additional mass
configurations were newly implemented and checked against SecDec. The increase
in two-loop mass configurations by more than a factor five with respect to the
previous setup in Ref. Borowka:2014wla  calls for a higher precision of the
integrals to avoid numerical instabilites due to cancellations. With the
integral reduction, unphysical thresholds can be introduced which cancel in
the sum of all contributing diagrams. Numerically, due to round-off errors,
the cancellation might however not always be exact, leading to numerical
instabilities. The latter are cured by introducing a small imaginary part to
the denominators of the coefficients arising from the integral reduction. We
have verified that the numerical dependence of the self-energies on this
technical regularization parameter is negligible.

The fact that we take a non-zero value of the bottom quark mass into account
leads to a large hierarchy among the different mass scales. Numerical
convergence at the desired accuracy is therefore difficult to accomplish. On
the other hand, we have analyzed the influence of the quark masses of the
first and second generation on the two-loop integrals in the self-energies.
For the second generation and $\tan{\beta}\gg 1$ a negative shift in the
Higgs-boson mass correction of only about $20$ MeV can be observed when
neglecting the light quark masses. The effect is even smaller for the quark
masses of the first generation. The terms which involve the light quark masses
in couplings are negligible, too. It is due to this reason that we will assume
the first and second generation quarks to be massless throughout the rest of
our numerical analysis. The numerical impact of the gauge contributions of the
light quarks will be discussed below.

In order to achieve a relative precision of at least $10^{-7}$ for each
integral, we use the deterministic integrator Cuhre included in the Cuba
library Agrawal:2011tm ; Hahn:2004fe  but have optimized the integration
parameters for each integral topology and mass configuration individually.

As a further crosscheck of our computation, we have compared the
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ contribution by the top–stop and
bottom–sbottom particles to the $Z$-boson self-energy which is required for
renormalization of the Higgs sector. Since Degrassi:2014pfa  uses massless
bottom quarks, we have reevaluated our result for the $Z$ self-energy in the
limit $m_{b}=0$. In order to avoid a dependence on the renormalization scheme
of the quark–squark sectors, the $Z$-boson self-energy has been evaluated in
the $\overline{\text{DR}}$ scheme by both groups for this comparison. Overall
we have found a very good agreement with discrepancies at the level of $0.3$
GeV2.

We find an overall uncertainty of the self-energies entering the light Higgs-
boson mass of maximally $0.2\%$ by adding all uncertainties on the numerical
evaluation of the two-loop integrals in quadrature. Given the resulting size
of our newly computed corrections analyzed in the next section, the absolute
uncertainty on the light Higgs boson mass is maximally $0.4$ MeV, which is
well below the shift coming from neglecting light quark masses from the first
and second generation.

The total of $513$ integrals have been computed numerically on the fly before
passing the resulting two-loop self-energies back to FeynHiggs, where they are
added to the corresponding matrix elements just before the determination of
the propagator poles.

##  5 Numerical results for the Higgs mass spectrum

In the following we analyze the numerical impact of the newly computed
corrections. We start with a comparison with earlier results in the literature
and then discuss our results in three different scenarios: an
$m_{h}^{\text{mod}}$-like scenario (based on Ref. Carena:2013ytb ), a scenario
with a particularly large value of $\tan\beta$ where contributions from the
bottom and sbottom sector are enhanced, and a low-$m_{H}$ scenario (inspired
by Refs. Carena:2013ytb ; Bechtle:2016kui ). For better readability of the
results, we define three different Higgs-boson masses resulting from different
higher-order contributions

$\displaystyle M_{h_{i}}^{\text{old}},$ $\displaystyle\text{ contains:
}\mathcal{O}(\alpha_{t}\alpha_{s})\rvert_{p^{2}=0}\text{ with complex
parameters},$ $\displaystyle\tilde{M}_{h_{i}}^{\text{old}},$
$\displaystyle\text{ contains: }\text{same as
}M_{h_{i}}^{\text{old}}+\left.\mathcal{O}(\alpha_{b}\alpha_{s})\right|_{p^{2}=0}\text{
with real parameters},$ $\displaystyle M_{h_{i}}^{\text{new}},$
$\displaystyle\text{ contains:
}\mathcal{O}(\alpha_{q}\alpha_{s}),\mathcal{O}(\alpha\alpha_{s}),\mathcal{O}(h_{q}h_{o}\alpha_{s})\text{
with non-zero }p^{2},$ $\displaystyle\text{with }i\in\\{1,2,3\\},\quad
q,o\in\\{b,t\\}.$ (56)

All the above results contain the full one-loop and leading
$\mathcal{O}(\alpha_{t}^{2})\rvert_{p^{2}=0}$ two-loop contributions, and the
$\tan\beta$-enhanced contributions to the relation between the bottom quark
mass and the bottom Yukawa coupling are resummed, see section 3.3. As
mentioned earlier, the quark masses and Yukawa couplings of the first and
second family are neglected. Thus, the first and second generation contributes
only at $\mathcal{O}(\alpha\alpha_{s})$ by D-term contributions of the
sfermions. We focus our numerical discussion on the fixed-order result up to
the two-loop level, i. e. no combination with resummed higher-order
logarithmic contributions as discussed in Refs. Hahn:2013ria ; Bahl:2016brp ;
Bahl:2017aev  is employed.

Using the definitions of Eq. (56), we assign

$\displaystyle\Delta
M_{h_{i}}=M_{h_{i}}^{\text{new}}-M_{h_{i}}^{\text{old}},\quad\Delta\tilde{M}_{h_{i}}=M_{h_{i}}^{\text{new}}-\tilde{M}_{h_{i}}^{\text{old}}\,.$
(57)

The size of the effects of our newly computed contributions is contained in
$\Delta M_{h_{i}}$, since all the previously known terms are subtracted. So
far, the two-loop terms of $\mathcal{O}(\alpha_{b}\alpha_{s})$ were only known
in the MSSM with real parameters and $m_{A}$ as input parameter.
$\Delta\tilde{M}_{h_{i}}$ shows our new contributions without these terms, if
$m_{A}$ is chosen as input parameter.

Below we will discuss our results for non-zero phases of complex parameters.
We investigate in particular the variation of the phases
$\phi_{M_{3}},\,\phi_{A_{t}}$ and $\phi_{A_{b}}$, which are much less
constrained by experimental bounds on EDMs than the phases of $\mu$, $M_{1}$
(in the usual convention where the parameter $M_{2}$ is chosen to be real) and
the phases of the trilinear couplings of the first and second generation. As
discussed e. g. in Ref. Arbey:2014msa , scenarios with relatively large phase
values are possible. In order to demonstrate the possible impact of the phase
variations on the Higgs spectrum, below we display the phase dependences over
the whole range $\left(-\pi,\pi\right]$.

###  5.1 Comparison with earlier results

In a first step, in Tab. 2 we show a comparison of the results for the light
Higgs-boson mass including our new contributions with the results of Ref.
Degrassi:2014pfa , where in the MSSM with real parameters the corrections of
$\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ and the full corrections of
$\mathcal{O}(\alpha\alpha_{s})$ have been evaluated, and with the results up
to $\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ in the MSSM with real parameters
from Ref. Borowka:2014wla . The comparison is carried out for the benchmark
scenarios $m_{h}^{\text{max}}$, $m_{h}^{\text{mod}+}$, $m_{h}^{\text{mod}-}$
defined in Ref. Carena:2013ytb  and for a modified light-stop scenario used in
Ref. Bagnaschi:2014zla . We find overall good agreement with the results of
Ref. Degrassi:2014pfa . The comparison of the corrections of
$\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ with the full corrections of
$\mathcal{O}(\alpha\alpha_{s})$ shows that the inclusion of momentum
dependence in the $\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ corrections yields
a downward shift in $M_{h}$ which is to a large extent compensated by the
further corrections of $\mathcal{O}(\alpha\alpha_{s})$ for the scenarios that
are considered here. The corrections beyond those of
$\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ yield an upward shift in $M_{h}$ of
$520$ MeV in the $m_{h}^{\text{mod}+}$ and more than $1$ GeV in the
$m_{h}^{\text{mod}-}$ scenario compared to the results of Ref. Borowka:2014wla
. The size of the corrections shows a significant dependence on the parameters
in the stop sector. The corrections are largest in the $m_{h}^{\text{mod}-}$
scenario, where the stop masses are near the SUSY scale and $A_{t}$ is
negative. In this case there is a large compensation between the downward
shift caused by the corrections of $\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$
and the upward shift caused by the further corrections of
$\mathcal{O}(\alpha\alpha_{s})$. On the other hand, the corrections are
smallest for the modified light-stop scenario, in which case we find that the
contributions beyond the ones of $\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ from
Ref. Borowka:2014wla  even yield a small downward shift. The numerical
differences between the results for the contributions of
$\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ from Ref. Degrassi:2014pfa  and Ref.
Borowka:2014wla , which amount up to 0.3 GeV for the examples considered here,
result from different renormalization scheme choices of
$\delta^{(2)}Z_{\mathcal{H}_{i}}$, see the discussion in Refs. Borowka:2014wla
; Degrassi:2014pfa ; Borowka:2015ura . Those differences in the
renormalization schemes also affect the comparison between our results for
$M_{h}^{\text{new}}$ and the results for
$M_{h}^{\text{old}}+\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})+\mathcal{O}(\alpha\alpha_{s})$
from Ref. Degrassi:2014pfa  in Tab. 2.

scenario | $m_{h}^{\text{max}}$ | $m_{h}^{\text{mod}+}$ | $m_{h}^{\text{mod}-}$ | modified light-stop
---|---|---|---|---
$M_{h}^{\text{old}}$ (GeV) | $128.31$ | $125.36$ | $124.84$ | $122.68$
$M_{h}^{\text{old}}$ (GeV)Degrassi:2014pfa  | $128.32$ | $125.36$ | $124.84$ | $122.67$
$M_{h}^{\text{old}}+\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ (GeV)Borowka:2014wla  | $128.25$ | $125.23$ | $123.83$ | $122.64$
$M_{h}^{\text{old}}+\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})$ (GeV)Degrassi:2014pfa  | $127.94$ | $124.98$ | $123.96$ | $122.33$
$M_{h}^{\text{old}}+\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})+\mathcal{O}(\alpha\alpha_{s})$ (GeV)Degrassi:2014pfa  | $128.38$ | $125.63$ | $124.90$ | $122.46$
$M_{h}^{\text{new}}$ (GeV) | $128.53$ | $125.75$ | $124.85$ | $122.61$

Table 1:  Comparison of the results for the light Higgs-boson mass with Ref.
Degrassi:2014pfa  and Ref. Borowka:2014wla  for four benchmark scenarios from
Refs. Carena:2013ytb  and Bagnaschi:2014zla  with $m_{A}=500$ GeV and
$\tan\beta=20$.

| this publication | Ref. Degrassi:2014pfa
---|---|---
| $\delta^{(2)}Z_{\mathcal{H}_{i}}^{\text{Ref.~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Borowka:2014wla}{\@@citephrase{(}}{\@@citephrase{)}}}}}$ | $\delta^{(2)}Z_{\mathcal{H}_{i}}^{\text{Ref.~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Borowka:2014wla}{\@@citephrase{(}}{\@@citephrase{)}}}}}$ | $\delta^{(2)}Z_{\mathcal{H}_{i}}^{\text{Ref.~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Degrassi:2014pfa}{\@@citephrase{(}}{\@@citephrase{)}}}}}$
$m_{h}^{\text{mod}+}$-like |  $M_{\text{SUSY}}=2$ TeV
$M_{h}^{\text{old}}$ (GeV) | $129.38$ | $129.38$ | $129.38$
$M_{h}^{\text{new}}$ (GeV) | $129.92$ | $129.92$ | $129.84$
|  $M_{\text{SUSY}}=3$ TeV
$M_{h}^{\text{old}}$ (GeV) | $128.63$ | $128.63$ | $128.63$
$M_{h}^{\text{new}}$ (GeV) | $129.62$ | $129.61$ | $129.59$
$m_{h}^{\text{mod}-}$-like |  $M_{\text{SUSY}}=2$ TeV
$M_{h}^{\text{old}}$ (GeV) | $126.92$ | $126.92$ | $126.92$
$M_{h}^{\text{new}}$ (GeV) | $127.34$ | $127.33$ | $127.44$
|  $M_{\text{SUSY}}=3$ TeV
$M_{h}^{\text{old}}$ (GeV) | $127.02$ | $127.02$ | $127.02$
$M_{h}^{\text{new}}$ (GeV) | $127.80$ | $127.80$ | $127.94$

Table 2:  Values for the lightest Higgs-boson mass in the
$m_{h}^{\text{mod}+}$-like and $m_{h}^{\text{mod}-}$-like scenarios of Ref.
Carena:2013ytb  using $M_{\text{SUSY}}=2,3$ TeV and $m_{A}=500$ GeV,
$\tan\beta=20$. The results are compared with those provided by the authors of
Ref. Degrassi:2014pfa  for two different wave-function renormalization
schemes.

The differences in the renormalization schemes and the dependence on the
parameters in the stop sector are further investigated in Tab. 2. Here the
shifts in the light Higgs-boson mass are shown for SUSY scales of $2$ TeV and
$3$ TeV, using otherwise the parameters of the $m_{h}^{\text{mod}+}$ and
$m_{h}^{\text{mod}-}$ scenarios. The results for
$M_{h}^{\text{old}}+\mathcal{O}(\alpha_{t}\alpha_{s}p^{2})+\mathcal{O}(\alpha\alpha_{s})$
from Ref. Degrassi:2014pfa , where the mass and Yukawa coupling of the bottom
quark have been neglected, are labelled as $M_{h}^{\text{new}}$ in Tab. 2. Two
versions of the results from Ref. Degrassi:2014pfa  are shown, one using the
renormalization scheme adopted in Ref. Degrassi:2014pfa  with
$\delta^{(2)}Z_{\mathcal{H}_{i}}=\delta^{(2)}Z_{\mathcal{H}_{i}}^{\text{Ref.~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Degrassi:2014pfa}{\@@citephrase{(}}{\@@citephrase{)}}}}}$,
and the other using the renormalization scheme of Ref. Borowka:2014wla , which
we have adopted in the present work, with
$\delta^{(2)}Z_{\mathcal{H}_{i}}=\delta^{(2)}Z_{\mathcal{H}_{i}}^{\text{Ref.~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Borowka:2014wla}{\@@citephrase{(}}{\@@citephrase{)}}}}}$.666We
are very grateful to S. de Vita for providing us with those results. It can be
seen in Tab. 2 that there is very good agreement, at the level of about 10
MeV, between our results and the results from Ref. Degrassi:2014pfa  using the
renormalization scheme of Ref. Borowka:2014wla . The different choices of
renormalization schemes in the result of Ref. Degrassi:2014pfa  amount to mass
shifts of up to 150 MeV for the displayed examples. The difference between
$M_{h}^{\text{new}}$ and $M_{h}^{\text{old}}$ increases with $M_{\text{SUSY}}$
and reaches up to $1$ GeV for the $m_{h}^{\text{mod}+}$-like scenario at 3
TeV.

###  5.2 Scenario 1: $m_{h}^{\text{mod}}$-like

In the following we further investigate the numerical impact of our results,
including the effect of non-zero phases of the complex parameters. We start
with an $m_{h}^{\text{mod}}$-like scenario. The MSSM model parameters in this
scenario are chosen as follows

$\displaystyle m_{H^{\pm}}$ $\displaystyle=1.5\,\text{TeV},$ $\displaystyle
M_{2}$ $\displaystyle=500\,\text{GeV},$ $\displaystyle\lvert M_{3}\rvert$
$\displaystyle=2.5\,\text{TeV},$ $\displaystyle
m_{\\{\tilde{t},\tilde{b}\\}_{\text{L}}}=m_{\tilde{Q}_{3}}$
$\displaystyle=2.1\,\text{TeV},$ $\displaystyle
m_{\\{\tilde{t},\tilde{b}\\}_{\text{R}}}$ $\displaystyle=2\,\text{TeV},$
$\displaystyle\lvert X_{t}\rvert$
$\displaystyle=1.3\,m_{\tilde{t}_{\text{R}}},$ $\displaystyle\lvert
A_{b}\rvert$ $\displaystyle=\lvert A_{t}\rvert,$ $\displaystyle
m_{\\{\tilde{q},\tilde{l}\\}_{\\{\text{L},\text{R}\\}}}$
$\displaystyle=2.5\,\text{TeV},$ $\displaystyle A_{\\{q,l\\}}$
$\displaystyle=0,$ $\displaystyle q$ $\displaystyle\in{u,d,s,c},$
$\displaystyle l$ $\displaystyle\in{e,\mu,\tau}\text{ .}$ (58)

Compared to the original $m_{h}^{\text{mod}}$ scenario we choose larger
bilinear soft-breaking parameters for the sfermions, and also larger absolute
values for $\mu$ (see below) and $M_{2}$. Thereby $m_{\tilde{Q}_{3}}$ is
slightly different from $m_{\\{\tilde{t},\tilde{b}\\}_{\text{R}}}$ in order to
avoid numerical instabilities by degeneracies. However, the general feature of
this scenario is kept: it allows for a wide range of
$X_{t}=A_{t}^{*}-\mu/\tan\beta$ to be in agreement with experimental bounds.
With our choice of parameters, $A_{t}$ and $A_{b}$ are not expected to be
affected by constraints from charge- and color-breaking minima Frere:1983ag ;
Gunion:1987qv ; Casas:1995pd ; Hisano:2006mj ; Hisano:2007cz ; Hisano:2008hn ;
Camargo-Molina:2013qva ; Hollik:2016dcm . As $A_{\tau}$ has negligible impact
on the Higgs mass prediction, we set it to zero.

Figure 7:  Prediction for the light Higgs-boson mass $M_{h_{1}}$ (left) and
the mass shifts $\Delta M_{h_{1}}$, $\Delta\tilde{M}_{h_{1}}$ (right, as
defined in Eq. (57)) as a function of $\text{tan}\beta$ using $m_{A}$ as input
mass for different values of $\mu$. Parameters are as described in (5.2).

Figure 8:  Prediction for the light Higgs-boson mass ${M}_{h_{1}}$ (left) and
the mass shifts $\Delta M_{h_{1}}$, $\Delta\tilde{M}_{h_{1}}$ (right, as
defined in Eq. (57)) as a function of $\text{tan}\beta$ using $m_{H^{\pm}}$ as
input mass for different values of $\mu$. The black lines show the results of
Fig. 8 for $\mu=-1500$ GeV. The results of Fig. 8 for $\mu=500$ GeV are
indicated by grey lines, which are underneath the blue lines. Parameters are
as described in Eq. (5.2).

First, the dependence of the lightest Higgs-boson mass $M_{h_{1}}$ on
$\tan\beta$ is analyzed for different values of the $\mu$ parameter. Setting
all phases of the parameters that can be complex to zero, our result can be
compared to previous ones in the MSSM with real parameters where the
corrections evaluated in the present paper were not included. In the
considered scenario, it is possible to choose either $m_{A}$ or
$m_{H^{\pm}}=\sqrt{m_{A}^{2}+M_{W}^{2}}$ as an input parameter which is
renormalized on-shell accordingly. The chosen input mass for Fig. 8 is
$m_{A}$. A comparison of the predicted mass from FeynHiggs-2.12.0, with
($M_{h_{1}}^{\text{new}}$) and without ($M_{h_{1}}^{\text{old}}$)
incorporating our new corrections is shown. Solid lines depict the new, dashed
lines the previous results. In order to illustrate the different relative
sizes of our new contributions, we further plot
$\tilde{M}_{h_{1}}^{\text{old}}$, where the FeynHiggs result for
$M_{h_{1}}^{\text{old}}$ is supplemented with the
$\mathcal{O}{\left(\alpha_{b}\alpha_{s}\right)}$ terms known in the MSSM with
real parameters (dotted lines). The prediction with $\mu=500$ GeV is shown in
blue, while the resulting Higgs-boson mass using $\mu=-1500$ GeV is shown in
red. The blue dashed and blue dotted lines are lying on top of each other
which means that the $\mathcal{O}{\left(\alpha_{b}\alpha_{s}\right)}$
corrections are negligible in this case. The red curves show that our new
corrections are significantly larger than the pure
$\mathcal{O}{\left(\alpha_{b}\alpha_{s}\right)}$ contributions and enter with
different sign. They therefore overcompensate the slight downward shift
induced by the pure $\mathcal{O}{\left(\alpha_{b}\alpha_{s}\right)}$
contributions. The differences $\Delta M_{h_{1}}$ and
$\Delta\tilde{M}_{h_{1}}$, as defined in Eq. (57), are plotted on the right-
hand side of Fig. 8. For low values of $\text{tan}\beta$ the new corrections
slightly increase and then stay constant over a wide range. Only for values
$\tan\beta>40$ and large negative $\mu$ they drop by about $20\%$. Values for
$\tan\beta$ above the depicted range and large negative $\mu$ lead to a
further rapid decrease of $M_{h_{1}}$, eventually yielding a tachyonic Higgs
boson. This is due to the large bottom Yukawa coupling with resummed
$\tan\beta$-enhanced terms which can become non-perturbative in that region of
the parameter space. The rise of the red dotted curve at large $\tan\beta$
reflects that this decrease happens for larger values of $\tan\beta$ once our
new corrections are taken into account.

In Fig. 8 the charged Higgs mass $m_{H^{\pm}}$ is used as an input parameter.
The latter implies the occurrence of terms of
$\mathcal{O}{\left(\sqrt{\alpha_{q}}\sqrt{\alpha_{o}}\alpha_{s}\right)}$ and
corresponds to the renormalization scheme compatible with both the MSSM with
real and complex parameters. On the left-hand side of Fig. 8, the blue
($\mu=500$ GeV) and red ($\mu=-1500$ GeV) lines show the prediction for the
lightest Higgs mass with (solid) and without (dashed) our new contributions.
In addition, the solid and dashed curves of Fig. 8 are indicated again as grey
($\mu=500$ GeV) and black ($\mu=-1500$ GeV) lines. In this way, the influence
of the two different renormalization schemes on the Higgs-mass prediction can
be seen. While the blue and grey lines lie on top of each other over the whole
range of $\tan\beta$, deviations of up to $1.5$ GeV can be observed between
the red and black curves in the region of large $\tan\beta$. Since the slope
of the red curves for large $\tan\beta$ is smaller than for the black curves,
the renormalization scheme with $m_{H^{\pm}}$ as input parameter is better
suited for this particular region in parameter space. On the right-hand side
of Fig. 8 the mass shifts $\Delta M_{h_{1}}$ and $\Delta\tilde{M}_{h_{1}}$
resulting from our new contributions are depicted. The color coding is the
same as described before. The size of the shifts is almost invariant under the
exchange of $m_{A}$ and $m_{H^{\pm}}$ as input parameter, since only small
differences between the two renormalization schemes can be noticed.

We note that setting $\mu=1500$ GeV and using $m_{H^{\pm}}$ as input, the same
qualitative behavior as for the lower positive $\mu$ value can be observed,
with the new contributions being of the same size as for $\mu=-1500$ GeV in
the low and intermediate $\tan\beta$ region. Furthermore, the size of the mass
shift $\Delta M_{h_{1}}$ in Figs. 8 and 8 shows a similar tendency with
respect to the chosen sfermion masses as depicted in Tab. 2, i. e. larger
scales increase the size of the new corrections. However, for stop- and
sbottom masses larger than $\approx 2$ TeV logarithmic contributions of higher
order also become important. Then, a resummation of these logarithms should be
taken into account for an accurate Higgs-mass prediction. The gluino mass can
have a sizable impact due to its appearance in the threshold correction of
$\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$.777In the scenario of Eq.
(5.2) the NNLL-resummation of logarithms with stop masses as implemented in
FeynHiggs gives rise to an upward shift of the lightest Higgs mass
$M_{h}^{\text{old}}$ by $\approx 1.5$ GeV over the whole range of $\tan\beta$;
for stop-mass scales at $3$ TeV this shift already amounts to $\approx 2.7$
GeV. Changing the gluino mass from $2.5$ TeV to $4$ TeV corresponds to a
downward shift of the SM-like Higgs mass due to the threshold corrections of
$\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$ of $\approx 0.5$ GeV over the
whole range of $\tan\beta$.

###  5.3 Scenario 2: large $\tan\beta$

Scenarios with large values of $\tan\beta$ are particularly interesting for
investigating effects of the new contributions in the bottom and sbottom
sector. In that parameter region, terms proportional to the bottom Yukawa
coupling can be as important as terms from the top sector. In the following,
we investigate the dependence of the new contributions on various parameters
at a fixed large $\tan\beta$ value. In order to be consistent with
experimental constraints by ATLAS and CMS we choose a sufficiently large value
of $m_{H^{\pm}}$ Khachatryan:2014wca ; ATLAS:2016fpj . If not stated
otherwise, the MSSM model parameters are

$\displaystyle\tan\beta$ $\displaystyle=50,$ $\displaystyle\mu$
$\displaystyle=-1.5\,\text{TeV},$ $\displaystyle m_{H^{\pm}}$
$\displaystyle=1.5\,\text{TeV},$ $\displaystyle M_{2}$
$\displaystyle=500\,\text{GeV},$ $\displaystyle|M_{3}|$
$\displaystyle=2.5\,\text{TeV},$ $\displaystyle
m_{\\{\tilde{t},\tilde{b}\\}_{\text{L}}}=m_{\tilde{Q}_{3}}$
$\displaystyle=2.1\,\text{TeV},$ $\displaystyle
m_{\\{\tilde{t},\tilde{b}\\}_{\text{R}}}$ $\displaystyle=2\,\text{TeV},$
$\displaystyle\lvert X_{t}\rvert$
$\displaystyle=1.3\,m_{\tilde{b}_{\text{R}}},$ $\displaystyle\lvert
A_{b}\rvert$ $\displaystyle=\lvert A_{t}\rvert,$ $\displaystyle
m_{\\{\tilde{q},\tilde{l}\\}_{\\{\text{L},\text{R}\\}}}$
$\displaystyle=2.5\,\text{TeV},$ $\displaystyle A_{\\{q,l\\}}$
$\displaystyle=0,$ $\displaystyle q$ $\displaystyle\in{u,d,s,c},$
$\displaystyle l$ $\displaystyle\in{e,\mu,\tau}\text{ .}$ (59)

In Fig. 10 the mass shift $\Delta M_{h_{1}}$ is displayed as a function of
$\mu$. Over a wide range the mass shift is nearly constant at about $\Delta
M_{h_{1}}\approx 0.85$ GeV. Only for large negative values $\mu\lesssim-1.8$
TeV, the correction to the lightest Higgs falls steeply indicating a parameter
region where the perturbative prediction for $M_{h_{1}}$ becomes unreliable
owing to the large value of the bottom Yukawa coupling. Thus, $\mu$ should be
kept above that value. The blue line shows the effect of only the third
generation quarks and squarks in our new contributions. The red line shows the
result where these contributions are supplemented with the corrections of the
first and second generation, neglecting the light quark masses and Yukawa
couplings of the first two generations, $m_{q}=0,\,q\in\\{c,s,u,d\\}$.
Accordingly, the difference between the two curves is given by the pure gauge
contributions of $\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ from the first
and second generation. They are rather small, amounting to about $30$ MeV.

Figure 9:  Variation of the mass shift $\Delta M_{h_{1}}$ with $\mu$. The blue
curve shows the result including contributions only from the 3rd generation.
The red line shows the result where also contributions of the 1st and 2nd
generation are included using the approximation $m_{q}=0,\,q\in\\{c,s,u,d\\}$.
Parameters are as described in Eq. (5.3).

Figure 10:  Variation of the mass shift $\Delta M_{h_{1}}$ with the absolute
value and phase of the gluino mass parameter $M_{3}=\lvert
M_{3}\rvert\exp{\left({i\mkern 1.0mu}\,\phi_{M_{3}}\right)}$. The vertical
dashed lines are at $\lvert M_{3}\rvert=1900$ GeV and $2500$ GeV. The
dependence on $\phi_{M_{3}}$ at those values of $\lvert M_{3}\rvert$ is
illustrated in Fig. 13. Parameters are as described in Eq. (5.3).

The variation of $\Delta M_{h_{1}}$ with the gluino-mass parameter
$M_{3}=|M_{3}|\exp{\left({i\mkern 1.0mu}\,\phi_{M_{3}}\right)}$ is shown in
Fig. 10. Close to $|M_{3}|\approx 1.9$ TeV, thresholds of the
gluino–fermion–sfermion system can be observed, which are introduced by one-
loop integrals entering via the subloop-renormalization and resummation of the
bottom Yukawa coupling. The effect of varying the absolute value of the
gluino-mass parameter $\lvert M_{3}\rvert$ on $\Delta M_{h_{1}}$ is strongest
for $\phi_{M_{3}}=0$ and successively weakened as $\phi_{M_{3}}$ approaches
$\pi$. The results for $\phi_{M_{3}}=\pm\frac{\pi}{2}$ almost lie on top of
each other.

Figure 11:  Variation of the light Higgs-boson mass $M_{h_{1}}$ (left) and the
mass shift $\Delta M_{h_{1}}$ (right) with the gluino phase $\phi_{M_{3}}$,
while all other phases are set to zero. Parameters are as described in Eq.
(5.3).

Figure 12:  Variation of the light Higgs-boson mass $M_{h_{1}}$ (left) and the
mass shift $\Delta M_{h_{1}}$ (right) with the phase $\phi_{A_{t}}$ for
different $\phi_{M_{3}}$ and $\phi_{A_{b}}=0$. Parameters are as described in
Eq. (5.3).

Figure 13:  Variation of the light Higgs-boson mass $M_{h_{1}}$ (left) and the
mass shift $\Delta M_{h_{1}}$ (right) with the phase $\phi_{A_{b}}$ for
different $\phi_{M_{3}}$ and $\phi_{A_{t}}=0$. The results for
$\phi_{M_{3}}=\pm\frac{\pi}{2}$ lie on top of each other. Parameters are as
described in Eq. (5.3).

In Figs. 13, 13 and 13 the dependence on the three phases
$\phi_{M_{3}},\,\phi_{A_{t}}$ and $\phi_{A_{b}}$ is displayed, respectively.
The impact of the new (solid) corrections in comparison with the ones
implemented so far in FeynHiggs (dashed) are shown for the lightest Higgs-
boson mass on the left-hand side of each figure, while the differences $\Delta
M_{h_{1}}$ are shown on the right-hand side. Comparing to the MSSM with real
parameters, where the phases are equal to zero or $\pi$, sizable differences
for the prediction of the lightest Higgs-boson mass are visible. Concerning
the total variation of $M_{h_{1}}$ including all now available corrections,
the impact of the phases $\phi_{A_{t}}$ and $\phi_{M_{3}}$ is seen to be
rather large with effects that can exceed 2 GeV, while varying the phase
$\phi_{A_{b}}$ yields only rather small shifts of $\approx 0.2\,\text{GeV}$.

The prediction for $M_{h_{1}}$ as function of $\phi_{M_{3}}$ shown in Fig. 13
is symmetric with respect to the sign of $\phi_{M_{3}}$. The variation of
$\Delta M_{h_{1}}$ with $\phi_{M_{3}}$ is shown on the right-hand side of Fig.
13. The pronounced dependence on the absolute value of $|M_{3}|$ seen in Fig.
10 can be observed again. The variation of $\phi_{M_{3}}$ changes $\Delta
M_{h_{1}}$ by up to $250$ MeV for an $|M_{3}|$ value around the
gluino–fermion–sfermion threshold, while for $|M_{3}|=2.5$ GeV $\Delta
M_{h_{1}}$ is shifted only by up to $70$ MeV.

The phase dependence of $\Delta M_{h_{1}}$ on $\phi_{A_{t}}$ and
$\phi_{A_{b}}$ is shown on the right-hand side of Fig. 13 and Fig. 13,
respectively. The variation of $\Delta M_{h_{1}}$ with $\phi_{A_{t}}$ and
$\phi_{A_{b}}$ is seen to be rather small. It reaches up to $150$ MeV for the
phase $\phi_{A_{t}}$ and up to $50$ MeV for $\phi_{A_{b}}$. It should be noted
that the results for $\phi_{M_{3}}=\pm\frac{\pi}{2}$ lie on top of each other
in Fig. 13. While the variation with $\phi_{A_{b}}$ is rather small for any
non-zero $\phi_{M_{3}}$, the variation with $\phi_{A_{t}}$ is minimal for
$\phi_{M_{3}}=0$ and maximal for $\phi_{M_{3}}=\pi$. Using different values of
$\phi_{A_{b}}$ (and keeping $\phi_{M_{3}}$ fixed) has only a small effect on
the variation of $\Delta M_{h_{1}}$ with $\phi_{A_{t}}$. The corresponding
plot is therefore not shown here.

###  5.4 Scenario 3: low $M_{H}$

In the low-$M_{H}$ scenario the observed SM-like Higgs boson with a mass of
about $125$ GeV can be identified with the next-to-lightest neutral $CP$-even
Higgs boson of the MSSM, see Ref. Bechtle:2016kui  for a recent update. We
choose the following MSSM model parameters,

$\displaystyle\tan\beta$ $\displaystyle=6.5,$ $\displaystyle\mu$
$\displaystyle=5\,\text{TeV},$ $\displaystyle M_{2}$
$\displaystyle=300\,\text{GeV},$ $\displaystyle\lvert M_{3}\rvert$
$\displaystyle=1.5\,\text{TeV},$ $\displaystyle
m_{\\{\tilde{t},\tilde{b}\\}_{\\{\text{L},\text{R}\\}}}$
$\displaystyle=750\,\text{GeV},$ $\displaystyle
m_{\tilde{\tau}_{\\{\text{L},\text{R}\\}}}$ $\displaystyle=500\,\text{GeV},$
$\displaystyle m_{\tilde{q}_{\\{\text{L},\text{R}\\}}}$
$\displaystyle=1.5\,\text{TeV},$ $\displaystyle
m_{\tilde{l}_{\\{\text{L},\text{R}\\}}}$ $\displaystyle=250\,\text{GeV},$
$\displaystyle A_{t}=A_{b}=A_{\tau}$ $\displaystyle=-70\,\text{GeV},$
$\displaystyle A_{\\{q,l\\}}$ $\displaystyle=0,$ $\displaystyle q$
$\displaystyle\in{u,d,s,c},$ $\displaystyle l$ $\displaystyle\in{e,\mu}\text{
.}$ (60)

Compared to the original scenario in Bechtle:2016kui  we had to choose a
smaller value of $\mu$ in order to avoid a tachyonic lightest Higgs boson for
a charged Higgs mass $m_{H^{\pm}}\approx 160$ GeV. Our value for $\tan\beta$
is chosen such that the scenario is valid according to Fig. 26 of
Bechtle:2016kui .

Figure 14:  Variation of the three neutral Higgs-boson masses $M_{h_{i}}$ with
the charged Higgs boson mass $m_{H^{\pm}}$. The results for
$M_{h_{i}}^{\text{new}}$ are shown as full lines and those for
$M_{h_{i}}^{\text{old}}$ as dotted lines. Parameters are as described in Eq.
(5.4). Figure 15: Variation of the three neutral Higgs-boson masses
$M_{h_{i}}$ with the charged Higgs boson mass $m_{H^{\pm}}$ for non-zero
phases $\phi_{M_{3}}$. The results for $M_{h_{i}}^{\text{new}}$ are shown as
full lines and those for $M_{h_{i}}^{\text{old}}$ as dotted lines. The results
of Fig. 14 with $\phi_{M_{3}}=0$ are depicted in grey for reference.
Parameters are as described in Eq. (5.4).

In Fig. 14 the three neutral Higgs-boson masses are depicted, varying the
charged Higgs-boson mass $m_{H^{\pm}}$ which is used as an input parameter.
The light green band illustrates the mass range of $125\pm 3$ GeV; it should
be interpreted as a rough indication of the mass range which is theoretically
in agreement with the detected Higgs boson. Up to $m_{H^{\pm}}\lesssim 188$
GeV the heavier Higgs $h_{2}$ could be associated with the discovered Higgs-
like particle; however, as can be seen in the low-$M_{H}^{\text{alt}+}$
scenario in Fig. 26 of Bechtle:2016kui , our choice of $\mu$ and $\tan\beta$
is already excluded for a charged Higgs mass $m_{H^{\pm}}=185$ GeV. Yet,
scenarios with values of $m_{H^{\pm}}$ closer to or below $m_{t}$ are still
allowed. In this region the new corrections presented here have a negligible
impact on $M_{h_{2}}$, but lead to a downward shift of about $1$ GeV for both
$M_{h_{1}}$ and $M_{h_{3}}$.

As shown in Fig. 15, using a non-zero value of the gluino phase of
$\phi_{M_{3}}=\pi/2$ or $\phi_{M_{3}}=\pi$ shifts all three neutral Higgs
masses to larger values as compared to the case $\phi_{M_{3}}=0$. For better
comparison, the results of Fig. 14 are underlaid in grey. The numerical impact
of the new contributions presented here rises with increasing $\phi_{M_{3}}$.
For $\phi_{M_{3}}=\pi$ all neutral Higgs masses can receive large corrections
of up to $5$ GeV.

##  6 Conclusions

We have computed the full two-loop QCD corrections to the lightest Higgs-boson
mass in the MSSM with complex parameters. Compared to previous works, this
primarily involves going beyond the gaugeless limit, and including a finite
bottom-quark mass; furthermore the momentum dependence of loop integrals is
taken into account. On the technical side, this involves the computation of
177 different mass topologies evaluated at different kinematical
configurations, amounting to a total of 513 two-point two-loop integrals with
up to five mass scales. These integrals have been computed numerically with
the program SecDec.

In the first part of our numerical analysis, we have compared our results with
earlier result in the literature taking the appropriate limit of real
parameters and / or vanishing external momentum of our results. We have found
very good agreement with the existing results in the approriate limit if the
same renormalization scheme is employed. The contributions evaluated in this
paper yield a shift in the lightest Higgs-boson mass at the level of 1 GeV,
where the impact has been seen to be more pronounced for an increasing mass
scale of the stops.

We have furthermore investigated the dependence of the new corrections on
$\tan\beta$ choosing different values of the $\mu$-parameter as well as
different renormalization schemes. For a large negative $\mu$ the corrections
are generally larger and amount to around $~{}0.9$ GeV in $M_{h_{1}}$. The
corrections are largest for $10<\tan\beta<30$, decrease by $3\%$ for lower
values and by about $20\%$ beyond $\tan\beta=30$.

We find non-vanishing mixed up- and down-type Yukawa corrections in the
charged Higgs-boson self-energy correction entering the mass predictions for
the neutral Higgs bosons as renormalization constant if the charged Higgs mass
$m_{H^{\pm}}$ instead of the neutral $\mathcal{CP}$-odd mass $m_{A}$ is chosen
as an input parameter. We have compared the mass prediction for the lightest
Higgs boson in both schemes and have found good agreement in general. However,
using the charged Higgs mass as an input parameter yields better numerical
stability at large $\tan\beta$ and large negative $\mu$.

The Yukawa contributions scale according to their Yukawa couplings, leading to
much smaller contributions from the first and second generation quarks and
squarks. The pure gauge terms of $\mathcal{O}{\left(\alpha\alpha_{s}\right)}$
in the limit of massless quarks are found to be of similar small size, below
$20$ MeV for one generation.

Analyzing the dependence on the gluino mass, we have found maximal shifts of
$\approx 900$ MeV in $M_{h_{1}}$. The corrections show a sensitive dependence
on the gluino–fermion–sfermion threshold, which enters via the counterterms of
our renormalization scheme, and the gluino phase. For the $\mu$-parameter a
mass shift of the lightest Higgs by $\approx 850$ MeV is found over large
regions of parameter space.

Concerning the impact of the three phases $\phi_{M_{3}}$, $\phi_{A_{t}}$ and
$\phi_{A_{b}}$, we find significant effects in our new corrections from
varying the gluino phase and the pase of $A_{t}$. For $\phi_{M_{3}}$ the phase
dependence becomes particularly pronounced in the threshold region of the
gluino–fermion–sfermion system, as mentioned above.

Besides scenarios where the lightest neutral Higgs boson in the spectrum of
the MSSM is the SM-like state that can be identified with the detected Higgs
signal, we have also analyzed the impact of the newly computed contributions
on the Higgs-mass predictions for the three neutral Higgs bosons within the
low-$M_{H}$ scenario for different values of the gluino phase $\phi_{M_{3}}$.
We have found mass corrections of $\approx 1$ GeV for $\phi_{M_{3}}=0$ and up
to $\approx 5$ GeV for $\phi_{M_{3}}=\pi$ in this case.

Accordingly, we have found that the subleading two-loop contributions that we
have evaluated in this paper yield a shift in the prediction for the mass of
the light SM-like Higgs boson of the MSSM of up to the level of 1 GeV. The
size of the correction sensitively depends on the mass scales of the stops and
sbottoms, on the absolute value and phase of the gluino mass parameter, as
well as on the absolute value and phase of the trilinear coupling in the stop
sector (and to a lesser extent on the trilinear coupling in the sbottom
sector). While these findings of course have an impact on the remaining
theoretical uncertainties from unknown higher-order corrections, we do not
attempt to provide an improved estimate of the remaining uncertainties here.
Such an improved estimate should be based on a combination of the fixed-order
result considered here with a resummation of higher-order logarithmic
contributions. We leave such an analysis to future work.

It should be noted in this context that our results for the corrections of
$\mathcal{O}{\left(\alpha\alpha_{s}\right)}$ beyond the gaugeless limit cannot
be used directly to infer the possible size of the corresponding contributions
of $\mathcal{O}{\left(\alpha^{2}\right)}$ to the Higgs-boson spectrum, which
are unknown up to now. This is due to the fact that the requirement of a
strong coupling in the corrections that we have evaluated significantly
constrains the structure of the contributing Feynman diagrams, while
additional classes of contributions will have to be taken into account for a
full calculation of the corrections of $\mathcal{O}{\left(\alpha^{2}\right)}$.

The new contributions evaluated in this paper will be made publicly available
in the program FeynHiggs.

## Acknowledgments

We are grateful to S. Di Vita for providing us with results for a detailed
comparison with Ref. Degrassi:2014pfa . We thank T. Hahn, S. Heinemeyer, W.
Hollik and P. Slavich for helpful discussions. S. B. gratefully acknowledges
financial support by the ERC Advanced Grant MC@NNLO (340983) and ERC Starting
Grant "MathAm" (39568) during different stages of this project. The work of S.
P. has been supported by the Collaborative Research Center SFB 676 of the DFG,
“Particles, Strings and the Early Universe”, and by the ANR grant
“HiggsAutomator” (ANR-15-CE31-0002) during different stages of the project.
The work of G. W. has been supported in part by the DFG through the SFB 676
“Particles, Strings and the Early Universe” and by the European Commission
through the “HiggsTools” Initial Training Network PITN-GA-2012-316704.

##  References

  * (1) ATLAS Collaboration,  Phys. Lett. B716, 1 (2012), arXiv:1207.7214.
  * (2) CMS Collaboration,  Phys. Lett. B716, 30 (2012), arXiv:1207.7235.
  * (3) ATLAS Collaboration, CMS Collaboration, D. Sperka,  53rd Rencontres de Moriond, March 2018, La Thuile, Italy, Measurements of the BEH scalar mass and other couplings in ATLAS and CMS.
  * (4) ATLAS Collaboration, CMS Collaboration,  JHEP 08, 045 (2016), arXiv:1606.02266.
  * (5) P. Bechtle, H. Haber, S. Heinemeyer, O. Stål, T. Stefaniak, G. Weiglein, L. Zeune,  Eur. Phys. J. C77, 67 (2017), arXiv:1608.00638.
  * (6) A. Pilaftsis,  Phys. Rev. D58, 096010 (1998), arXiv:hep-ph/9803297.
  * (7) D. Demir,  Phys. Rev. D60, 055006 (1999), arXiv:hep-ph/9901389.
  * (8) A. Pilaftsis, C. Wagner,  Nucl. Phys. B553, 3 (1999), arXiv:hep-ph/9902371.
  * (9) S. Heinemeyer,  Eur. Phys. J. C22, 521 (2001), arXiv:hep-ph/0108059.
  * (10) M. Frank, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein,  JHEP 0702, 047 (2007), arXiv:hep-ph/0611326.
  * (11) G. Degrassi, S. Heinemeyer, W. Hollik, P. Slavich, G. Weiglein,  Eur. Phys. J. C28, 133 (2003), arXiv:hep-ph/0212020.
  * (12) S. Heinemeyer, W. Hollik, G. Weiglein,  Eur. Phys. J. C9, 343 (1999), arXiv:hep-ph/9812472.
  * (13) S. Heinemeyer, W. Hollik, G. Weiglein,  Phys. Rev. D58, 091701 (1998), arXiv:hep-ph/9803277.
  * (14) S. Heinemeyer, W. Hollik, G. Weiglein,  Phys. Lett. B455, 179 (1999), arXiv:hep-ph/9903404.
  * (15) H. Haber, R. Hempfling,  Phys. Rev. Lett. 66, 1815 (1991).
  * (16) J. Ellis, G. Ridolfi, F. Zwirner,  Phys. Lett. B257, 83 (1991).
  * (17) Y. Okada, M. Yamaguchi, T. Yanagida,  Prog. Theor. Phys. 85, 1 (1991).
  * (18) Y. Okada, M. Yamaguchi, T. Yanagida,  Phys. Lett. B262, 54 (1991).
  * (19) J. Ellis, G. Ridolfi, F. Zwirner,  Phys. Lett. B262, 477 (1991).
  * (20) K. Sasaki, M. Carena, C. Wagner,  Nucl. Phys. B381, 66 (1992).
  * (21) P. Chankowski, S. Pokorski, J. Rosiek,  Phys. Lett. B274, 191 (1992).
  * (22) A. Brignole,  Phys. Lett. B281, 284 (1992).
  * (23) R. Hempfling, A. Hoang,  Phys. Lett. B331, 99 (1994), arXiv:hep-ph/9401219.
  * (24) J. Casas, J. Espinosa, M. Quiros, A. Riotto,  Nucl. Phys. B436, 3 (1995), arXiv:hep-ph/9407389,  [Erratum: Nucl. Phys. B439, 466 (1995)].
  * (25) A. Dabelstein,  Z. Phys. C67, 495 (1995), arXiv:hep-ph/9409375.
  * (26) M. Carena, J. Espinosa, M. Quiros, C. Wagner,  Phys. Lett. B355, 209 (1995), arXiv:hep-ph/9504316.
  * (27) M. Carena, M. Quiros, C. Wagner,  Nucl. Phys. B461, 407 (1996), arXiv:hep-ph/9508343.
  * (28) D. Pierce, J. Bagger, K. Matchev, R. Zhang,  Nucl. Phys. B491, 3 (1997), arXiv:hep-ph/9606211.
  * (29) H. Haber, R. Hempfling, A. Hoang,  Z. Phys. C75, 539 (1997), arXiv:hep-ph/9609331.
  * (30) S. Heinemeyer, W. Hollik, G. Weiglein,  Phys. Lett. B440, 296 (1998), arXiv:hep-ph/9807423.
  * (31) R. Zhang,  Phys. Lett. B447, 89 (1999), arXiv:hep-ph/9808299.
  * (32) J. Espinosa, R. Zhang,  JHEP 03, 026 (2000), arXiv:hep-ph/9912236.
  * (33) M. Carena, H. Haber, S. Heinemeyer, W. Hollik, C. Wagner, G. Weiglein,  Nucl. Phys. B580, 29 (2000), arXiv:hep-ph/0001002.
  * (34) J. Espinosa, R. Zhang,  Nucl. Phys. B586, 3 (2000), arXiv:hep-ph/0003246.
  * (35) J. Espinosa, I. Navarro,  Nucl. Phys. B615, 82 (2001), arXiv:hep-ph/0104047.
  * (36) G. Degrassi, P. Slavich, F. Zwirner,  Nucl. Phys. B611, 403 (2001), arXiv:hep-ph/0105096.
  * (37) S. Martin,  Phys. Rev. D65, 116003 (2002), arXiv:hep-ph/0111209.
  * (38) A. Brignole, G. Degrassi, P. Slavich, F. Zwirner,  Nucl. Phys. B631, 195 (2002), arXiv:hep-ph/0112177.
  * (39) A. Dedes, P. Slavich,  Nucl. Phys. B657, 333 (2003), arXiv:hep-ph/0212132.
  * (40) S. Martin,  Phys. Rev. D66, 096001 (2002), arXiv:hep-ph/0206136.
  * (41) A. Dedes, G. Degrassi, P. Slavich,  Nucl. Phys. B672, 144 (2003), arXiv:hep-ph/0305127.
  * (42) S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein,  Eur. Phys. J. C39, 465 (2005), arXiv:hep-ph/0411114.
  * (43) S. Martin,  Phys. Rev. D68, 075002 (2003), arXiv:hep-ph/0307101.
  * (44) S. Martin,  Phys. Rev. D70, 016005 (2004), arXiv:hep-ph/0312092.
  * (45) B. Allanach, A. Djouadi, J. Kneur, W. Porod, P. Slavich,  JHEP 09, 044 (2004), arXiv:hep-ph/0406166.
  * (46) S. Heinemeyer, W. Hollik, G. Weiglein,  Phys. Rept. 425, 265 (2006), arXiv:hep-ph/0412214.
  * (47) S. Martin,  Phys. Rev. D71, 116004 (2005), arXiv:hep-ph/0502168.
  * (48) S. Martin, D. Robertson,  Comput. Phys. Commun. 174, 133 (2006), arXiv:hep-ph/0501132.
  * (49) R. Harlander, P. Kant, L. Mihaila, M. Steinhauser,  Phys. Rev. Lett. 100, 191602 (2008), arXiv:0803.0672.
  * (50) P. Kant, R. Harlander, L. Mihaila, M. Steinhauser,  JHEP 1008, 104 (2010), arXiv:1005.5709.
  * (51) A. Brignole, G. Degrassi, P. Slavich, F. Zwirner,  Nucl. Phys. B643, 79 (2002), arXiv:hep-ph/0206101.
  * (52) R. Harlander, J. Klappert, A. Voigt,  Eur. Phys. J. C77, 814 (2017), arXiv:1708.05720.
  * (53) S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, W. Hollik,  Eur. Phys. J. C74, 2994 (2014), arXiv:1404.7074.
  * (54) G. Degrassi, S. Di Vita, P. Slavich,  Eur. Phys. J. C75, 61 (2015), arXiv:1410.3432.
  * (55) S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, W. Hollik,  Eur. Phys. J. C75, 424 (2015), arXiv:1505.03133.
  * (56) P. Draper, G. Lee, C. Wagner,  Phys. Rev. D89, 055023 (2014), arXiv:1312.5743.
  * (57) J. Pardo Vega, G. Villadoro,  JHEP 07, 159 (2015), arXiv:1504.05200.
  * (58) G. Lee, C. Wagner,  (2015), arXiv:1508.00576.
  * (59) T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein,  Phys. Rev. Lett. 112, 141801 (2014), arXiv:1312.4937.
  * (60) H. Bahl, W. Hollik,  Eur. Phys. J. C76, 499 (2016), arXiv:1608.01880.
  * (61) H. Bahl, S. Heinemeyer, W. Hollik, G. Weiglein,  (2017), arXiv:1706.00346.
  * (62) P. Athron, J. Park, T. Steudtner, D. Stöckinger, A. Voigt,  (2016), arXiv:1609.00371.
  * (63) E. Bagnaschi, G. F. Giudice, P. Slavich, A. Strumia,  JHEP 09, 092 (2014), arXiv:1407.4081.
  * (64) E. Bagnaschi, J. Pardo Vega, P. Slavich,  Eur. Phys. J. C77, 334 (2017), arXiv:1703.08166.
  * (65) S. Martin,  Phys. Rev. D67, 095012 (2003), arXiv:hep-ph/0211366.
  * (66) S. Martin,  Phys. Rev. D71, 016012 (2005), arXiv:hep-ph/0405022.
  * (67) S. Martin,  Phys. Rev. D75, 055005 (2007), arXiv:hep-ph/0701051.
  * (68) S. Choi, M. Drees, J. Lee,  Phys. Lett. B481, 57 (2000), arXiv:hep-ph/0002287.
  * (69) T. Ibrahim, P. Nath,  Phys. Rev. D63, 035009 (2001), arXiv:hep-ph/0008237.
  * (70) T. Ibrahim, P. Nath,  Phys. Rev. D66, 015005 (2002), arXiv:hep-ph/0204092.
  * (71) M. Carena, J. Ellis, A. Pilaftsis, C. Wagner,  Nucl. Phys. B586, 92 (2000), arXiv:hep-ph/0003180.
  * (72) S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein,  Phys. Lett. B652, 300 (2007), arXiv:0705.0746.
  * (73) W. Hollik, S. Paßehr,  Phys. Lett. B733, 144 (2014), arXiv:1401.8275.
  * (74) W. Hollik, S. Paßehr,  JHEP 10, 171 (2014), arXiv:1409.1687.
  * (75) S. Paßehr, G. Weiglein,  (2017), arXiv:1705.07909.
  * (76) M. Goodsell, F. Staub,  Eur. Phys. J. C77, 46 (2017), arXiv:1604.05335.
  * (77) S. Liebler, S. Patel, G. Weiglein,  Eur. Phys. J. C77, 305 (2017), arXiv:1611.09308.
  * (78) S. Heinemeyer, W. Hollik, G. Weiglein,  Comput. Phys. Commun. 124, 76 (2000), arXiv:hep-ph/9812320.
  * (79) J. Lee, A. Pilaftsis, M. Carena, S. Choi, M. Drees, J. Ellis, C. Wagner,  Comput. Phys. Commun. 156, 283 (2004), arXiv:hep-ph/0307377.
  * (80) K. Williams, H. Rzehak, G. Weiglein,  Eur. Phys. J. C71, 1669 (2011), arXiv:1103.1335.
  * (81) ACME Collaboration,  Science 343, 269 (2014), arXiv:1310.7534.
  * (82) M. Pospelov, A. Ritz,  Annals Phys. 318, 119 (2005), arXiv:hep-ph/0504231.
  * (83) J. Pendlebury, _et al._ ,  Phys. Rev. D92, 092003 (2015), arXiv:1509.04411.
  * (84) C. Baker, _et al._ ,  Phys. Rev. Lett. 97, 131801 (2006), arXiv:hep-ex/0602020.
  * (85) C. Baker, _et al._ ,  Phys. Rev. Lett. 98, 149102 (2007), arXiv:0704.1354.
  * (86) A. Serebrov, _et al._ ,  Phys. Rev. C92, 055501 (2015).
  * (87) B. Graner, Y. Chen, E. Lindahl, B. Heckel,  Phys. Rev. Lett. 116, 161601 (2016), arXiv:1601.04339,  [Erratum: Phys. Rev. Lett.119,no.11,119901(2017)].
  * (88) N. Yamanaka, B. Sahoo, N. Yoshinaga, T. Sato, K. Asahi, B. Das,  Eur. Phys. J. A53, 54 (2017), arXiv:1703.01570.
  * (89) P. Drechsel, L. Galeta, S. Heinemeyer, G. Weiglein,  Eur. Phys. J. C77, 42 (2017), arXiv:1601.08100.
  * (90) P. Drechsel, R. Gröber, S. Heinemeyer, M. Mühlleitner, H. Rzehak, G. Weiglein,  Eur. Phys. J. C77, 366 (2017), arXiv:1612.07681.
  * (91) F. Domingo, P. Drechsel, S. Paßehr,  Eur. Phys. J. C77, 562 (2017), arXiv:1706.00437.
  * (92) G. Weiglein, R. Scharf, M. Böhm,  Nucl. Phys. B416, 606 (1994), arXiv:hep-ph/9310358.
  * (93) T. Hahn, S. Paßehr,  Comput. Phys. Commun. 214, 91 (2017), arXiv:1508.00562.
  * (94) J. Carter, G. Heinrich,  Comput. Phys. Commun. 182, 1566 (2011), arXiv:1011.5493.
  * (95) S. Borowka, J. Carter, G. Heinrich,  Comput. Phys. Commun. 184, 396 (2013), arXiv:1204.4152.
  * (96) S. Borowka, G. Heinrich,  Comput. Phys. Commun. 184, 2552 (2013), arXiv:1303.1157.
  * (97) T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein,  Nucl. Phys. Proc. Suppl. 205-206, 152 (2010), arXiv:1007.0956.
  * (98) R. D. Peccei, H. R. Quinn,  Phys. Rev. Lett. 38, 1440 (1977).
  * (99) S. Dimopoulos, S. D. Thomas,  Nucl. Phys. B465, 23 (1996), arXiv:hep-ph/9510220.
  * (100) E. Fuchs, G. Weiglein,  JHEP 09, 079 (2017), arXiv:1610.06193.
  * (101) J. Küblbeck, M. Böhm, A. Denner,  Comput. Phys. Commun. 60, 165 (1990).
  * (102) T. Hahn,  Comput. Phys. Commun. 140, 418 (2001), arXiv:hep-ph/0012260.
  * (103) A. von Manteuffel, C. Studerus,  (2012), arXiv:1201.4330.
  * (104) T. Hahn, M. Pérez-Victoria,  Comput. Phys. Commun. 118, 153 (1999), arXiv:hep-ph/9807565.
  * (105) A. Freitas, D. Stockinger,  Phys. Rev. D66, 095014 (2002), arXiv:hep-ph/0205281.
  * (106) M. Sperling, D. Stöckinger, A. Voigt,  JHEP 07, 132 (2013), arXiv:1305.1548.
  * (107) M. Sperling, D. Stöckinger, A. Voigt,  JHEP 01, 068 (2014), arXiv:1310.7629.
  * (108) N. Baro, F. Boudjema, A. Semenov,  Phys. Rev. D78, 115003 (2008), arXiv:0807.4668.
  * (109) S. Heinemeyer, H. Rzehak, C. Schappacher,  Phys. Rev. D82, 075010 (2010), arXiv:1007.0689.
  * (110) T. Banks,  Nucl. Phys. B303, 172 (1988).
  * (111) L. Hall, R. Rattazzi, U. Sarid,  Phys. Rev. D50, 7048 (1994), arXiv:hep-ph/9306309.
  * (112) R. Hempfling,  Phys. Rev. D49, 6168 (1994).
  * (113) M. Carena, M. Olechowski, S. Pokorski, C. Wagner,  Nucl. Phys. B426, 269 (1994), arXiv:hep-ph/9402253.
  * (114) M. Carena, D. Garcia, U. Nierste, C. Wagner,  Nucl. Phys. B577, 88 (2000), arXiv:hep-ph/9912516.
  * (115) H. Eberl, K. Hidaka, S. Kraml, W. Majerotto, Y. Yamada,  Phys. Rev. D62, 055006 (2000), arXiv:hep-ph/9912463.
  * (116) T. Takagi,  Japanese J. Math. 1, 83 (1927).
  * (117) S. Bloch, P. Vanhove,  J. Number Theor. 148, 328 (2015), arXiv:1309.5865.
  * (118) S. Bloch, M. Kerr, P. Vanhove,  Compos. Math. 151, 2329 (2015), arXiv:1406.2664.
  * (119) L. Adams, C. Bogner, S. Weinzierl,  J. Math. Phys. 54, 052303 (2013), arXiv:1302.7004.
  * (120) E. Remiddi, L. Tancredi,  Nucl. Phys. B880, 343 (2014), arXiv:1311.3342.
  * (121) L. Adams, C. Bogner, S. Weinzierl,  J. Math. Phys. 56, 072303 (2015), arXiv:1504.03255.
  * (122) L. Adams, C. Bogner, S. Weinzierl,  J. Math. Phys. 57, 032304 (2016), arXiv:1512.05630.
  * (123) S. Bloch, M. Kerr, P. Vanhove,  (2016), arXiv:1601.08181.
  * (124) E. Remiddi, L. Tancredi,  Nucl. Phys. B907, 400 (2016), arXiv:1602.01481.
  * (125) L. Adams, C. Bogner, A. Schweitzer, S. Weinzierl,  J. Math. Phys. 57, 122302 (2016), arXiv:1607.01571.
  * (126) J. Brödel, C. Duhr, F. Dulat, L. Tancredi,  (2017), arXiv:1712.07095.
  * (127) S. Agrawal, T. Hahn, E. Mirabella,  J. Phys. Conf. Ser. 368, 012054 (2012), arXiv:1112.0124.
  * (128) T. Hahn,  Comput. Phys. Commun. 168, 78 (2005), arXiv:hep-ph/0404043.
  * (129) M. Carena, S. Heinemeyer, O. Stål, C. E. M. Wagner, G. Weiglein,  Eur. Phys. J. C73, 2552 (2013), arXiv:1302.7033.
  * (130) A. Arbey, J. Ellis, R. Godbole, F. Mahmoudi,  Eur. Phys. J. C75, 85 (2015), arXiv:1410.4824.
  * (131) E. Bagnaschi, R. Harlander, S. Liebler, H. Mantler, P. Slavich, A. Vicini,  JHEP 06, 167 (2014), arXiv:1404.0327.
  * (132) J. Frere, T. Jones, S. Raby,  Nucl. Phys. B222, 11 (1983).
  * (133) J. Gunion, H. Haber, M. Sher,  Nucl. Phys. B306, 1 (1988).
  * (134) J. Casas, A. Lleyda, C. Muñoz,  Nucl. Phys. B471, 3 (1996), arXiv:hep-ph/9507294.
  * (135) J. Hisano, M. Nagai, P. Paradisi,  Phys. Lett. B642, 510 (2006), arXiv:hep-ph/0606322.
  * (136) J. Hisano, M. Nagai, P. Paradisi,  Phys. Rev. D78, 075019 (2008), arXiv:0712.1285.
  * (137) J. Hisano, M. Nagai, P. Paradisi,  Phys. Rev. D80, 095014 (2009), arXiv:0812.4283.
  * (138) J. Camargo-Molina, B. O’Leary, W. Porod, F. Staub,  Eur. Phys. J. C73, 2588 (2013), arXiv:1307.1477.
  * (139) W. G. Hollik,  JHEP 08, 126 (2016), arXiv:1606.08356.
  * (140) CMS Collaboration,  JHEP 10, 160 (2014), arXiv:1408.3316.
  * (141) ATLAS Collaboration,  (2016), CDS 2206278.