text:
\subsection{Higgs sector at the tree level}

The Higgs sector of the MSSM consists out of two Higgs doublets $\mathcal{H}_1$ and $\mathcal{H}_2$ with opposite hypercharge. They are given by
\begin{align}\label{Eq:HiggsExpansion}
\mathcal{H}_1= \begin{pmatrix} v_1 + \frac{1}{\sqrt{2}}(\phi_1-i \chi_1)\\ -\phi_1^- \end{pmatrix},\hspace{2cm} \mathcal{H}_2= \begin{pmatrix} \phi_2^+\\ v_2 + \frac{1}{\sqrt{2}}(\phi_2+i \chi_2) \end{pmatrix},
\end{align}
where $v_1$ and $v_2$ are the vacuum expectation values (vevs) of the doublets. The ratio of the two vevs is denoted by
\begin{align}
\tan\beta = \frac{v_2}{v_1}.
\end{align}
The mass matrices of the component fields are diagonalized by the means of unitary transformations yielding the mass eigenstates. The transformations can be written as
\begin{align}\label{eq:HiggsMassStateTrafo}
\begin{pmatrix} h \\ H \\ A \\ G \end{pmatrix} = \textbf{U}_{n} \begin{pmatrix} \phi_1 \\ \phi_2 \\ \chi_1 \\ \chi_2 \end{pmatrix}, \hspace{0.5cm}
\begin{pmatrix} H^\pm \\ G^\pm \end{pmatrix} = \textbf{U}_{c} \begin{pmatrix} \phi_1^\pm \\ \phi_2^\pm \end{pmatrix}.
\end{align}
The unitary matrices $\textbf{U}_n$ and $\textbf{U}_c$ can be parametrized using the three mixing angles $\alpha$, $\beta_n$ and $\beta_c$,
\begin{align}
\textbf{U}_n = \begin{pmatrix}
-s_\alpha & c_\alpha & 0      & 0     \\
 c_\alpha & s_\alpha & 0      & 0     \\
0     & 0    & -s_{\beta_n} & c_{\beta_n} \\
0     & 0    &  c_{\beta_n} & s_{\beta_n}
\end{pmatrix}, \hspace{0.5cm}
\textbf{U}_c = \begin{pmatrix}
-s_{\beta_c} & c_{\beta_c} \\
c_{\beta_c} & s_{\beta_c}
\end{pmatrix}.
\end{align}
where we introduced the short-hand notation
\begin{align}
c_\gamma = \cos\gamma, \hspace{.3cm} s_\gamma = \sin\gamma, \hspace{.3cm} t_\gamma = \tan\gamma
\end{align}
for a generic angle $\gamma$ (the abbreviation for the tangent is used later).

By fixing
\begin{align}
\tan 2\alpha &= \frac{M_A^2+M_Z^2}{M_A^2-M_Z^2}\tan 2\beta,
\end{align}
where  $-\pi/2 \le \alpha \le 0$ is chosen by convention, and
\begin{align}
\beta_c &= \beta_n = \beta
\end{align}
the resulting mass matrices are diagonal. $M_Z$ is the mass of the $Z$~boson. $M_A$ denotes the mass of the $A$ boson, which is typically used as input parameter in the \ensuremath{{\cal CP}}-conserving case (i.e., if all parameters are real). In the \ensuremath{{\cal CP}}-violating case, the mass of the ${H^\pm}$ boson is used as input parameter.

The resulting tree-level masses are,
\begin{subequations}\label{eq:tree-level_masses}
\begin{align}
m_{H^\pm}^2  &= M_A^2 + M_W^2 \\
m_{h,H}^2 &= \frac{1}{2}\left(M_A^2 + M_Z^2 \mp \sqrt{(M_A^2 + M_Z^2)^2 - 4 M_A^2 M_Z^2 c_{2\beta}^2}\right), \\
m_G^2     &= m_{G^\pm}^2 = 0,
\end{align}
\end{subequations}
where $M_W$ is the mass of the $W$~boson.

The irreducible two-point vertex functions at lowest order, marked by the superscript ``0'', are diagonal and given by
\begin{subequations}
\begin{align}
\mathbf{\Gamma}^{(0)}_{hHAG}(p^2)     = i\left(p^2 \mathbf{1}_{4\times 4} - \mathbf{M}_{hHAG}\right),\\
\mathbf{\Gamma}^{(0)}_{{H^\pm}{G^\pm}}(p^2) = i\left(p^2 \mathbf{1}_{2\times 2} - \mathbf{M}_{{H^\pm}{G^\pm}}\right),
\end{align}
\end{subequations}
where $\mathbf{M}_{hHAG}$ and $\mathbf{M}_{{H^\pm}{G^\pm}}$ are the diagonal tree-level mass matrices with the entries given in \EqsSub{eq:tree-level_masses}.

\subsection{Renormalization of the Higgs sector}

At higher orders, loop corrections modify the simple tree-level expressions presented above,
\begin{subequations}
\begin{align}
i\left(p^2 \mathbf{1}_{4\times 4} - \mathbf{M}_{hHAG}\right) \rightarrow i\left(p^2 \mathbf{1}_{4\times 4} - \mathbf{M}_{hHAG} + \mathbf{\hat\Sigma}_{hHAG}(p^2)\right), \\
i\left(p^2 \mathbf{1}_{2\times 2} - \mathbf{M}_{{H^\pm}{G^\pm}}\right) \rightarrow i\left(p^2 \mathbf{1}_{2\times 2} - \mathbf{M}_{{H^\pm}{G^\pm}} + \mathbf{\hat\Sigma}_{{H^\pm}{G^\pm}}(p^2)\right),
\end{align}
\end{subequations}
where the matrices $\mathbf{\hat\Sigma}$ mark the renormalized loop corrections to the Higgs-boson propagators, called self-energies. Note that those matrices are not diagonal. Consequently, the loop corrections lead to mixing of the tree-level mass eigenstates.

For a description of the mass and tadpole renormalization, we refer to the discussions in e.g.~\cite{Frank:2006yh,Hollik:2014bua}. Here, we instead focus on the Higgs field renormalization.

\subsubsection{Higgs field renormalization}
\label{subsubsec:fieldren}

To render the results for physical observables UV finite, mass and tadpole renormalization is sufficient. If Green's functions involving external Higgs bosons are required to be UV finite, also the Higgs fields have to be renormalized.

We introduce the field renormalization in the original gauge eigenstate basis (up to the two-loop order),
\begin{align} \label{eq:introfieldren}
\begin{pmatrix}
\mathcal{H}_1 \\ \mathcal{H}_2
\end{pmatrix}
\rightarrow
\mathbf{Z_\mathcal{H}}
\begin{pmatrix}
\mathcal{H}_1 \\ \mathcal{H}_2
\end{pmatrix}
=
\begin{pmatrix}
1 + \frac{1}{2}\dZOL{11} + \frac{1}{2}\DZ{11} & \frac{1}{2}\dZOL{12} + \frac{1}{2}\DZ{12} \\
\frac{1}{2}\dZOL{21} + \frac{1}{2}\DZ{21} & 1 + \frac{1}{2}\dZOL{22} + \frac{1}{2}\DZ{22}
\end{pmatrix}
\begin{pmatrix}
\mathcal{H}_1 \\ \mathcal{H}_2
\end{pmatrix}
\end{align}
with
\begin{align}\label{eq:dZ_def}
\DZ{ij} = \dZTL{ij} - \frac{1}{4}\left(\dZOL{ij}\right)^2\hspace{.5cm} (i,j = 1,2).
\end{align}
We label the components of $\mathbf{Z_\mathcal{H}}$ by $Z_{ij}$ ($i,j=1,2$). Note that in contrast to e.g.~\cite{Hollik:2014bua}, we also introduced an off-diagonal field renormalization. This off-diagonal field renormalization is not required in order to cancel UV divergencies. We will, however, use this renormalization prescription in \Sec{sec:04_poledet} to prevent the introduction of unwanted higher-order terms through the determination of the Higgs-boson propagator poles. Note also that the field renormalization as introduced in \Eq{eq:introfieldren} goes beyond the one employed in \cite{Bahl:2018jom}: we do not assume $Z_{12} = Z_{21}$ here.

For simplification, we restrict us here to the case of an hermitian field renormalization matrix,
\begin{align}\label{eq:hermitianZH}
\mathbf{Z_\mathcal{H}} = \mathbf{Z_\mathcal{H}}^\dagger.
\end{align}
Consequently, $Z_{11}$ and $Z_{22}$ have to be real and $Z_{21}$ has to be the complex conjugated of $Z_{12}$.

The field renormalization constants in the mass eigenstate basis are obtained straightforwardly by the rotation introduced in \Eq{eq:HiggsMassStateTrafo},
\begin{align}
\begin{pmatrix}h \\ H \\ A \\ G \end{pmatrix} &\rightarrow
\mathbf{Z}_{hHAG}
\begin{pmatrix}h \\ H \\ A \\ G \end{pmatrix}, \hspace{.5cm}
\begin{pmatrix} H^\pm \\ G^\pm \end{pmatrix}\rightarrow
\mathbf{Z}_{H^\pm G^\pm}
\begin{pmatrix} H^\pm \\ G^\pm \end{pmatrix},
\end{align}
with
\begin{align}
\mathbf{Z}_{hHAG} &=
\begin{pmatrix}
1 + \frac{1}{2}\delta Z_{hh} &     \frac{1}{2}\delta Z_{hH} &     \frac{1}{2}\delta Z_{hA} &     \frac{1}{2}\delta Z_{hG} \\
    \frac{1}{2}\delta Z_{hH} & 1 + \frac{1}{2}\delta Z_{HH} &     \frac{1}{2}\delta Z_{HA} &     \frac{1}{2}\delta Z_{HG} \\
    \frac{1}{2}\delta Z_{hA} &     \frac{1}{2}\delta Z_{HA} & 1 + \frac{1}{2}\delta Z_{AA} &     \frac{1}{2}\delta Z_{AG} \\
    \frac{1}{2}\delta Z_{hG} &     \frac{1}{2}\delta Z_{HG} &     \frac{1}{2}\delta Z_{AG} & 1 + \frac{1}{2}\delta Z_{GG}
\end{pmatrix}, \label{eq:ZhHAG}\\
\mathbf{Z}_{H^\pm G^\pm} &=
\begin{pmatrix}
1 + \frac{1}{2}\delta Z_{H^\pm H^\pm} &     \frac{1}{2}\delta Z_{H^\pm G^\pm} \\
    \frac{1}{2}\delta Z_{G^\pm H^\pm} & 1 + \frac{1}{2}\delta Z_{G^\pm G^\pm}
\end{pmatrix},
\end{align}
where $\delta Z_{ij} = \dZOL{ij}+\dZTL{ij}$ ($i,j=h,H,A,G$). Explicit expressions for these one- and two-loop counterterms can be found in \App{app:07_ren}.

As a consequence of \Eq{eq:hermitianZH}, all field renormalization constants in the \ensuremath{{\cal CP}}-even subsector, involving $h$ and $H$, get contributions only from the real parts of $\mathbf{Z_\mathcal{H}}$. The same holds true for the $\ensuremath{{\cal CP}}$-odd subsector, involving $A$ and $G$, and the diagonal field renormalization constants in the charged Higgs sector. In contrast, the field renormalization constants connecting the $\ensuremath{{\cal CP}}$-even and the $\ensuremath{{\cal CP}}$-odd subsectors (e.g.\ $\delta Z_{hA}$) receive contributions only from the imaginary parts of $\mathbf{Z_\mathcal{H}}$. The off-diagonal field renormalization constants in the charged Higgs sector get contributions from both the real and imaginary parts of $\mathbf{Z_\mathcal{H}}$.

\subsubsection[Renormalization of \texorpdfstring{$\tan\beta$}{tanb}]{Renormalization of \texorpdfstring{$\boldsymbol{\tan}\boldsymbol{\beta}$}{tanb}}
\label{sec:tanb_ren}

The vevs $v_1$ and $v_2$ are renormalized according to
\begin{align}\label{eq:vevren}
\begin{pmatrix}
v_1 \\
v_2
\end{pmatrix}
\rightarrow&
\begin{pmatrix}
1 + \frac{1}{2}\dZOL{11} + \frac{1}{2}\DZ{11} & \frac{1}{2}\Re\left(\dZOL{12} + \DZ{12}\right) \\
\frac{1}{2}\Re\left(\dZOL{21} + \DZ{21}\right) & 1 + \frac{1}{2}\dZOL{22} + \frac{1}{2}\DZ{22}
\end{pmatrix}
\begin{pmatrix}
v_1 + \dL{1} v_1 + \dL{2} v_1 \\
v_2 + \dL{1} v_2 + \dL{2} v_2
\end{pmatrix}
\end{align}
The appearance of the field renormalizaton constants follows from the renormalization prescription in \Eq{eq:introfieldren}.

The difference of the one-loop vev counterterms is UV~finite~\cite{Chankowski:1992er,Dabelstein:1994hb}. Therefore, we set
\begin{align}\label{eq:dv1eqdv2OL}
\frac{\dL{1} v_1}{v_1} = \frac{\dL{1} v_2}{v_2}.
\end{align}
In the approximation of vanishing electroweak gauge couplings, which we apply at the two-loop level, also the difference of the two-loop vev counterterms is UV~finite~\cite{Sperling:2013xqa}. So, we set
\begin{align}\label{eq:dv1eqdv2TL}
\frac{\dL{2} v_1}{v_1}\bigg|_\text{gl} = \frac{\dL{2} v_2}{v_2}\bigg|_\text{gl},
\end{align}
where the subscript ``gl'' is used to indicate the limit of vanishing electroweak gauge couplings, which is also often called gaugeless limit.

Applying these simplifications, the resulting one- and two-loop counterterms of $\tan\beta$, introduced via the renormalization transformation
\begin{align}
t_\beta\rightarrow t_\beta + \dOLt_\beta + \dTLt_\beta,
\end{align}
are then given by
\begin{subequations}
\begin{align}
\dOLt_\beta ={}& \frac{1}{2}t_\beta\left(\dZOL{22} - \dZOL{11}\right) + \frac{1}{4}\left(1-t_\beta^2)(\dZOL{12}+\dZOL{21}\right), \\
\dTLt_\beta\big|_\text{gl} ={}&\bigg[
                 \frac{1}{2}t_\beta\left(\dZTL{22}-\dZTL{11}\right) +\frac{1}{4}\left(1-t_\beta^2\right)\left(\dZTL{12}+\dZTL{21}\right) \nonumber\\
&\hspace{.1cm} - \frac{1}{2}t_\beta\left(\dZOL{22}-\dZOL{11}\right)\left(\frac{1}{4}\left(\dZOL{22}-\dZOL{11}\right)+\dZOL{11}\right) \nonumber\\
&\hspace{.1cm} - \frac{1}{8}\left((1-2t_\beta^2)\dZOL{11} + t_\beta^2\dZOL{22}\right)\left(\dZOL{12}+\dZOL{21}\right) \nonumber\\
&\hspace{.1cm} - \frac{1}{16}(1-t_\beta^2)\left((1+t_\beta)\left((\dZOL{12})^2 + (\dZOL{21})^2\right) + 2t_\beta\dZOL{12}\dZOL{21}\right)
\bigg]_\text{gl}.
\end{align}
\end{subequations}

\subsubsection{Renormalization conditions for the Higgs fields}

The UV divergent part of the field renormalization constants is fixed by demanding that the derivatives of the Higgs self-energies with respect to the external momentum are UV finite,
\begin{subequations}
\begin{align}\label{eq:dZUV}
\dZOL{11}\Big|_\text{div} &= - \Re\left[\Sigma_{11}^{(1)\prime}\right]_\text{div}, \\
\dZOL{22}\Big|_\text{div} &= - \Re\left[\Sigma_{22}^{(1)\prime}\right]_\text{div}, \\
\dZOL{12}\Big|_\text{div} &=  \dZOL{21}\Big|_\text{div} = 0,
\end{align}
\end{subequations}
where we used the subscript ``div'' to denote that only the UV divergent part is taken into account. We do not specify the momentum at which the derivatives of the self-energies are evaluated, since their UV divergent part does not depend on the external momentum.

As noted already above, the off-diagonal field renormalization constants are not needed to render Green's functions with external Higgs bosons UV finite. In the addition to the UV divergent part, we also allow for a UV finite contribution to the field renormalization constants. This part will be specified in \Sec{sec:04_poledet}. Note that the specification of the field renormalization constants also fixes the renormalization and therefore definition of $\tan\beta$.

In the approximation of vanishing electroweak gauge couplings and vanishing external momentum, as applied at the two-loop level, all two-loop field renormalization constants drop out of the calculation~\cite{Hollik:2014bua}. They also do not appear at higher-orders, since we work in the approximation of vanishing external momentum at the two-loop level. Therefore, we do not impose any condition on them in this work.

\subsubsection{Renormalized self-energies}

At the one-loop level, the renormalized self-energies are given in terms of the unrenormalized self-energies and counterterms by
\begin{subequations}
\begin{align}\label{eq:ren_seOL}
\hat\Sigma^{(1)}_{hh}(p^2) &= \Sigma^{(1)}_{hh}(p^2) - \dmOL{h}  + \frac{1}{2}\dZOL{hh}\left(p^2-m_h^2\right),\\
\hat\Sigma^{(1)}_{HH}(p^2) &= \Sigma^{(1)}_{HH}(p^2) - \dmOL{H}  + \frac{1}{2}\dZOL{HH}\left(p^2-m_H^2\right),\\
\hat\Sigma^{(1)}_{AA}(p^2) &= \Sigma^{(1)}_{AA}(p^2) - \dmOL{A}  + \frac{1}{2}\dZOL{AA}\left(p^2-M_A^2\right),\\
\hat\Sigma^{(1)}_{hH}(p^2) &= \Sigma^{(1)}_{hH}(p^2) - \dmOL{hH} + \frac{1}{2}\dZOL{hH}\left(p^2-\frac{1}{2}(m_h^2+m_H^2)\right),\\
\hat\Sigma^{(1)}_{hA}(p^2) &= \Sigma^{(1)}_{hA}(p^2) - \dmOL{hA} + \frac{1}{2}\dZOL{hA}\left(p^2-\frac{1}{2}(m_h^2+M_A^2)\right),\\
\hat\Sigma^{(1)}_{HA}(p^2) &= \Sigma^{(1)}_{HA}(p^2) - \dmOL{HA} + \frac{1}{2}\dZOL{HA}\left(p^2-\frac{1}{2}(m_H^2+M_A^2)\right), \\
\hat\Sigma^{(1)}_{{H^\pm}\Hpm}(p^2) &= \Sigma^{(1)}_{{H^\pm}\Hpm}(p^2) - \dmOL{{H^\pm}}  + \dZOL{{H^\pm}\Hpm}\left(p^2-m_{{H^\pm}}^2\right),
\end{align}
\end{subequations}
where the various $\dL{1} m^2$ denote the one-loop mass counterterms, which are given explicitly in \App{app:07_ren}. The renormalized self-energies involving Goldstone bosons are obtained analogously.

The renormalized two-loop self-energies are given by
\begin{align}\label{eq:renTLse}
\hat\Sigma^{(2)}_{ij}(0) &= \Sigma^{(2)}_{ij}(0) - \delta^{(2)}m^\mathbf{Z}_{ij},\hspace{.5cm}(i,j = h, H, A, G, {H^\pm}, {G^\pm}),
\end{align}
where $\delta^{(2)}m^\mathbf{Z}_{ij}$ is a combination of mass counterterms and field renormalization constants. Explicit expressions are listed in \App{app:07_ren}. Note that the unrenormalized two-loop self-energies contain subloop renormalization contributions.

\subsection{Numerical pole determination with {\ensuremath{\overline{\text{DR}}}}\xspace field renormalization}
\label{subsec:num_poledet}

The most direct approach to determine the Higgs-boson propagator poles is to solve \EqsSub{eq:poleeqs} numerically. This can be done employing different schemes for the Higgs field renormalization. First, we discuss the case of {\ensuremath{\overline{\text{DR}}}}\xspace Higgs field renormalization, which was used as default scheme in \mbox{{\tt FeynHiggs}}\xspace until version \texttt{2.14.2}. This means that no UV finite parts are added to the UV divergent parts as specified in \EqsSub{eq:dZUV}.

In~\cite{Bahl:2017aev}, this procedure was closely investigated in the limit of $M_A \gg M_Z$. In this limit, \mbox{$h$-$H$}~mixing effects in the \ensuremath{{\cal CP}}-even Higgs sector are suppressed by powers of $M_A$ and the mass of the $h$ boson can be obtained by solving the simpler equation
\begin{align}
p^2 - m_h^2 + \hat\Sigma_{hh}(p^2)= 0
\end{align}
up to corrections suppressed by $M_A$. Solving this equation up to the two-loop order yields
\begin{align}\label{Mh2_TL_DecouplingLimit_Eq}
M_h^2 ={}& m_h^2 - \hat\Sigma_{hh}^{(1)}(m_h^2) - \hat\Sigma_{hh}^{(2)}(0)\big|_\text{gl} + \hat\Sigma_{hh}^{(1)\prime}(m_h^2) \hat\Sigma_{hh}^{(1)}(m_h^2) = \nonumber\\
={}& m_h^2 - \hat\Sigma_{hh}^{(1)}(m_h^2) - \hat\Sigma_{hh}^{(2)}(0)\big|_\text{gl} \nonumber\\
& + \left(\hat\Sigma_{hh}^{{\text{heavy}},(1)\prime}(m_h^2)+\hat\Sigma_{hh}^{{\text{light}},(1)\prime}(m_h^2)\right) \hat\Sigma_{hh}^{(1)}(m_h^2),
\end{align}
where the number in the superscript is used to indicate the loop order. Again, the subscript ``gl'' is used to indicate the approximation of vanishing electroweak gauge couplings at the two-loop level. In the last step, we split up the derivative of the $hh$ self-energy into a part containing all contributions involving heavy particles (omitting terms suppressed by the heavy mass scale), labelled by the superscript ``heavy'' particles, and a part containing all contributions from light particles (including terms suppressed by a heavy mass scale), labelled by the superscript ``light''. Note that, in order to simplify the notation in this Section, the labels have a slightly different meaning in this Section as compared to \Sec{sec:03_p2heavy}. Here, the light contribution also includes terms suppressed by a large mass, whereas these terms are omitted in the heavy contribution. In~\cite{Bahl:2017aev}, a scenario with all non-SM particles sharing one common mass scale was assumed. Therefore, $\hat\Sigma_{hh}^{{\text{light}}}$ contains contributions from SM particles as well as terms suppressed by a heavy non-SM mass.

The analysis in~\cite{Bahl:2017aev} showed that the term labelled by ``heavy'' in \Eq{Mh2_TL_DecouplingLimit_Eq} is cancelled by parts of the subloop renormalization contained in the two-loop self-energy. The discussion in \Sec{sec:03_p2heavy} clearly shows the reason for this cancellation: The momentum dependence of heavy contributions to two-point vertex functions, i.e.\ the term proportional to $\hat\Sigma_{hh}^{{\text{heavy}},(1)\prime}(m_h^2)$ in \Eq{Mh2_TL_DecouplingLimit_Eq}, could be absorbed into a field renormalization constant and therefore drops out order by order in the perturbative expansion.\footnote{Note that we included terms suppressed by a heavy mass, which can contain additional momentum-dependent terms, into the ``light'' contribution in this Section.} If the Higgs field renormalization is fixed in the {\ensuremath{\overline{\text{DR}}}}\xspace scheme, it is therefore important for a proper cancellation to ensure that all corrections at the two-loop level are included at the same level of accuracy and that no higher order terms proportional to $\hat\Sigma_{hh}^{{\text{heavy}},(1)\prime}(m_h^2)$ are induced. Determining the Higgs-boson propagator poles numerically, however, does not allow to easily control the inclusion or exclusion of certain terms.

Concerning the corrections included in \mbox{{\tt FeynHiggs}}\xspace, this means that in order to achieve a complete cancellation at the two-loop level only \order{\al_t^2,\al_t\al_b,\al_b^2} corrections should be taken into account for the term $\hat\Sigma_{hh}^{{\text{heavy}},(1)\prime}(m_h^2) \hat\Sigma_{hh}^{(1)}(m_h^2)$, since the two-loop self-energy is included up to this order.

\subsection{Fixed-order pole determination}
\label{subsec:fo_poledet}

As an approach to solve this problem, it was suggested in \cite{Bahl:2017aev} to simply employ \Eq{Mh2_TL_DecouplingLimit_Eq} and to ensure that the last term in the equation is computed at the same level of accuracy as the two-loop self-energy. This simple procedure is, however, only applicable in case of high $M_A$ for which \mbox{$h$-$H$}~mixing effects are negligible. Here, we work out a method based upon the same idea -- solving \EqsSub{eq:poleeqs} at a fixed-order -- applicable also for low $M_A$. This method was used in the \mbox{{\tt FeynHiggs}}\xspace versions~\texttt{2.14.0} through~\texttt{2.14.2}.

A proper cancellation is guaranteed by including all two-loop corrections at the same level of accuracy. In case of low $M_A$, we cannot neglect the off-diagonal elements of the Higgs two-point vertex function, since their inclusion is crucial for a precise prediction. Instead, we consider each element of the two-point vertex function separately.

The complex pole corresponding to the $i$-th boson ($\mbox{i = 1,2,3}$) is given at the one-loop level by
\begin{align}\label{eq:oneloopsol}
\left(\mathcal{M}_{i}^{(1)}\right)^2 = m_{i}^2 - \hat\Sigma_{ii}^{(1)}(m_{i}^2),
\end{align}
where the superscript indicates the loop-order. We choose the tree-level masses to be
\begin{align}
m_{1}^2 &= m_h^2,\quad m_{2}^2 = m_H^2,\quad m_{3}^2 = M_A^2
\end{align}
and the one-loop self-energies to be
\begin{align}
\hat\Sigma_{{1}{1}}^{(1)} &= \hat\Sigma_{hh}^{(1)},\quad \hat\Sigma_{{2}{2}}^{(1)} = \hat\Sigma_{HH}^{(1)},\quad \hat\Sigma_{{3}{3}}^{(1)} = \hat\Sigma_{AA}^{(1)}.
\end{align}
Expanding each of the self-energies appearing in the two-point vertex function matrix (see \Eq{eq:propmatrix}) around the one-loop solution of \Eq{eq:oneloopsol} up to the two-loop level yields the expanded two-point vertex function matrix $\hat\Gamma_{hHA}^{i-\text{exp}}$ with the matrix elements given by
\begin{align}
 \left(\hat\Gamma_{hHA}^{i-\text{exp}}(p^2)\right)_{jk} &= D_{jk}(p^2) + \hat\Sigma_{jk}^{(1)}(m_{i}^2) + \hat\Sigma_{jk}^{(2)}(0)\big|_\text{gl}  - \left[\hat\Sigma_{jk}^{(1)\prime}(m_{i}^2)\hat\Sigma_{{i}{i}}^{(1)}(m_{i}^2)\right]_\text{gl}
\end{align}
with $j,k = 1,2,3$ and
\begin{align}
D_{jk}(p^2) = (p^2 - m_j^2)\delta_{jk}.
\end{align}
The subscript of the reducible terms $\hat\Sigma^{(1)\prime}\hat\Sigma^{(1)}$, ``gl'', indicates that we calculate these in the gaugeless limit (i.e.\ vanishing electroweak gauge couplings corresponding to the approximation used at the two-loop level). Consequently, these terms are of \order{\al_t^2,\al_t\al_b,\al_b^2} and therefore of the same order as the corresponding parts of the two-loop self-energies $\hat\Sigma^{(2)}$. Thereby, it is guaranteed that cancellations between the various two-loop contributions are correctly taken into account.

The next step is to find the zeroes of $\hat\Gamma_{hHA}^{i-\text{exp}}(p^2)$. We pick the pole with the largest overlap of the associated physical state with the tree-level mass eigenstate. This pole corresponds to the one-loop solution given in \Eq{eq:oneloopsol}. The physical mass squared is then given by taking the real part of this pole. The mass of the ${H^\pm}$ boson is determined accordingly. This procedure successfully avoids inducing higher order terms which would cancel in a more complete calculation.

However, as we will argue in \Sec{subsec:poledet_comp} and explicitly show in \Sec{sec:05_results}, the approach presented in this Section yields unsatisfactory results close to crossings points where two mixing Higgs bosons are almost mass degenerate. This issue does not appear if the propagator poles are determined numerically.

\subsection{Numerical pole determination with heavy-OS field renormalization}
\label{subsec:hOS_poledet}

Last, we present a third approach to determine the Higgs-boson propagator poles avoiding the problems of the approaches presented above. As in \Sec{subsec:num_poledet}, the Higgs-boson propagator poles are determined numerically. However, to avoid the problem with uncancelled terms beyond the order of the fixed-order calculation, we do not employ the {\ensuremath{\overline{\text{DR}}}}\xspace scheme for the Higgs field renormalization.

As shown in \Sec{sec:03_p2heavy}, the Higgs field renormalization can be used to absorb the momentum dependence of contributions from heavy particles. As discussed in \Sec{subsec:num_poledet}, their momentum dependence is the origin of the uncancelled higher-order terms. Consequently, the issue is resolved by fixing the Higgs field renormalization according to the heavy-OS scheme (as presented in \Sec{sec:03_p2heavy}). This means that the one-loop Higgs field renormalization constants (introduced in \Sec{sec:02_MSSMintro}) are fixed by
\begin{align}\label{eq:MSSMhOS}
\delta^{(1)}Z_{ij}\Big|_\text{hOS} &= \delta^{(1)}Z_{ij}\Big|_{\ensuremath{\overline{\text{DR}}}}\xspace - \Sigma_{ij}^{{\text{heavy}},(1)\prime}\Big|_\text{fin} - \gamma_{ij}^{{\text{light}},(1)}\ln\frac{\mu_R^2}{Q_\text{pd}^2}\hspace{.3cm}(i,j=h,H,A).
\end{align}
Note that we did not specify at which external momentum the second term is evaluated, since it does not depend on the external momentum. In addition to absorbing the momentum dependence of the heavy contributions, we also add a term proportional to the Higgs anomalous dimension involving light fields only,\footnote{Terms suppressed by a heavy scale do not contribute to the Higgs anomalous dimensions.}
\begin{align}\label{eq:anomal_dim_light}
\gamma_{ij}^{{\text{light}},(1)} = \frac{\partial}{\partial \ln \mu_R^2} \Sigma_{ij}^{{\text{light}},(1)\prime}\Big|_{\mu_R = Q_\text{pd}},
\end{align}
which is used to absorb the dependence of Higgs self-energies' derivatives on the renormalization scale $\mu_R$ (see end of \Sec{sec:03_p2heavy}). In this way, we separate the scale dependence associated with the pole determination from the renormalization scale dependence of the fixed-order corrections. The corresponding reference scale at which the poles are determined is called~$Q_\text{pd}$. In order to avoid inducing large logarithms via the pole determination we choose it equal to the OS top mass, $Q_\text{pd}=M_t$.

Explicit expressions for the one-loop renormalization constants in the gauge eigenstate basis are listed in \App{app:09_finfieldren}.

Since the fixed-order calculation used by default in \mbox{{\tt FeynHiggs}}\xspace applies the approximation of vanishing external momentum at the two-loop level, no finite pieces have to be added to the two-loop field renormalization constants. Going beyond this approximation (see~\cite{Borowka:2014wla,Borowka:2015ura,Borowka:2018anu}) implies that also the two-loop field renormalization constants should be fixed in the heavy-OS scheme.

Having calculated the Higgs self-energies using the field renormalization scheme specified in \Eq{eq:MSSMhOS}, the Higgs-boson propagator poles are then determined by solving \EqsSub{eq:poleeqs} numerically.

Note that employing the full OS scheme, meaning that the full derivatives of the Higgs self-energies are renormalized to be zero, is not a good alternative. In the OS scheme, unphysical threshold effects coming from the ``light'' contributions are induced (see e.g.~\cite{Frank:2002qf}).

It is important to note that the renormalization scheme choice made in \Eq{eq:MSSMhOS} affects the definition of $\tan\beta$ (see \Sec{sec:02_MSSMintro}). In particular $\tan\beta$ is not anymore a MSSM quantity but rather a Two-Higgs-Doublet Model (THDM) one. This is explained in detail in \cite{Bahl:2018jom}. By choosing the field renormalization as in \Eq{eq:MSSMhOS}, $\tan\beta$ is defined at the scale $Q_\text{pd}$, which is set to the OS top mass. In situations where the THDM scale $M_A$ is much larger than the top-quark mass, this definition is not well-suited. Furthermore, it would be preferable to have an MSSM $\tan\beta$ as input for a MSSM calculation. We present two methods to address these shortcomings in \App{app:08_tbrepara}.

In addition to the Higgs poles, also the $\mathbf{Z}$-matrix relating the tree-level mass eigenstates and the external physical states in scattering amplitudes is calculated~\cite{Chankowski:1992er,Heinemeyer:2001iy,Frank:2006yh,
Fuchs:2016swt,Fuchs:2017wkq,Domingo:2017rhb}. This matrix depends on the choice of the Higgs field renormalization. For other calculations (e.g.~decay rates or productions cross-sections), the heavy-OS scheme is, however, not needed to avoid uncancelled higher order terms, since these processes involve more than two external legs. In order to ease the use of the output $\mathbf{Z}$-matrix, we therefore discuss how it can be transferred to the {\ensuremath{\overline{\text{DR}}}}\xspace scheme, which is more widely used, in \App{app:10_Zmatrix}.

\subsection{Comparison of the different methods}
\label{subsec:poledet_comp}

To further illustrate the difference between the various approaches, we investigate the corrections to the mass of the $h$ boson up to the two-loop level.

In the first approach, numerical pole determination with {\ensuremath{\overline{\text{DR}}}}\xspace Higgs field renormalization, they are given by
\begin{align}\label{eq:Mh_pd1}
\left(M_h^2\right)_\text{num,{\ensuremath{\overline{\text{DR}}}}\xspace} ={}& m_h^2 - \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} - \hat\Sigma_{hh}^{(2)}(0)\big|_{\text{gl},\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} + \hat\Sigma_{hh}^{(1)\prime}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} \nonumber\\
&+ \frac{\left(\hat\Sigma_{hH}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace}\right)^2}{m_h^2 - m_H^2} + \frac{\left(\hat\Sigma_{hA}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace}\right)^2}{m_h^2 - M_A^2},
\end{align}
where we added the subscript ``$\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace$'' to clarify that the field renormalization constants are fixed in the {\ensuremath{\overline{\text{DR}}}}\xspace scheme.

In the fixed-order approach, we obtain
\begin{align}\label{eq:Mh_pd2}
\left(M_h^2\right)_\text{FO} ={}& m_h^2 - \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} - \hat\Sigma_{hh}^{(2)}(0)\big|_{\text{gl},\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} + \hat\Sigma_{hh}^{(1)\prime}(m_h^2)\big|_{\text{gl},\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\text{gl},\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace} \nonumber\\
&+ \frac{\left(\hat\Sigma_{hH}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace}\right)^2}{m_h^2 - m_H^2} + \frac{\left(\hat\Sigma_{hA}^{(1)}(m_h^2)\big|_{\delta Z|_{\ensuremath{\overline{\text{DR}}}}\xspace}\right)^2}{m_h^2 - M_A^2}.
\end{align}
We see that the fixed-order pole determination avoids including the terms identified in \Sec{subsec:num_poledet}, which would cancel in a more complete calculation.

Using numerical pole determination with OS pole determination, we get
\begin{align}\label{eq:Mh_pd3}
\left(M_h^2\right)_\text{num,hOS} ={}& m_h^2 - \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\delta Z_\text{hOS}} - \hat\Sigma_{hh}^{(2)}(0)\big|_{\text{gl},\delta Z_\text{hOS}} + \hat\Sigma_{hh}^{(1)\prime}(m_h^2)\big|_{\delta Z_\text{hOS}} \hat\Sigma_{hh}^{(1)}(m_h^2)\big|_{\delta Z_\text{hOS}} \nonumber\\
&+ \frac{\left(\hat\Sigma_{hH}^{(1)}(m_h^2)\big|_{\delta Z_\text{hOS}}\right)^2}{m_h^2 - m_H^2} + \frac{\left(\hat\Sigma_{hA}^{(1)}(m_h^2)\big|_{\delta Z_\text{hOS}}\right)^2}{m_h^2 - M_A^2},
\end{align}
where we added the subscript ``$\delta Z_\text{hOS}$'' to clarify that the field renormalization constants are fixed in the heavy-OS scheme. Again, terms which would cancel in a more complete calculation are not included.

The difference between the numerical pole determination employing heavy-OS Higgs field renormalization and the fixed-order pole determination is given by\footnote{Note that $\hat\Sigma_{hh}^{(1)}(m_h^2)$ does not depend explicitly on the scheme chosen for the Higgs field renormalization. Assuming the same definition of $\tan\beta$, $\hat\Sigma_{hh}^{(1)}(m_h^2)$ is therefore equal for {\ensuremath{\overline{\text{DR}}}}\xspace Higgs field renormalization and for heavy-OS Higgs field renormalization.}
\begin{align}\label{eq:diffpd12}
\left(M_h^2\right)_\text{num,hOS} - \left(M_h^2\right)_\text{FO} ={}& \Sigma_{hh}^{{\text{light}},(1)\prime}(m_h^2)\big|_\text{fin}\cdot\hat\Sigma_{hh}^{(1)}(m_h^2) \nonumber\\
&- \left[\Sigma_{hh}^{{\text{light}},(1)\prime}(m_h^2)\big|_\text{fin}\cdot\hat\Sigma_{hh}^{(1)}(m_h^2)\right]_\text{gl} .
\end{align}
We see that the numerical pole determination employing heavy-OS renormalization includes additional terms. The term containing the light contribution to the derivative of the $hh$~self-energy and the $hh$~self-energy is included fully and not only in the limit of vanishing electroweak gauge couplings. The electroweak next-to-leading logarithms which are contained in this term are of the same order as some of the logarithms resummed by the means of the EFT calculation, which are added on top of the fixed-order result. Therefore, this term, which was missed in \cite{Bahl:2017aev}, should be included fully and not only in the limit of vanishing electroweak gauge couplings.

\medskip

In addition to this differences at the two-loop level, there are also differences at higher orders. In the case of fixed-order pole determination, unphysical discontinuities in the prediction of the Higgs boson masses close to crossings points, where two mixing Higgs bosons are almost mass degenerate, appear (numerical results displaying this behavior will be shown in \Sec{sec:05_results}). At these crossing points, a specific property of the Higgs-boson propagator matrix becomes important: Every complex pole is not only the solution of \EqsSub{eq:poleeqs} but also a pole of every element of the Higgs-boson propagator matrix~\cite{Fuchs:2016swt}. This property is fulfilled in case of numerical pole determination. In case of the fixed-order pole determination, this property, however, holds only up to the level of the truncation, i.e., the two-loop level. The violation at higher orders is the origin of the discontinuities in the Higgs mass predictions.

\subsection{One-loop renormalization}

The one-loop field renormalization constants of the neutral Higgs bosons in the mass basis are given by
\begin{subequations}
\begin{align}
\dZOL{hh} &= s_\alpha^2 \dZOL{11} - s_\alphac_\alpha \left(\dZOL{12} + \dZOL{21}\right) + c_\alpha^2 \dZOL{22}, \\
\dZOL{HH} &= c_\alpha^2 \dZOL{11} + s_\alphac_\alpha \left(\dZOL{12} + \dZOL{21}\right) + s_\alpha^2 \dZOL{22}, \\
\dZOL{hH} &= s_\alphac_\alpha \left(\dZOL{22} - \dZOL{11}\right) + \frac{1}{2}c_{2\alpha}\left(\dZOL{12}+\dZOL{21}\right), \\
\dZOL{AA} &= s_{\beta_n}^2 \dZOL{11} - s_{\beta_n}c_{\beta_n} \left(\dZOL{12} + \dZOL{21}\right) + c_{\beta_n}^2 \dZOL{22}, \\
\dZOL{GG} &= c_{\beta_n}^2 \dZOL{11} + s_{\beta_n}c_{\beta_n} \left(\dZOL{12} + \dZOL{21}\right) + s_{\beta_n}^2 \dZOL{22}, \\
\dZOL{AG} &= s_{\beta_n}c_{\beta_n} \left(\dZOL{22} - \dZOL{11}\right) + \frac{1}{2}c_{2\beta_n}\left(\dZOL{12}+\dZOL{21}\right), \\
\dZOL{hA} &= \frac{i}{2}s_{\beta_n-\alpha}\left(\dZOL{12}-\dZOL{21}\right), \\
\dZOL{hG} &= -\frac{i}{2}c_{\beta_n-\alpha}\left(\dZOL{12}-\dZOL{21}\right), \\
\dZOL{HA} &= \frac{i}{2}c_{\beta_n-\alpha}\left(\dZOL{12}-\dZOL{21}\right), \\
\dZOL{HG} &= \frac{i}{2}s_{\beta_n-\alpha}\left(\dZOL{12}-\dZOL{21}\right).
\end{align}
\end{subequations}
The corresponding constants in the charged Higgs sector read
\begin{subequations}
\begin{align}
\dZOL{{H^\pm}\Hpm} &= s_{\beta_c}^2 \dZOL{11} - s_{\beta_c}c_{\beta_c} \left(\dZOL{12} + \dZOL{21}\right) + c_{\beta_c}^2 \dZOL{22}, \\
\dZOL{{G^\pm}\Gpm} &= c_{\beta_c}^2 \dZOL{11} + s_{\beta_c}c_{\beta_c} \left(\dZOL{12} + \dZOL{21}\right) + s_{\beta_c}^2 \dZOL{22}, \\
\dZOL{{H^\pm}{G^\pm}} &= s_{\beta_c}c_{\beta_c} \left(\dZOL{22} - \dZOL{11}\right) + c_{\beta_c}^2\dZOL{21} - s_{\beta_c}^2\dZOL{12}, \\
\dZOL{{G^\pm}{H^\pm}} &= s_{\beta_c}c_{\beta_c} \left(\dZOL{22} - \dZOL{11}\right) + c_{\beta_c}^2\dZOL{12} - s_{\beta_c}^2\dZOL{21}.
\end{align}
\end{subequations}
The field renormalization counterterms in the gaugeless limit can be obtained by setting $\alpha = \beta - \pi/2$ in the expressions listed above.

The one-loop mass counterterms of the neutral Higgs sector are given by
\begin{subequations}
\begin{align}
\dL{1} m_h^2 ={}& \dL{1} M_A^2 c_{\alpha-\beta}^2 + \dL{1} M_Z^2 s_{\alpha+\beta}^2 \nonumber\\
&+ \frac{e}{2 M_W s_W c_W}\left(\dL{1} T_H c_{\alpha-\beta}s_{\alpha-\beta}^2 + \dL{1} T_h s_{\alpha-\beta}(1+c_{\alpha-\beta}^2)\right) \nonumber\\
&+ c_\beta^2\dOLt_\beta (M_Z^2 s_{2(\alpha+\beta)} + M_A^2 s_{2(\alpha-\beta)}) \nonumber\\
&+ c_\alpha^2\dOLt_\alpha (M_Z^2 s_{2(\alpha+\beta)}-M_A^2 s_{2(\alpha-\beta)}), \\
\dL{1} m_{hH}^2 ={}& \frac{1}{2}\left(\dL{1} M_A^2 s_{2(\alpha-\beta)} - \dL{1} M_Z^2 s_{2(\alpha+\beta)}\right) \nonumber\\
&+ \frac{e}{2 M_W s_W c_W}\left(\dL{1} T_H s_{\alpha-\beta}^3 - \dL{1} T_h c_{\alpha-\beta}^3\right) \nonumber\\
&- c_\beta^2\dOLt_\beta (M_Z^2 c_{2(\alpha+\beta)} + M_A^2 c_{2(\alpha-\beta)}) \nonumber\\
&- c_\alpha^2\dOLt_\alpha (M_Z^2 c_{2(\alpha+\beta)}-M_A^2 c_{2(\alpha-\beta)}), \label{eq:dmhH}\\
\dL{1} m_{H}^2 ={}& \dL{1} M_A^2 s_{\alpha-\beta}^2 + \dL{1} M_Z^2 c_{\alpha+\beta}^2 \nonumber\\
&- \frac{e}{2 M_W s_W c_W}\left(\dL{1} T_H c_{\alpha-\beta}(1+s_{\alpha-\beta}^2) + \dL{1} T_h s_{\alpha-\beta}c_{\alpha-\beta}^2\right) \nonumber\\
&- c_\beta^2\dOLt_\beta (M_Z^2 s_{2(\alpha+\beta)} + M_A^2 s_{2(\alpha-\beta)}) \nonumber\\
&- c_\alpha^2\dOLt_\alpha (M_Z^2 s_{2(\alpha+\beta)}-M_A^2 s_{2(\alpha-\beta)}),\\
\dL{1} m_{AG}^2 ={}& \frac{e}{2 M_W s_W c_W}\left(-\dL{1} T_H s_{\alpha-\beta} - \dL{1} T_h c_{\alpha-\beta}\right) \nonumber\\
&- M_A^2 c_\beta^2\left(\dOLt_\beta - \dOLt_{\beta_n}\right),\\
\dL{1} m_{G}^2 ={}& \frac{e}{2 M_W s_W c_W}\left(-\dL{1} T_H c_{\alpha-\beta} + \dL{1} T_h s_{\alpha-\beta}\right)\\
\dL{1} m_{hA}^2 ={}& \frac{e}{2 M_W s_W c_W}\dL{1} T_A s_{\alpha-\beta},\\
\dL{1} m_{HA}^2 ={}& - \dL{1} m_{hG}^2, \\
\dL{1} m_{HG}^2 ={}& \dL{1} m_{hA}^2
\end{align}
\end{subequations}
and for the charged Higgs sector by
\begin{subequations}
\begin{align}
\dL{1} m_{{H^\pm}}^2 ={}& \dL{1} M_A^2 - \dL{1} M_W^2, \\
\dL{1} m_{{H^\pm}{G^\pm}}^2 ={}& \frac{e}{2 M_W s_W c_W}\left(-\dL{1} T_H s_{\alpha-\beta} - \dL{1} T_h c_{\alpha-\beta} - i \dL{1} T_A \right) \nonumber\\
&- m_{H^\pm}^2 c_\beta^2 \left(\dOLt_\beta - \dOLt_{\beta_c}\right) \\
\dL{1} m_{{G^\pm}{H^\pm}}^2 ={}& \left(\dL{1} m_{H^\pm G^\pm}^2\right)^* \\
\dL{1} m_{{G^\pm}}^2 ={}& \frac{e}{2 M_W s_W c_W}\left(-\dL{1} T_H c_{\alpha-\beta} + \dL{1} T_h s_{\alpha-\beta}\right).
\end{align}
\end{subequations}

\subsection{Two-loop renormalization}

The relations between the two-loop field renormalization constants of the neutral Higgs sector in the mass basis and in the original gauge basis are given by
\begin{subequations}
\begin{align}
\dZTL{hh} &= s_\alpha^2 \DZ{11} - s_\alphac_\alpha \left(\DZ{12} + \DZ{21}\right) + c_\alpha^2 \DZ{22}, \\
\dZTL{HH} &= c_\alpha^2 \DZ{11} + s_\alphac_\alpha \left(\DZ{12} + \DZ{21}\right) + s_\alpha^2 \DZ{22}, \\
\dZTL{hH} &= s_\alphac_\alpha \left(\DZ{22} - \DZ{11}\right) + \frac{1}{2}c_{2\alpha}\left(\DZ{12}+\DZ{21}\right), \\
\dZTL{AA} &= s_{\beta_n}^2 \DZ{11} - s_{\beta_n}c_{\beta_n} \left(\DZ{12} + \DZ{21}\right) + c_{\beta_n}^2 \DZ{22}, \\
\dZTL{AG} &= s_{\beta_n}c_{\beta_n} \left(\DZ{22} - \DZ{11}\right) + \frac{1}{2}c_{2\beta_n}\left(\DZ{12}+\DZ{21}\right), \\
\dZTL{hA} &= \frac{i}{2}s_{\beta_n-\alpha}\left(\DZ{12}-\DZ{21}\right), \\
\dZTL{HA} &= \frac{i}{2}c_{\beta_n-\alpha}\left(\DZ{12}-\DZ{21}\right),
\end{align}
\end{subequations}
where $\DZ{ij} = \dZTL{ij} - \frac{1}{4}(\dZOL{ij})^2$. The corresponding expressions for the charged Higgs sector read
\begin{subequations}
\begin{align}
\dZTL{{H^\pm}\Hpm} &= s_{\beta_c}^2 \DZ{11} - s_{\beta_c}c_{\beta_c} \left(\DZ{12} + \DZ{21}\right) + c_{\beta_c}^2 \DZ{22}, \\
\dZTL{{G^\pm}\Gpm} &= c_{\beta_c}^2 \DZ{11} + s_{\beta_c}c_{\beta_c} \left(\DZ{12} + \DZ{21}\right) + s_{\beta_c}^2 \DZ{22}, \\
\dZTL{{H^\pm}{G^\pm}} &= s_{\beta_c}c_{\beta_c} \left(\DZ{22} - \DZ{11}\right) + c_{\beta_c}^2\DZ{21} - s_{\beta_c}^2\DZ{12}, \\
\dZTL{{G^\pm}{H^\pm}} &= s_{\beta_c}c_{\beta_c} \left(\DZ{22} - \DZ{11}\right) + c_{\beta_c}^2\DZ{12} - s_{\beta_c}^2\DZ{21}.
\end{align}
\end{subequations}
The field renormalization counterterms in the gaugeless limit can be obtained by setting $\alpha = \beta - \pi/2$.

The two-loop mass counterterms in the neutral Higgs sector are given by
\begin{subequations}
\begin{align}
\dL{2} m_h^2 ={}& M_A^2 \left(c_\beta^2\dOLt_\beta - s_\beta^2\dOLt_\alpha\right)^2 \nonumber\\
&- \frac{e}{2 M_W s_W}\left[\delta^{(2)}T_h + \dL{1} T_h\delta^{(1)}Z_W + \dL{1} T_H\left(c_\beta^2\dOLt_\beta - s_\beta^2\dOLt_\alpha\right) \right],\\
\dL{2} m_{hH}^2 ={}& M_A^2\left[c_\beta^2\dTLt_\beta - s_\beta^2\dTLt_\alpha - s_\betac_\beta^3(\dOLt_\beta)^2 - s_\beta^3c_\beta(\dOLt_\alpha)^2\right] \nonumber\\
& + \dL{1} M_A^2\left(c_\beta^2\dOLt_\beta - s_\beta^2\dOLt_\alpha\right) \nonumber\\
&- \frac{e}{2 M_W s_W}\left[\delta^{(2)}T_H+\delta^{(1)}T_H\delta^{(1)}Z_W \right],\\
\dL{2} m_H^2 =& \delta^{(2)}M_A^2 + \frac{e}{M_W s_W}\dL{1} T_H \left(c_\beta^2\dOLt_{\beta_n} - s_\beta^2\dOLt_\alpha\right) \nonumber\\
&+ M_A^2\left[c_\beta^4(\dOLt_{\beta_n})^2 - s_\beta^4(\dOLt_\alpha)^2 - 2c_\beta^2\dOLt_\beta\left(c_\beta^2\dOLt_{\beta_n} - s_\beta^2\dOLt_\alpha\right)\right],\\
\dL{2} m_{A G}^2 ={}& M_A^2c_\beta^2 \left[\dTLt_{\beta_n} - \dTLt_\beta + s_\betac_\beta\left((\dOLt_\beta)^2-(\dOLt_{\beta_n})^2\right)\right] \nonumber\\
& + \dL{1} M_A^2c_\beta^2\left(\dOLt_{\beta_n} - \dOLt_\beta\right) \nonumber\\
& + \frac{e}{2 M_W s_W}\left[\dL{2} T_H + \dL{1} T_H\delta^{(1)}Z_W + \dL{1} T_h \left(c_\beta^2\dOLt_{\beta_n} - s_\beta^2\dOLt_\alpha\right) \right],\\
\dL{2} m_G^2 ={}& M_A^2c_\beta^4 \left(\dOLt_{\beta_n} - \dOLt_\beta\right)^2 \nonumber\\
& - \frac{e}{2 M_W s_W}\left[\dL{2} T_h  + \dL{1} T_h\delta^{(1)}Z_W  \right.\nonumber\\
&\left.\hspace{2.3cm}- \dL{1} T_H \left(2c_\beta^2\dOLt_{\beta_n} - c_\beta^2\dOLt_\beta - s_\beta^2\dOLt_\alpha\right) \right],\\
\dL{2} m_{hA}^2 ={}& - \frac{e}{2 M_W s_W} \left(\dL{2} T_A  + \dL{1} T_A\delta^{(1)}Z_W\right), \\
\dL{2} m_{hG}^2 ={}& - \frac{e}{2 M_W s_W} \dL{1} T_A \left(c_\beta^2\dOLt_{\beta_n} - s_\beta^2\dOLt_\alpha\right), \\
\dL{2} m_{HA}^2 ={}& - \dL{2} m_{hG}^2, \\
\dL{2} m_{HG}^2 ={}& \dL{2} m_{hA}^2
\end{align}
\end{subequations}
and for the charged Higgs sector by
\begin{subequations}
\begin{align}
\dL{2} m_{{H^\pm}{G^\pm}}^2 ={}& m_{H^\pm}^2c_\beta^2 \left[\dTLt_{\beta_c} - \dTLt_\beta + s_\betac_\beta\left((\dOLt_\beta)^2-(\dOLt_{\beta_c})^2\right)\right] \nonumber\\
& + \dL{1} m_{H^\pm}^2c_\beta^2\left(\dOLt_{\beta_c} - \dOLt_\beta\right) \nonumber\\
& + \frac{e}{2 M_W s_W}\left[\dL{2} T_H - i \dL{2} T_A + \dL{1} T_h \left(c_\beta^2\dOLt_{\beta_c} - s_\beta^2\dOLt_\alpha\right) \right],\\
\dL{2} m_{{G^\pm}{H^\pm}}^2 ={}& (\dL{2} m_{H^\pm G^\pm}^2)^*,\\
\dL{2} m_{{G^\pm}}^2 ={}& m_{H^\pm}^2c_\beta^4 \left(\dOLt_{\beta_c} - \dOLt_\beta\right)^2 \nonumber\\
& - \frac{e}{2 M_W s_W}\left[\dL{2} T_h - \dL{1} T_H \left(2c_\beta^2\dOLt_{\beta_c} - c_\beta^2\dOLt_\beta - s_\beta^2\dOLt_\alpha\right) \right].
\end{align}
\end{subequations}
Note that the expressions are only valid in the limit of vanishing electroweak gauge couplings. Correspondingly, also the counterterm of $\alpha$ has to be inserted applying the limit of vanishing electroweak gauge couplings.

The remaining two-loop mass counterterms, $\dmTL{A}$ and $\dmTL{{H^\pm}}$, are fixed by employing OS conditions for either the mass of the $A$ boson or the mass of the ${H^\pm}$ boson (and by using $\dmTL{A} = \dmTL{{H^\pm}}$, which holds also only in the limit of vanishing electroweak gauge couplings).

In order to obtain the renormalized two-loop self-energies (see \Eq{eq:renTLse}) the following combinations of mass and field counterterms are needed in the neutral Higgs sector
\begin{subequations}\label{Eq:TLseCTs}
\begin{align}
\delta^{(2)}m_h^{\textbf{Z}} ={}
&   \frac{1}{4}M_A^2\left[\left(\dZOL{hH}\right)^2+\left(\dZOL{hA}\right)^2\right] \nonumber\\
& + \dZOL{hh}\dmOL{h} + \dZOL{hH}\dmOL{hH} \nonumber\\
& + \dZOL{hA}\dmOL{hA} + \dZOL{hG}\dmOL{hG} \nonumber\\
& + \dmTL{h}, \\
\delta^{(2)}m_H^{\textbf{Z}} ={}
&   M_A^2\left[\dZTL{HH}+\frac{1}{4}\left(\dZOL{HH}\right)^2+\frac{1}{4}\left(\dZOL{HA}\right)^2\right]\nonumber\\
& + \dZOL{HH}\dmOL{H}+\dZOL{hH}\dmOL{hH} \nonumber\\
& + \dZOL{HA}\dmOL{HA}+\dZOL{HG}\dmOL{HG} \nonumber\\
& + \dmTL{H},\\
\delta^{(2)}m_A^{\textbf{Z}} ={}
&   M_A^2\left[\dZTL{AA}+\frac{1}{4}\left(\dZOL{AA}\right)^2+\frac{1}{4}\left(\dZOL{HA}\right)^2\right]\nonumber\\
& + \dZOL{AA}\dmOL{A}+\dZOL{AG}\dmOL{AG} \nonumber\\
& + \dZOL{hA}\dmOL{hA}+\dZOL{HA}\dmOL{HA} \nonumber\\
& + \dmTL{A},\\
\delta^{(2)}m_{hH}^{\textbf{Z}} ={}
&   \frac{1}{2}M_A^2\left[\delta^{(2)}Z_{hH} + \frac{1}{2}\dZOL{HH}\dZOL{hH} + \frac{1}{2}\dZOL{hA}\dZOL{HA}\right] \nonumber\\
& + \frac{1}{2}\left[\left(\dZOL{hh}+\dZOL{HH}\right)\dmOL{hH}+\dZOL{hH}\left(\dmOL{h}+\dmOL{H}\right)\right] \nonumber\\
& + \frac{1}{2}\left[\dZOL{hA}\dmOL{HA}+\dZOL{HA}\dmOL{hA}+\dZOL{hG}\dmOL{HG}+\dZOL{HG}\dmOL{hG}\right] \nonumber\\
& + \dmTL{hH},\\
\delta^{(2)}m_{hA}^{\textbf{Z}} ={}
&   \frac{1}{2}M_A^2\left[\delta^{(2)}Z_{hA} + \frac{1}{2}\dZOL{hA}\dZOL{AA} + \frac{1}{2}\dZOL{hH}\dZOL{HA}\right] \nonumber\\
& + \frac{1}{2}\left[\left(\dZOL{hh}+\dZOL{AA}\right)\dmOL{hA}+\dZOL{hH}\dmOL{HA}+\dZOL{AG}\dmOL{hG}\right] \nonumber\\
& + \frac{1}{2}\left[\dZOL{hA}\left(\dmOL{hh}+\dmOL{AA}\right)+\dZOL{hG}\dmOL{AG}+\dZOL{HA}\dmOL{hH}\right] \nonumber\\
& + \dmTL{hA},\\
\delta^{(2)}m_{HA}^{\textbf{Z}} ={}
&   M_A^2\left[\delta^{(2)}Z_{HA} + \frac{1}{4}\dZOL{HH}\dZOL{HA} + \frac{1}{4}\dZOL{HA}\dZOL{AA}\right] \nonumber\\
& + \frac{1}{2}\left[\left(\dZOL{HH}+\dZOL{AA}\right)\dmOL{HA}+\dZOL{hH}\dmOL{hA}+\dZOL{AG}\dmOL{HG}\right] \nonumber\\
& + \frac{1}{2}\left[\dZOL{HA}\left(\dmOL{HH}+\dmOL{AA}\right)+\dZOL{HG}\dmOL{AG}+\dZOL{hA}\dmOL{hH}\right] \nonumber\\
& + \dmTL{HA}
\end{align}
\end{subequations}
and the charged Higgs sector
\begin{align}
\delta^{(2)}m_{{H^\pm}}^{\textbf{Z}} ={}
& m_{{H^\pm}}^2\left[\dZTL{{H^\pm}\Hpm}+\frac{1}{4}\left(\dZOL{{H^\pm}\Hpm}\right)^2\right] \nonumber\\
& + \dZOL{{H^\pm}\Hpm}\dmOL{{H^\pm}} \nonumber\\
& + \frac{1}{2}\left(\dZOL{{H^\pm}{G^\pm}}\dmOL{{G^\pm}{H^\pm}}+\dZOL{{G^\pm}{H^\pm}}\dmOL{{H^\pm}{G^\pm}}\right)\nonumber\\
& + \dmTL{{H^\pm}}.
\end{align}
These expressions are again only valid in the gaugeless limit. Correspondingly, this limit has also to be applied for the appearing field renormalization constants.

The two-loop tadpole counterterms are obtained by
\begin{samepage}
\begin{subequations}
\begin{align}
\dL{2} T_h ={}& - \frac{1}{2}\left(\delta^{(1)}Z_{hh}\delta^{(1)}T_h+\delta^{(1)}Z_{hH}\delta^{(1)}T_H \right.\nonumber\\
&\left.\hspace{.9cm}+\delta^{(1)}Z_{hA}\delta^{(1)}T_A+\delta^{(1)}Z_{hG}\delta^{(1)}T_G\right), \\
\dL{2} T_H ={}& -\frac{1}{2}\left(\delta^{(1)}Z_{HH}\delta^{(1)}T_H+\delta^{(1)}Z_{hH}\delta^{(1)}T_h \right.\nonumber\\
&\left.\hspace{.9cm}+\delta^{(1)}Z_{HA}\delta^{(1)}T_A+\delta^{(1)}Z_{HG}\delta^{(1)}T_G\right), \\
\dL{2} T_A ={}& - \frac{1}{2}\left(\delta^{(1)}Z_{AA}\delta^{(1)}T_A+\delta^{(1)}Z_{AG}\delta^{(1)}T_G \right.\nonumber\\
&\left.\hspace{.9cm}+\delta^{(1)}Z_{hA}\delta^{(1)}T_h+\delta^{(1)}Z_{HA}\delta^{(1)}T_H\right).
\end{align}
\end{subequations}
\end{samepage}

\section{Introduction}
\label{sec:01_intro}

\input{01_intro.tex}

\section{The MSSM Higgs sector and its renormalization}
\label{sec:02_MSSMintro}

\input{02_MSSMintro.tex}

\section{Momentum dependence of heavy contributions}
\label{sec:03_p2heavy}

\input{03_p2heavy.tex}

\section{Higgs pole determination in the MSSM}
\label{sec:04_poledet}

\input{04_poledet.tex}

\section{Numerical results}
\label{sec:05_results}

\input{05_results.tex}

\section{Conclusions}
\label{sec:06_conclusion}

\input{06_conclusion.tex}

\section*{Acknowledgments}
\sloppy{
We thank Wolfgang Hollik, Sven Heinemeyer, Georg Weiglein for helpful comments and carefully reading the manuscript. We thank Thomas Hahn, Stefan Liebler, Sebastian Paßehr, Heidi Rzehak, Pietro Slavich, Ivan Sobolev as well as Tim Stefaniak for useful discussions.
}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Computer Vision is one of the advanced AI domains that allows computers to see and understand visual data (images/videos) in a way similar to how humans do. Advancements in computer vision research has revolutionized the usefulness of AI applications in daily lives. Also, with an increase in computational capacity along with research for developing neural networks with low latency, deep learning models can be deployed on mobile devices- making them accessible and affordable to the common people. One such field which has reaped significant benefits from these advancements is assistive technology for the visually impaired people.

A major problem faced by visually impaired people all over the world is crossing roads. White canes or walking sticks can help them scan their surroundings to detect any obstacles and/or orientation marks. But road crossing is a very sensitive activity, and requires a quick decision from pedestrian which is beyond the scope of normal walking sticks or white canes.

Companies like eSight \cite{eSight}, Seeing AI \cite{SeeingAI}, and the Sound of Vision System \cite{Caraiman_2017_ICCV} have developed some useful products for visually impaired people over the years; but they are not specifically designed for helping them cross roads. Also, these products are very costly so they are not feasible for use in developing countries. Low cost is a very essential factor to create a large scale impact in developing countries like India. Also, one major constraint while crossing roads in India (and other developing countries) is that there are no traffic lights or pedestrian crosswalks in small towns and small internal roads. So, our work in this paper is an attempt to develop a model to predict road crossing safety on Indian roads.

To summarize, we make the following contributions:
\begin{enumerate}
    \item Data: This paper introduces INDRA (INdian Dataset for RoAd crossing), a first of its kind dataset capturing Indian traffic scenarios (from pedestrian point-of-view) labelled with road crossing safety information.
    \item Prediction models: This paper presents analysis of rigorous training experiments performed to develop a classifier to predict road crossing safety.
    \item Working prototype: This paper also presents a low-cost and feasible working prototype of road crossing assistant which can potentially be used by visually impaired people to cross roads in real-time.
\end{enumerate}

\section{Related Work}

To the best of our knowledge, all the existing attempts to help visually impaired people cross roads are based on detection and localization of pedestrian crosswalk signs and traffic lights (hence applicable in developed countries only). In this section, we summarize some of these attempts.

One attempt to help visually impaired people cross roads has been the use of acoustic pedestrian lights \cite{harkey2007accessible}. These traffic lights use audible tones, verbal messages and/or vibrating surfaces to notify the pedestrians about useful information such as the color of traffic light and road crossing safety. Though this concept appears to be very useful and reliable, it is implemented only in urban areas. Also, there is almost no presence of acoustic lights in developing countries like India.

Recent years have witnessed significant improvements in the task of detecting traffic lights for autonomous driving \cite{almeida2018traffic, li2021improved, gong2010recognition}. These advancements have been extended to improve the detection of pedestrian traffic lights for navigation support of visually impaired people \cite{martinez2017using}. But in developing countries like India, crossing roads is something much beyond just detecting pedestrian crosswalks and traffic lights. Initiatives like NavAssistAI \cite{NavAssistAI} and XIMIRA \cite{XIMIRA} are noteworthy when reviewing the existing solutions to help blind people cross roads. But again, they are not feasible for Indian traffic scenarios.

There are India road datasets available for autonomous navigation (see Table 1 for comparison), but they are not recorded from a pedestrian's point of view. Also, to the best of our knowledge, no existing Indian road dataset is labelled for road crossing safety.

\section{Data}

Any dataset suitable for the task of predicting road crossing safety on Indian roads did not exist. There are  datasets available for semantic segmentation and object detection (see Table 1) for autonomous driving, but they are not recorded from a pedestrian's point of view, and neither are they labelled for road crossing safety. So we introduce our own dataset: \textit{\textbf{INDRA (INdian Dataset for RoAd crossing)}}. Our dataset consists of 104 videos ($>$26K frames) capturing various traffic patterns from roads of a town of Gujarat, India. We collected the videos using head-mounted SJCam SJ4000 12MP action camera (Figure 2). Video specifications: MOV format, 1920*1080 resolution, 30 fps.

\begin{figure}[htp]
    \centering
    \includegraphics[width=0.5\columnwidth]{images/hardware-action-camera.jpeg}
    \caption{SJCam SJ4000 12MP action camera}

\end{figure}

We also manually annotated the videos with timestamps for durations safe for crossing (example images in Figure 1). For annotations, we had made a simple key-press based utility software: the annotator had to press pre-defined keys to store the safe-start and safe-end durations (in terms of frame number) in a csv file. The software also had provision to store multiple safe durations from a given video. We also release vehicle bounding boxes detected by state of the art RetinaNet for each frame of INDRA as auxiliary data for the community to use in further research.

\section{Problem Statement}

The task is to train models to predict road crossing safety using the dataset mentioned above. More precisely, given a sequence of k frames (say, $f_{n}$ to $f_{n-k}$), the goal is to predict whether $f_{n+1}$ is safe (here $f_{n}$ represents $n^{th}$ frame of the sequence). Precision-recall characteristics can be used to evaluate such classifiers \cite{hossin2015review, goutte2005probabilistic}. To elaborate:
\begin{enumerate}
    \item \textbf{Precision score} represents the proportion of positive predictions made by the model which are actually positive (i.e., high precision implies low number of false positives). It is extremely important to have high precision because the safety of the user depends on a low number of false positives.

    \item \textbf{Recall score} represents the proportion of positive predictions identified correctly by the model (i.e., high recall implies low number of false negatives). While high recall score is not the foremost priority in this application, higher recall is better because in the greed of getting a high precision score, we do not want the model to predict every frame as unsafe for crossing the road.

\end{enumerate}

\section{Machine Learning Methods}

In this section, we present the details of support vector machines (SVMs) trained on hand crafted features extracted from video frames of INDRA. We used 80 videos for training and remaining 24 videos for testing the trained SVMs.
\subsection{Single frame SVM}
It is our simplest approach where we extracted simple per-frame features capturing number, location, size, speed and direction of vehicles. Further, we trained SVM classifier to predict whether a particular frame is safe for crossing road. For extracting features from a particular frame, we divided the video frame into 6 regions as shown in Figure 3.
\begin{figure}[htp]
    \centering
    \includegraphics[width=7.5cm]{images/single_frame_svm.png}
    \caption{Feature extraction for single frame SVM model (Note: s1 = speed of truck \& s2 = speed of bike)}

\end{figure}

 It is very intuitive because even we as humans, when crossing roads, consciously/unconsciously consider information about vehicles in different regions of our field of vision. We extract following features from each of the 6 regions:

\begin{enumerate}
    \item Number of vehicles in the region.
    \item Total area covered by vehicles in the region
    \item Distance (from origin) of the vehicle closest to the origin. Here, the origin approximately represents the position of user with head mounted camera.
    \item Maximum speed of any vehicle in the region.
\end{enumerate}
Hence, a total of 6 * 4 = 24 features are extracted from each frame.
\begin{figure}[htp]
    \centering
    \includegraphics[width=\columnwidth]{images/directional_filtering.png}
    \caption{Directional vehicle filtering using KCF tracker}

\end{figure}

\begin{figure}[htp]
    \centering
    \includegraphics[width=\columnwidth]{images/speed.png}
    \caption{Vehicle speed detection using KCF tracker}

\end{figure}

We also applied directional vehicle filtering on these features. As seen in Figure 4, for crossing road, the information of vehicles on the other side of the divider (i.e., the vehicles marked in red) is not required. So, we did not consider those vehicles during feature extraction.
To calculate the number, area and distance of vehicles, we used ImageAI library for video object detection using RetinaNet \cite{lin2017focal} to generate and save the vehicle bounding boxes. To get direction and speed (see Figure 5) of these vehicle bounding boxes, we used OpenCV's KCF (Kernelized Correlation Filters) \cite{henriques2014high} tracker. We applied minmax scaling to the extracted features and trained an SVM on the train set.

\subsection{Multi frame SVM}
As it is obvious that even we as humans do not decide whether it is safe to cross a road by just having one glance at the road, we extended the above approach and extracted multi-frame features instead of per-frame features.
\begin{figure}[htp]
    \centering
    \includegraphics[width=\columnwidth]{images/multi_frame_svm.jpg}
    \caption{Feature extraction for multi frame SVM model}

\end{figure}

\begin{figure*}
    \centering
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/tp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/tp2.jpg}} \hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/tn1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/tn2.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/fp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/ml1/fn1.jpg}}

\end{figure*}

\begin{figure*}
    \centering
    \subfloat[True positive]{\includegraphics[width=0.15\textwidth]{images/ml2/tp1.jpg}}\hfil
    \subfloat[True positive]{\includegraphics[width=0.15\textwidth]{images/ml2/tp2.jpg}} \hfil
    \subfloat[True negative]{\includegraphics[width=0.15\textwidth]{images/ml2/tn1.jpg}}\hfil
    \subfloat[True negative]{\includegraphics[width=0.15\textwidth]{images/ml2/tn2.jpg}}\hfil
    \subfloat[False positive]{\includegraphics[width=0.15\textwidth]{images/ml2/fp1.jpg}}\hfil
    \subfloat[False negative]{\includegraphics[width=0.15\textwidth]{images/ml2/fn1.jpg}}

    \caption{From top to bottom: sample prediction outputs of single frame SVM and multi frame SVM respectively.}
\end{figure*}

As shown in Figure 6 , we used a sliding window based approach for feature extraction. For predicting label of a particular frame, we considered an entire window of frames (i.e., the current frame plus past k frames). The feature vector for a particular frame consists of 24 features of that frame (using the per frame features explained in single frame SVM section), and it additionally consists of predicted labels for previous k-1 frames using saved model of single frame SVM. So if we take window size = 10, the length of feature vector will be 24 + 9 = 33. We applied minmax scaling to the extracted features and trained an SVM on the train set.

\subsection{Results}

\begin{table}[htp]
  \caption{ML methods: Results on test data}
  \begin{tabular}{ccl}
    \toprule
    Method & Precision & Recall\\
    \midrule
    Single frame SVM & 0.68 & 0.87 \\
    Multi frame SVM & 0.79 & 0.84 \\
  \bottomrule
\end{tabular}
\end{table}
Table 2 shows the precision and recall scores (using threshold of 0.5 to predict the frame as safe/unsafe) on test split of INDRA achieved by the aforementioned ML methods. Figure 7 shows sample prediction outputs for the ML methods. While we feel that it was a good idea to start with simple SVMs as we were working on a completely new dataset (without any prior baseline performance for this task), instead of trying to improve these methods, we switched to more advanced deep learning (DL) methods because of the following reasons:
\begin{enumerate}
    \item The precision scores of 0.68 and 0.79 are very poor given that road crossing is a very sensitive application.
    \item For the training experiments, we saved the vehicle bounding boxes generated using RetinaNet beforehand. The entire inference pipeline starting from detecting vehicles, extracting features and predicting the road crossing safety using pre-trained SVM would have been very slow in real-time, thus making it infeasible for deployment.
\end{enumerate}

\section{Deep Learning Methods}

Deep learning (DL) has shown to learn highly effective features from image and video data, yielding high accuracy in many tasks. Hence, we performed many training experiments using convolutional neural networks (CNNs) on individual video frames of INDRA. We performed these experiments on AWS g4dn.xlarge instance (with Tesla T4 GPU). The data preparation and CNN architecture details are explained in this section.

\subsection{Data Preparation}

Firstly, we performed train-test-validation split on the dataset. We used 66 videos in train set, 22 in test set and 16 in validation set.

Then, we created an input data pipeline using Tensorflow's Dataset API \cite{DataPipeline}. For the same, we stored all the individual frames from the videos in a separate folder (as jpg images) and made a list containing all the filenames. Next, we performed the following steps to create input pipeline:

\begin{enumerate}
    \item Created tensorflow dataset instance from slices of the filenames and labels.
    \item Shuffled the data with a buffer size equal to length of the dataset.
    \item Parsed the images from filename to pixel values (used multiple threads to improve speed). Also, we resized the images to the required input size of the models.
    \item Applied data augmentation by using random brightness, contrast and saturation (used multiple threads to improve the speed).
    \item Batched the images (batch size = 16).
    \item Prefetched one batch to make sure that a batch is ready to be served at all time.
\end{enumerate}

\subsection{Class weights}

Using class weights was essential during model training since the dataset is highly unbalanced (obviously there are far more unsafe frames than safe frames in traffic scenarios). Hence we used class weights as: 1 for unsafe class, 1.92 for safe class in each of the following DL methods. We obtained the value of 1.92 by calculating the ratio of number of unsafe frames to the number of safe frames from train set.

\subsection{Baseline Model - MobileNetV2}

As a baseline model, we used the MobileNetV2 architecture \cite{sandler2018mobilenetv2} with additional dense layers at the top. We used MobileNetV2 architecture because it is very memory efficient and lightweight, making it highly efficient for mobile and embedded vision applications. The key properties of MobileNetV2 architecture that makes it lightweight and efficient are the use of depthwise separable convolutions and an inverted residual structure.

The model architecture is as shown in Figure 8. On top of the last convolutional layer of MobileNetV2, we added an average pooling layer for downsampling, followed by a dropout layer and a fully connected layer with 1 unit (i.e., the output layer giving the probability of a frame being safe).

\begin{figure}[htp]

            \includegraphics[width=\columnwidth]{images/a1.png}
            \caption{Baseline Model - MobileNetV2}

\end{figure}

We trained the model for 150 epochs minimizing the binary cross entropy loss. We used dropout rate of 0.4 for the additional layer, and used Adam's optimizer with learning rate of 5e-4.

\subsection{RoadCrossNet}

\begin{figure}[htp]
            \centering
            \includegraphics[width=\columnwidth]{images/approach_3.2.png}
            \caption{RoadCrossNet}

\end{figure}

We observed that MobileNetV2 exhibited a clear case of overfitting, so we attempted to develop a new lightweight architecture: \textbf{\textit{RoadCrossNet}} from scratch. The model architecture is as shown in Figure 9.

Similar to the trend in many neural network architectures, in RoadCrossNet, we increased the number of filters and decreased the filter size while going deeper in the network. As shown in the architecture diagram, we also used multiple batch-normalization and dropout layers to prevent overfitting. Aditionally, for downsampling, we used max pooling layers.

We trained the model for 150 epochs minimizing binary cross entropy loss. We used dropout rate of 0.3 for first three dropout layers and used dropout rate of 0.4 for subsequent dropout layers. Also, we used Adam's optimizer with learning rate of 5e-4 for training, and used ReLU activation for the layers.

\subsection{DilatedRoadCrossNet}

\begin{figure}[htp]
            \centering
            \includegraphics[width=\columnwidth]{images/approach_3.3.png}
            \caption{DilatedRoadCrossNet}

\end{figure}

Since it is difficult to optimize convolution layers with larger kernel size (optimization required for real-time deployment), we improved the RoadCrossNet architecture by reducing the kernel sizes of initial layers (which were 7*7 and 5*5). Also, to compensate for the reduced receptive field caused by decrease in kernel size, we used dilated convolutions instead of simple convolutions \cite{luo2016understanding}. Intuitively speaking, receptive field plays an important role in this application since for crossing roads, one needs to consider his/her entire field of vision to make the decision, and not just a particular region from the field of vision. The model architecture is as shown in Figure 10. Apart from using dilated convolutions, the architecture is the same as RoadCrossNet (containing batch normalization, dropout and max pooling layers, with 3 dense layers at the top).

We trained the model for 150 epochs minimizing binary cross entropy loss. We used dropout rate of 0.3 in first three dropout layers and used dropout rate of 0.4 in later two layers. Also, we used Adam's optimizer with learning rate of 5e-4 for training.
\begin{figure*}
    \centering
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/tp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/tp2.jpg}} \hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/tn1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/tn2.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/fp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl1/fn1.jpg}}

\end{figure*}

\begin{figure*}
    \centering
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/tp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/tp2.jpg}} \hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/tn1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/tn2.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/fp1.jpg}}\hfil
    \subfloat{\includegraphics[width=0.15\textwidth]{images/dl2/fn1.jpg}}

\end{figure*}

\begin{figure*}
    \centering
    \subfloat[True positive]{\includegraphics[width=0.15\textwidth]{images/dl3/tp1.jpg}}\hfil
    \subfloat[True positive]{\includegraphics[width=0.15\textwidth]{images/dl3/tp2.jpg}} \hfil
    \subfloat[True negative]{\includegraphics[width=0.15\textwidth]{images/dl3/tn1.jpg}}\hfil
    \subfloat[True negative]{\includegraphics[width=0.15\textwidth]{images/dl3/tn2.jpg}}\hfil
    \subfloat[False positive]{\includegraphics[width=0.15\textwidth]{images/dl3/fp1.jpg}}\hfil
    \subfloat[False negative]{\includegraphics[width=0.15\textwidth]{images/dl3/fn1.jpg}}

    \caption{From top to bottom: sample prediction outputs of MobileNetv2, RoadCrossNet and DilatedRoadCrossNet respectively.}
\end{figure*}

\subsection{Results}

\begin{table}[htp]
  \caption{DL methods: Results on test data}
  \begin{tabular}{c c c c}
    \toprule
    Method & Precision & Recall & Throughput\\
    \midrule
    MobileNetv2 & 0.90 & 0.60 & 15 fps\\
    RoadCrossNet & 0.90 & 0.72 & 3 fps \\
    DilatedRoadCrossNet & 0.90 & 0.77 & 8 fps \\
  \bottomrule
\end{tabular}
\end{table}
Table 3 shows the precision and recall scores on test set of INDRA achieved by the DL methods as explained in detail above. The table also shows inference throughput obtained on Nvidia Jetson Nano in real-time (deployment details explained in the next section). Sample prediction outputs of the DL methods are shown in Figure 11. We chose DilatedRoadCrossNet for deployment on our working prototype as it exhibited a high precision score, along with decent enough recall score and inference throughput.

\section{Deployment}

\begin{figure*}[htp]
            \centering
            \includegraphics[width=0.8\textwidth]{images/setup1.jpg}
            \caption{Working prototype of our Road Crossing Assistant}

\end{figure*}

We used Nvidia Jetson Nano and few other supporting components (powerbank, headphones, camera, USB cables \& adapters)  to develop a wearable road-crossing assistant for use in real-time. Low cost is essential for a large-scale impact in developing countries like India. Therefore, we chose to work with a separate camera and mobile computer (Nvidia Jetson Nano) rather than developing a smartphone app. Smartphones with the required deep learning inference hardware are currently significantly more expensive than our discrete combination of components. The workflow of our prototype is as shown in Figure 13.

Firstly we converted our pretrained DilatedRoadCrossNet model into a TensorRT model and deployed it on Jetson Nano. Conversion to TensorRT was essential to optimize the processing power and achieve low latency with high throughput on Jetson Nano. We then wrote a script that is supposed to execute at the start up of Jetson Nano. This script ran the model inference on real-time video stream coming from the camera attached to the Jetson Nano, and also communicated with the user via headphones. Specifically, on startup, our assistant acknowledges the activation through simple audio code and guides the user to look in the direction of oncoming traffic. Next, it runs the deployed pre-trained model on the video stream and communicates with the user when it is safe to cross the road (more precisely, we set a threshold of 0.85, and notify the user to start crossing only when atleast past 5 frames have been predicted as safe). Figure 12 shows the setup of our working prototype. A video demonstrating the use of this assistant in real-time can be found in the supplementary material.
\begin{figure}[htp]
            \centering
            \includegraphics[width=\columnwidth]{images/deployment_workflow.jpg}
            \caption{Workflow of our prototype}

\end{figure}

\section{Conclusion and Future Work}

We presented INDRA, the first dataset capturing pedestrian point-of-view traffic patterns of Indian traffic scenarios. We also proposed DilatedRoadCrossNet, which showed promising results for predicting road crossing safety in real-time. Considering that this is the first attempt to train models for predicting road crossing safety from pedestrian point-of-view video stream, our work can serve as a useful baseline for upcoming research in this area. We feel that there are still a lot of ideas worth exploring in order to build a more reliable prediction model. We list some of our future work ideas below:
\begin{enumerate}
    \item Although there was diversity in terms of locations, times, traffic patterns in the dataset, all the videos of this dataset were collected from one city only. We feel that collecting more videos from different cities, different type of roads and over different seasons will add much needed diversity in the dataset.

    \item Multi-frame features: We think that using multi-frame ConvNets will significantly improve the model performance. Recurrent networks appear to be a tempting starting point for most sequence modelling projects, but after a rigorous literature study, we feel that convolutional models can prove out be a better option \cite{bai2018empirical}. Also, since predicting the road crossing safety is an application that requires a lightweight model in order to be deployed for use in real-time, we feel that using convolutional models (capturing spatio-temporal features) over recurrent models (e.g., LSTM, GRU which are comparatively heavier and computationally expensive) is the way to go. We have shortlisted the following approaches which may prove out to be useful for upcoming research in this area:

    \begin{enumerate}
        \item Collaborative spatiotemporal feature learning \cite{li2019collaborative}: This approach involves performing 2D convolution along all three orthogonal views of volumetric video data (which learns spatial appearance and temporal motion cues respectively). Also, by sharing the convolution kernels of different views, spatial and temporal features are collaboratively learned and thus benefit from each other; additionally the model doesn't turn out to be very heavy.

        \item Using raw frames and optical flow between adjacent frames as two input streams \cite{simonyan2014two, feichtenhofer2017spatiotemporal}: This is a very intuitive approach in video prediction; video can naturally be decomposed into spatial and temporal components. The individual frames carry the spatial information, while the optical flow (between adjacent frames) carries the temporal information. Both the streams are implemented using a convolutional network, and are subsequently combined to predict the class scores.

    \end{enumerate}

\end{enumerate}