title.1title.1FeynHiggs and a template for two-loop self-energiesFeynHiggs and
a template for two-loop self-energiestitle.1

MPP

DESY

FeynHiggs and a template for two-loop self-energies

Thomas Hahn1***email<EMAIL_ADDRESS>and Sebastian Paßehr2†††email:
<EMAIL_ADDRESS>
1Max-Planck-Institut für Physik
(Werner-Heisenberg-Institut),
Föhringer Ring 6, D–80805 München, Germany

2Deutsches Elektronen-Synchrotron DESY,
Notkestraße 85, D–22607 Hamburg, Germany

###### Abstract

Recently, the leading two-loop corrections of
$\mathcal{O}{\left(\alpha_{t}^{2}\right)}$ to the Higgs-boson spectrum of the
complex MSSM have been implemented into the public code FeynHiggs. For this
purpose the whole calculation of these contributions has been entirely
revised, automatized and split up into several tasks which are explained in
detail in this letter. Although this procedure is specified for the
$\mathcal{O}{\left(\alpha_{t}^{2}\right)}$ self-energy corrections of the
Higgs bosons, it can be utilized as a template and adapted to similar
calculations.

##  1 Introduction

The Higgs-like boson discovered at LHC Aad:2012tfa ; Chatrchyan:2012ufa
initiated great efforts to pin it down as a particle responsible for
electroweak symmetry breaking. Although the measured properties of this new
boson are consistent with the expectations for the Standard Model (SM) Higgs
boson Moriond2015 , a variety of beyond-SM interpretations is still possible.
One such possibility is the interpretation as a light state within a richer
spectrum of scalar particles in the theoretically well motivated Minimal
Supersymmetric Standard Model (MSSM). Its Higgs sector consists of two complex
scalar doublets leading to five physical Higgs bosons and three (pseudo-)
Goldstone bosons. The physical mass eigenstates at the tree level are the
neutral $CP$-even $h$, $H$, the $CP$-odd $A$, and the charged $H^{\pm}$
bosons. Their masses can be parameterized by the charged (or $CP$-odd) Higgs-
boson mass and the ratio of the two vacuum expectation values,
$\tan\beta=v_{2}/v_{1}$. $CP$-violation in the Higgs sector is induced by
complex parameters in other sectors of the MSSM via loop corrections, leading
to mixing between $h$, $H$, and $A$ in the mass eigenstates Pilaftsis:1998pe ;
Demir:1999hj ; Pilaftsis:1999qt ; Heinemeyer:2001qd .

Since the masses of the neutral Higgs bosons are strongly affected by loop
contributions, a lot of work has been invested into higher-order calculations
of the mass spectrum from the SUSY parameters, in the case of the real MSSM
Heinemeyer:1998jw ; Heinemeyer:1998np ; Heinemeyer:1999be ; Heinemeyer:2004xw
; Borowka:2014wla ; Degrassi:2014pfa ; Borowka:2015ura ; mhiggsFD3l ;
Zhang:1998bm ; Espinosa:2000df ; Brignole:2001jy ; Casas:1994us ;
Degrassi:2002fi ; Heinemeyer:2004gx ; Allanach:2004rh ; Martin:2001vx  as well
as for the MSSM with complex parameters (cMSSM) Demir:1999hj ; Choi:2000wz ;
Pilaftsis:1999qt ; Heinemeyer:2001qd ; Carena:2000yi ; Frank:2006yh ;
Heinemeyer:2007aq ; Hollik:2014wea ; Hollik:2014bua . At one-loop order, the
dominant contributions arise from the Yukawa sector with the large top-Yukawa
coupling $h_{t}$, or $\alpha_{t}=\left.h_{t}^{2}\middle/(4\pi)\right.$. The
class of leading two-loop Yukawa-type corrections of
$\mathcal{O}(\alpha_{t}^{2})$ has recently been calculated for the case of
complex parameters Hollik:2014wea ; Hollik:2014bua . Together with the full
one-loop result Frank:2006yh  and the leading
$\mathcal{O}{\left(\alpha_{t}\alpha_{s}\right)}$ terms Heinemeyer:2007aq  it
is now available in the public program FeynHiggs Heinemeyer:1998np ;
Degrassi:2002fi ; Frank:2006yh ; Heinemeyer:1998yj ; Hahn:2009zz ; Hahn:2010te
.

The $\mathcal{O}(\alpha_{t}^{2})$ contributions in the cMSSM have been
evaluated in the Feynman-diagrammatic approach, extending the on-shell
renormalization scheme of Ref. Frank:2006yh  to the two-loop level. This
ensures that the obtained analytical results for the renormalized two-loop
self-energies can consistently be incorporated in FeynHiggs. For
implementation in FeynHiggs the whole calculation of the
$\mathcal{O}(\alpha_{t}^{2})$ contributions was reworked and broken up into
several working steps. Together they provide a semi-automated framework for a
two-loop calculation in a model with non-trivial renormalization. Although the
code is specific to the evaluation of the mentioned corrections, it can be
utilized as a template for similar calculations. The aim of this paper is to
explain in detail how the code works and what is done at each step, and where
process-specific adjustments are necessary.

The paper is organized as follows: Sect. 2 provides the theoretical foundation
of the computation of the Higgs-mass corrections in the cMSSM. The
approximations used for the evaluation of the $\mathcal{O}(\alpha_{t}^{2})$
contributions to the Higgs masses are explained in Sect. 3. A detailed
description of the implementation of the $\mathcal{O}(\alpha_{t}^{2})$
corrections in FeynHiggs is given in Sect. 4.

##  2 Masses of the Higgs bosons in the Complex MSSM

Since the calculation of the $\mathcal{O}(\alpha_{t}^{2})$ corrections has
already been presented in great detail elsewhere Hollik:2014wea ;
Hollik:2014bua  we shall recap only the parts necessary for the implementation
in FeynHiggs here.

In the basis of the tree-level Higgs states the loop-corrected mass matrices
are given by Frank:2006yh

$\displaystyle\mathbf{M}_{0}^{2}$
$\displaystyle=\begin{pmatrix}m_{h}^{2}-\hat{\Sigma}_{hh}&-\hat{\Sigma}_{hH}&-\hat{\Sigma}_{hA}&\greyed{-\hat{\Sigma}_{hG}}\\\
-\hat{\Sigma}_{hH}&m_{H}^{2}-\hat{\Sigma}_{HH}&-\hat{\Sigma}_{HA}&\greyed{-\hat{\Sigma}_{HG}}\\\
-\hat{\Sigma}_{hA}&-\hat{\Sigma}_{HA}&m_{A}^{2}-\hat{\Sigma}_{AA}&\greyed{-\hat{\Sigma}_{AG}}\\\
\greyed{-\hat{\Sigma}_{hG}}&\greyed{-\hat{\Sigma}_{HG}}&\greyed{-\hat{\Sigma}_{AG}}&\greyed{m_{G}^{2}-\hat{\Sigma}_{GG}}\\\
\end{pmatrix},$ (1) $\displaystyle\mathbf{M}_{\pm}^{2}$
$\displaystyle=\begin{pmatrix}m_{H^{\pm}}^{2}-\hat{\Sigma}_{H^{-}H^{+}}&\greyed{-\hat{\Sigma}_{H^{-}G^{+}}}\\\
\greyed{-\hat{\Sigma}_{G^{-}H^{+}}}&\greyed{m_{G^{\pm}}^{2}-\hat{\Sigma}_{G^{-}G^{+}}}\end{pmatrix}.$
(2)

where the $\hat{\Sigma}$ denote the renormalized self-energies. In the
following we shall neglect the Goldstone matrix elements (greyed out above)
since their numerical effect is tiny. Mixing with Goldstone bosons is taken
into account inside the loop diagrams and for a consistent renormalization, of
course. The physical masses are obtained as the real parts of the poles of the
propagator matrix

$\displaystyle\Delta(p^{2})=-{i\mkern 1.0mu}\left[p^{2}\mathrm{\rm
1\mskip-4.5mul}-\mathbf{M}^{2}\right]^{-1}.$ (3)

This is tantamount to finding the zeroes of the determinant of
$p^{2}\mathrm{\rm 1\mskip-4.5mul}-\mathbf{M}^{2}$, which is core functionality
of FeynHiggs.

For the cMSSM, the following contributions to Eq. (1) are available in
FeynHiggs: the one-loop self-energies with full $p^{2}$ dependence
Frank:2006yh , the leading momentum-independent two-loop
$\mathcal{O}(\alpha_{s}\alpha_{t})$ self-energies Heinemeyer:2007aq , and now
also the $\mathcal{O}(\alpha_{t}^{2})$ Hollik:2014bua  self-energies. The
subleading two-loop $\mathcal{O}(\alpha_{s}\alpha_{b})$,
$\mathcal{O}(\alpha_{t}\alpha_{b})$ parts are implemented only in the rMSSM
yet, and are interpolated in the phases when complex parameters are chosen
Frank:2006yh .

For the latter corrections the following components have to be computed:

  * •

The unrenormalized genuine two-loop self-energies $\Sigma_{hh}^{(2)}$,
$\Sigma_{hH}^{(2)}$, $\Sigma_{hA}^{(2)}$, $\Sigma_{HH}^{(2)}$,
$\Sigma_{HA}^{(2)}$, $\Sigma_{AA}^{(2)}$, $\Sigma_{H^{+}H^{-}}^{(2)}$ at
$p^{2}=0$ in $\mathcal{O}(\alpha_{t}^{2})$ approximation.

  * •

The one-loop diagrams with insertions of one-loop counterterms. For the
$\mathcal{O}{\left(\alpha_{t}^{2}\right)}$ corrections the couplings and
masses of the colored sector and chargino–neutralino sector are affected by
this subrenormalization.

  * •

The two-loop counterterms for these self-energies. The counterterms, the
imposed renormalization conditions, and all renormalization constants are
given in Ref. Hollik:2014bua .

  * •

The two-loop tadpoles $T_{h}^{(2)}$, $T_{H}^{(2)}$, $T_{A}^{(2)}$ in
$\mathcal{O}(\alpha_{t}^{2})$ approximation, which enter the two-loop
counterterms.

##  3 Approximations in the
$\mathcal{O}{\scalebox{1.2}[1.3]{(}\alpha_{t}^{2}\scalebox{1.2}[1.3]{)}}$
contributions

The dominant terms of the $\mathcal{O}{(\alpha_{t}^{2})}$ contributions are
enhanced by an additional factor $m_{t}^{2}$. The following approximations are
applied to the renormalized two-loop self-energies in Eq. (1) to yield only
this leading part.

###  3.1 Gaugeless limit

The gauge couplings $g_{1}$ and $g_{2}$ are set to zero, and also the strong
coupling $g_{\text{s}}$ is discarded. Consequently, also the gauge-boson
masses $M_{W}$ and $M_{Z}$ vanish, while the weak mixing angle
$\theta_{\text{w}}$ retains its original value. These choices cannot naively
be substituted in the original MSSM model file of FeynArts Hahn:2001rv  where
the couplings are not yet expressed in terms of the top-Yukawa coupling
$h_{t}=e\,m_{t}/(\sqrt{2}\,s_{\beta}\,s_{\text{w}}\,M_{W})$. Furthermore, the
non-zero ratios $\delta M_{W}/M_{W}$ and $\delta M_{Z}/M_{Z}$ appear in the
renormalization of $h_{t}$, while the gauge-boson-mass renormalization
constants $\delta M_{W}$ and $\delta M_{Z}$ themselves vanish. The technical
implementation of these issues is described in Sect. 4.

###  3.2 Vanishing external momentum

The external momentum is set to zero, i. e.the integrals that appear in the
two-loop self-energies and the renormalization constants are calculated at
$p^{2}=0$. As an immediate consequence all appearing two-loop integrals are
independent of $p^{2}$, leading to vacuum diagrams that are known
analytically.

###  3.3 Bottom mass equal to zero

For the $\mathcal{O}(\alpha_{t}^{2})$ contributions to form a supersymmetric
and gauge-invariant class of two-loop corrections to the Higgs-boson masses,
the bottom-quark mass $m_{b}$ has to be set to zero.

###  3.4 Consequences of the approximations

The approximations mentioned above modify the MSSM parameters that enter the
calculation as follows:

  * •

The Higgs- and Goldstone-boson masses at the tree level become
$m_{h}^{2}=m_{G}^{2}=m_{G^{\pm}}^{2}=0$ and
$m_{H}^{2}=m_{A}^{2}=m_{H^{\pm}}^{2}$.

  * •

The Higgs-sector mixing angles $\alpha$ and $\beta$ are related through
$\alpha=\beta-\pi/2$.

  * •

The stop masses are derived as the eigenvalues of

$\displaystyle\begin{pmatrix}m_{\tilde{t}_{\text{L}}}^{2}+m_{t}^{2}&m_{t}X_{t}^{*}\\\
m_{t}X_{t}&m_{\tilde{t}_{\text{R}}}^{2}+m_{t}^{2}\end{pmatrix},$ (4)

with $X_{t}=A_{t}-\mu^{*}/t_{\beta}$, where the soft-SUSY-breaking parameters
$m_{\tilde{t}_{\text{L,R}}}^{2}$ and $A_{t}$ enter.

  * •

Only one sbottom remains in the Feynman diagrams with
$m_{\tilde{b}}^{2}=m_{\tilde{t}_{\text{L}}}^{2}$. The sbottom mixing matrix is
set to the unity matrix.

  * •

Of the neutralinos and charginos only the original Higgsinos
$\tilde{\chi}_{3,4}^{0},\,\tilde{\chi}_{2}^{\pm}$ occur in the Feynman
diagrams and their degenerate masses are equal to $|\mu|$. The mixing matrices
$\mathbf{N}$ for the neutralinos and $\mathbf{U}$, $\mathbf{V}$ for the
charginos simplify to

$\displaystyle\mathbf{N}$
$\displaystyle=\begin{pmatrix}\begin{pmatrix}\mathrm{e}^{\frac{{i\mkern
1.0mu}}{2}\phi_{M_{1}}}&0\\\ 0&\mathrm{e}^{\frac{{i\mkern
1.0mu}}{2}\phi_{M_{2}}}\end{pmatrix}&\mathbf{0}\\\
\mathbf{0}&\frac{1}{\sqrt{2}}\mathrm{e}^{\frac{{i\mkern
1.0mu}}{2}\phi_{\mu}}\begin{pmatrix}1&-1\\\ {i\mkern 1.0mu}&{i\mkern
1.0mu}\end{pmatrix}\end{pmatrix},$ $\displaystyle\mathbf{U}$
$\displaystyle=\begin{pmatrix}\mathrm{e}^{{i\mkern 1.0mu}\phi_{M_{2}}}&0\\\
0&\mathrm{e}^{{i\mkern 1.0mu}\phi_{\mu}}\end{pmatrix},$
$\displaystyle\mathbf{V}$ $\displaystyle=\mathrm{\rm 1\mskip-4.5mul}\,.$ (5)

  * •

Many couplings are simplified. Technically this is taken care of by a modified
version of the FeynArts model file for the complex MSSM, generated as
described in Sect. 4. Explicit results can also be found in Ref.
Hollik:2014bua .

##  4 Implementation in FeynHiggs

The original calculation of the renormalized two-loop self-energies at
$\mathcal{O}(\alpha_{t}^{2})$ Hollik:2014bua  was overhauled entirely for
inclusion in FeynHiggs. The new code is included in the FeynHiggs tarball and,
while still specific to the corrections mentioned, is very adaptable and may
serve as a template for similar calculations, which is the reason we describe
it here in some detail.

The calculation has been divided into seven working steps implemented by the
scripts described in the following sections. Additionally, several packages
that are required for specific purposes are accessed during each step. A
makefile coordinates the entire calculation. The total running time is about
15–20 minutes and the final output is a highly optimized and very compact
Fortran code with a file size of about 350 kBytes. The result can directly be
moved into the FeynHiggs source tree.

The advantages of this setup are twofold: First, the compartmentalization
takes care of ‘concurrent’ software packages (which cannot be loaded into the
same Mathematica session, typically because of symbol conflicts). Likewise, if
an extension of this code required a different package for a particular step
(say, tensor reduction), only this one step would need to be modified, and the
additional package could be an arbitrary one, not just a Mathematica package,
since it is invoked from a shell script. Second, the makefile sees to it that
only the necessary parts of the calculation are redone in case of a change,
which was a particular time-saver during the debugging phase.

Some details of the software setup:

  * •

The code described below is included in FeynHiggs 2.11.0 or later, in the
gen/tlsp tree.

  * •

The shell scripts in the ‘scripts’ subdirectory perform the actual
calculation. They are executed directly from the command-line (i. e.not loaded
within Mathematica), usually by the makefile. Internally they run the
Mathematica Kernel as described in Ref. Hahn:2006ig .

  * •

The Mathematica programs in the ‘packages’ subdirectory provide functions for
specific purposes and are read by the scripts.

  * •

As in the FeynArts model file for the MSSM Hahn:2001rv  we choose the particle
labels

$h=\texttt{h0},\quad H=\texttt{HH},\quad A=\texttt{A0},\quad
H^{+}=\texttt{Hp},\quad H^{-}=\texttt{Hm}.$

  * •

Except for the last step (Code Generation) each script operates on a single
self-energy or tadpole only, which is referred to on the command-line and
stored under one of the following labels:

self-energies: $\displaystyle\texttt{h0h0},\ \texttt{h0HH},\ \texttt{h0A0},\
\texttt{HHHH},\ \texttt{HHA0},\texttt{A0A0},\ \texttt{HmHp},$ tadpoles:
$\displaystyle\texttt{h0},\ \texttt{HH},\ \texttt{A0}.$

  * •

Symbolic expressions are generally output to the ‘m’ directory, Fortran code
to the ‘f’ directory.

  * •

Input and output files are defined in statements ‘in=...’ and ‘out=...’ in the
first few lines of each script. For debugging, a script’s session dialog is
written to the output file with .log.gz appended. For example,
m/h0h0/6-comb.log.gz is the log of the run that made m/h0h0/6-comb.

###  4.1 Model-file preparation, Gaugeless Limit

Our calculation uses the FeynArts model file for the MSSM with counterterms
Fritzsche:2013fta . Since the calculation is carried out in the gaugeless
limit (cf. Sect. 3), it speeds up calculations to prepare a variant of the
model file where this limit is taken already at the level of the Feynman
rules. At this stage we also introduce $h_{t}$ (where
$h_{t}^{2}=4\pi\alpha_{t}$) and substitute $A_{t}$ by
$X_{t}+\mu^{*}/t_{\beta}$.

The original MSSMCT.mod model file contains both the counterterm vertices and
the definitions of the (one-loop) renormalization constants. We moved the
latter to model/MSSMCT.rc1 to make modification of the Feynman rules easier.

For each step we summarize usage as follows:

  * •

Command line: model/0-glmod

  * •

Main packages used: packages/Gaugeless.m, packages/XtSimplify.m

  * •

Input: model/MSSMCT.mod.in

  * •

Output: model/MSSMCT.mod, model/MSSMCTgl.mod

###  4.2 Step 1: Diagram generation

Diagrams are generated with FeynArts Hahn:2000kx  using the MSSMCTgl.mod model
file prepared above. Customized wrapper functions for the FeynArts functions
are provided by packages/FASettings.m.

Diagram selection in particular is simplified by filters defined in
FASettings.m. Defining a diagram filter which works for an arbitrary topology
is not entirely straightforward at the two-loop level but since only a few
topologies contribute we have chosen a per-topology approach where particle
choices can be made for each propagator.

For example, the two-loop self-energy diagrams for the neutral Higgs bosons
are picked with the following statement:

       sel[0][S[_] -> S[_]] = {
         t[3] && htb[6],
         t[3] && tb[6],
         t[3] && tb[6],
         t[3] && t[4] && htb[5],
         t[3] && htb[5|6],
         t[3] && htb[5],
         t[3] && t[5],
         t[5] && ht[3|4],
         t[3|4|5] && ht[3|4|5] }

The nine elements of the list correspond to the nine two-loop topologies of
Fig. 1. An entry like t[3] mandates that propagator #3 must carry a top or
stop. Likewise, htb[5|6] enforces a Higgs/Higgsino, top/stop, or
bottom/sbottom on propagators #5 and/or #6.

Figure 1: The 1PI two-point two-loop topologies with propagator numbers.

  * •

Command line: scripts/1-amps $\mathit{se}$ $\mathit{cto}$

$\mathit{se}$ is one of the self-energy/tadpole names given above (h0h0, HHHH,
etc.) and $\mathit{cto}$ is the counterterm order, i. e.0 for the virtual two-
loop diagrams and 1 for the one-loop plus counterterm diagrams.

  * •

Main packages used: packages/FASettings.m, packages/Gaugeless.m

  * •

Input: (none)

  * •

Output: m/$\mathit{se}$/1-amps.$\mathit{cto}$ (amplitude),
m/$\mathit{se}$/1-amps.$\mathit{cto}$.ps.gz (diagrams)

###  4.3 Step 2: Preparation for Tensor Reduction

In preparation for Step 3 (Tensor Reduction) several simplifications are
carried out: Firstly, the limit $p^{2}\to 0$ is applied. This limit is part of
the approximation employed here and moreover admits closed-form expressions
(logs and dilogs) for the two-loop integrals. Secondly, simplifications are
done in particular on the ubiquitous sfermion mixing matrix elements111The USf
actually carry two more indices for sfermion type and generation:
USf[$i$,$j$,$t$,$g$], which we suppress here.
$U_{ij}\equiv\texttt{USf[$i$,$j$]}$ $(i,j=1,2)$. Products are substituted as

UCSf[$i$,$j$]
$\displaystyle=\texttt{USf[$i$,$j$]}\,\texttt{USf[$i$,$j$]}^{*},$ (6)
UCSf[$i$,3] $\displaystyle=\texttt{USf[$i$,1]}\,\texttt{USf[$i$,2]}^{*},$
UCSf[3,3] $\displaystyle=\texttt{USf[1,1]}\,\texttt{USf[2,2]}^{*},$ (7)
UCSf[3,$j$] $\displaystyle=\texttt{USf[1,$j$]}\,\texttt{USf[2,$j$]}^{*},$
UCSf[3,4] $\displaystyle=\texttt{USf[1,2]}\,\texttt{USf[2,1]}^{*}$ (8)

which makes it a lot simpler to exploit unitarity: the relation
$UU^{\dagger}=\mathrm{\rm 1\mskip-4.5mul}$, for example, directly leads to
$\texttt{UCSf[1,1]}=\texttt{UCSf[2,2]}$ and
$\texttt{UCSf[1,2]}=\texttt{UCSf[2,1]}$. A subtle technical aspect is that not
all unitarity relations could even be formulated as definitions for the USf
(being too ‘deep’ for Mathematica’s UpValues) and applying them as rules would
not work reliably since Mathematica has no particular incentive to arrange the
USf in just the way to make the match.

After applying unitarity it is observed that almost all UCSf in the amplitude
appear in the form of only four linear combinations:

U2s1[$x$]
$\displaystyle=(\texttt{UCSf[1,3]}\,x^{*}+\texttt{UCSf[1,3]}^{*}\,x)/2\,,$ (9)
U2s2[$x$]
$\displaystyle=(\texttt{UCSf[1,3]}\,x^{*}-\texttt{UCSf[1,3]}^{*}\,x)/(2{i\mkern
1.0mu})\,,$ (10) U2c1[$x$]
$\displaystyle=(\texttt{UCSf[3,3]}\,x^{*}+\texttt{UCSf[3,4]}^{*}\,x)/2\,,$
(11) U2c2[$x$]
$\displaystyle=(\texttt{UCSf[3,3]}\,x^{*}-\texttt{UCSf[3,4]}^{*}\,x)/(2{i\mkern
1.0mu})\,.$ (12)

Substituting these combinations indeed shortens the amplitude by about 1/3.
Again, this choice is empirical and may well be specific to the corrections
computed here. All relations for substitution and simplification are contained
in packages/U2Simplify.m. Note in particular that, despite the suggestive
notation, Eqs. (11) and (12) are complex quantities since the two summands are
not conjugates of each other.

  * •

Command line: scripts/2-prep $\mathit{se}$ $\mathit{cto}$

$\mathit{se}$ and $\mathit{cto}$ as in Step 1.

  * •

Main packages used: packages/USfSimplify.m, packages/U2Simplify.m,
packages/SimplificationDefinitions.m

  * •

Input: m/$\mathit{se}$/1-amps.$\mathit{cto}$

  * •

Output: m/$\mathit{se}$/2-prep.$\mathit{cto}$

###  4.4 Step 3: Tensor Reduction

The tensor reduction of the loop integrals is performed with the packages
TwoCalc Weiglein:1993hd ; Weiglein:1995qs  (virtual two-loop diagrams,
$\mathit{cto}$ = 0) and FormCalc Hahn:1998yk  (one-loop plus counterterm
diagrams, $\mathit{cto}$ = 1).

One feature of this script is that, because TwoCalc and FormCalc cannot be
loaded into one Mathematica session due to symbol conflicts, the shell script
invokes not one but two different Mathematica sessions depending on the value
of $\mathit{cto}$.

  * •

Command line: scripts/3-calc $\mathit{se}$ $\mathit{cto}$

$\mathit{se}$ and $\mathit{cto}$ as in Step 1.

  * •

Main packages used: TwoCalc, FormCalc, packages/FCSettings.m,
packages/UseSimplePackage.m

  * •

Input: m/$\mathit{se}$/2-prep.$\mathit{cto}$

  * •

Output: m/$\mathit{se}$/3-calc.$\mathit{cto}$

###  4.5 Step 4: Simplification after Tensor Reduction

The tensor reduction is traditionally the step that increases the number of
terms in the amplitude most. Step 4 reduces the size of each individual
amplitude before they are combined in Step 6.

The key to decent performance comes from a ‘tag’ function222The DiagMark
function was originally developed for Refs. Heinemeyer:2003dq ;
Heinemeyer:2004yq . DiagMark[$m_{i}$] which the FeynArts wrapper functions of
Step 1 (packages/FASettings.m) inserted for each diagram, where $m_{i}$ are
the masses that run in the loop. This DiagMark function induces a
partitioning, i. e.the entire amplitude can be written as the sum of pieces,
each of which is multiplied by a different DiagMark. Typically only few
simplifications can be made across these pieces, hence one can restrict
application of a simplification function $\mathit{simp}$ to each piece, as in

> Collect[$\mathit{amp}$, _DiagMark, $\mathit{simp}$]

This is of course much faster than applying $\mathit{simp}$ to the entire
expression $\mathit{amp}$.

  * •

Command line: scripts/4-simp $\mathit{se}$ $\mathit{cto}$

$\mathit{se}$ and $\mathit{cto}$ as in Step 1.

  * •

Main packages used: packages/SimplificationDefinitions.m,
packages/FCSettings.m

  * •

Input: m/$\mathit{se}$/3-calc.$\mathit{cto}$

  * •

Output: m/$\mathit{se}$/4-simp.$\mathit{cto}$

###  4.6 Step 5: Calculation of the Renormalization Constants

The one-loop Renormalization Constants (RCs) are defined in MSSMCT.rc1 (see
Sect. 4.1), where the $\delta M_{V}^{2}$ ($V=W,Z$) require special attention,
however. As already pointed out in Sect. 3, they are proportional to
$M_{V}^{2}$ and would naively vanish in the gaugeless limit, while $\delta
M_{V}^{2}/M_{V}^{2}$ must not vanish. Technically, we accomplish this by
explicitly writing out the proportionality to $M_{V}^{2}$:
$\texttt{dMVsq1}\to\texttt{dMVsq1MV2*MV2}$, such that a possible MV2 in the
denominator of the coefficient of dMVsq1 can cancel.

For the two-loop RCs of the cMSSM no complete set of definitions is available
as a model file yet.333A complete list for the counterterms of the Higgs
sector in the complex MSSM can be found in Ref. Hollik:2014bua . We added the
necessary ones in the approximation relevant here in model/MSSMCT.rc2. For
better extendability we tried to preserve the mathematical notation as much as
possible, e. g.the two-loop wavefunction renormalization is defined via (the
mass counterterm is defined separately)

       dMHmat[Msq_, dZ1_, dZ2_, dM1_] :=
         1/2 (Msq.dZ2 + ConjugateTranspose[dZ2].Msq) +
         1/2 (ConjugateTranspose[dZ1].dM1 + ConjugateTranspose[dM1].dZ1) +
         1/4 ConjugateTranspose[dZ1].Msq.dZ1

which preserves the matrix notation of Eq. (2.30) of Hollik:2014bua .

The RCs are further expanded in the dimension parameter $\varepsilon=(D-4)/2$.
For example, the one-loop top-mass RC, dMf1[3,3], is decomposed as

       dMf1[3,3] = RC[-1, dMf1[-1,3,3]] + RC[0, dMf1[0,3,3]]

where the two terms on the r.h.s. represent the coefficients of
$\varepsilon^{-1}$ and $\varepsilon^{0}$, respectively. The loop integrals are
expanded with the ExpandDel function and powers of $\varepsilon$ are collected
with DelSeries, both defined in packages/ExpandDel.m. The notation allows to
easily extract the terms actually needed in the final result in Step 6.

The $\mathcal{O}(\varepsilon^{1})$ parts of the one-loop integrals often
cancel in the final result. This cancellation should (and in the present case
has been) checked, yet it needs some effort to see the coefficients of the
$\varepsilon^{1}$-terms actually vanish. Once this is established, however, it
is useful to set the $\varepsilon^{1}$-terms to zero manually, which is faster
and gives a more compact result. In packages/ExpandDel.m this is done through
rules defined in AssumeVanishingDel. To re-enable checking, set
AssumeVanishingDel = .

The actual calculation of the RCs is carried out with FormCalc, for each self-
energy and tadpole. So as not to re-compute RCs many times, they are cached in
the m/rc subdirectory. Take care that this cache is not controlled by the
makefile.

  * •

Command line: scripts/5-rc $\mathit{se}$

$\mathit{se}$ as in Step 1.

  * •

Main packages used: packages/RenConst.m, packages/ExpandDel.m

  * •

Input: m/$\mathit{se}$/4-simp.1

  * •

Output: m/$\mathit{se}$/5-rc, m/rc/*

###  4.7 Step 6: Combining Results

The final analytical step for each self-energy/tadpole joins the amplitudes
with the renormalization constants and, unless a debug flag is set, retains
only the coefficient of $\varepsilon^{0}$, i. e.the finite part. It also
applies a (mostly empirical) blend of simplification functions to get the
compact expression for code generation.

Two of these simplifications merit a closer look, for they affect also the
numerical precision of the final result. The difference and sum of the
sfermion mass squares not only appear frequently in the amplitude, they are
also computed as a by-product of diagonalization anyway, with effectively no
loss in precision even for degenerate sfermion masses. We thus substitute

> MSf2[3,$t$,$g$] = MSf2[2,$t$,$g$] - MSf2[1,$t$,$g$],
> MSf2[4,$t$,$g$] = MSf2[1,$t$,$g$] + MSf2[2,$t$,$g$].

Considering that in the gaugeless limit with vanishing bottom mass the
remaining sbottom mass square is equal to the stop-sector breaking parameter,
$m_{\tilde{b}}^{2}=m_{\tilde{t}_{\text{L}}}^{2}$, the difference of stop and
sbottom mass squares can be computed with better precision, too, which is
particularly relevant for the $H^{\pm}$ self-energy where this difference
appears in denominators. Hence we also introduce

> MSq2Diff[$i$,$j$] = MSf2[$i$,4,3] - MSf2[$j$,3,3].

  * •

Command line: scripts/6-comb $\mathit{se}$

$\mathit{se}$ as in Step 1.

  * •

Main packages used: packages/FinalSimp.m, packages/ExpandDel.m

  * •

Input: m/$\mathit{se}$/4-simp.0, m/$\mathit{se}$/4-simp.1, m/rc/*

  * •

Output: m/$\mathit{se}$/6-comb

###  4.8 Step 7: Generating Code

The last step writes out optimized Fortran code using FormCalc’s code-
generation functions Hahn:2010zi . Self-energies and tadpoles are arranged in
three groups, such that FeynHiggs can also perform a partial evaluation
depending on its flags (e. g.if $2\times 2$ mixing only is requested):

> 1\. A0A0, HmHp,
> 2\. h0h0, HHHH, h0HH, h0, HH,
> 3\. h0A0, HHA0, A0.

Abbreviations are introduced for both loop integrals and common
subexpressions, which shortens the entire expression significantly. Some
variables need to be disambiguated, for example the sfermion masses in the
gaugeless limit must not use the same identifiers as other sfermion masses in
FeynHiggs.

Lastly, the handwritten part of the code from the ‘static’ subdirectory is
copied to the Fortran output directory ‘f’, which is then complete and
replaces the entire src/TwoLoop/TLsp branch of the FeynHiggs source tree.

  * •

Command line: scripts/7-code

  * •

Main packages used: FormCalc, packages/FCSettings.m

  * •

Input: m/*/6-comb, static/*

  * •

Output: f/*

##  5 Conclusions

We have described the implementation of the $\mathcal{O}(\alpha_{t}^{2})$
Higgs-mass corrections for the complex MSSM in FeynHiggs. The original
calculation of Ref. Hollik:2014bua  was restructured entirely for this purpose
and the new code is included in the FeynHiggs source code, available at
http://feynhiggs.de. The calculation is split into seven tasks which were
explained in detail. While each step is specific to the corrections computed
here the code can be used as a template due to its modular structure and
adapted with little effort for similar calculations.

## Acknowledgments

We thank Stefano Di Vita, Thomas Hahn, Sven Heinemeyer, Heidi Rzehak, Pietro
Slavich, Dominik St ckinger, Alexander Voigt, and Georg Weiglein for helpful
discussions.

refs.1refs.1ReferencesReferencesrefs.1

## References

  * (1) G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 716 (2012) 1 [arXiv:1207.7214].
  * (2) S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 716 (2012) 30 [arXiv:1207.7235].
  * (3) G. Piacquadio [ATLAS Collaboration], J. Bendavid [CMS Collaboration], M. Dührssen [ATLAS+CMS Collaborations], talks given at the 50th Rencontres de Moriond, La Thuile, Aosta valley, Italy, March 2015.
  * (4) A. Pilaftsis, Phys. Rev. D 58 (1998) 096010 [hep-ph/9803297], Phys. Lett. B 435 (1998)88 [hep-ph/9805373].
  * (5) D.A. Demir, Phys. Rev. D 60 (1999) 055006 [hep-ph/9901389].
  * (6) A. Pilaftsis, C.E.M. Wagner, Nucl. Phys. B 553 (1999) 3 [hep-ph/9902371].
  * (7) S. Heinemeyer, Eur. Phys. J. C 22 (2001) 521 [hep-ph/0108059].
  * (8) J.A. Casas, J.R. Espinosa, M. Quiros, A. Riotto, Nucl. Phys. B 436 (1995) 3 [Erratum-ibid. B 439 (1995) 466] [hep-ph/9407389]. M.S. Carena, J.R. Espinosa, M. Quiros, C.E.M. Wagner, Phys. Lett. B 355 (1995) 209 [hep-ph/9504316].
  * (9) S. Heinemeyer, W. Hollik, G. Weiglein, Phys. Rev. D 58 (1998) 091701 [hep-ph/9803277], Phys. Lett. B 440 (1998) 296 [hep-ph/9807423].
  * (10) S. Heinemeyer, W. Hollik, G. Weiglein, Eur. Phys. J. C 9 (1999) 343 [hep-ph/9812472].
  * (11) S. Heinemeyer, W. Hollik, G. Weiglein, Phys. Lett. B 455 (1999) 179 [hep-ph/9903404]. M.S. Carena, H.E. Haber, S. Heinemeyer, W. Hollik, C.E.M. Wagner, G. Weiglein, Nucl. Phys. B 580 (2000) 29 [hep-ph/0001002].
  * (12) S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Eur. Phys. J. C 39 (2005) 465 [hep-ph/0411114].
  * (13) S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, W. Hollik, Eur. Phys. J. C 74 (2014) 2994 [arXiv:1404.7074].
  * (14) G. Degrassi, S. Di Vita, P. Slavich, Eur. Phys. J. C 75 (2015) 2, 61 [arXiv:1410.3432 [hep-ph]].
  * (15) S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, W. Hollik, arXiv:1505.03133 [hep-ph].
  * (16) R. Harlander, P. Kant, L. Mihaila, M. Steinhauser, Phys. Rev. Lett. 100 (2008) 191602; _ibid._ 101 (2008) 039901 [arXiv:0803.0672], JHEP 1008 (2010) 104 [arXiv:1005.5709].
  * (17) R.-J. Zhang, Phys. Lett. B 447 (1999) 89 [hep-ph/9808299]. J.R. Espinosa, R.-J. Zhang, Nucl. Phys. B 586 (2000) 3 [hep-ph/0003246], JHEP 0003 (2000) 026 [hep-ph/9912236]. J.R. Espinosa, I. Navarro, Nucl. Phys. B 615 (2001) 82 [hep-ph/0104047]. G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 611 (2001) 403 [hep-ph/0105096]. R. Hempfling, A.H. Hoang, Phys. Lett. B 331 (1994) 99 [hep-ph/9401219]. A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 643 (2002) 79 [hep-ph/0206101]. A. Dedes, G. Degrassi, P. Slavich, Nucl. Phys. B 672 (2003) 144 [hep-ph/0305127].
  * (18) J.R. Espinosa, R.-J. Zhang, Nucl. Phys. B 586 (2000) 3 [hep-ph/0003246].
  * (19) A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 631 (2002) 195 [hep-ph/0112177].
  * (20) G. Degrassi, S. Heinemeyer, W. Hollik, P. Slavich, G. Weiglein, Eur. Phys. J. C 28 (2003) 133 [hep-ph/0212020].
  * (21) S. Heinemeyer, W. Hollik, G. Weiglein, Phys. Rept. 425, 265 (2006) [hep-ph/0412214].
  * (22) B.C. Allanach, A. Djouadi, J. L. Kneur, W. Porod, P. Slavich, JHEP 0409 (2004) 044 [hep-ph/0406166].
  * (23) S.P. Martin, Phys. Rev. D 65 (2002) 116003 [hep-ph/0111209], Phys. Rev. D 66 (2002) 096001 [hep-ph/0206136], Phys. Rev. D 67 (2003) 095012 [hep-ph/0211366], Phys. Rev. D 68 (2003) 075002 [hep-ph/0307101], Phys. Rev. D 70 (2004) 016005 [hep-ph/0312092], Phys. Rev. D 71 (2005) 016012 [hep-ph/0405022], Phys. Rev. D 71 (2005) 116004 [hep-ph/0502168]. S.P. Martin, D.G. Robertson, Comput. Phys. Commun. 174 (2006) 133 [hep-ph/0501132].
  * (24) S.Y. Choi, M. Drees, J.S. Lee, Phys. Lett. B 481 (2000) 57 [hep-ph/0002287]. T. Ibrahim, P. Nath, Phys. Rev. D 63 (2001) 035009 [hep-ph/0008237], Phys. Rev. D 66 (2002) 015005 [hep-ph/0204092].
  * (25) M.S. Carena, J.R. Ellis, A. Pilaftsis, C.E.M. Wagner, Nucl. Phys. B 586 (2000) 92 [hep-ph/0003180].
  * (26) M. Frank, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, JHEP 0702 (2007) 047 [hep-ph/0611326].
  * (27) S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Phys. Lett. B 652 (2007) 300 [arXiv:0705.0746].
  * (28) W. Hollik, S. Paßehr, Phys. Lett. B 733, 144 (2014) [arXiv:1401.8275].
  * (29) W. Hollik, S. Paßehr, JHEP 1410 (2014) 171 [arXiv:1409.1687].
  * (30) S. Heinemeyer, W. Hollik, G. Weiglein, Comput. Phys. Commun. 124 (2000) 76 [hep-ph/9812320].
  * (31) T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Comput. Phys. Commun. 180 (2009) 1426.
  * (32) T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Nucl. Phys. Proc. Suppl. 205-206 (2010) 152 [arXiv:1007.0956].
  * (33) T. Hahn, C. Schappacher, Comput. Phys. Commun. 143 (2002) 54 [hep-ph/0105349].
  * (34) T. Hahn, J.I. Illana, Nucl. Phys. Proc. Suppl. 160 (2006) 101 [hep-ph/0607049].
  * (35) T. Fritzsche, T. Hahn, S. Heinemeyer, F. von der Pahlen, H. Rzehak, C. Schappacher, Comput. Phys. Commun. 185 (2014) 1529 [arXiv:1309.1692].
  * (36) T. Hahn, Comput. Phys. Commun. 140 (2001) 418 [hep-ph/0012260].
  * (37) G. Weiglein, R. Scharf, M. Böhm, Nucl. Phys. B 416 (1994) 606 [hep-ph/9310358].
  * (38) G. Weiglein, R. Mertig, R. Scharf, M. Böhm, in: New computing techniques in physics research II, 617–623, La Londe-les-Maures 1992.
  * (39) T. Hahn, M. Perez-Victoria, Comput. Phys. Commun. 118 (1999) 153 [hep-ph/9807565].
  * (40) S. Heinemeyer, D. Stöckinger, G. Weiglein, Nucl. Phys. B 690 (2004) 62 [hep-ph/0312264].
  * (41) S. Heinemeyer, D. Stöckinger, G. Weiglein, Nucl. Phys. B 699 (2004) 103 [hep-ph/0405255].
  * (42) T. Hahn, PoS ACAT 2010 (2010) 078 [arXiv:1006.2231].