---
abstract: 'For systems of unstable particles that mix with each other, an approximation of the fully momentum-dependent propagator matrix is presented in terms of a sum of simple Breit–Wigner propagators that are multiplied with finite on-shell wave function normalisation factors. The latter are evaluated at the complex poles of the propagators. The pole structure of general propagator matrices is carefully analysed, and it is demonstrated that in the proposed approximation imaginary parts arising from absorptive parts of loop integrals are properly taken into account. Applying the formalism to the neutral MSSM Higgs sector with complex parameters, very good numerical agreement is found between cross sections based on the full propagators and the corresponding cross sections based on the described approximation. The proposed approach does not only technically simplify the treatment of propagators with non-vanishing off-diagonal contributions, it is shown that it can also facilitate an improved theoretical prediction of the considered observables via a more precise implementation of the total widths of the involved particles. It is also well-suited for the incorporation of interference effects arising from overlapping resonances.'
author:
- |
    <span style="font-variant:small-caps;">Elina Fuchs</span>$^{a,}$^1^,   <span style="font-variant:small-caps;">Georg Weiglein</span>$^{b,}$^2^\
    *(a) Department of Particle Physics and Astrophysics,*\
    *Weizmann Institute of Science, Rehovot 76100, Israel*\
    *(b) DESY, Deutsches Elektronen-Synchrotron, Notkestr. 85,*\
    *D-22607 Hamburg, Germany*
bibliography:
- 'mthesis\_literature.bib'
- 'gNWALoop.bib'
- 'PhDthesis.bib'
- 'bbH.bib'
- 'propagator.bib'
title: '\'
---

Introduction {#chap:intro}
============

As a consequence of electroweak symmetry breaking, particles with the same quantum numbers of the electric charge and colour can mix with each other. This is reflected in their propagators contributing to production or decay processes in particle physics. Beyond lowest order those propagators receive contributions from all possible mixings with other particles, such that the propagators of the system of particles that can mix with each other are of matrix-type. Since the involved particles are in general unstable, a proper treatment of imaginary parts is necessary. The mass of an unstable particle is determined from the real part of the (gauge-invariant) complex pole [@Willenbrock:1991hu; @Sirlin:1991fd; @Stuart:1991xk; @Veltman:1992tm], while the imaginary part yields the total decay width of the particle. The eigenstates associated with the masses obtained from the complex poles of the propagator matrix in general differ from the interaction eigenstates that contribute at lowest order.

As an example, in the Standard Model (SM) the neutral interaction eigenstates arising from the ${\rm SU}(2)_{\rm I}$ and ${\rm U}(1)_{\rm Y}$ gauge groups mix with each other to form the mass eigenstates photon, $\gamma$, and $Z$. Higher-order contributions lead to a mixing between these states, giving rise to a $\gamma-Z$ propagator matrix. In a renormalisable gauge, the (would-be) Goldstone boson as an unphysical degree of freedom does not vanish, and mixing contributions with this unphysical scalar particle need to be taken into account as well. In the quark sector of the SM the relation between the interaction and the mass eigenstates is encoded in terms of the CKM matrix [@Cabibbo:1963yz; @Kobayashi:1973fv]. The additional states in models of physics beyond the SM (BSM) usually follow the same pattern, i.e. interaction eigenstates mix with each other to form mass eigenstates, and the propagators of the latter are of matrix-type as a consequence of mixing contributions at higher orders (in some cases it can also be convenient to derive loop-corrected mass eigenstates directly from lowest-order interaction eigenstates, avoiding the introduction of lowest-order mass eigenstates). In particular, extended Higgs sectors usually contain several neutral Higgs states that can mix with each other. If ${\mathcal{CP}}$-conservation is assumed, ${\mathcal{CP}}$-even and ${\mathcal{CP}}$-odd states can only mix among themselves. In the general case where the possibility of ${\mathcal{CP}}$-violating interactions affecting the Higgs sector is taken into account, non-vanishing mixing contributions can occur between all neutral Higgs bosons of the model. Examples are general two-Higgs-doublet models, the Higgs sector of the minimal supersymmetric extension of the SM (MSSM), as well as the Higgs sectors of non-minimal supersymmetric models. Similarly, propagator mixing between vector resonances can occur for instance for the Kaluza–Klein excitations of the $Z$-boson and the photon ($Z'$ and $A'$, respectively), which can have similar masses and have a large mixing with each other[@Cacciapaglia:2009ic]. New, possibly nearly degenerate vector bosons can also arise in extended gauge symmetries and theories with a strongly coupled sector, e.g. composite Higgs models. For the case of the mixing of two $Z'$-like vectors, see e.g. Ref.[@deBlas:2012qp]. An example involving mixing with a new fermion is the mixing between the top quark $t$ and a top partner $t'$[@Pilaftsis:1997dr].

In order to properly treat external particles in physical processes, the correct on-shell properties of in- and outgoing particles have to be ensured. While the S-matrix theory is formulated only for stable in- and outgoing states [@Eden:AnaSmatrix], in practical applications one often needs to deal with the case of unstable external particles. The proper normalisation of external particles requires in particular that all external particles are on their mass shell, that their residues are equal to one, and that the mixing contributions with all other particles vanish on-shell at the considered order. This requirement holds independently from the specific renormalisation scheme that has been adopted for a certain calculation. In case the renormalisation conditions of the chosen scheme do not impose the proper normalisation of external particles, the correct on-shell properties need to be established via UV-finite wave function normalisation factors (“Z-factors”). In the case of unstable particles, finite-width effects have to be taken into account, and more generally imaginary parts arising from the absorptive parts of loop integrals need to be properly treated. A further complication arises if ${\mathcal{CP}}$-violating effects are taken into account, since in this case the presence of complex parameters yields an additional source of imaginary parts. In the SM the different sources of imaginary parts affect for instance the renormalisation of the CKM matrix from the two-loop level onwards, see e.g. the discussion in Ref. [@Espriu:2002xv] and references therein. In the MSSM with complex parameters the corresponding effects arising from absorptive parts of loop integrals and complex model parameters enter predictions for physical observables in the chargino and neutralino sector already at the one-loop level, see Refs. [@Fowler:2009ay; @Fowler:2010eba; @Bharucha:2012nx].

In the present paper we demonstrate that on-shell wave function normalisation factors, evaluated at the complex poles, are well suited for approximating the full mixing propagators also off-shell. In order to achieve this, the elements of the propagator matrix are expressed in terms of simple Breit–Wigner propagators that are multiplied with wave function normalisation factors. We show for various examples that this approximation works very well in practice. In this context we discuss in detail the uncertainties that are associated with the neglected terms in the approximation and provide estimates of their impact, which we find to be in good agreement with the observed numerical deviations between the full result and the approximation. We furthermore demonstrate that besides very significant technical simplifications our approach also has conceptual advantages.

Concerning technical simplifications, the evaluation of the propagator matrix for a given value of the external momentum squared requires the full momentum dependence of all contributing diagonal and off-diagonal self-energies. In our approximation based on wave function normalisation factors that are evaluated at the complex poles, on the other hand, the only momentum dependence that needs to be taken into account is the one of the simple Breit–Wigner propagators. We emphasize that our approach goes significantly beyond a simple “effective coupling”-type approximation where all external momenta at the loop level are neglected. While the latter approximation would disregard the important absorptive parts of the loop functions, those contributions are taken account in our approach through the evaluation of the Z-factors at the complex poles. We will demonstrate in particular that expanding the full propagators around all of their complex poles indeed results in the sum of Breit–Wigner propagators in combination with on-shell Z-factors. Approximating the full propagator matrix in this simple and convenient way has several conceptual advantages. In particular, the formulation in terms of Breit–Wigner propagators and Z-factors facilitates the implementation of a more precise total width than the width obtained from the imaginary part of the complex pole with self-energies evaluated at the same order. Furthermore, for states with masses that are nearly degenerate (i.e. the difference between the masses of the two states is of the same order as the sum of the total widths of the two states) the mixing becomes resonantly enhanced and the resonances overlap, see e.g. Refs. [@Pilaftsis:1997dr; @Frank:2006yh; @Cacciapaglia:2009ic; @deBlas:2012qp; @Fuchs:2014ola; @Gomez-Bock:2015sia]. In such a case a single-pole approximation is not applicable [@Cacciapaglia:2009ic; @deBlas:2012qp; @Fuchs:2014ola]. Instead, the contributions from the relevant complex poles of the propagators have to be taken into account. The approximation of the propagator matrix in terms of on-shell wave function normalisation factors evaluated at the complex poles and Breit–Wigner propagators allows the implementation of interference contributions using the prescription of Ref. [@Fuchs:2014ola].

We perform in this paper a careful treatment of imaginary parts and analyse the pole structure of matrix-type propagators in detail. While in the no-mixing case every element of an $(n \times n)$ propagator matrix has a single pole, we will demonstrate that in general every single element of an $(n \times n)$ propagator matrix (diagonal and off-diagonal) by itself has $n$ poles. This fact makes the association of the (loop-corrected) mass eigenstates with the propagator poles non-trivial. We will show that in principle all permutations of the $n$ propagator poles with the $n$ mass eigenstates are physically equivalent. For practical purposes, however, it is important to make a choice that leads to numerically stable and well-behaved results.

To be specific, in the following we use the language of the Higgs sector of the MSSM with complex parameters. In this case ${\mathcal{CP}}$-violating effects entering at the loop level give rise to a mixing between the lowest-order mass eigenstates $h, H, A$, which are ${\mathcal{CP}}$-eigenstates ($h, H$ are ${\mathcal{CP}}$-even states, while $A$ is a ${\mathcal{CP}}$-odd state). The $(3 \times 3)$ mixing between the neutral Higgs states is the smallest type of propagator mixing that exhibits the qualitative features of general $(n \times n)$ propagator matrices. Furthermore, the MSSM example is well-suited to emphasise the distinction between interaction eigenstates (in this case the component field of the two Higgs doublets, e.g. $\phi_1$ and $\phi_2$), lowest-order mass eigenstates ($h, H, A$) and loop-corrected mass eigenstates, which in the MSSM case are called $h_1, h_2, h_3$. The MSSM case also illustrates the possibility of mixing with other spin states ($\gamma, Z$) and with unphysical degrees of freedom ($G$). An important feature of models like the MSSM is the fact that the masses are [*predicted*]{} from the input parameters of the model, which makes the need for an appropriate assignment of mass eigenstates with propagator poles particularly apparent. The generalisation of our results to other models, cases with $(n \times n)$ mixing for $n > 3$ and other spin states should be straightforward.

The paper is organised as follows. After reviewing aspects of higher-order contributions and mixing effects in a system of three scalar states in Section\[sect:LoopMix\], we derive the pole structure of the propagator matrix and resulting relations between the wave function normalisation factors in Section\[sect:poles\]. An analytical derivation of the multi-Breit–Wigner approximation of the full propagators including an uncertainty estimate is given in Section\[sect:fullBWZ\]. We perform a detailed numerical comparison of the two approaches in Section\[sect:NumBWfull\] before we conclude in Section\[sect:conclusions\].

Loop-level mixing of scalar propagators {#sect:LoopMix}
=======================================

In a system of $n$ particles which have the same conserved quantum numbers in general a non-trivial mixing between the different states will occur. While at lowest order the mass matrix given in terms of the interaction eigenstates can be diagonalised in order to obtain the mass eigenstates, at the loop level momentum-dependent self-energy corrections enter. These contributions give rise to a momentum-dependent $(n \times n)$ propagator matrix which is in general non-diagonal. The physical masses that are associated with the loop-corrected propagator matrix need to be determined from the (complex) poles of the propagator matrix.

While the case of $(2\times 2)$ mixing provides the minimal setup for studying such a propagator structure, the issue of how to associate the mass eigenstates with the $n$ poles of the propagator matrix becomes fully apparent only from the $(3\times 3)$ case onwards. The case of $(3\times 3)$ mixing captures all relevant aspects of the pole structure and permutation of the states so that the conclusions can be generalised for larger mixing systems.

For illustration and definiteness we use in this paper the language and the formalism of the neutral Higgs sector of the MSSM with ${\mathcal{CP}}$-violating mixing, i.e. for the three lowest-order mass eigenstates $i,j,k = h,H,A$ and their mixing into the loop-corrected mass eigenstates $h_a,~a=1,2,3$. Since the analytic discussions in Chapters \[sect:LoopMix\]-\[sect:fullBWZ\] do not rely on any model-dependent relations of the MSSM, the results can be transferred also to scalar sectors of other models by replacing $h,H,A$ by a different set of $i,j,k$, and analogously $h_a$ by a general $X_a$, where the generalisation to the case of mixing among more than three particles should be straightforward. Upon incorporation of the appropriate spin structure, the results can also be generalised to the case of mixing propagators of vector bosons or fermions.

The Higgs propagators in the MSSM in fact receive contributions from the mixing with the longitudinal components of the neutral gauge bosons $Z$ and $\gamma$. While the $Z$-boson propagator contributes to the diagonal elements of the Higgs propagator matrix at the two-loop level, the photon enters in those elements only at the four-loop level. Besides the mixing with the longitudinal component of the $Z$ boson, also the mixing with the unphysical (would-be) Goldstone boson $G$ needs to be taken into account. Instead of considering a $(6\times6)$ propagator matrix of the states $\left\lbrace h,H,A,G,Z,\gamma\right\rbrace$, we focus on the full mixing contributions of the physical Higgs fields, i.e. we treat the case of a $(3
\times 3)$ propagator matrix (which would generalise to $(n \times n)$ for the case of $n$ physical scalar fields). For the $(3 \times 3)$ propagator matrix of physical Higgs fields, higher-order effects from the inversion of the matrix of the irreducible 2-point vertex functions are taken into account. In contrast, we treat the mixing contributions with the gauge bosons and the (would-be) Goldstone boson in a strict perturbative expansion up to the desired order. It should be noted that Higgs–$G/Z$ mixing contributions already appear at the one-loop level in processes with external Higgs bosons, see e.g.  Refs. [@Williams:2007dc; @Fowler:2009ay; @Williams:2011bu; @Bharucha:2012nx], and therefore need to be included in order to obtain a complete one-loop result for processes of this kind. For contributions to the propagator matrix with an incoming and outgoing Higgs boson the mixing contributions with $G$ and $Z$ enter from the order of $({\hat{\Sigma}}_{i,G/Z})^{2}$ onwards, which is of subleading two-loop order. For our numerical analysis below, see Sect.\[sect:NumBWfull\], we will use use the program `FeynHiggs`[@Heinemeyer:1998np; @Heinemeyer:1998yj; @Degrassi:2002fi; @Heinemeyer:2007aq], where dominant two-loop corrections and leading higher-order contributions are incorporated in the irreducible self-energies. Since sub-leading two-loop contributions that are of the same order as the mixing contributions with $G$ and $Z$[@Frank:2006yh] are neglected in this analysis, we include no further corrections beyond the $(3 \times 3)$ propagator matrix of the physical Higgs fields in this case.

In this section we introduce the relevant quantities and fix our notation for the propagator matrix and the wave function normalisation factors, see Refs.[@Dabelstein:1994hb; @Dabelstein:1995js; @Heinemeyer:2000fa; @Heinemeyer:2001iy; @Hahn:2002gm; @Frank:2006yh; @Hahn:2006np; @Williams:2007dc; @Fowler:2010eba; @Williams:2011bu]. In Sect.\[sect:poles\] we will analyse the pole structure of propagator matrices of unstable particles.

Propagator matrix and the effective self-energy
-----------------------------------------------

As explained above, we use the case of the MSSM with complex parameters in order to illustrate propagator mixing between three physical scalar fields. If ${\mathcal{CP}}$ were assumed to be conserved in the MSSM, the ${\mathcal{CP}}$-violating self-energies would vanish, ${\hat{\Sigma}}_{hA}={\hat{\Sigma}}_{HA}=0$, so that only the two ${\mathcal{CP}}$-even states $h$ and $H$ would mix with each other. Our treatment corresponds to the case where non-zero phases from complex parameters are taken into account. Hence all renormalised self-energies ${\hat{\Sigma}_{ij}}{(p^{2})}$ of the Higgs bosons $i,j = h,H,A$ are in general non-vanishing, so that the matrix ${\textbf{M}}$ of mass squares consists of the tree-level masses $m_i^{2}$ on the diagonal and renormalised self-energies on the diagonal and off-diagonal entries. Expressed in terms of the lowest-order mass eigenstates $h, H, A$, which are also ${\mathcal{CP}}$ eigenstates, the matrix takes the form $${\textbf{M}}({p^{2}}) = {\begin{pmatrix}}m_h^{2}-\hat{\Sigma}_{hh}({p^{2}}) & -\hat{\Sigma}_{hH}({p^{2}}) &-\hat{\Sigma}_{hA}({p^{2}})\\
  -\hat{\Sigma}_{Hh}({p^{2}}) & m_H^{2}-\hat{\Sigma}_{HH}({p^{2}}) &-\hat{\Sigma}_{HA}({p^{2}})\\
  -\hat{\Sigma}_{Ah}({p^{2}}) & -\hat{\Sigma}_{AH}({p^{2}}) &m_A^{2}-\hat{\Sigma}_{AA}({p^{2}})
{\end{pmatrix}}. \label{eq:HiggsMassMatrix}$$ The renormalised irreducible 2-point vertex functions $$\hat{\Gamma}_{ij}{(p^{2})}= i\left[({p^{2}}-m_i^{2})\delta_{ij}+{\hat{\Sigma}}_{ij}{(p^{2})}\right] \label{eq:IrepVert}$$ can be collected in the $3\times3$ matrix ${\hat{\boldsymbol{\Gamma}}}_{hHA}$ in terms of ${\textbf{M}}$ as $${\hat{\boldsymbol{\Gamma}}}_{hHA}({p^{2}}) = i\left[{p^{2}}\textbf{1} -{\textbf{M}}({p^{2}}) \right] \label{eq:Gamma}.$$ Finally, the propagator matrix ${\boldsymbol{\Delta}}_{hHA}$ equals, up to the sign, the inverse of ${\hat{\boldsymbol{\Gamma}}}_{hHA}$, $${\boldsymbol{\Delta}}_{hHA}({p^{2}}) = - \left[{\hat{\boldsymbol{\Gamma}}}_{hHA}({p^{2}}) \right]^{-1}. \label{eq:DeltaGamInv}$$ Accordingly, the matrix inversion yields the individual propagators $\Delta_{ij}{(p^{2})}$ as the the $(ij)$ elements of the $3\times3$ matrix ${\boldsymbol{\Delta}}_{hHA}({p^{2}})$, $${\boldsymbol{\Delta}}_{hHA} = {\begin{pmatrix}}\Delta_{hh} &\Delta_{hH}    &\Delta_{hA}\\
            \Delta_{Hh} &\Delta_{HH}    &\Delta_{HA}\\
            \Delta_{Ah} &\Delta_{AH}    &\Delta_{AA} 
 {\end{pmatrix}}\label{eq:DeltaMat}.$$ The off-diagonal entries (for $i\neq j$) result in: $$\Delta_{ij}({p^{2}}) = \frac{{\hat{\Gamma}_}{ij}{\hat{\Gamma}_}{kk}- {\hat{\Gamma}_}{jk}{\hat{\Gamma}_}{ki}}{{\hat{\Gamma}_}{ii}{\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{kk} +2{\hat{\Gamma}_}{ij}{\hat{\Gamma}_}{jk}{\hat{\Gamma}_}{ki}-{\hat{\Gamma}_}{ii}{\hat{\Gamma}_}{jk}^{2} -{\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{ki}^{2}-{\hat{\Gamma}_}{kk}{\hat{\Gamma}_}{ij}^{2}} \label{eq:Dij}.$$ All 2-point vertex functions $\hat{\Gamma}{(p^{2})}$ depend on ${p^{2}}$ via Eq.(\[eq:IrepVert\]). Here we do not write the ${p^{2}}$-dependence explicitly for the purpose of a simpler notation, but the full ${p^{2}}$-dependence is implied also below. The solutions of the diagonal propagators, $\Delta_{ii}$, can be expressed in the following compact way: $$\begin{aligned}
 \Delta_{ii}{(p^{2})}&= \frac{{\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{kk}-{\hat{\Gamma}_}{jk}^{2}}{-{\hat{\Gamma}_}{ii}{\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{kk} + {\hat{\Gamma}_}{ii}{\hat{\Gamma}_}{jk}^{2} - 2{\hat{\Gamma}_}{ij}{\hat{\Gamma}_}{jk}{\hat{\Gamma}_}{ki} + {\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{ki}^{2} +{\hat{\Gamma}_}{kk}{\hat{\Gamma}_}{ij}^{2}}\label{eq:Dii}\\
&= \frac{i}{{p^{2}}-m_i^{2} + {\hat{\Sigma}^{\rm{eff}}_{ii}}{(p^{2})}},\label{eq:Diieff}\end{aligned}$$ where the effective self-energy is introduced, $$\begin{aligned}
   {\hat{\Sigma}^{\rm{eff}}_{ii}}(p^{2}) &= \hat{\Sigma}_{ii}(p^{2}) - i \frac{2{\hat{\Gamma}_}{ij}(p^{2}){\hat{\Gamma}_}{jk}(p^{2}){\hat{\Gamma}_}{ki}(p^{2}) - {\hat{\Gamma}_}{ki}^{2}({p^{2}}){\hat{\Gamma}_}{jj}({p^{2}})-{\hat{\Gamma}_}{ij}^{2}({p^{2}}){\hat{\Gamma}_}{kk}({p^{2}})}{{\hat{\Gamma}_}{jj}({p^{2}}){\hat{\Gamma}_}{kk}({p^{2}})-{\hat{\Gamma}_}{jk}^{2}({p^{2}})} \label{eq:seff}.\end{aligned}$$ It contains the diagonal self-energy, ${\hat{\Sigma}_{ii}}$ (which exists already at 1-loop order), and the mixing 2-point functions (whose products only contribute to ${\hat{\Sigma}^{\rm{eff}}_{ii}}$ from 2-loop order onwards). Hence, replacing the pure self-energy $\hat{\Sigma}_{ii}$ by the effective one, ${\hat{\Sigma}^{\rm{eff}}_{ii}}$, includes also the $3\times 3$ mixing contributions to the diagonal propagator in Eq.(\[eq:Diieff\]) while preserving formally the structure of the propagator as in the unmixed case. In the limit of no mixing, the second term in Eq.(\[eq:seff\]) vanishes.

The composition of ${\hat{\Sigma}^{\rm{eff}}_{ii}}$ in terms of the unmixed and the mixing contributions in Eq.(\[eq:seff\]) can be written in the alternative way [@KarinaPhD; @Fowler:2010eba; @Fuchs:2015jwa] $$\begin{aligned}
 {\hat{\Sigma}^{\rm{eff}}_{ii}}{(p^{2})}&= {\hat{\Sigma}_{ii}}{(p^{2})}+ \frac{\Delta_{ij}{(p^{2})}}{\Delta_{ii}{(p^{2})}}{\hat{\Sigma}_{ij}}{(p^{2})}+ \frac{\Delta_{ik}{(p^{2})}}{\Delta_{ii}{(p^{2})}}{\hat{\Sigma}_{ik}}{(p^{2})}\label{eq:seffratio}\end{aligned}$$ for $i,j,k$ all different. Eq.(\[eq:seffratio\]) represents the sum of the diagonal and off-diagonal self-energies involving a Higgs boson $i$ where the off-diagonal contributions are weighted by the ratios of the respective off-diagonal and the diagonal propagators. This expression follows from Eq.(\[eq:seff\]) with the help of the equality $$\begin{aligned}
 \frac{\Delta_{ij}}{\Delta_{ii}} &= -\frac{{\hat{\Gamma}_}{ij}{\hat{\Gamma}_}{kk}-{\hat{\Gamma}_}{jk}{\hat{\Gamma}_}{ki}}{{\hat{\Gamma}_}{jj}{\hat{\Gamma}_}{kk}-{\hat{\Gamma}_}{jk}^{2}}
 = -\frac{{\hat{\Sigma}_{ij}}\,(D_k+{\hat{\Sigma}_{kk}})-{\hat{\Sigma}_{jk}}{\hat{\Sigma}_{ki}}}{(D_j+{\hat{\Sigma}_{jj}})\,(D_k+{\hat{\Sigma}_{kk}})-{\hat{\Sigma}_{jk}}^{2}}
 ,\label{eq:ratioijii}\end{aligned}$$ where the shorthand $$\begin{aligned}
 D_i{(p^{2})}&= {p^{2}}- m_i^{2} \label{eq:Di}\end{aligned}$$ has been used, analogously for $j\leftrightarrow k$, and ${\hat{\Sigma}}_{ij}=-i{\hat{\Gamma}_}{ij}$ from Eq.(\[eq:IrepVert\]) with $i\neq j$ (for the special case of $(2\times 2)$ mixing, this relation can be directly read off from Eqs.(\[eq:Dii2x2\]-\[eq:seff2x2\]) given below).

On-shell wave function normalisation factors {#sect:onshellZ}
--------------------------------------------

### Complex poles

The complex poles ${\mathcal{M}^{2}}_{a}$, where $a = 1, 2, 3$, are determined as solutions of the equations $$\begin{aligned}
{p^{2}}-m_i^{2} + {\hat{\Sigma}^{\rm{eff}}_{ii}}{(p^{2})}= 0 ,\label{eq:complexpole}\end{aligned}$$ with $i = h, H, A$, see Sect. \[sect:poles\] below. While the propagators have poles for ${p^{2}}={\mathcal{M}^{2}}_a$, it is interesting to note that the ratios of off-diagonal and diagonal propagators $$R_{ij}^{(a)} := \frac{\Delta_{ij}{(p^{2})}}{\Delta_{ii}{(p^{2})}}\bigg\rvert_{{p^{2}}={\mathcal{M}^{2}}_a},\label{eq:Rija}$$ stay finite at the complex poles ${\mathcal{M}^{2}}_a,\,~a=1,2,3$.

### Z-matrix for on-shell properties of external particles {#sect:Zmatrix}

Higgs bosons appearing as *external* particles in a process need the appropriate on-shell properties for a correct normalisation of the S-matrix. In the hybrid on-shell/$\overline{\rm{DR}}$ renormalisation scheme [@Frank:2006yh], the masses are renormalised on-shell, but the $\overline{\textrm{DR}}$ renormalisation conditions for the fields and $\tan\beta$ [@Frank:2006yh], $$\begin{aligned}
 \delta Z_{\mathcal{H}_1}^{\overline{\rm{DR}}} &=& -\text{Re}\left[\Sigma^{'(\rm{div})}_{HH}(m_{H}^{2})\right]_{\alpha=0}\label{eq:dZH1},\\
 \delta Z_{\mathcal{H}_2}^{\overline{\rm{DR}}} &=& -\text{Re}\left[\Sigma^{'(\rm{div})}_{hh}(m_{h}^{2})\right]_{\alpha=0}\label{eq:dZH2},\\
 \delta \tan\beta^{\overline{\rm{DR}}} &=& \frac{1}{2}(\delta Z_{\mathcal{H}_2}^{\overline{\rm{DR}}} - \delta Z_{\mathcal{H}_1}^{\overline{\rm{DR}}})\label{eq:dtanb},\end{aligned}$$ where $\alpha$ is the mixing angle between $h$ and $H$, see Eq.(\[eq:ncHmix\]) below, do not ensure proper on-shell properties of the Higgs bosons. In fact, the loop-corrected mass eigenstates $h_1,h_2,h_3$, which occur as external, on-shell, particles e.g. in decay processes, are a mixture of the lowest order states $h,H,A$. Thus, finite wave function normalisation factors need to be introduced in order to guarantee that the mixing vanishes on-shell and that the propagators of the external particles have unit residue.

These so-called $Z$-factors for a neutral Higgs boson $i=h,H,A$ on an external line are obtained from the residue of the propagators at the complex pole ${\mathcal{M}^{2}}_a,\,a=1,2,3$,[@Dabelstein:1995js; @Heinemeyer:2001iy] $$\begin{aligned}
 {\hat{Z}}_i^{a} &:= \textrm{Res}_{{\mathcal{M}^{2}}_a}\left\lbrace\Delta_{ii}{(p^{2})}\right\rbrace   \label{eq:ZhaRes}.
 \end{aligned}$$ Expanding ${\hat{\Sigma}^{\rm{eff}}_{ii}}({p^{2}})$ around the complex pole ${p^{2}}={\mathcal{M}^{2}}_a$, one obtains the diagonal propagator $$\begin{aligned}
  \Delta_{ii}(p^{2}) &=\frac{i}{p^{2}-m_i^{2}+{\hat{\Sigma}^{\rm{eff}}_{ii}}(p^{2})}
=\frac{i}{p^{2}-{\mathcal{M}^{2}}_a}\cdot
\frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a)+ \mathcal{O}({p^{2}}- {\mathcal{M}^{2}}_a)}
\label{eq:DiiApprdSeff}\,.\end{aligned}$$ Performing the limit ${p^{2}}\rightarrow {\mathcal{M}^{2}}_a$ yields the residue of Eq.(\[eq:ZhaRes\]), $$\begin{aligned}
 {\hat{Z}}_i^{a}
 &= \frac{1}{\frac{\partial}{\partial {p^{2}}}\frac{i}{\Delta_{ii}{(p^{2})}}}\bigg\rvert_{{p^{2}}={\mathcal{M}^{2}}_a}
 =\frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a)} \label{eq:Zha}\,.\end{aligned}$$ Considering a diagram with the Higgs boson $i$ on an external line, whose propagator has three poles (see Sect. \[sect:poles\]), there are three possibilities which residue to compute. If the amputated Green’s function is evaluated at ${\mathcal{M}^{2}}_a$, the external $i$-line has to be multiplied by $\sqrt{{\hat{Z}}_i^{a}}$ for the correct S-matrix normalisation[^1]. So the resulting mass eigenstate as an outgoing particle is $h_a$. Alternatively, if the Green’s function is evaluated at ${\mathcal{M}^{2}}_b$, it has to be normalised by $$\begin{aligned}
\sqrt{{\hat{Z}}_i^{b}}=\frac{1}{\sqrt{1+{\hat{\Sigma}^{\rm{eff'}}}_{ii}({\mathcal{M}^{2}}_b)}}\end{aligned}$$ to achieve the correct S-matrix element. For this choice, the external mass eigenstate is $h_b$. Moreover, the wave function normalisation factor for $i-j$ mixing on an external on-shell line at ${\mathcal{M}^{2}}_a$ is composed of the overall normalisation factor $\sqrt{{\hat{Z}}_i^{a}}$ times the on-shell transition ratio $${\hat{Z}}_{ij}^{a} \equiv R_{ij}^{(a)} = \frac{\Delta_{ij}{(p^{2})}}{\Delta_{ii}{(p^{2})}}\bigg\rvert_{{p^{2}}={\mathcal{M}^{2}}_a} \label{eq:transij}\,.$$ Since $\Delta_{ii}$ and $\Delta_{ij}$ have – in the case of $3\times 3$ mixing – 3 complex poles (see Sect. \[sect:poles\]), any of them can be chosen for the evaluation of ${\hat{Z}}_i$, ${\hat{Z}}_{ij}$ and ${\hat{Z}}_{ik}$; in the considered example this is ${\mathcal{M}^{2}}_a$. Correspondingly, ${\hat{Z}}_j, {\hat{Z}}_{ji}$ and ${\hat{Z}}_{jk}$ will be evaluated at ${\mathcal{M}^{2}}_b$ and ${\hat{Z}}_k, {\hat{Z}}_{ki}$ and ${\hat{Z}}_{kj}$ at ${\mathcal{M}^{2}}_c$ where $a,b,c$ are a permutation of 1,2,3, and $i,j,k$ a permutation of $h,H,A$[@Williams:2007dc; @Fuchs:2015jwa]. For $2\times 2$ mixing, of course only the two indices involved in the mixing can be permuted.

All choices allowed by the mixing structure are generally possible because of the pole structure of each propagator. However, they might not be equally numerically stable. If the loop-corrected mass eigenstate $h_a$ contains only a small admixture of the lowest-order state $i$, the propagator still has a pole at ${\mathcal{M}^{2}}_a$, but the contribution of $i$ to $h_a$ is suppressed at ${p^{2}}\neq {\mathcal{M}^{2}}_a$.

In order to be definite, it is in any case necessary to define at which pole to evaluate which normalisation and mixing ${\hat{Z}}$-factor. This choice corresponds to fixing an assignment between the indices $i,j,k$ of the lowest-order states and the indices $a,b,c$ of the higher-order mixed states and then using it consistently. The assignment $(i,a), (j,b), (k,c)$, which we label as scheme $I$, prescribes to evaluate ${\hat{Z}}_i$, ${\hat{Z}}_{ij}$ and ${\hat{Z}}_{ik}$ at ${\mathcal{M}^{2}}_a$. Once the indices have been assigned we can clear up the notation by writing $${\hat{Z}}_a\big\rvert_I := {\hat{Z}}_{i}^{a},~~ {\hat{Z}}_{aj}\big\rvert_I := {\hat{Z}}_{ij}^{a},~~ {\hat{Z}}_{bi}\big\rvert_I := {\hat{Z}}_{ji}^{b}\,,\label{eq:Zidentification}$$ and accordingly for the other indices such that the first index always refers to a loop-corrected mass eigenstate ($a,b,c\in\{1,2,3\}$)[^2] and the second index to a lowest-order state ($i,j,k\in\{h,H,A\}$). Note that ${\hat{Z}}_{ai}={\hat{Z}}_{bj}={\hat{Z}}_{ck}\equiv1$ in the index scheme $I$ defined above. Once the index scheme has been specified, one can leave out the subscript $I$. For the scheme-independence of physical results, see Sect.\[sect:schemeIndep\].

Furthermore, it is convenient[@Frank:2006yh] to arrange the products of the normalisation factors $\sqrt{{\hat{Z}}_a}$ and transition ratios ${\hat{Z}}_{aj}$ as $$\begin{aligned}
{\hat{\textbf{Z}}}_{aj} &= \sqrt{{\hat{Z}}_a}{\hat{Z}}_{aj}
\label{eq:ZiZij}\end{aligned}$$ (note the difference between ${\hat{Z}}_{aj}$ and ${\hat{\textbf{Z}}}_{aj}$) into a non-unitary matrix: $$\begin{aligned}
 \hat{\textbf{Z}} &= {\begin{pmatrix}}\sqrt{\hat{Z}_1}{\hat{Z}}_{1h} & \sqrt{\hat{Z}_1}\hat{Z}_{1H} &\sqrt{\hat{Z}_1}\hat{Z}_{1A}\\
\sqrt{\hat{Z}_2}\hat{Z}_{2h} & \sqrt{\hat{Z}_2}{\hat{Z}}_{2H} &\sqrt{\hat{Z}_2}\hat{Z}_{2A}\\
\sqrt{\hat{Z}_3}\hat{Z}_{3h} & \sqrt{\hat{Z}_3}\hat{Z}_{3H} &\sqrt{\hat{Z}_3}{\hat{Z}}_{3A}
{\end{pmatrix}}\label{eq:Zmatrix}.\end{aligned}$$ The ${\hat{\textbf{Z}}}$-matrix defined in Eq.(\[eq:Zmatrix\]) fulfils the on-shell conditions (unit residue and vanishing mixing), which can be written in the following compact form[@Williams:2007dc; @KarinaPhD; @Fowler:2010eba; @Williams:2011bu]: $$\begin{aligned}
 \lim_{{p^{2}}\rightarrow {\mathcal{M}^{2}}_a} -\frac{i}{{p^{2}}-{\mathcal{M}^{2}}_a}\left({\hat{\textbf{Z}}}\cdot 
{\hat{\boldsymbol{\Gamma}}}_{hHA}\cdot {\hat{\textbf{Z}}}^{T} \right)_{hh} ~&= 1
\label{eq:Reshha},\\
 \lim_{{p^{2}}\rightarrow {\mathcal{M}^{2}}_b} -\frac{i}{{p^{2}}-{\mathcal{M}^{2}}_b}\left({\hat{\textbf{Z}}}\cdot 
{\hat{\boldsymbol{\Gamma}}}_{hHA}\cdot 
{\hat{\textbf{Z}}}^{T} \right)_{HH} &= 1
\label{eq:ResHHb},\\
 \lim_{{p^{2}}\rightarrow {\mathcal{M}^{2}}_c} -\frac{i}{{p^{2}}-{\mathcal{M}^{2}}_c}\left({\hat{\textbf{Z}}}\cdot 
{\hat{\boldsymbol{\Gamma}}}_{hHA}\cdot {\hat{\textbf{Z}}}^{T} \right)_{AA} &= 1
\label{eq:ResAAc},\end{aligned}$$ with vanishing off-diagonal entries. It is equally possible to begin with these equations (\[eq:Reshha\]-\[eq:ResAAc\]), i.e. to require unit residues and to demand vanishing mixing on-shell, to derive the elements of the ${\hat{\textbf{Z}}}$-matrix whose solutions are given in Eqs.(\[eq:Zha\]) and (\[eq:transij\]).

In Ref.[@Williams:2007dc] the evaluation at the full complex poles and the inclusion of imaginary parts were introduced (see Eq.(\[eq:SigijReIm\]) below for the treatment of loop functions of complex arguments). The evaluation at the complex poles leads to numerically more stable results than the evaluation at real ${p^{2}}=M_{h_a}^{2}$, as well as to the physically equivalent choices of index assignments. The calculation of the ${\hat{\textbf{Z}}}$-factors of MSSM Higgs bosons can be performed with the program `FeynHiggs`.

Use of Z-factors for external states {#sect:ExtH}
------------------------------------

The ${\hat{\textbf{Z}}}$-factors have been introduced for the correct normalisation of matrix elements with external (on-shell) Higgs bosons $h_a,~{p^{2}}={\mathcal{M}^{2}}_a$. It should be noted that ${\hat{\textbf{Z}}}$ does not provide a unitary transformation between the basis of lowest-order states and the basis of loop-corrected mass eigenstates. The fact that ${\hat{\textbf{Z}}}$ is a non-unitary matrix is related to the imaginary parts appearing in the propagators of unstable particles. Using the ${\hat{\textbf{Z}}}$-matrix, one can express the one-particle irreducible (1PI) vertex functions ${\hat{\Gamma}_}{h_a}$ involving a loop-corrected mass eigenstate $h_1,h_2,h_3$ as an external particle as a linear combination of the 1PI vertex functions of the lowest-order states, ${\hat{\Gamma}_}i$: $$\begin{aligned}
 {\hat{\Gamma}_}{h_a}&={\hat{\textbf{Z}}}_{ah}{\hat{\Gamma}_}h+{\hat{\textbf{Z}}}_{aH}{\hat{\Gamma}_}H+{\hat{\textbf{Z}}}_{aA}{\hat{\Gamma}_}A + ...\label{eq:Zextsum}\\
 &=\sqrt{{\hat{Z}}_a}\,\left({\hat{Z}}_{ah}{\hat{\Gamma}_}h+{\hat{Z}}_{aH}{\hat{\Gamma}_}H+{\hat{Z}}_{aA}{\hat{\Gamma}_}A\right) + \dots \label{eq:Vertexha},\end{aligned}$$ where the ellipsis refers to additional terms arising from the mixing with Goldstone and vector bosons, which are not described by the ${\hat{\textbf{Z}}}$-matrix. Thus, the overall normalisation factor $\sqrt{{\hat{Z}}_a}$ accounts for the particle $h_a$ appearing at an external line. In addition, the factors ${\hat{Z}}_{ai}$ given in Eqs.(\[eq:transij\]) and (\[eq:Zidentification\]) as ratios of propagators at ${p^{2}}={\mathcal{M}^{2}}_a$ describe the transition between the states $h_a$ and $i$. The transition factor ${\hat{Z}}_{ai}$ occurs in a diagram where $h_a$ is the external particle, but $i$ directly couples to the vertex. All possibilities for $i=h,H,A$ need to be included for each $h_a$, hence the sum arises. This is depicted in Fig.\[fig:Zext\] (cf. also Refs.[@Hahn:2007fq; @KarinaPhD]).

at (2.5,-0.6) [${p^{2}}={\mathcal{M}^{2}}_a$]{}; (1,0)–(4,0); at (2.5,0.6) [$\boldsymbol{h_a}$]{}; (4.3,0) circle (0.3); (5.6,1)–(4.6,0.1); (5.6,-1)–(4.6,-0.1); at (5.5,0) [$\hat\Gamma_{{\color{blue} h_a}}$]{}; at (7.2,0) [= $\sqrt{{\hat{Z}}_{{\color{blue} a}}}\Bigg($]{};

at (2.5,-0.8) [${\hat{Z}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; (1,0)–(2.2,0); (2.8,0)–(4,0); at (1.6,0.6) [$\boldsymbol{h_a}$]{}; at (3.4,0.6) [$\boldsymbol{h}$]{}; (2.5,0) circle (0.3); (5,1)–(4,0)–(5,-1); at (5,0) [$\hat\Gamma_{{\color{OliveGreen} h}}$]{}; at (6,0) [+]{};

at (2.5,-0.8) [${\hat{Z}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; (1,0)–(2.2,0); (2.8,0)–(4,0); at (1.6,0.6) [$\boldsymbol{h_a}$]{}; at (3.4,0.6) [$\boldsymbol{H}$]{}; (2.5,0) circle (0.3); (5,1)–(4,0)–(5,-1); at (5,0) [$\hat\Gamma_{{\color{OliveGreen} H}}$]{}; at (6,0) [+]{};

at (2.5,-0.8) [${\hat{Z}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; (1,0)–(2.2,0); (2.8,0)–(4,0); at (1.6,0.6) [$\boldsymbol{h_a}$]{}; at (3.4,0.6) [$\boldsymbol{A}$]{}; (2.5,0) circle (0.3); (5,1)–(4,0)–(5,-1); at (5,0) [$\hat\Gamma_{{\color{OliveGreen} A}}$]{};

at (25,0) [$\Bigg)_{{p^{2}}={\mathcal{M}^{2}}_a}$]{}; at (26.4,0) [$+\dots$]{};

Conveniently, Eq.(\[eq:Vertexha\]) can be written in matrix form for all $h_1,h_2,h_3$ as $${\begin{pmatrix}}{\hat{\Gamma}_}{h_1}\\{\hat{\Gamma}_}{h_2}\\{\hat{\Gamma}_}{h_3}{\end{pmatrix}}= {\hat{\textbf{Z}}}\cdot {\begin{pmatrix}}{\hat{\Gamma}_}{h}\\{\hat{\Gamma}_}{H}\\{\hat{\Gamma}_}{A}{\end{pmatrix}}+\dots\label{eq:vertZ}~.$$ In this way, propagator corrections at external legs are effectively absorbed into the vertices of neutral Higgs bosons. If ${\hat{\textbf{Z}}}$-factors are applied to supplement the Born result, only non-Higgs propagator type corrections (such as mixing with the Goldstone and $Z$-bosons) as well as vertex, box and real corrections need to be calculated individually.

As it is the case for the usual (UV-divergent) on-shell field renormalisation constants, the finite wave function normalisation factors are in general gauge-dependent quantities. While the gauge-parameter dependence drops out at a fixed order in perturbation theory, the inclusion of higher-order contributions into the ${\hat{\textbf{Z}}}$-factors described above can give rise to a residual gauge-parameter dependence, which is formally of sub-leading higher order. This is caused by the fact that the ${\hat{\textbf{Z}}}$-factors contain not only gauge-parameter independent leading higher-order contributions but also sub-leading higher-order contributions that are in general gauge-parameter dependent. Those sub-leading higher-order contributions would cancel with corresponding sub-leading higher-order contributions from other Green functions such as vertex and box contributions.

Effective couplings {#sect:Ueff}
-------------------

Since the ${\hat{\textbf{Z}}}$-matrix is not unitary, it does not represent a unitary transformation between the $\{h,H,A\}$ and the $\{h_1,h_2,h_3\}$ basis. However, it is not necessary to diagonalise the mass matrix for the determination of the poles of the propagators. Hence there is a priori no need to introduce a unitary transformation. Though, if a unitary matrix ${\textbf{U}}$ is desired for the definition of effective couplings, an approximation of the momentum dependence of ${\hat{\textbf{Z}}}$ is required. There is no unique prescription of how to achieve a unitary mixing matrix as an approximation of the ${\hat{\textbf{Z}}}$-matrix, but a possible choice is the ${p^{2}}=0$ approximation[@Frank:2006yh; @Hahn:2007fq]. As in the effective potential approach, the external momentum ${p^{2}}$ is set to zero in the renormalised self-energies ${\hat{\Sigma}_{ij}}({p^{2}})\rightarrow {\hat{\Sigma}_{ij}}(0)$ so that they become real. Then ${\textbf{U}}$ diagonalises the real matrix ${\textbf{M}}(0)$. ${\textbf{U}}$ can be chosen real and it transforms the ${\mathcal{CP}}$-eigenstates into the mass eigenstates, $$\begin{aligned}
 {\begin{pmatrix}}h_1\\ h_2\\ h_3 {\end{pmatrix}}= {\textbf{U}}{\begin{pmatrix}}h\\ H\\ A {\end{pmatrix}},~~~~
 {\textbf{U}}= {\begin{pmatrix}}{\textbf{U}}_{1h}& {\textbf{U}}_{1H}& {\textbf{U}}_{1A}\\
        {\textbf{U}}_{2h}& {\textbf{U}}_{2H}& {\textbf{U}}_{2A}\\
        {\textbf{U}}_{3h}& {\textbf{U}}_{3H}& {\textbf{U}}_{3A}
    {\end{pmatrix}}\label{eq:Urotate},\end{aligned}$$ so that $U_{aA}^{2}$ quantifies the admixture of a ${\mathcal{CP}}$-odd component inside $h_a$[@Frank:2006yh]. The elements of ${\textbf{U}}$ can then be used to introduce effective couplings of the loop-corrected states $h_a$ to any other particles $X$ in terms of the couplings of the unmixed states $i$ by the relation $$C_{h_a X}^{U} = \sum_{i=h,H,A} {\textbf{U}}_{ai} C_{i X} \label{eq:EffCoup}.$$ absorbing some higher-order corrections, but neglecting imaginary parts and the full momentum dependence of the self-energies. Hence, the application of ${\textbf{U}}$ resembles the use of ${\hat{\textbf{Z}}}$-factors in Eq.(\[eq:Zextsum\]) for the purpose of implementing partial higher-order effects into an improved Born result. Yet, the rotation matrix ${\textbf{U}}$ introduced for effective couplings as a unitary approximation is conceptually different from the ${\hat{\textbf{Z}}}$-matrix arising from propagator corrections and introduced for the correct normalisation of the $S$-matrix.

Technical treatment of imaginary parts
--------------------------------------

In order to account for complex momenta and imaginary parts of self-energies, we employ an expansion of the self-energies around the real part of the complex momentum, $$\begin{aligned}
 {p^{2}}&\equiv {p^{2}}_r + i {p^{2}}_i \label{eq:psqReIm},\\
 {\hat{\Sigma}_{ij}}({p^{2}}) &\simeq {\hat{\Sigma}_{ij}}({p^{2}}_r) + i {p^{2}}_i {\hat{\Sigma}^{'}}_{ij}({p^{2}}_r)\label{eq:SigijReIm},\end{aligned}$$ where ${\hat{\Sigma}^{'}}_{ij}{(p^{2})}\equiv \frac{d{\hat{\Sigma}}_{ij}{(p^{2})}}{d{p^{2}}}$. This expansion enables the calculation of self-energies at complex momentum in terms of self-energies evaluated at real momentum (e.g. by `FeynHiggs`). For the inclusion of all products of real and imaginary parts, we do not expand the effective self-energy from Eq.(\[eq:seff\]) directly according to Eq.(\[eq:psqReIm\]). Instead, in the same way as in Refs.[@KarinaPhD; @Fowler:2010eba; @Fuchs:2015jwa], we expand all $\hat{\Gamma}_{ij}{(p^{2})}$ in Eq.(\[eq:seff\]) individually before combining them into ${\hat{\Sigma}^{\rm{eff}}_{ii}}$.

Pole structure of propagator matrices {#sect:poles}
=====================================

In this section, we will first discuss in detail the solutions of the poles of the propagator matrix depending on the mixing. Subsequently we will prove the equivalence of different assignments between the lowest-order states and the loop-corrected mass eigenstates for the physical ${\hat{\textbf{Z}}}$-matrix.

Determination of the poles
--------------------------

Imaginary parts of the self-energies lead to propagator poles in the complex momentum plane. The Higgs masses are determined from the real part of the complex poles ${\mathcal{M}^{2}}$ of the diagonal propagators. Equivalently, the complex poles are obtained as the zeros of the inverse diagonal propagators. Due to $\boldsymbol{A}^{-1}\propto \det\left[\boldsymbol{A}\right]^{-1}$, $\det\left[\boldsymbol{A}^{-1}\right]=\det\left[\boldsymbol{A}\right]^{-1}$ and $\det\left[-\boldsymbol{A}\right]=(-1)^n\,\det\left[\boldsymbol{A}\right]$ for any $n\times n$-matrix $\boldsymbol{A}$, finding the roots of ${\boldsymbol{\Delta}}^{-1}{(p^{2})}$ is equivalent to solving $$\begin{aligned}
-\frac{1}{\det\left[{\boldsymbol{\Delta}}_{hHA}{(p^{2})}\right]}
&=\det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}\right]   \stackrel{!}{=} 0,\label{eq:det0}
 \end{aligned}$$ because of the relation between ${\boldsymbol{\Delta}}_{hHA}$ and ${\hat{\boldsymbol{\Gamma}}}_{hHA}$ from Eq.(\[eq:DeltaGamInv\])[^3]. Then the loop-corrected masses $M$ are obtained from the real parts of the complex poles and the total widths $\Gamma$ from the imaginary parts via $${\mathcal{M}^{2}}= M^{2} - i M\Gamma. \label{eq:MGamma}$$ In the following, the impact of higher-order and mixing contributions on the pole structure of the propagators will be discussed.

#### Lowest order

At lowest order, the self-energy contributions in Eq.(\[eq:IrepVert\]) are absent and the matrix ${\hat{\boldsymbol{\Gamma}}}$ simply reads $$\begin{aligned}
 {\hat{\boldsymbol{\Gamma}}}^{(0)}_{hHA}{(p^{2})}&= i\,\textrm{diag}\left\lbrace D_h{(p^{2})}, D_H{(p^{2})},
D_A{(p^{2})}\right\rbrace \label{eq:GamLowestOrder},\end{aligned}$$ where the shorthand of Eq.(\[eq:Di\]) has been used. The solutions of Eq.(\[eq:det0\]) are the three tree-level masses $m_i^{2}$.

#### Higher order without mixing

Beyond tree-level, the self-energies are added at the available order. Restricting them to the unmixed case, ${\hat{\Sigma}}_{ij}=0$ for $i\neq j$, leads to $$\begin{aligned}
 {\hat{\boldsymbol{\Gamma}}}^{(\textrm{no mix})}_{hHA}{(p^{2})}&= i\,\textrm{diag}\left\lbrace D_h{(p^{2})}+ {\hat{\Sigma}}_{hh}{(p^{2})}, D_H{(p^{2})}+{\hat{\Sigma}}_{HH}{(p^{2})}, D_A{(p^{2})}+{\hat{\Sigma}}_{AA}{(p^{2})}\right\rbrace \label{eq:GamNoMix},\end{aligned}$$ so that $$\det\left[{\hat{\boldsymbol{\Gamma}}}^{(\textrm{no mix})}_{hHA}{(p^{2})}\right]=\prod_{i=h,H,A}\left(D_i{(p^{2})}+{\hat{\Sigma}}_{ii}{(p^{2})}\right)=0 \label{eq:DetUnmix}$$ is achieved if ${p^{2}}$ fulfils the following on-shell relation $${p^{2}}- m_i^{2} + {\hat{\Sigma}}_{ii}{(p^{2})}= 0 \label{eq:0unmixed}$$ for any $i=h,H,A$. Thus, the full propagator matrix ${\boldsymbol{\Delta}}$ has three poles and each propagator $\Delta_{ii}{(p^{2})}$ has exactly one pole ${p^{2}}={\mathcal{M}^{2}}_i$ that solves Eq.(\[eq:0unmixed\]) in this case.

#### Higher order with 2x2 mixing

If now the mixing between $h$ and $H$ is taken into account, corresponding to the ${\mathcal{CP}}$-conserving case, the matrix ${\hat{\boldsymbol{\Gamma}}}$ becomes block-diagonal with the $2\times 2$ matrix ${\hat{\boldsymbol{\Gamma}}}_{hH}$ and the 2-point vertex function of $A$, which does not mix with the other states: $$\begin{aligned}
 {\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}&= {\begin{pmatrix}}{\hat{\boldsymbol{\Gamma}}}_{hH}{(p^{2})}&0\\
            0       &\hat{\Gamma}_A{(p^{2})}{\end{pmatrix}}, \label{eq:blockdiag}\\
 \det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}\right] &= \det\left[{\hat{\boldsymbol{\Gamma}}}_{hH}{(p^{2})}\right]\cdot\hat\Gamma_{A}{(p^{2})}\,.\end{aligned}$$ For a closer look at the relation between the roots of the determinant and the roots of the inverse propagator, we write down the propagators and the effective self-energy of the $\left\lbrace h,H\right\rbrace$ system explicitly. They follow from Eqs.(\[eq:Dij\]), (\[eq:Dii\]) and (\[eq:seff\]) by setting ${\hat{\Sigma}}_{hA}={\hat{\Sigma}}_{HA}=0$ or equivalently from the inversion of the $2\times 2$ submatrix ${\hat{\boldsymbol{\Gamma}}}_{hH}$: $$\begin{aligned}
\Delta_{ii}({p^{2}}) &= \frac{i\left[D_j{(p^{2})}+ {\hat{\Sigma}_{jj}}({p^{2}})\right]}{\left[D_i{(p^{2})}+{\hat{\Sigma}_{ii}}({p^{2}})\right] \left[D_j{(p^{2})}+{\hat{\Sigma}_{jj}}({p^{2}})\right]-{\hat{\Sigma}_{ij}}^{2}({p^{2}})}
=\frac{i}{{p^{2}}-m_i^{2}+{\hat{\Sigma}^{\rm{eff}}_{ii}}({p^{2}})},\label{eq:Dii2x2}\\
 \Delta_{ij}({p^{2}}) 
&= \frac{-i{\hat{\Sigma}_{ij}}({p^{2}})}{\left[D_i{(p^{2})}+{\hat{\Sigma}_{ii}}({p^{2}})\right]\,\left[D_j{(p^{2})}+{\hat{\Sigma}_{jj}}({p^{2}})\right]-{\hat{\Sigma}_{ij}}^{2}({p^{2}})}\label{2x2eq:Dij},\\
{\hat{\Sigma}^{\rm{eff}}_{ii}}({p^{2}}) &= {\hat{\Sigma}_{ii}}({p^{2}}) -\frac{{\hat{\Sigma}_{ij}}^{2}({p^{2}})}{D_j({p^{2}})+{\hat{\Sigma}_{jj}}({p^{2}})}\label{eq:seff2x2}.\end{aligned}$$ Comparing the inverse diagonal propagators with the determinant of the submatrix ${\hat{\boldsymbol{\Gamma}}}_{hH}$, we find for $i,j\in\left\lbrace h,H \right\rbrace,~i\neq j,$ $$\begin{aligned}
 \frac{1}{\Delta_{ii}{(p^{2})}} = \frac{i}{D_j{(p^{2})}+{\hat{\Sigma}}_{jj}{(p^{2})}}\,\det\left[{\hat{\boldsymbol{\Gamma}}}_{hH}{(p^{2})}\right]. \label{eq:InvPropDet2x2}\end{aligned}$$ Eq.(\[eq:InvPropDet2x2\]) reveals that both inverse diagonal propagators, $1/\Delta_{hh}$ and $1/\Delta_{HH}$, are proportional to the determinant of ${\hat{\boldsymbol{\Gamma}}}_{hH}$, which has two zeros. As opposed to the unmixed case, both zeros of $\det\left[{\hat{\boldsymbol{\Gamma}}}_{hH}{(p^{2})}\right]$ are poles of *each* of the diagonal propagators $\Delta_{hh},\Delta_{HH}$. The lowest-order states $h$ and $H$ are mixed into the loop-corrected mass eigenstates $h_1$ and $h_2$ with the loop-corrected masses $M_{h_1}, M_{h_2}$. The corresponding poles ${p^{2}}={\mathcal{M}^{2}}_{h_1}, {\mathcal{M}^{2}}_{h_2}$ solve $${p^{2}}- m_i^{2} + \hat{\Sigma}^{\textrm{eff}}_{ii}({p^{2}}) = 0 \label{eq:0mixed}$$ for ${p^{2}}={\mathcal{M}^{2}}_{h_a}$ in any combination of $i=h,H$ and $a=1,2$, where the effective self-energy is given in Eq.(\[eq:seff2x2\]). In the $2\times 2$ mixing system, it is convenient to choose $M_{h_1}\leq M_{h_2}$. As for the nomenclature in the $2\times 2$ case, the lighter mass eigenstate $h_1$ is often denoted as $h$ and the heavier one as $H$ because both are ${\mathcal{CP}}$-even states. It should be noted that both roots of ${\hat{\boldsymbol{\Gamma}}}_{hH}$, ${\mathcal{M}^{2}}_{h_1}$ and ${\mathcal{M}^{2}}_{h_2}$, are also complex poles of the off-diagonal propagators $\Delta_{hH}{(p^{2})}\equiv\Delta_{Hh}{(p^{2})}$ due to $$\begin{aligned}
  \frac{1}{\Delta_{ij}{(p^{2})}} = \frac{-i}{{\hat{\Sigma}}_{ij}{(p^{2})}}\,\det\left[{\hat{\boldsymbol{\Gamma}}}_{hH}{(p^{2})}\right]. \label{eq:InvPropijDet2x2}\end{aligned}$$ Since in this case $A$ does not mix with $h$ and $H$, the third pole ${\mathcal{M}^{2}}_A$ solely solves $${\mathcal{M}^{2}}_{A} - m_A^{2} + {\hat{\Sigma}}_{AA}({\mathcal{M}^{2}}_{A}) = 0 \label{eq:solveA},$$ but no other combination of $A$ and $h_a$ satisfies the on-shell condition. $M_A$ is the loop-corrected mass of the (mass and interaction) eigenstate $A$.

#### Higher order with 3x3 mixing

Now we turn to the case where complex MSSM parameters lead to ${\mathcal{CP}}$-violating self-energies ${\hat{\Sigma}}_{hA}, {\hat{\Sigma}}_{HA}$. Thus, all three neutral Higgs lowest-order and ${\mathcal{CP}}$-eigenstates $h,H,A$ mix into the loop-corrected mass eigenstates $h_1, h_2, h_3$, which have no longer well-defined ${\mathcal{CP}}$ quantum numbers, but are admixtures of ${\mathcal{CP}}$-even and ${\mathcal{CP}}$-odd components. In this framework, ${\hat{\boldsymbol{\Gamma}}}_{hHA}$ is a full $3\times 3$ matrix with the determinant $$\begin{aligned}
 \det[{\hat{\boldsymbol{\Gamma}}}_{hHA}] &= -i\left[(D_{h}+{\hat{\Sigma}}_{hh})(D_{H}+{\hat{\Sigma}}_{HH})(D_{A}+{\hat{\Sigma}}_{AA})+2{\hat{\Sigma}}_{hH}{\hat{\Sigma}}_{HA}{\hat{\Sigma}}{hA}\right.\nonumber\\
 &\left.~~~~~~~~-(D_{h}+{\hat{\Sigma}}_{hh}){\hat{\Sigma}}_{HA}^2-(D_{H}+{\hat{\Sigma}}_{HH}){\hat{\Sigma}}_{hA}^2-(D_{A}+{\hat{\Sigma}}_{AA}){\hat{\Sigma}}_{hH}^2 \right] \label{eq:Det3mix},\end{aligned}$$ where we dropped the explicit ${p^{2}}$-dependence of each term for an ease of notation. Comparing Eq.(\[eq:Det3mix\]) with the diagonal and off-diagonal propagators from Eqs.(\[eq:Dii\]) and (\[eq:Dij\]), respectively, we see that their inverse is proportional to the determinant of ${\hat{\boldsymbol{\Gamma}}}_{hHA}$: $$\begin{aligned}
 \frac{1}{\Delta_{ii}} &= \frac{\det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}\right]}{(D_{j}+{\hat{\Sigma}}_{jj})(D_{k}+{\hat{\Sigma}}_{kk})-{\hat{\Sigma}}_{jk}^{2}} \label{eq:InvPropDet},\\
 \frac{1}{\Delta_{ij}} &= \frac{\det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}\right]}{{\hat{\Sigma}}_{jk}{\hat{\Sigma}}_{ki}-{\hat{\Sigma}}_{ij}(D_k+{\hat{\Sigma}}_{kk}) } \label{eq:InvPropijDet}\,,\end{aligned}$$ where $i\neq j\neq k\neq i$ and without summation over the indices. From Eq.(\[eq:InvPropDet\]) we conclude that all three roots ${p^{2}}={\mathcal{M}^{2}}_{h_a},\,a=1,2,3,$ of $\det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}\right]$ are complex poles of *each* of the three diagonal propagators $\Delta_{ii},\,i=h,H,A$. This means that $${\mathcal{M}^{2}}_{h_a} - m_i^{2}+\hat{\Sigma}^{\rm{eff}}_{ii}({\mathcal{M}^{2}}_{h_a}) =0\label{eq:SolvePole}$$ holds for any combination of $i$ and $a$ in the presence of $3\times 3$ mixing. Moreover, Eq.(\[eq:InvPropijDet\]) implies that also the off-diagonal propagators have as many poles as the determinant has zeros, namely three in the case of ${\mathcal{CP}}$-violating mixing. In the unmixed case, ${\hat{\Sigma}^{\rm{eff}}_{ii}}= {\hat{\Sigma}_{ii}}$ and each propagator has exactly one pole so that there is a unique mapping between $i$ and $a$, see Eq.(\[eq:0unmixed\]). On the other hand, for the general mixing case it is not unique how to relate the mass eigenstates to the interaction eigenstates. An assignment will be needed for the definition of on-shell wave function normalisation factors in Sect.\[sect:onshellZ\].

Index scheme independence of the Z-matrix {#sect:schemeIndep}
-----------------------------------------

As discussed above, the pole structure of the full propagators provides the freedom at which pole to evaluate which ${\hat{Z}}$-factor, or equivalently which loop-corrected mass eigenstate index ($a,b,c$) to assign to a lowest-order mass eigenstate index ($i,j,k$). We denote two choices for such index schemes by $$\begin{aligned}
 I &\leftrightarrow  (i,a), (j,b), (k,c)\label{eq:schemeI},\\
II &\leftrightarrow  (j,a), (i,b), (k,c)\label{eq:schemeII}.\end{aligned}$$ This initial ambiguity, however, results in physically equivalent results. Using properties of ratios of diagonal and off-diagonal propagators and exploiting relations at a complex pole (for details, see Ref.[@Fuchs:2015jwa]), we are able to show that $$\begin{aligned}
 \frac{1+{\hat{\Sigma}^{\rm{eff'}}}_{jj}({\mathcal{M}^{2}}_a)}{1+{\hat{\Sigma}^{\rm{eff'}}}_{ii}({\mathcal{M}^{2}}_a)}
=\left(\frac{\Delta_{ji}{(p^{2})}}{\Delta_{jj}{(p^{2})}}\right)^{2}_{{p^{2}}={\mathcal{M}^{2}}_a}.\end{aligned}$$ This equality provides a transformation between scheme $I$ (where $i$ and $a$ are associated, hence ${\hat{Z}}_{ai}\equiv1$) and scheme $II$ (where $j$ and $a$ are matched): $$\begin{aligned}
 {\hat{\textbf{Z}}}_{ai}|_{I} = \left(\sqrt{{\hat{Z}}_a}{\hat{Z}}_{ai}\right)_{I} &= \frac{1}{\sqrt{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a)}}\nonumber\\
 &= \frac{1}{\sqrt{1+{\hat{\Sigma}^{\rm{eff'}}}_{jj}({\mathcal{M}^{2}}_a)}}\frac{\Delta_{ji}}{\Delta_{jj}}\bigg\rvert_{{p^{2}}={\mathcal{M}^{2}}_a} 
 = \left(\sqrt{{\hat{Z}}_a} {\hat{Z}}_{ai} \right)_{II} = {\hat{\textbf{Z}}}_{ai}\bigg\rvert_{II}\label{eq:TransformScheme}.
 \end{aligned}$$ While the values of ${\hat{Z}}_a$ and ${\hat{Z}}_{ai}$ depend on the choice of the index mapping, Eq.(\[eq:TransformScheme\]) ensures that the elements of the ${\hat{\textbf{Z}}}$-matrix, which appear in physical processes, are invariant under the choice of $a,b,c$ as a permutation of ${1,2,3}$. We also tested this relation numerically for various parameter points and always found agreement within the numerical precision.

Breit–Wigner approximation of the full propagators {#sect:fullBWZ}
==================================================

The propagator matrix depends on the momentum ${p^{2}}$ in a twofold way. On the one hand, the tree-level propagator factors $D_i{(p^{2})}={p^{2}}-m_i^{2}$ give rise to an explicit ${p^{2}}$-dependence. On the other hand, the self-energies ${\hat{\Sigma}_{ij}}{(p^{2})}$ depend on the momentum as well, but away from thresholds their ${p^{2}}$ dependence is not particularly pronounced. In this chapter, we will develop an approximation of the full mixing propagators with the aim to maintain the leading momentum dependence, but to greatly simplify the mixing contributions by making use of the on-shell ${\hat{\textbf{Z}}}$-factors.

Unstable particles and the total decay width
--------------------------------------------

In the context of determining complex poles of propagators, we now briefly discuss resonances and unstable particles, see e.g.Refs. [@Wackeroth:1996hz; @Freitas:2002ja; @Grunewald:2000ju; @Argyres:1995ym; @Beenakker:1996kn]. Stable particles are associated with a real pole of the S-matrix, whereas the self-energies of unstable particles develop an imaginary part, so that the pole of the propagator is located within the complex momentum plane off the real axis. For a single pole, the scattering amplitude $\mathcal{A}$ can be schematically written near the complex pole ${\mathcal{M}^{2}}_a$ in a gauge-invariant way as $$\mathcal{A}(s) = \frac{R}{s-{\mathcal{M}^{2}}_a}+F(s) \label{eq:gaugeinv},$$ where $s$ is the squared centre-of-mass energy, $R$ denotes the residue, while $F$ represents non-resonant contributions. The mass $M_{h_a}$ of the unstable particle $h_a$ is obtained from the real part of the complex pole ${\mathcal{M}^{2}}_a=M_{h_a}^{2}-iM_{h_a}\Gamma_{h_a}$, while the imaginary part gives rise to the total width. Accordingly, the expansion around the complex pole ${\mathcal{M}^{2}}_a$ leads to a Breit–Wigner propagator with a constant width, $${\Delta^{\textrm{BW}}}_a{(p^{2})}:=\frac{i}{{p^{2}}- {\mathcal{M}^{2}}_{a}} = \frac{i}{{p^{2}}- M_{h_a}^2 + i M_{h_a}\Gamma_{h_a}} . \label{eq:BWdef}$$ In the following, we will use a Breit–Wigner propagator of this form to describe the contribution of the unstable scalar $h_a$ with mass $M_{h_a}$ and total width $\Gamma_{h_a}$ in the resonance region.

Expansion of the full propagators around the complex poles
----------------------------------------------------------

Eqs.(\[eq:InvPropDet\]) and (\[eq:InvPropijDet\]) imply for $3\times
3$ mixing that each propagator $\Delta_{ii}, \Delta_{ij}$ has a pole at ${\mathcal{M}^{2}}_1,{\mathcal{M}^{2}}_2$ and ${\mathcal{M}^{2}}_3$. Because of this structure, an expansion of the full propagators near one single pole is not expected to yield a sufficient approximation. Instead, we will perform an expansion of the full propagators around all of their complex poles. The final expression obtained from combining the contributions from the different poles will constitute a main result of the present paper.

### Expansion of the diagonal propagators

We begin with an expansion of $\Delta_{ii}{(p^{2})}$ in the vicinity of ${\mathcal{M}^{2}}_a$ as in Eq.(\[eq:DiiApprdSeff\]), where the first factor equals the definition of the Breit–Wigner propagator of the state $h_a$, and the second factor corresponds to ${\hat{Z}}_a$ in scheme $I$ where $i$ and $a$ are associated indices. On top of that, ${\hat{Z}}_a\big\rvert_I={\hat{\textbf{Z}}}_{ai}^{2}$ as defined in Eq.(\[eq:ZiZij\]), and the elements of the ${\hat{\textbf{Z}}}$-matrix are independent of the index scheme (see Eq.(\[eq:TransformScheme\])). Thus, the following scheme-independent approximation holds for ${p^{2}}\simeq {\mathcal{M}^{2}}_a$: $$\begin{aligned}
 \Delta_{ii}{(p^{2})}&\simeq {\Delta^{\textrm{BW}}}_a{(p^{2})}\,{\hat{Z}}_a\big\rvert_I = {\Delta^{\textrm{BW}}}_a{(p^{2})}\,{\hat{\textbf{Z}}}_{ai}^{2}
\label{eq:DeltaiiBWa}.\end{aligned}$$ In this approach, the mixing contributions are summarised in the on-shell $Z$-factor evaluated at ${\mathcal{M}^{2}}_a$. In contrast, the leading momentum dependence is contained in the Breit–Wigner propagator parametrised by the loop-corrected mass $M_{h_a}$ and the total width $\Gamma_{h_a}$ from the complex pole. In addition, $\Delta_{ii}({p^{2}})$ has a second pole at ${\mathcal{M}^{2}}_b$ because ${p^{2}}-m_i^{2}+{\hat{\Sigma}^{\rm{eff}}_{ii}}(p^{2})=0$ holds also at ${p^{2}}={\mathcal{M}^{2}}_b$. Analogously, we can expand ${\hat{\Sigma}^{\rm{eff}}_{ii}}$ around ${\mathcal{M}^{2}}_b$ and obtain for the diagonal propagator $$\begin{aligned}
  \Delta_{ii}(p^{2}) &=\frac{i}{p^{2}-m_i^{2}+{\hat{\Sigma}^{\rm{eff}}_{ii}}(p^{2})}\\
&\simeq \frac{i}{p^{2}-m_i^{2} + {\hat{\Sigma}^{\rm{eff}}_{ii}}({\mathcal{M}^{2}}_b) +(p^{2}-{\mathcal{M}^{2}}_b)\cdot {\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_b)}\\
&=\frac{i}{(p^{2}-{\mathcal{M}^{2}}_b)\cdot\left[1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_b)\right]}\label{eq:DeltaiiMb}.\end{aligned}$$ Formally, $\frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_b)}$ has the structure of the definition of a ${\hat{Z}}$-factor from Eq.(\[eq:Zha\]), but in the index scheme $II$ where $b$ is assigned to $i$, whereas Eq.(\[eq:DeltaiiBWa\]) has been obtained in scheme $I$ with the $(i,a)$ assignment. Using the relation (\[eq:TransformScheme\]), we can rewrite Eq.(\[eq:DeltaiiMb\]) as $$\begin{aligned}
 \Delta_{ii}({p^{2}}) &\simeq {\Delta^{\textrm{BW}}}_b({p^{2}}) \cdot \frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_b)}\\
 &= {\Delta^{\textrm{BW}}}_b({p^{2}})\cdot \frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}}_{jj}({\mathcal{M}^{2}}_b)} \left(\frac{\Delta_{ji}}{\Delta_{jj}}\right)^{2}_{{p^{2}}={\mathcal{M}^{2}}_b} \label{eq:BWbinDeltaii}\\
 &= {\Delta^{\textrm{BW}}}_b({p^{2}})\cdot \left({\hat{Z}}_b\, {\hat{Z}}_{bi}^{2}\right)\bigg\rvert_{I}\label{eq:BWbSchemeia}\\
 &= {\Delta^{\textrm{BW}}}_b({p^{2}})\cdot {\hat{\textbf{Z}}}_{bi}^{2} \label{eq:BWbZbi},\end{aligned}$$ where the ${\hat{Z}}$-factors in Eq.(\[eq:BWbSchemeia\]) are expressed in the same scheme as in Eq.(\[eq:DeltaiiBWa\]). Hence, in the vicinity of ${p^{2}}\simeq {\mathcal{M}^{2}}_b$, the diagonal propagator $\Delta_{ii}$ can be approximated by the Breit–Wigner propagator of $h_b$ weighted by the square of ${\hat{\textbf{Z}}}_{bi}$ that ensures the coupling to the incoming fields as Higgs boson $i$, propagation as the mass eigenstate $h_b$ and the coupling to the outgoing fields again as Higgs boson $i$. In the same manner, $\Delta_{ii}$ can be expanded around the third complex pole, ${\mathcal{M}^{2}}_c$, yielding $$\begin{aligned}
 \Delta_{ii}({p^{2}}) &\simeq \frac{i}{(p^{2}-{\mathcal{M}^{2}}_c)\cdot\left[1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_c)\right]}\label{eq:DeltaiiMc}\\
 &\simeq {\Delta^{\textrm{BW}}}_c({p^{2}})\cdot \frac{1}{1+{\hat{\Sigma}^{\rm{eff'}}}_{kk}({\mathcal{M}^{2}}_c)} \left(\frac{\Delta_{ki}}{\Delta_{kk}}\right)^{2}\bigg\rvert_{{p^{2}}={\mathcal{M}^{2}}_c} \label{eq:BWcinDeltaii}\\
 &= {\Delta^{\textrm{BW}}}_c({p^{2}})\cdot {\hat{\textbf{Z}}}_{ci}^{2}\label{eq:BWcZci}.\end{aligned}$$ Thus, close to one of the complex poles (e.g. ${\mathcal{M}^{2}}_a$), the dominant contribution to the full propagator $\Delta_{ii}$ can be approximated by the corresponding Breit–Wigner propagator (${\Delta^{\textrm{BW}}}_a$) multiplied by the square of the respective ${\hat{\textbf{Z}}}$-factor, (${\hat{\textbf{Z}}}_{ai}^{2}$). However, close-by poles may cause overlapping resonances. In order to include this possibility and to extend the range of validity of the Breit–Wigner approximation to a more general case, we take the sum of all three Breit–Wigner contributions into account: $$\begin{aligned}
 \Delta_{ii}({p^{2}})\simeq {\Delta^{\textrm{BW}}}_a({p^{2}}) \,{\hat{\textbf{Z}}}_{ai}^{2} + {\Delta^{\textrm{BW}}}_b({p^{2}}) \,{\hat{\textbf{Z}}}_{bi}^{2} + {\Delta^{\textrm{BW}}}_c({p^{2}}) \,{\hat{\textbf{Z}}}_{ci}^{2} 
 = \sum_{a=1}^{3}{\Delta^{\textrm{BW}}}_a({p^{2}})\,{\hat{\textbf{Z}}}_{ai}^{2}
 \label{eq:iiBWsum}.\end{aligned}$$

### Expansion of the off-diagonal propagators {#sec:ExpansionOffD}

We proceed similarly for the off-diagonal propagators, which also have three complex poles so that we can expand the propagators around them. Note that ${\hat{\textbf{Z}}}_{ai}=\sqrt{{\hat{Z}}_a}$ and ${\hat{\textbf{Z}}}_{aj}=\sqrt{{\hat{Z}}_a}{\hat{Z}}_{aj}$ as defined in Eq.(\[eq:ZiZij\]). Starting at ${p^{2}}\simeq{\mathcal{M}^{2}}_a$, we express the ${\hat{Z}}$-factors in scheme $I$, $$\begin{aligned}
 \Delta_{ij}{(p^{2})}&= \frac{\Delta_{ij}{(p^{2})}}{\Delta_{ii}{(p^{2})}} \Delta_{ii}{(p^{2})}\simeq {\hat{Z}}_{aj} {\hat{\textbf{Z}}}_{ai}^{2}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}= {\hat{\textbf{Z}}}_{aj}{\hat{\textbf{Z}}}_{ai}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}\label{eq:DeltaijBWa},\end{aligned}$$ Next, we approximate $\Delta_{ij}$ near ${p^{2}}={\mathcal{M}^{2}}_b$: $$\begin{aligned}
 \Delta_{ij}{(p^{2})}&= \frac{\Delta_{ji}{(p^{2})}}{\Delta_{jj}{(p^{2})}} \Delta_{jj}{(p^{2})}\simeq {\hat{Z}}_{bi} {\hat{\textbf{Z}}}_{bj}^{2}\,{\Delta^{\textrm{BW}}}_b{(p^{2})}= {\hat{\textbf{Z}}}_{bi}{\hat{\textbf{Z}}}_{bj}\,{\Delta^{\textrm{BW}}}_b{(p^{2})}\label{eq:DeltaijBWb}.\end{aligned}$$ For ${p^{2}}\simeq{\mathcal{M}^{2}}_c$, we switch to a scheme where the indices $i$ and $c$ belong together. Thereby we can write $$\begin{aligned}
 \Delta_{ij}{(p^{2})}&= \frac{\Delta_{ij}{(p^{2})}}{\Delta_{ii}{(p^{2})}} \Delta_{ii}{(p^{2})}\simeq {\hat{Z}}_{cj} {\hat{\textbf{Z}}}_{ci}^{2}\,{\Delta^{\textrm{BW}}}_c{(p^{2})}= {\hat{\textbf{Z}}}_{cj}{\hat{\textbf{Z}}}_{ci}\,{\Delta^{\textrm{BW}}}_c{(p^{2})}\label{eq:DeltaijBWc},\end{aligned}$$ which is expressed in terms of scheme-invariant ${\hat{\textbf{Z}}}$-factors. Finally, we take the sum of Eqs.(\[eq:DeltaijBWa\])-(\[eq:DeltaijBWc\]) to obtain $$\begin{aligned}
 \Delta_{ij}({p^{2}}) &\simeq \sum_{a=1}^{3}{\hat{\textbf{Z}}}_{ai}\,{\Delta^{\textrm{BW}}}_a({p^{2}})\,{\hat{\textbf{Z}}}_{aj} \label{eq:ijBWsum}.\end{aligned}$$ This sum is illustrated diagrammatically in Fig.\[fig:fullvsBWZdiagram\].

(0,0)–(2,0); at (1,0.6) [$\boldsymbol{i}$]{}; (2.5,0) circle (0.5); (3,0)–(5,0); at (4,0.6) [$\boldsymbol{j}$]{}; at (5.6,0) [$\simeq$]{};

(0,0)–(1,0); at (0.5,0.6) [$\boldsymbol{i}$]{}; (1.3,0) circle (0.3); (1.6,0)–(3,0); at (2.3,0.6) [$\boldsymbol{h_1}$]{}; (3.3,0) circle (0.3); (3.6,0)–(4.6,0); at (4.1,0.6) [$\boldsymbol{j}$]{}; at (5.1,0) [$+$]{}; at (1.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}1}{\color{OliveGreen}i } }$]{}; at (3.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}1}{\color{OliveGreen}j } }$]{};

(0,0)–(1,0); at (0.5,0.6) [$\boldsymbol{i}$]{}; (1.3,0) circle (0.3); (1.6,0)–(3,0); at (2.3,0.6) [$\boldsymbol{h_2}$]{}; (3.3,0) circle (0.3); (3.6,0)–(4.6,0); at (4.1,0.6) [$\boldsymbol{j}$]{}; at (5.1,0) [$+$]{}; at (1.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}2}{\color{OliveGreen}i } }$]{}; at (3.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}2}{\color{OliveGreen}j } }$]{};

(0,0)–(1,0); at (0.5,0.6) [$\boldsymbol{i}$]{}; (1.3,0) circle (0.3); (1.6,0)–(3,0); at (2.3,0.6) [$\boldsymbol{h_3}$]{}; (3.3,0) circle (0.3); (3.6,0)–(4.6,0); at (4.1,0.6) [$\boldsymbol{j}$]{}; at (1.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}3}{\color{OliveGreen}i } }$]{}; at (3.3,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}3}{\color{OliveGreen}j } }$]{};

Eq.(\[eq:ijBWsum\]) represents the central result of this section, covering also the diagonal propagators in the special case of $i=j$. It shows how the full propagator can be approximated by the contributions of the three resonance regions, expressed by the Breit–Wigner propagators $\Delta_a{(p^{2})},\,a=1,2,3$, reflecting the main momentum dependence. The mixing among the Higgs bosons is comprised in the ${\hat{\textbf{Z}}}$-factors which are evaluated on-shell. Nonetheless, even a part of the momentum dependence of the self-energies is accounted for because the derivation of Eq.(\[eq:iiBWsum\]) is based on a first-order expansion of the momentum-dependent effective self-energies. Furthermore, the ${\hat{\textbf{Z}}}$-factors serve as transition factors between the loop-corrected mass eigenstates $h_a$ and the lowest-order states $i$ (although ${\hat{\textbf{Z}}}$ is not a unitary matrix transforming the states into each other). Pictorially, $\Delta_{ij}$ is a propagator that begins on the state $i$ and ends on $j$ while mixing occurs in between. Also on the RHS of Fig.\[fig:fullvsBWZdiagram\] and Eq.(\[eq:ijBWsum\]) the propagator in the $h_a$-basis begins with $i$ and ends on $j$. Thus, the coupling to the rest of any diagram connected to the propagator is well-defined. In between, each of the $h_a$ can propagate, and the correct transition is ensured by ${\hat{\textbf{Z}}}_{ai}$ and ${\hat{\textbf{Z}}}_{aj}$. All three combinations are visualised in Fig.\[fig:fullvsBWZdiagram\].

If ${\mathcal{CP}}$ is conserved and only $h$ and $H$ mix, or if two states are nearly degenerate and their resonances widely separated from the remaining complex pole, the full $3\times 3$ mixing is (exactly or approximately) reduced to the $2\times 2$ mixing case. Then the off-diagonal ${\hat{\textbf{Z}}}$-factors involving the unmixed state vanish or become negligible so that some terms in Eq.(\[eq:ijBWsum\]) approach zero.

Beyond that, if no mixing occurs among the neutral Higgs bosons, all off-diagonal full propagators as well as the off-diagonal ${\hat{\textbf{Z}}}$-factors vanish and each diagonal full propagator consists of only a single Breit–Wigner term where the ${\hat{\textbf{Z}}}$-factor is based on the diagonal self-energy instead of the effective self-energy. Thus, Eq.(\[eq:ijBWsum\]) covers all special cases of the a priori $3\times 3$ mixing among the neutral Higgs bosons.

### Uncertainty estimate {#sec:uncertainty}

As for the usual Breit–Wigner approximation in the case of a single pole, the approximation of Eq.(\[eq:ijBWsum\]) contains all resonant contributions for the case with mixing, while the non-resonant contributions that are neglected in this approximation are suppressed by powers of $\Gamma_{h_a}/M_{h_a}$. While a narrow width enhances the resonance contribution near its pole, it also causes it to fall rapidly for momenta away from the pole and thereby suppresses the non-resonant, far off-shell propagators, as well as their possible interference with a resonant term. The approximation in the off-shell regime is therefore of a comparable quality, determined by the ratios $\Gamma_{h_a}/M_{h_a}$ for the different resonances, as the standard narrow-width approximation without mixing. Interference effects between different resonances are accounted for by employing the approximation of Eq.(\[eq:ijBWsum\]). As we will discuss in more detail below, interference effects can be large for resonances that are close in mass, while on the other hand interferences between widely separated resonances are suppressed. The effects of the interference of resonances with a non-resonant background continuum can be inferred from the case without mixing.

While the on-shell approximation of Eq.(\[eq:ijBWsum\]) is based on the expansion around the complex pole, the description of the momentum dependence can be improved by incorporating higher powers of $({p^{2}}- {\mathcal{M}^{2}}_a)$. It should be noted in this context that apart from thresholds the self-energies generally depend only rather mildly on the incoming momentum.

In the following, we will assess the main sources of uncertainties that apply in the vicinity of the complex poles. The approximation given in Eq.(\[eq:ijBWsum\]) of the fully momentum-dependent propagators in terms of Breit–Wigner propagators and on-shell ${\hat{\textbf{Z}}}$-factors introduces an uncertainty by neglecting higher orders in the expansion of the effective self-energies and in the evaluation of the self-energies depending on complex momenta in Eq.(\[eq:SigijReIm\]). Besides, the pole condition is numerically not exactly fulfilled, which may represent an additional source of uncertainty.

#### Expansion of self-energies around the real part of a complex momentum

As defined in Eq.(\[eq:SigijReIm\]), we expand all self-energies as functions of a complex squared momentum ${p^{2}}$ around the real part of ${p^{2}}$ up to the first order in the imaginary part $p_i^2\equiv{\rm Im}\,{p^{2}}$. Furthermore, higher powers of $p_i^2$ can enter, $$\begin{aligned}
 {\hat{\Sigma}_{ij}}({p^{2}}) &\simeq {\hat{\Sigma}_{ij}}({p^{2}}_r) + i {p^{2}}_i\, {\hat{\Sigma}^{'}}_{ij}({p^{2}}_r)
 + \frac{1}{2}\left(i {p^{2}}_i\right)^2\, {\hat{\Sigma}^{''}}_{ij}({p^{2}}_r)\,.\end{aligned}$$ In order to expand the propagators in terms of the second-order contribution we define $$\begin{aligned}
 \epsilon^{{\rm Im}}_{ij} := \frac{\frac{1}{2}\left(i {p^{2}}_i\right)^2\, {\hat{\Sigma}^{''}}_{ij}({p^{2}}_r)}{{\hat{\Sigma}_{ij}}({p^{2}}_r) + i {p^{2}}_i\, {\hat{\Sigma}^{'}}_{ij}({p^{2}}_r)}\,,\end{aligned}$$ and calculate $\Delta_{ii}$ to first order in $ \epsilon^{{\rm Im}}_{ij}$. One should keep in mind, though, that we also evaluate the full propagators using the above-mentioned expansion of self-energies around the real part of the momentum. Thus, the comparison between the full and the approximated propagators is not directly affected by omitting terms of $\mathcal{O}(\epsilon^{{\rm Im}}_{ij})$, but it represents a source of uncertainty in both methods.

#### Expansion of effective self-energies around a complex pole

Regarding the effective self-energies, the expansion of ${\hat{\Sigma}^{\rm{eff}}_{ii}}{(p^{2})}$ around ${\mathcal{M}^{2}}_a$ in Eq.(\[eq:DiiApprdSeff\]) includes the term of $\mathcal{O}\left(p^2-{\mathcal{M}^{2}}_a\right)$. Additional terms in the expansion up to $\mathcal{O}\left((p^2-{\mathcal{M}^{2}}_a)^2\right)$, $$\begin{aligned}
 p^{2}-m_i^{2}+{\hat{\Sigma}^{\rm{eff}}_{ii}}(p^{2})
\simeq\left(p^{2}-{\mathcal{M}^{2}}_a\right)\cdot
\left(1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a)+\frac{1}{2}({p^{2}}- {\mathcal{M}^{2}}_a){\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a)\right)
\label{eq:d2Seff}\,,\end{aligned}$$ contribute if the term involving the second derivative of the effective self-energy at the complex pole, ${\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a)$, is non-negligible. We define $$\begin{aligned}
 \epsilon^p_{a,i} &= \frac{\frac{1}{2}({p^{2}}- {\mathcal{M}^{2}}_a){\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a)}{1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a)}
 \simeq \frac{1}{2} (x-1) {\mathcal{M}^{2}}_a\, {\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a)\,,
 \label{eq:epsp}\end{aligned}$$ where the approximation holds for ${\rm Re}\,{\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a),\,{\rm Im}\,{\hat{\Sigma}^{\rm{eff''}}_{ii}}({\mathcal{M}^{2}}_a) \ll 1$ and the variation of ${p^{2}}= x{\mathcal{M}^{2}}_a$ around ${\mathcal{M}^{2}}_a$. The impact of the neglected terms on the propagator $\Delta_{ii}$ reads at first order in $\epsilon^p_{a,i}$: $$\begin{aligned}
 \Delta_{ii}{(p^{2})}\simeq \frac{i}{({p^{2}}-{\mathcal{M}^{2}}_a)\, (1+{\hat{\Sigma}^{\rm{eff'}}_{ii}}({\mathcal{M}^{2}}_a))}\, \cdot (1-\epsilon^p_{a,i})\,.\end{aligned}$$

The higher powers of $({p^{2}}- {\mathcal{M}^{2}}_a)$ directly translate to an expansion in terms of $\Gamma_{h_a}/M_{h_a}$ via ${\rm Im}\,{\mathcal{M}^{2}}_a = -iM_{h_a} \Gamma_{h_a}$, analogously to the unmixed case. Expanding the diagonal propagator $\Delta_{ii}{(p^{2})}$ around the complex pole ${p^{2}}\sim {\mathcal{M}^{2}}_a$ yields $$\begin{aligned}
 \Delta_{ii}{(p^{2})}\simeq &\frac{i}{{p^{2}}- M_{h_a}^2 + iM_{h_a}\Gamma_{h_a}  }
  \cdot \left[ 1-\sum_{n=1}^{\infty} \frac{1}{n!} \left(\frac{{p^{2}}}{M_{h_a}^2 }-1 +
i\frac{\Gamma_{h_a} }{M_{h_a} }\right)^{n-1}\,M_{h_a}^{2(n-1)}\,
\hat\Sigma_{ii}^{\textrm{eff},(n)}\left({\mathcal{M}^{2}}_a\right)
 \right]\,, \label{eq:expandGammaOverM}\end{aligned}$$ where $\hat\Sigma_{ii}^{\textrm{eff},(n)}\left({\mathcal{M}^{2}}_a\right) \equiv \left.\frac{\partial^n {\hat{\Sigma}^{\rm{eff}}_{ii}}{(p^{2})}}{\partial ({p^{2}})^n}\right|_{{\mathcal{M}^{2}}_a}$ denotes the $n$th derivative of the effective self-energy evaluated at the complex pole ${\mathcal{M}^{2}}_a$. The expression $\frac{{p^{2}}}{M_{h_a}^2 }-1$ approaches zero at the resonance for ${p^{2}}\rightarrow M_{h_a}$, both in the Breit–Wigner propagator in front of the square bracket and in the correction term. Beyond that, the correction scales with powers of $\Gamma_{h_a}/M_{h_a}$ multiplied by derivatives of ${\hat{\Sigma}^{\rm{eff}}_{ii}}$. While the term of the first derivative giving rise to the diagonal ${\hat{\textbf{Z}}}$-factor $Z_{ai}^2$ is included in the approximation in Eq.(\[eq:iiBWsum\]), the term linear in $\Gamma_{h_a}/M_{h_a}$ times the second derivative, and higher-order terms, are omitted and therefore represent a source of uncertainty. These higher terms exist at the real part of the resonance, ${p^{2}}=M_{h_a}^2$, but vanish at the complex pole, ${p^{2}}= {\mathcal{M}^{2}}_a$. In any case, they are not only suppressed by a narrow width compared to the mass, but also by the smallness of the higher derivatives at the pole. We confirmed numerically that already the second derivatives are suppressed by several orders compared to the first derivatives.

The expansion in Eq.(\[eq:expandGammaOverM\]) indeed corresponds exactly to the well-known unmixed case for the simplification of $M_{h_a}\rightarrow
M, \Gamma_{h_a} \rightarrow \Gamma$ and ${\hat{\Sigma}^{\rm{eff}}_{ii}}\rightarrow \Sigma$. In contrast, the presence of mixing necessitates the expansion around all complex poles. The mixing is accounted for by the effective self-energies, but the structure of the dependence on $\Gamma_{h_a}/M_{h_a}$ remains in analogy to the unmixed case.

Even in the case of no mixing, the inclusion of the diagonal ${\hat{\textbf{Z}}}$-factor does not only represent a conceptual, but also a numerical improvement of the agreement between the full and the approximate propagator compared to the pure Breit–Wigner propagator.

#### Pole condition

As an additional source of uncertainties, Eq.(\[eq:0mixed\]) is numerically not exactly fulfilled. Instead, $${\mathcal{M}^{2}}_a - m_i^{2} + \hat{\Sigma}^{\textrm{eff}}_{ii}({\mathcal{M}^{2}}_a) = o_{a,i} \label{eq:PoleNumericalZero}$$ sums up to a small number $o_{a,i}$. However, evaluating the dimensionless ratio $\epsilon^o_{a,i}:= \frac{o_{a,i}}{m_i^2}$ for all combinations of $a$ and $i$ and its impact on the propagators, this source of uncertainty is found to be negligible for the scenarios analysed in Sections \[sect:prop3light\] and \[sect:propMixextreme\].

### Amplitude with mixing based on full or Breit–Wigner propagators

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{h}$]{}; at (4,0.6) [$\boldsymbol{h}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{hh}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{h}$]{}; at (4,0.6) [$\boldsymbol{H}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{hH}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{h}$]{}; at (4,0.6) [$\boldsymbol{A}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{hA}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{H}$]{}; at (4,0.6) [$\boldsymbol{h}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{Hh}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{H}$]{}; at (4,0.6) [$\boldsymbol{H}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{HH}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{H}$]{}; at (4,0.6) [$\boldsymbol{A}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{HA}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{A}$]{}; at (4,0.6) [$\boldsymbol{h}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{Ah}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{A}$]{}; at (4,0.6) [$\boldsymbol{H}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{AH}}$]{};

(-1,1)–(0,0)–(-1,-1); (0,0)–(2,0); (2.5,0) circle (0.5); (3,0)–(5,0); (6,1)–(5,0)–(6,-1); at (1,0.6) [$\boldsymbol{A}$]{}; at (4,0.6) [$\boldsymbol{A}$]{}; at (2.5,-1) [$\boldsymbol{\Delta_{AA}}$]{};

In a physical process where neutral Higgs bosons can appear as intermediate particles, all of them need to be included in the prediction, see Fig.\[fig:FullPropAmpl\] and Ref.[@Fowler:2010eba]. The Higgs part of the amplitude then contains a sum over the irreducible vertex functions ${\hat{\Gamma}_}i^{X}$ (for a coupling of Higgs $i$ at the first vertex $X$) and ${\hat{\Gamma}_}j^{Y}$ (for a coupling of Higgs $j$ at the second vertex $Y$) times the fully momentum-dependent mixing propagators, $$\begin{aligned}
 \mathcal{A} = \sum_{i,j=h,H,A} {\hat{\Gamma}_}i^{X} \,\Delta_{ij}{(p^{2})}\, {\hat{\Gamma}_}j^{Y} \label{eq:FullSum}.\end{aligned}$$ Applying Eq.(\[eq:ijBWsum\]), the amplitude in Eq.(\[eq:FullSum\]) can be approximated by the sum over Breit–Wigner propagators multiplied by on-shell $Z$-factors, in agreement with Ref.[@Fowler:2010eba], $$\begin{aligned}
 \mathcal{A} &\simeq \sum_{i,j=h,H,A} {\hat{\Gamma}_}i^{X} \,\left[\sum_{a=1}^{3}{\hat{\textbf{Z}}}_{ai}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}\,{\hat{\textbf{Z}}}_{aj}\right]\, {\hat{\Gamma}_}j^{Y} \label{eq:BWsum}\\
 &= \sum_{a=1}^{3} \left(\sum_{i=h,H,A}{\hat{\textbf{Z}}}_{ai}{\hat{\Gamma}_}i^{X}\right)\,{\Delta^{\textrm{BW}}}_a{(p^{2})}\,\left(\sum_{j=h,H,A}{\hat{\textbf{Z}}}_{aj}{\hat{\Gamma}_}j^{Y}\right)\label{eq:BWsums}\\
 &=\sum_{a=1}^{3} {\hat{\Gamma}_}{h_a}^{X}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}\,{\hat{\Gamma}_}{h_a}^{Y}\label{eq:VerthaBW}.\end{aligned}$$ The first bracket in Eq.(\[eq:BWsums\]) represents ${\hat{\Gamma}_}{h_a}^{X}$, i.e., the vertex $X$ connected to the mass eigenstate $h_a$ as for an external Higgs boson in Eq.(\[eq:Vertexha\]). Subsequently, the second bracket is equal to the coupling of $h_a$ at vertex $Y$, ${\hat{\Gamma}_}{h_a}^{Y}$. As opposed to Sect.\[sect:ExtH\], the $h_a$ is not on-shell here, but a propagator with momentum ${p^{2}}$ between the vertices $X$ and $Y$, represented by the Breit–Wigner propagator ${\Delta^{\textrm{BW}}}_a{(p^{2})}$. So the ${\hat{\textbf{Z}}}$-factors are not only useful for the on-shell properties of external Higgs bosons, but they can also be used as an on-shell approximation of the mixing between Higgs propagators. This will be investigated numerically in Sect.\[sect:NumBWfull\].

Concerning the propagators of unstable particles, the long-standing problems related to the treatment of unstable particles in quantum field theory should be kept in mind. While the introduction of a non-zero width into the propagator of an unstable particle is indispensable for treating the resonance region where the unstable particle is close to its mass shell, such a modification of the propagator in general mixes orders of perturbation theory. This can give rise to a violation of unitarity and / or to gauge-dependent results, see e.g. Ref.[@Denner:2014zga] for a recent discussion of this issue and references therein. The approximation of the full off-shell propagators, where the effective self-energy appears in the denominator of the diagonal propagator according to Eq.(\[eq:Diieff\]), in terms of on-shell ${\hat{\textbf{Z}}}$-factors and Breit–Wigner propagators is not meant to provide a solution to this long-standing problem. It should be noted, however, that the expression of the off-shell propagator in terms of on-shell ${\hat{\textbf{Z}}}$-factors and simple Breit–Wigner factors significantly remedies the gauge dependence of the off-shell propagator. This is due to the fact that the approximation of Eq. (\[eq:DeltaiiBWa\]) and Eq. (\[eq:ijBWsum\]) is on the one hand numerically close to the full propagator (in the considered kinematic region and applying the Feynman gauge) but on the other hand does not contain gauge-dependent contributions that depend on the squared momentum $p^2$ of the propagator. While those $p^2$-dependent contributions of the gauge part of a Dyson-resummed self-energy in the denominator of a propagator often lead to violations of unitarity for high values of $p^2$ (see e.g. Ref.[@Denner:1996gb]), the Breit–Wigner factors multiplied with ${\hat{\textbf{Z}}}$-factors that are evaluated on-shell are much better behaved in the limit where $p^2$ gets large. We defer a more detailed discussion of this issue to future work.

Calculation of interference terms in the Breit–Wigner formulation {#sect:IntBW}
-----------------------------------------------------------------

In Eq.(\[eq:ijBWsum\]), the Breit–Wigner propagators are combined such that they approximate a given full propagator. Conversely, we will now separate the $h_a$ part from the contribution of the other mass eigenstates in the amplitude with Higgs exchange between the vertices $X$ and $Y$: $$\begin{aligned}
 \mathcal{A}_{h_a} \equiv {\hat{\Gamma}_}{h_a}^{X}\, {\Delta^{\textrm{BW}}}_a{(p^{2})}\, {\hat{\Gamma}_}{h_a}^{Y}
 = \sum_{i,j=h,H,A} {\hat{\Gamma}_}i^{X}\, {\hat{\textbf{Z}}}_{ai}\, {\Delta^{\textrm{BW}}}_a{(p^{2})}\, {\hat{\textbf{Z}}}_{aj}\, {\hat{\Gamma}_}j^{Y} \label{eq:amplha},\end{aligned}$$ i.e. the exchange of the state $h_a$ coupling with the mixed vertices ${\hat{\Gamma}_}{h_a}$ from Eq.(\[eq:Vertexha\]) as for an external Higgs.

at (0.3,0) [$\hat\Gamma_{{\color{blue} h_a}}^{X}$]{}; (0.4,1)–(1.4,0.1); (0.4,-1)–(1.4,-0.1); (1.7,0) circle (0.3); (2,0)–(4,0); at (3,0.6) [$\boldsymbol{h_a}$]{}; (4.3,0) circle (0.3); (5.6,1)–(4.6,0.1); (5.6,-1)–(4.6,-0.1); at (5.5,0) [$\hat\Gamma_{{\color{blue} h_a}}^{Y}$]{}; at (-0.5,-3) [=]{};

at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{X}$]{}; (0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; (3.2,0)–(4.7,0); (5,0) circle (0.3); at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; (5.3,0)–(6.5,0); at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2,0.6) [$\boldsymbol{h}$]{}; at (5.9,0.6) [$\boldsymbol{h}$]{}; (7.5,1)–(6.5,0)–(7.5,-1); at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{Y}$]{}; at (8.4,0) [+]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (8.4,0) [+]{}; at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (2,0.6) [$\boldsymbol{h}$]{}; at (5.9,0.6) [$\boldsymbol{H}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{Y}$]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (2,0.6) [$\boldsymbol{h}$]{}; at (5.9,0.6) [$\boldsymbol{A}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{Y}$]{};

at (-0.6,0) [+]{}; (0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (8.4,0) [+]{}; at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; at (2,0.6) [$\boldsymbol{H}$]{}; at (5.9,0.6) [$\boldsymbol{h}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{Y}$]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (8.4,0) [+]{}; at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (2,0.6) [$\boldsymbol{H}$]{}; at (5.9,0.6) [$\boldsymbol{H}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{Y}$]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (2,0.6) [$\boldsymbol{H}$]{}; at (5.9,0.6) [$\boldsymbol{A}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{Y}$]{};

at (-0.6,0) [+]{}; (0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (8.4,0) [+]{}; at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}h } }$]{}; at (2,0.6) [$\boldsymbol{A}$]{}; at (5.9,0.6) [$\boldsymbol{h}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} h}}^{Y}$]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (8.4,0) [+]{}; at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}H } }$]{}; at (2,0.6) [$\boldsymbol{A}$]{}; at (5.9,0.6) [$\boldsymbol{H}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} H}}^{Y}$]{};

(0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); (3.2,0)–(4.7,0); (5,0) circle (0.3); (5.3,0)–(6.5,0); (7.5,1)–(6.5,0)–(7.5,-1); at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}A } }$]{}; at (2,0.6) [$\boldsymbol{A}$]{}; at (5.9,0.6) [$\boldsymbol{A}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{X}$]{}; at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} A}}^{Y}$]{};

at (-2,0) [$=\displaystyle{\sum\limits_{{\color{OliveGreen}\boldsymbol{i,j=h,H,A}}}}$]{}; at (0.3,0) [$\hat\Gamma_{{\color{OliveGreen} \boldsymbol{i} }}^{X}$]{}; (0.4,1)–(1.4,0)–(0.4,-1); (1.4,0)–(2.6,0); (2.9,0) circle (0.3); at (2.9,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}i } }$]{}; (3.2,0)–(4.7,0); (5,0) circle (0.3); at (5,-0.8) [${\hat{\textbf{Z}}}_{{\color{blue}a}{\color{OliveGreen}j } }$]{}; (5.3,0)–(6.5,0); at (3.95,0.6) [$\boldsymbol{h_a}$]{}; at (2,0.6) [$\boldsymbol{i}$]{}; at (5.9,0.6) [$\boldsymbol{j}$]{}; (7.5,1)–(6.5,0)–(7.5,-1); at (7.4,0) [$\hat\Gamma_{{\color{OliveGreen} \boldsymbol{j}}}^{Y}$]{};

In order to calculate the squared amplitude as a *coherent* sum, all contributions of $h_1, h_2, h_3$ are summed up first before taking the absolute square, $$\begin{aligned}
 |\mathcal{A}|^{2}_{\textrm{coh}} =  \bigg\lvert \sum_{a=1}^{3} \mathcal{A}_{h_a} \bigg\rvert^{2} \label{eq:amplcoh}\,.\end{aligned}$$ On the contrary, the *incoherent* sum is the sum of the squared individual amplitudes, which misses the interference contribution, $$\begin{aligned}
 |\mathcal{A}|^{2}_{\textrm{incoh}} =   \sum_{a=1}^{3} \bigg\lvert\mathcal{A}_{h_a} \bigg\rvert^{2} \label{eq:amplincoh}\,.\end{aligned}$$ Thus, an advantage of the Breit–Wigner propagators is also the possibility to conveniently discern the interference of several resonances from their individual contributions in a squared amplitude $$\begin{aligned}
 |\mathcal{A}|^{2}_{\textrm{int}} = |\mathcal{A}|^{2}_{\textrm{coh}} -|\mathcal{A}|^{2}_{\textrm{inccoh}}
  = \sum_{a<b}2\,\textrm{Re}\left[\mathcal{A}_{h_a} \mathcal{A}_{h_b}^{*}\right].\end{aligned}$$ In contrast, the squared amplitude based on the full propagators $$\begin{aligned}
 \lvert \mathcal{A}_{\textrm{full}}\rvert^{2} &= \Big\lvert\sum_{i,j=h,H,A} {\hat{\Gamma}_}i^{X} \,\Delta_{ij}{(p^{2})}\, {\hat{\Gamma}_}j^{Y} \Big\rvert^{2} \label{eq:fullA}\end{aligned}$$ does not allow for a straightforward determination of the pure interference term.

Numerical comparison in the MSSM with complex parameters {#sect:NumBWfull}
========================================================

After the analytical considerations so far, we will now compare the full propagators with their approximation as a combination of Breit–Wigner propagators and ${\hat{\textbf{Z}}}$-factors numerically. For definiteness, we will evaluate the propagators in the MSSM with complex parameters. We use `FeynHiggs-2.10.2` [^4][@Heinemeyer:1998np; @Heinemeyer:1998yj; @Degrassi:2002fi; @Heinemeyer:2007aq] for the numerical evaluation of Higgs masses, widths and mixing properties. In order to investigate the applicability of the expansion of the full propagators in one or all three resonance regions, we will first use as input a complex squared momentum around the three complex poles. For the later application to physical processes where the squared momentum equals the centre-of-mass energy $s$, we will also evaluate the propagators at ${p^{2}}=s$ near the real parts of the complex poles. In Sect.\[sect:prop3light\] we choose a scenario where all three Higgs bosons are relatively light so that we can study their mutual overlap. As a test of the ${\hat{\textbf{Z}}}$-factor approximation, we work in a scenario with large mixing between $H$ and $A$ in Sect.\[sect:propMixextreme\].

The MSSM with complex parameters at tree level {#chap:MSSM}
----------------------------------------------

In this section we fix the notation for the different particle sectors of the MSSM with complex parameters, following Ref.[@Frank:2006yh].

#### Sfermion sector {#sect:sfermion}

The mixing of sfermions $\tilde{f}_L, \tilde{f}_R$ within one generation into mass eigenstates $\tilde{f}_1, \tilde{f}_2$ is parametrised by the matrix $$\begin{aligned}
 M_{\tilde{f}}^{2}&=\left(
 \begin{matrix}
      M_{\tilde{f}_L}^{2}+m_f^{2}+M_Z^{2}\cos 2 \beta (I_f^{3}-Q_f{s_{W}}^{2}) & m_f X_f^{*}\\
      m_f X_f & M_{\tilde{f}_R}^{2}+m_{f}^{2}+M_Z^{2}\cos2\beta Q_f{s_{W}}^{2}
                        \end{matrix} \right),\label{eq:Msf}\\
 X_f &:= A_f - \mu^{*}\cdot \left\{ \begin{matrix}
                               \cot \beta,~f = \text{up-type}~~~~\\\tan \beta,~f = \text{down-type}.
                               \end{matrix} \right.\end{aligned}$$ The trilinear couplings $A_f=|A_f|e^{i\phi_{A_f}}$, as well as $\mu=|\mu|e^{i\phi_{\mu}}$, can be complex. These phases enter the Higgs sector via sfermion loops starting at one-loop order. Diagonalising $M_{\tilde{f}}^{2}$ for all $\tilde{f}$ separately, one obtains the sfermion masses $m_{\tilde{f}_1}\leq m_{\tilde{f}_2}$.

#### Gluino sector

The gluino $\tilde{g}^{a},~a=1,2,3,$ has a mass of $ m_{\tilde{g}}=|M_3|$, where $M_3=|M_3|\,e^{i\phi_{M_3}}$ is the possibly complex gluino mass parameter. Since the gluino does not directly couple to Higgs bosons, the phase $\phi_{M_3}$ enters the Higgs sector only at the two-loop level, but has an impact for example on the bottom Yukawa coupling already at one-loop order.

#### Neutralino and chargino sector {#sect:neucha}

At tree-level, mixing in the chargino sector is governed by the higgsino and wino mass parameters $\mu$ and $M_2$, respectively. In the neutralino sector it additionally depends on the bino mass parameter $M_1$. The charginos $\widetilde{\chi}_i^{\pm},~i=1,2$, as mass eigenstates are superpositions of the charged winos $\widetilde{W}^{\pm}$ and higgsinos $\widetilde{H}^{\pm}$, with the mass matrix, $$X =  {\begin{pmatrix}}M_2 & \sqrt{2}M_W{s_{\beta}}\\ \sqrt{2}M_W {c_{\beta}}& \mu {\end{pmatrix}}.
  \label{eq:X}$$ Likewise in the neutralino sector, the neutral electroweak gauginos $\widetilde{B},\,\widetilde{W}^{3}$ and the neutral Higgsinos $\tilde{h}_d^0,\, \tilde{h}_u^0$ mix into the mass eigenstates ${\widetilde{\chi}_i^{0}},\,i=1,...,4$. The mixing is encoded in the gaugino mass matrix $Y$, $$Y = \left(\begin{matrix}
    M_1 & 0 & -M_Z {c_{\beta}}s_W & M_Z {s_{\beta}}s_W\\
    0 & M_2 & M_Z {c_{\beta}}c_W & -M_Z {s_{\beta}}c_W \\
    -M_Z {c_{\beta}}s_W & M_Z {c_{\beta}}c_W & 0 &-\mu\\
    M_Z {s_{\beta}}s_W & -M_Z {s_{\beta}}c_W &  -\mu &0\\
         \end{matrix} \right) .
\label{eq:neutmix}$$ The gaugino mass parameters $M_1$ and $M_2$ as well as the higgsino mass parameter can in principle be complex. However, only two of the phases are independent, and a frequently used convention is to set $\phi_{M_2}=0$.

#### Higgs sector {#sect:Higgstree}

The MSSM requires two complex scalar Higgs doublets with opposite hypercharge $Y_{\mathcal{H}_{1,2}} = \pm 1$, $$\begin{aligned}
 \mathcal{H}_1 &= \begin{pmatrix}h_d^{0}\\h_d^{-}\end{pmatrix}
                  ={\begin{pmatrix}}v_d+\frac{1}{\sqrt{2}}(\phi_1^{0}-i\chi_1^{0})\\ -\phi_1^{-}{\end{pmatrix}}\label{eq:H1doublet}\\
 \mathcal{H}_2 &= \begin{pmatrix}h_u^{+}\\h_u^{0}\end{pmatrix}
                  = {\begin{pmatrix}}\phi_2^{+}\\v_u+\frac{1}{\sqrt{2}}(\phi_2^{0}+i\chi_2^{0}){\end{pmatrix}}. \label{eq:H2doublet}\end{aligned}$$ A relative phase between both Higgs doublets vanishes at the minimum of the Higgs potential, and the possible phase of the coefficient of the bilinear term in the Higgs potential can be rotated away. Hence, the Higgs sector conserves ${\mathcal{CP}}$ at lowest order, and the $4\times 4$ mass matrix of the neutral states $\mathbf{M}_{\phi\phi \chi\chi}$ becomes block-diagonal. The neutral tree level mass eigenstates are obtained by a diagonalisation of $\mathbf{M}_{\phi\phi \chi\chi}$: $$\begin{aligned}
 {\begin{pmatrix}}h\\H\\A\\G {\end{pmatrix}}= {\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}}{\begin{pmatrix}}\phi_1^{0}\\ \phi_2^{0}\\ \chi_1^{0}\\ \chi_2^{0}{\end{pmatrix}},\hspace*{1cm}
 {\begin{pmatrix}}H^{\pm}\\G^{\pm}{\end{pmatrix}}= {\begin{pmatrix}}-{s_{\beta_c}}& {c_{\beta_c}}\\ {c_{\beta_c}}&{s_{\beta_c}}{\end{pmatrix}}{\begin{pmatrix}}\phi_1^{\pm}\\ \phi_2^{\pm}{\end{pmatrix}}\label{eq:ncHmix},\end{aligned}$$ where we introduced the short-hand notation $s_x\equiv\sin x,~c_x\equiv\cos x$. For later use we define ${t_{\beta}}\equiv \tan\beta$. The mixing angle $\alpha$ is applied for the $\mathcal{CP}$-even Higgs bosons $h,H$; $\beta_n$ for the neutral $\mathcal{CP}$-odd Higgs $A$ and Goldstone boson $G$, and $\beta_c$ for the charged Higgs $H^{\pm}$ and the charged Goldstone boson $G^{\pm}$. The minimum conditions for $V_H$ lead to $\beta = \beta_n = \beta_c$ at tree level. At higher orders, however, $\tan\beta$ is renormalised whereas the mixing angles $\alpha, \beta_n$ and $\beta_c$ are not renormalised. The masses of the $\mathcal{CP}$-odd and the charged Higgs bosons are at tree level related by $$\begin{aligned}
 m_{H^{\pm}}^{2} &=& m_A^{2} + M_W^{2}\label{eq:mH+}.\end{aligned}$$ At lowest order, the Higgs sector is fully determined by the two SUSY input parameters (in addition to SM masses and gauge couplings) $\tan\beta$ and $m_{H^{\pm}}$ (or, for conserved $\mathcal{CP}$, equivalently $m_A$).

#### Higher-order corrections

Higher-order corrections in the MSSM Higgs sector are very relevant. Particles from other sectors contribute via loop diagrams to Higgs observables, in particular the trilinear couplings $A_f$, the stop and sbottom masses, and in the sub-leading terms the higgsino mass parameter $\mu$. Thus, beyond the lowest order, the Higgs sector is influenced by more parameters than only $M_{H^{\pm}}$ (or $M_A$) and $\tan\beta$. We adopt the hybrid on-shell and $\overline{\rm{DR}}$-renormalisation scheme defined in Ref.[@Frank:2006yh].

Scenario with three light Higgs bosons {#sect:prop3light}
--------------------------------------

For the numerical evaluation of the propagators, we work in the ${M_h^{\rm{mod+}}}$ scenario[@Carena:2013qia]. In this example, we fix the variable parameters $$\begin{aligned}
\mu &= 200{\,\textrm{GeV}},\nonumber\\
 {M_{H^{\pm}}}&= 160{\,\textrm{GeV}},\nonumber\\
 \tan\beta &= 50,\label{eq:choiceTbMHp}\end{aligned}$$ and choose for the complex phase ${\phi_{A_t}}=\pi/4$ to allow for ${\mathcal{CP}}$-violating mixing. This parameter choice is not meant to be the experimentally viable; in fact it has been experimentally ruled out by searches for $H,A\rightarrow \tau\tau$. The purpose is rather to provide a setting for illustration with nearby, but resolvable resonances so that the test of the Breit–Wigner approximation is not limited to well separated poles. These parameter values result in the following complex poles: $$\begin{aligned}
 {\mathcal{M}^{2}}_1 &=(15791-70i){\,\textrm{GeV}}^{2},~~~{\mathcal{M}^{2}}_2=(16202-525i){\,\textrm{GeV}}^{2},~~~{\mathcal{M}^{2}}_3=(17388-385i){\,\textrm{GeV}}^{2}\label{eq:polevalues}.
 \end{aligned}$$ All of the loop-corrected masses obtained from the real parts of the complex poles listed above are relatively light: $$\begin{aligned}
 M_{h_1}&=125.7{\,\textrm{GeV}},~~~M_{h_2}=127.3{\,\textrm{GeV}},~~~M_{h_3}=131.9{\,\textrm{GeV}}\label{eq:massvalues}
 \end{aligned}$$ so that the mass differences are of the order of – but not smaller than – the total widths from the imaginary parts of the complex poles, $\Gamma_{h_1}=0.6{\,\textrm{GeV}},~\Gamma_{h_2}=4.1{\,\textrm{GeV}},~\Gamma_{h_3}=2.9{\,\textrm{GeV}}$. The phase of $A_t$ induces ${\mathcal{CP}}$-violating mixing, and the on-shell mixing properties are reflected by the ${\hat{\textbf{Z}}}$-matrix obtained with `FeynHiggs`, $$\begin{aligned}
 {\hat{\textbf{Z}}}&= {\begin{pmatrix}}0.95-0.04i     & 0.34+0.09i    &-0.05-0.05i\\
         0.05-0.05i & 0.02+0.03i    & 0.99-0.02i\\
        -0.35-0.09i & 0.94-0.05i    &-0.006-0.003i  
 {\end{pmatrix}}\label{eq:Zvalues},\end{aligned}$$ which indicates that $h_1$ couples mostly $h$-like, $h_2$ mostly $A$-like and $h_3$ mostly $H$-like. Thus, the contribution of $h_a$ to $\Delta_{ij}$, i.e. $$\Delta_{ij}\Big\rvert_{h_a}{(p^{2})}={\hat{\textbf{Z}}}_{ai}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}\,{\hat{\textbf{Z}}}_{aj}\label{eq:hacontrib},$$ is only significant if the product ${\hat{\textbf{Z}}}_{ai}{\hat{\textbf{Z}}}_{aj}$ is not suppressed. In this case, we can already estimate that, for example, $h_3$ hardly contributes to $\Delta_{AA}$.

We have analysed all propagators $\Delta_{ij}$ around ${\mathcal{M}^{2}}_1,\,{\mathcal{M}^{2}}_2$ and ${\mathcal{M}^{2}}_3$ as well as for real momenta. In the following, we show and discuss a selection of these cases.

### Propagators depending on complex momenta

The analytical derivation of Eq.(\[eq:ijBWsum\]) builds on the expansion of the full propagators around the complex poles, and the on-shell condition in Eq.(\[eq:SolvePole\]) holds exactly only at complex momentum. Therefore, we evaluate the self-energies and propagators around the complex poles. Fig.\[fig:diag1\] displays $\Delta_{hh}{(p^{2})}$ for ${p^{2}}= 
0.5\,{\mathcal{M}^{2}}_1~...~1.5\,{\mathcal{M}^{2}}_1$. In particular, Fig.\[fig:diag1Re\] shows $\textrm{Re}\left[\Delta_{hh}\right]$, and Fig.\[fig:diag1Im\] shows $\textrm{Im}\left[\Delta_{hh}\right]$ versus the ratio $x_1={p^{2}}/{\mathcal{M}^{2}}_1$ such that $x_1=1$ corresponds to the complex pole ${p^{2}}={\mathcal{M}^{2}}_1$. The black line (labelled by $\Delta$ full) represents the fully momentum-dependent mixing propagator from Eq.(\[eq:Diieff\]). Since the three poles do not have the same ratio between the real and imaginary parts, scaling $x_1$ does not run into ${\mathcal{M}^{2}}_2$ and ${\mathcal{M}^{2}}_3$. $\Delta_{hh}$ has a pole at $x=1$ and a second peak at $x\simeq1.1$ which is close to the real part of ${\mathcal{M}^{2}}_3$. This structure is very precisely reproduced (for a discussion of the uncertainties, see below) by the sum $\sum_{a=1}^{3}{\hat{\textbf{Z}}}_{ah}^{2}\,{\Delta^{\textrm{BW}}}_a{(p^{2})}$ according to Eq.(\[eq:iiBWsum\]) – as can be seen from the red dotted line (labelled as $\sum\textrm{BW}\cdot Z$), which lies directly on top of the black solid line.

In order to understand which of the Breit–Wigner propagators and ${\hat{\textbf{Z}}}$-factors dominate at which momentum, the individual contributions are shown by the dashed curves. The blue line (labelled by $h_1$) represents the contribution of $h_1$ to $\Delta_{hh}$, i.e., ${\hat{\textbf{Z}}}_{1h}^{2}\,{\Delta^{\textrm{BW}}}_1{(p^{2})}$. It clearly displays the pole at $x=1$, but strongly deviates from the full propagator at momenta away from ${\mathcal{M}^{2}}_1$. The orange line (labelled by $h_2$) represents ${\hat{\textbf{Z}}}_{2h}^{2}\,{\Delta^{\textrm{BW}}}_2{(p^{2})}$. Since ${\hat{\textbf{Z}}}_{2h}$ is small in this scenario, the contribution of $h_2$ to $\Delta_{hh}$ is numerically suppressed, but a tiny share is visible. The green line (labelled by $h_3$) stands for ${\hat{\textbf{Z}}}_{3h}^{2}\,{\Delta^{\textrm{BW}}}_3{(p^{2})}$ and it contributes significantly to $\Delta_{hh}$ near ${\mathcal{M}^{2}}_3$ because ${\hat{\textbf{Z}}}_{3h}=-0.35-0.09i$ is sizeable in this scenario. So we notice that none of the individual Breit–Wigner propagators multiplied by the appropriate ${\hat{\textbf{Z}}}$-factors suffices to approximate the full propagator, which has three complex poles. On the other hand, the sum of all three Breit–Wigner propagators times ${\hat{\textbf{Z}}}$-factors yields an accurate approach to the full mixing. This holds for the real and the imaginary part.

Having discussed the example of a diagonal propagator, we will now assess whether the ${\hat{\textbf{Z}}}$-factor approximation succeeds also for off-diagonal propagators. For instance, Fig.\[fig:offdiag2\] depicts $\Delta_{HA}$ versus $x_2={p^{2}}/{\mathcal{M}^{2}}_2$ such that $x=1$ matches ${p^{2}}={\mathcal{M}^{2}}_2$ where the propagator diverges. As above, owing to the different ratio between the real and imaginary part of each complex pole, scaling $x_2$ does not run into ${\mathcal{M}^{2}}_1$ and ${\mathcal{M}^{2}}_3$, but $\Delta_{HA}$ peaks close to their real parts. As in Fig.\[fig:diag1\], the black line representing the full propagator and the red, dotted line representing the sum of Breit–Wigner propagators according to Eq.(\[eq:ijBWsum\]) agree very well. Additionally, one can see the individual Breit–Wigner shapes. Because the products of the relevant ${\hat{\textbf{Z}}}$-factors, here ${\hat{\textbf{Z}}}_{aH}{\hat{\textbf{Z}}}_{aA}$, are non-negligible for all $a=1,2,3$, each Breit–Wigner propagator is important in the approximation of both the real part (Fig.\[fig:HARe\]) and the imaginary part (Fig.\[fig:HAIm\]) of $\Delta_{HA}$. The other diagonal and off-diagonal propagators in this scenario which are not displayed here have an equally good agreement between the full calculation and the approximation.

In order to quantify this agreement, we compute the relative deviation of the approximated propagators from the full ones, for each diagonal and off-diagonal element of the propagator matrix, depending on the momentum around each of the complex poles. For ${p^{2}}=x\cdot {\mathcal{M}^{2}}_1$, $0.5\leq x\leq 1.5$ as in Fig.\[fig:diag1\], the relative deviation of the approximated propagators from the full ones is below $4\cdot 10^{-3}$ for $\Delta_{hh}$, $\Delta_{HH}$ and below $9\cdot 10^{-3}$ for $\Delta_{AA}$ (only $\Delta_{hh}$ is displayed here), apart from the regions where the propagator itself goes through zero and changes sign. A similar agreement holds both for the diagonal and off-diagonal propagators such as in Fig.\[fig:offdiag2\], and for their real and imaginary parts. As expected, the relative deviation is minimal at the complex pole ($x=1$). The largest deviation between the full and the approximated propagators occurs (apart from the sign flips) at momenta far away from the pole where the full propagators approach zero and in some cases the individual propagators of the states $h_a$ contribute with opposite signs to sum.

The relative deviation for the displayed momenta should be compared to the expected uncertainty based on the estimates in Sect.\[sec:uncertainty\]. The higher powers of the imaginary part of the momentum that are neglected in the approximate evaluation of the self-energies at complex momenta give rise to only a small effect ($\epsilon^{\rm Im}_{ij}\sim \mathcal{O}(10^{-7})$. This induces an effect on the propagator of $\mathcal{O}(10^{-9}-10^{-6})$), which does not account for the actual difference (as explained above, both the full and the approximated propagators are subject to this kind of uncertainties arising from the evaluation at complex momenta). In contrast, the relative effect of the expansion of the momentum-dependent effective self-energy around the complex pole has a more significant effect of up to $\epsilon^p_{a,i}\lesssim 6\cdot 10^{-3}$ for the analysed momentum range, where the largest deviations occur furthest away from the poles. As a conclusion, the uncertainty arising from the approximate evaluation of the propagators at complex momenta is found to be negligible, while the expected uncertainty from the expansion around a complex pole is very close to the observed deviation.

### Propagators depending on the real momentum ${p^{2}}=s$

The calculation of the propagators at and around the complex poles together with the evaluation of the self-energies at complex momenta according to Eq.(\[eq:SigijReIm\]) was needed to fulfil the assumptions of the approximation. However, in collider processes, the Higgs propagator might appear for example in the s-channel of a $2\rightarrow 2$ scattering process where the squared momentum equals the square of the centre-of-mass energy $s$. Therefore here we will check the Breit–Wigner approximation around the real parts of the complex poles.

Fig.\[fig:hhsRe\] shows $\textrm{Re}\left[\Delta_{hh}\right]$ in the range $\sqrt{{p^{2}}}\simeq M_{h_1}, M_{h_2}, M_{h_3}$ of the three loop-corrected masses given in Eq.(\[eq:massvalues\]). The propagator has a pronounced peak around $M_{h_1}$ and a smaller and broader one at $M_{h_3}$. Again, the approximation (red, dotted) defined in Eq.(\[eq:iiBWsum\]) meets the full propagator (black) very precisely. The contribution of $h_1$ multiplied by ${\hat{\textbf{Z}}}_{1h}^{2}$ of $\mathcal{O}(1)$ dominates near $M_{h_1}$. At $M_{h_3}$ the Breit–Wigner shape of $h_3$ is dominant although multiplied only by ${\hat{\textbf{Z}}}_{3h}^{2}=(-0.35-0.09i)^{2}$, but also the tail of ${\Delta^{\textrm{BW}}}_1$ is relevant. The resonance of $h_2$ is strongly suppressed by the small ${\hat{\textbf{Z}}}_{2h}$. The described behaviour is analogous in $\textrm{Re}\left[\Delta_{HH}\right]$ (not shown here), where the strongest peak around $M_{h_3}$ is dominated by the $h_3$-contribution and the second peak near $M_{h_1}$ by the $h_1$-contribution, whereas the $h_2$ component is negligible. Also in this case, the tails of the propagators of the relevant mass eigenstates extend also to the other resonance regions. Fig.\[fig:AAsRe\] visualises $\textrm{Re}\left[\Delta_{AA}\right]$ with a broad peak at $M_{h_2}$. The black curve of the full propagator is again directly underneath the red, dotted curve of the Breit–Wigner approximation, which in this case stems nearly entirely from $h_2$ because of ${\hat{\textbf{Z}}}_{2A}\simeq 1$. Within $\Delta_{AA}$, the contribution of $h_1$ only has a minor impact, which can be seen as a small kink in Fig.\[fig:AAsRe\]. Although ${\Delta^{\textrm{BW}}}_1{(p^{2})}$ gets close to its pole, the resonance of $h_1$ is strongly suppressed by the small ${\hat{\textbf{Z}}}$-factor ${\hat{\textbf{Z}}}_{1A}=0.05 (1+i)$. As we anticipated above from the structure of ${\hat{\textbf{Z}}}$ in Eq.(\[eq:Zvalues\]), ${\Delta^{\textrm{BW}}}_{3}$ is a negligible component of $\Delta_{AA}$ for this parameter point.\

So far we have seen that the Breit–Wigner formulation combined with on-shell ${\hat{\textbf{Z}}}$-factors accurately reproduces (up to a relative deviation of ${\mathcal{O}}(10^{-3}- 10^{-2})$) the momentum dependence of the full diagonal and off-diagonal propagators by adding the contributions from all resonance regions. If all resonances are sufficiently separated (not shown in this example) or if all but one contributing products of ${\hat{\textbf{Z}}}$-factors are negligible, a single Breit–Wigner term is enough to approximate the full propagator in one of the resonance regions. In the general case, however, all Breit–Wigner terms need to be included. Even if the peaks are not located very close to each other compared to their widths, the tail of one resonance supported by a substantial product of ${\hat{\textbf{Z}}}$-factors can leak into another resonance region.

Scenario with large mixing {#sect:propMixextreme}
--------------------------

While the scenario in the previous section was characterised by three relatively similar masses, we now choose a setting with quasi degenerate heavy states $h_2$ and $h_3$. In Sect.\[sect:prop3light\] we considered the $M_h^{\rm{mod}+}$-scenario with the standard value of $\mu= 200{\,\textrm{GeV}}$ in combination with the phase $\phi_{A_t}=\pi/4$, leading to a moderate mixing predominantly between $h$ and $A$. In Ref.[@Carena:2013qia] it was suggested to choose also different values, $\mu=\pm200,\pm500, \pm1000{\,\textrm{GeV}}$. So in addition to the choice above, we now apply the following modification of the parameters in Eq.(\[eq:choiceTbMHp\]): $$\begin{aligned}
 \mu &= 1000{\,\textrm{GeV}},\nonumber\\
 {M_{H^{\pm}}}&= 650{\,\textrm{GeV}},\nonumber\\
 \tan\beta &= 20.\label{eq:choicemue1000}\end{aligned}$$ This results in the complex poles $${\mathcal{M}^{2}}_1 = (15797-0.2i){\,\textrm{GeV}}^{2},~~{\mathcal{M}^{2}}_2 = (415336-1673i){\,\textrm{GeV}}^{2},~~{\mathcal{M}^{2}}_3 = (415554-1857i){\,\textrm{GeV}}^{2},$$ and therefore in similar masses of the heavy Higgs bosons, $$\begin{aligned}
 M_{h_1} &= 125.33\,{\,\textrm{GeV}},~~~M_{h_2}= 644.47{\,\textrm{GeV}},~~~M_{h_3}= 644.63{\,\textrm{GeV}}\,.\end{aligned}$$ There is a large mixing between $H$ and $A$, visible in the ${\hat{\textbf{Z}}}$-factors evaluated with `FeynHiggs`, $$\begin{aligned}
 {\hat{\textbf{Z}}}&\simeq {\begin{pmatrix}}1.01   &0      &0\\
 0  &1.15-0.27i &-0.47-0.66i\\
 0  &0.49+0.65i &1.13-0.28i
 {\end{pmatrix}}\,.\end{aligned}$$ In the scenario in Sect.[\[sect:prop3light\]]{}, ${\hat{\textbf{Z}}}$ is approximately unitary, $$\begin{aligned}
 {\hat{\textbf{Z}}}\cdot{\hat{\textbf{Z}}}^{\dagger}\Big\rvert_{\textrm{Eq.}~(\ref{eq:choiceTbMHp})} &\simeq {\begin{pmatrix}}1   &-0.1i  &0.2i\\
    0.1i    &1  &0\\
    -0.2i   &0  &1
 {\end{pmatrix}}\simeq \mathbb{1}.\end{aligned}$$ On the contrary, in the present scenario with $\mu=1000{\,\textrm{GeV}}$, the product ${\hat{\textbf{Z}}}\cdot{\hat{\textbf{Z}}}^{\dagger}$ deviates strongly from $\mathbb{1}$: $$\begin{aligned}
 {\hat{\textbf{Z}}}\cdot{\hat{\textbf{Z}}}^{\dagger}\Big\rvert_{\textrm{Eq.}~(\ref{eq:choicemue1000})}&\simeq {\begin{pmatrix}}1  &0  &0\\
 0  &2  &-1.7i\\
 0  &1.7i   &2 \label{eq:Znonunitary}
 {\end{pmatrix}}.\end{aligned}$$ Hence it is useful to examine whether the Breit–Wigner propagators with ${\hat{\textbf{Z}}}$-factors still yield a viable approximation of the fully momentum-dependent propagators in this scenario with large mixing.

### Propagators depending on complex momenta

Due to the large difference between $M_{h_1}$ and $M_{h_2}\simeq M_{h_3}$, the propagators $\Delta_{hh},\,\Delta_{hH}$ and $\Delta_{hA}$ are strongly suppressed around ${\mathcal{M}^{2}}_2$. Fig.\[fig:mue1000ReHH\] shows $\textrm{Re}\left[\Delta_{HH}\right]$ for complex momentum around ${\mathcal{M}^{2}}_2$. The full calculation (black) and the Breit–Wigner approximation (red, dotted) are in relatively good agreement with each other although the curves do not lie directly on top of one another as in the scenario of Sect.\[sect:prop3light\]. The two heavy states $h_2$ and $h_3$ with a mass difference of less than $0.2{\,\textrm{GeV}}$ and total widths of $\Gamma_{h_2}=2.6{\,\textrm{GeV}}$ and $\Gamma_{h_3}=2.9{\,\textrm{GeV}}$ are too close to be resolved. ${\hat{\textbf{Z}}}_{2H}^{2}\Delta_{h_2}$ (orange) and ${\hat{\textbf{Z}}}_{3H}^{2}\Delta_{h_3}$ (green) contribute with similar magnitude, but opposite signs so that the result differs strongly from the single terms. Indeed, only their sum provides a reliable approximation. A comparable situation is shown in Fig.\[fig:mue1000ReHA\] for $\textrm{Re}\left[\Delta_{HA}\right]$.

In this scenario of highly admixed, almost mass-degenerate states $h_2,\,h_3$, the relative deviation of the approximated propagators from the full ones is larger, in particular in the region where the propagator itself is close to zero. The uncertainty of the evaluation of complex momenta is negligible here, as in the scenario of Sect.\[sect:prop3light\]. Instead, the reasons for the resulting deviation are the following. The neglected terms from the expansion around the complex pole (see Eq.(\[eq:epsp\])) amount to $\epsilon^p_{a,i}\sim{\mathcal{O}}(10^{-2}-10^{-1})$ for the states $i=H,\,A$ involved in the mixing around the nearby complex poles ${\mathcal{M}^{2}}_a,\,a=2,3$. The relative effect of the deviation is enhanced by the fact that the full propagator is close to zero. The approximation consists of the sum of two resonant contributions ($h_2,\,h_3$), which are in this case several orders of magnitude larger than the full propagator and occur with opposite signs. Hence, despite a cancellation of several orders of magnitude the approximation yields the correct sign of the propagator and a good description of its behaviour and size. In contrast, a single-resonance approach would result in a prediction that would be off by several orders of magnitude. As a result, although slight deviations between the the approximated propagator and the full one are visible in this scenario that is characterised by large mixing effects, the approximation based on Eq.(\[eq:ijBWsum\]) still represents a very significant improvement compared to a single-resonance treatment.

### Propagators depending on the real momentum ${p^{2}}=s$

Fig.\[fig:mue1000S\] shows the same selection of propagators as in Fig.\[fig:mue1000Mm2\], but in this case evaluated at real momentum. The approximation from Eq.(\[eq:ijBWsum\]) leads to a good agreement between the propagators in the full mixing calculation (black) and the Breit–Wigner formulation (red, dotted), displayed for $\textrm{Re}[\Delta_{HH}]$ in Fig.\[fig:mue1000sReHH\] and for $\textrm{Re}[\Delta_{HA}]$ in Fig.\[fig:mue1000sReHA\]. While the $h_1$-part (blue) is negligible due to the much lower mass $M_{h_1}$, $h_2$ (orange) and $h_3$ (green) both contribute substantially because their complex poles are very close to each other. For the evaluation at real momentum, the diagonal propagators $\Delta_{HH}, \Delta_{AA}$ are reproduced with an accuracy of ${\mathcal{O}}(10^{-3})$ by the approximation while the relative deviation is ${\mathcal{O}}(10^{-1})$ for $\Delta_{HA}$ where a cancellation of larger terms occurs. These comparisons show that the Breit–Wigner approximation is also applicable in scenarios of quasi-degenerate states and a strong resonance-enhanced mixing. We note, however, that the agreement between the full propagators and those with on-shell mixing factors is less accurate here than in the scenario with moderate, nearly unitary mixing.

### Comparison of ${\hat{\textbf{Z}}}$-factors with effective couplings {#sect:UZ}

The effective coupling approach mentioned in Sect.\[sect:Ueff\] makes use of the unitary, real ${\textbf{U}}$-matrix instead of the ${\hat{\textbf{Z}}}$-matrix. However, ${\textbf{U}}$ is not evaluated at the complex pole, but at ${p^{2}}=0$ (see Sect.\[sect:Ueff\]) and it does not comprise the imaginary parts of the self-energies. We compare the effective coupling approach based on the matrix ${\textbf{U}}$ with the ${\hat{\textbf{Z}}}$-factor approach which does take the imaginary parts into account, but cannot be directly interpreted as a unitary transformation between the states of the different bases.

Based on ${\hat{\textbf{Z}}}$- and ${\textbf{U}}$-factors from `FeynHiggs`, Fig.\[fig:mue1000UZ\] displays in the scenario of Eq.(\[eq:choicemue1000\]) real parts of $\Delta_{HH}$ and $\Delta_{HA}$ at real momentum around $M_{h_2}\simeq M_{h_3}$, calculated as a fully momentum-dependent mixing propagator (black), using the ${\hat{\textbf{Z}}}$-matrix approach (red, dotted) defined in Eq.(\[eq:ijBWsum\]), i.e.$\Delta_{ij}^{Z}({p^{2}})\simeq \sum_{a=1}^{3}{\hat{\textbf{Z}}}_{ai}\,{\Delta^{\textrm{BW}}}_a({p^{2}})\,{\hat{\textbf{Z}}}_{aj}$, and the ${\textbf{U}}$-matrix variant (grey, dashed) in the ${p^{2}}=0$ approximation as $$\Delta_{ij}^{U}({p^{2}})\simeq \sum_{a=1}^{3}{\textbf{U}}_{ai}({p^{2}}=0)\,{\Delta^{\textrm{BW}}}_a({p^{2}})\,{\textbf{U}}_{aj}({p^{2}}=0) \label{eq:UijBWsum}.$$

\
The black curves in Fig.\[fig:mue1000UZ\] representing the full propagators are identical to those shown in Fig.\[fig:mue1000S\]. While in Fig.\[fig:UZReHH\] the curve representing the approximation in terms of ${\hat{\textbf{Z}}}$-factors is almost identical to the curve of the full $\Delta_{HH}$ (with a relative deviation at the peak of $0.8\%$), the approximation in terms of ${\textbf{U}}$-factors differs from the full result by up to $14\%$. In Fig.\[fig:UZReHA\], not even the shape of $\Delta_{HA}$ is correctly approximated by the approach using the ${\textbf{U}}$-matrix whereas the approach using the ${\hat{\textbf{Z}}}$-matrix comes close to the full $\Delta_{HA}$ (up to the small deviations that can be seen in the plots).

This analysis indicates that the ${\hat{\textbf{Z}}}$-factors combined with Breit–Wigner propagators are well-suited to describe the Higgs propagators including their mixing also in scenarios with close-by resonances and strong mixing. This approach captures the leading momentum dependence and adequately accounts for the imaginary parts. In contrast, the combination of ${\textbf{U}}$-factors and Breit–Wigner propagators is – despite its unitary nature – incomplete with respect to the mixing effects in the resonance region and regarding the significance of imaginary parts.

Breit–Wigner and full propagators in cross sections {#sect:fullBWxsec}
---------------------------------------------------

As an application of the derivations above, we calculate a cross section with Higgs exchange. We study the example process $b\overline{b} \rightarrow h, H, A \rightarrow \tau^{+}\tau^{-}$, where the intermediate Higgs bosons are once represented by the full mixing propagators $\Delta_{ij}$ and once by Breit–Wigner propagators multiplied by ${\hat{\textbf{Z}}}$-factors. In order to disentangle this investigation from other higher-order effects, we restrict the Higgs-fermion-fermion vertices to the tree-level and do not include the emission of real particles in the initial or final state, but focus on the propagator corrections.

For the implementation of the full propagator method, we extended a `FeynArts` model file. New scalars $ij$ are introduced that correspond to the full propagator $\Delta_{ij}{(p^{2})}$ and couple to the first vertex as the interaction eigenstate $i$ and to the second vertex as $j$. Those propagators are used in the `FormCalc` calculation supplemented by self-energies from `FeynHiggs` incorporating corrections up to the two-loop level and the full momentum dependence at the one-loop level.

Considering only mixing between $h$ and $H$ at this point, we choose as a ${\mathcal{CP}}$-conserving scenario the $M_{h}^{\textrm{max}}$-scenario[@Carena:1999xa; @Carena:2002qg] with $\tan\beta = 50,\, {M_{H^{\pm}}}= 153\,$GeV, but we modify it by setting $A_{f_3}=2504\,$GeV. As before, this scenario has been selected for illustration purposes and it is not meant to be phenomenologically viable. An outcome of this parameter choice are large off-diagonal $Z$-factors $\hat{\textbf{Z}}_{12}\simeq 0.65+0.29i,\, \hat{\textbf{Z}}_{21}\simeq 
-0.64-0.29i$, and $\hat{\textbf{Z}}_{11}\simeq 0.85-0.22i,\, 
\hat{\textbf{Z}}_{22}\simeq 0.84-0.23i$. The masses of the $\mathcal{CP}$-even Higgs bosons are very close to each other, $M_{h_1}=126.20\,$GeV and $M_{h_2}=127.55$GeV, while the widths obtained from the imaginary part of the complex poles are $\Gamma_{h_1}=0.94\,$GeV and $\Gamma_{h_2}=1.21\,$GeV. Despite its large width of $\Gamma_A=3.58\,$GeV, the third neutral Higgs boson does not overlap significantly with the other two resonances due to the mass of $M_{h_3}=119.91\,$GeV, and no mixing with the other two states occurs because we are considering here the ${\mathcal{CP}}$-conserving case of real parameters.

at (0,0) [![The partonic cross section $\hat\sigma(b\overline{b}\rightarrow 
\tau^{+}\tau^{-})$ in a modified $M_{h}^{\textrm{max}}$-scenario with $\tan\beta 
= 50$ and ${M_{H^{\pm}}}= 153\,$GeV. The cross section is calculated with the full mixing propagators (blue, solid), approximated by the coherent sum of Breit–Wigner propagators times ${\hat{\textbf{Z}}}$-factors with the interference term (red, dashed) and the incoherent sum without the interference term (grey, dot-dashed). The individual contributions mediated by $h_1$ (light blue), $h_2$ (green) and $h_3$ (purple) are shown as dotted lines.[]{data-label="fig:DeltaBWMhmax"}](figures/bbtautau_Mhmax_DeltaBW_hHA_p2 "fig:"){width="80.00000%"}]{}; at (-67pt,112pt) [**\^**]{};

Fig.\[fig:DeltaBWMhmax\] shows the partonic cross section $\hat\sigma(b\overline{b} \rightarrow h, H, A \rightarrow \tau^{+}\tau^{-})$ as a function of the centre-of-mass energy $\sqrt{\hat s}$, where $\hat s =
(p_b+p_{\overline{b}})^{2}$ is the squared sum of the momenta of the $b$- and $\overline{b}$-quarks in the initial state. The calculation based on the full propagators (represented by the blue, solid line) is in very good agreement with the cross section based on the coherent sum of the $h_1, h_2,
h_3$ contributions of the Breit–Wigner propagators multiplied with ${\hat{\textbf{Z}}}$-factors (red, dashed) according to Eq.(\[eq:amplcoh\]). Both curves lie on top of each other and contain two peaks originating from $h_1$ (light blue, dotted) and $h_2$ (green, dotted). The resonances of $h_1$ and $h_2$ partly overlap as the mass difference is of the order of the total widths, but the two peaks can still be distinguished. The $h_3$ contribution peaks at a lower mass in this scenario, but for completeness it is also shown (purple, dotted). The incoherent sum $|h_1|^{2}+|h_2|^{2}+|h_3|^{2}$ (grey, dash-dotted) from Eq.(\[eq:amplincoh\]) clearly overestimates the full cross section as a consequence of the missing interference term that turns out to be destructive in this case. It is taken into account in the full calculation and in the coherent sum of Breit–Wigner propagators with ${\hat{\textbf{Z}}}$-factors. For the efficient calculation of interference terms of quasi degenerate resonances, see e.g. Refs.[@Fuchs:2014ola; @Fuchs:2014zra].

While the comparison in Fig.\[fig:DeltaBWMhmax\] is restricted to the case of $2\times 2$-mixing due to a scenario with real parameters, the agreement of the cross section $\hat\sigma(b\overline{b} \rightarrow h, H, A \rightarrow \tau^{+}\tau^{-})$ calculated with the full or Breit–Wigner propagators can be seen also in a scenario with complex parameters in Fig.\[fig:epsilon\]. For the numerical evaluation, we choose a scenario with large mixing defined in Eq.(\[eq:choicemue1000\]), i.e. the ${M_h^{\rm{mod+}}}$ scenario with $\mu=1000{\,\textrm{GeV}}$ and the phase ${\phi_{A_t}}=\pi/4$ where we vary ${M_{H^{\pm}}}$ and ${t_{\beta}}$.

We investigate the relative deviation $\epsilon$ between the cross section $\sigma_{\rm{full}}$ based on the full propagators and the cross section $\sigma_{\rm{coh}}^{\rm{BW}{\hat{\textbf{Z}}}}$ based on the coherent sum of Breit–Wigner propagators with ${\hat{\textbf{Z}}}$-factors, where the total widths are obtained from the imaginary parts of the complex poles, $\Gamma^{\textrm{Im}}$, defined in Eq.(\[eq:GammaIm\]), $$\begin{aligned}
 \epsilon &= \frac{\sigma_{\rm{coh}}^{\rm{BW}{\hat{\textbf{Z}}}}({\phi_{A_t}}) 
}{\sigma_{\rm{full}}({\phi_{A_t}}) }-1\label{eq:epsilon}.\end{aligned}$$ Fig.\[fig:epsilon\] reveals that both methods agree very well, with a maximum deviation of $\pm2\%$ around ${M_{H^{\pm}}}=500{\,\textrm{GeV}},~\tan\beta=28$ and of about $0.8\%$ along the green band from this parameter point to larger values of ${M_{H^{\pm}}}$ and lower values of ${t_{\beta}}$. Otherwise the two calculations lead to the same results within $0.1\%$. Hence the use of Breit–Wigner propagators is suitable also for phenomenological applications such as the calculation of cross sections even in a scenario with complex parameters and large, ${\mathcal{CP}}$-violating mixing. A phenomenological investigation of such a scenario will be addressed in a forthcoming publication[@Fuchs:2017poi; @Fuchs:2017wkq; @IntCalc:InProgress].

Impact of the total width {#sect:Gtot}
-------------------------

This section addresses the impact of the precise value of the total width. So far, we have obtained the Higgs widths from the imaginary part of the complex poles as in Eq.(\[eq:MGamma\]) in order to consistently compare with the full propagator mixing. If the self-energies in ${\hat{\Sigma}^{\rm{eff}}_{ii}}$ are calculated at the one-loop level, the total width extracted from a complex pole of $\Delta_{ii}$ is then a tree-level width. Correspondingly, partial two-loop contributions to the imaginary parts of the self-energies give rise to partial one-loop corrections of the decay width. However, two-loop self-energies evaluated at ${p^{2}}=0$, as they are approximated in `FeynHiggs`, do not contribute to the imaginary part of the pole so that the width determined from the imaginary part of the complex pole remains at its tree-level value. Corrections to Higgs boson decays in the MSSM at and beyond the one-loop level are known and have been found to be important, see e.g. Refs.[@Dabelstein:1995js; @Djouadi:1995gt; @Williams:2011bu; @Heinemeyer:2015pfa]. Thus, the sum of the partial decay widths into any final state $X$ of a Higgs boson $h_a$, $$\Gamma_{h_a}^{\textrm{tot}} = \sum\limits_X \Gamma(h_a \rightarrow X),\label{eq:GammaTot}$$ leads to a more accurate result for its total width than from the imaginary part of the corresponding complex pole, $$\Gamma_{h_a}^{\textrm{Im}}=-\textrm{Im}[{\mathcal{M}^{2}}_a]/M_{h_a}, \label{eq:GammaIm}$$ even if ${\mathcal{M}^{2}}_a$ is based on self-energies at the same order as used for the calculation of $\Gamma_{h_a}^{\textrm{tot}} $ as in Eq.(\[eq:GammaTot\]). `FeynHiggs` contains the partial Higgs decay widths and their sum at the leading two-loop order. Having checked the compelling agreement between the full propagators and the Breit–Wigner propagators with the width from the imaginary part of the complex pole in the previous sections, now we implement the total width from `FeynHiggs` into the Breit–Wigner propagators in order to obtain the most precise phenomenological prediction.

at (0,0) [![Effect of the total width as an input for Breit–Wigner propagators: The partonic cross section $\hat\sigma(b\overline{b}\rightarrow \tau^{+}\tau^{-})$ in the same modified $M_{h}^{\textrm{max}}$-scenario as in Fig.\[fig:DeltaBWMhmax\] with $\tan\beta = 50$ and $M_{H^{+}} = 153\,$GeV. The Breit–Wigner propagators with the total widths from the imaginary part of the complex pole including (red, dashed) and excluding (grey, dash-dotted) the interference term (as in Fig.\[fig:DeltaBWMhmax\]) are compared with the Breit–Wigner propagators where the total widths are obtained from `FeynHiggs`. The corresponding results are shown including (black, solid) and excluding (black, dotted) the interference term. []{data-label="fig:GammaImTotMhmax"}](figures/bbtautau_Mhmax_BW_Im_Gtot "fig:"){width="80.00000%"}]{}; at (-60.8pt,110.5pt) [**\^**]{}; at (-42pt,114.5pt) [-]{}; at (-176pt,-1pt) [**\^**]{}; at (13pt,-109pt) [**\^**]{};

In the modified ${M_h^{\rm{max}}}$ scenario, the higher-order corrections have a significant impact on the Higgs decay widths so that $\Gamma_{h_1}^{\textrm{tot}}=2.55\,$GeV and $\Gamma_{h_2}^{\textrm{tot}}=3.24\,$GeV are much larger than the widths obtained from the imaginary part of the complex pole. This affects the order of magnitude of the cross-section $\hat\sigma\left(b\overline{b}\rightarrow \tau^{+}\tau^{-}\right)$ and the structure of the resonances, as can be seen in Fig.\[fig:GammaImTotMhmax\]. The coherent sum of Breit–Wigner propagators including the interference term (red, dashed) and the incoherent sum without the interference term (grey, dash-dotted) using $\Gamma^{\textrm{Im}}$ from the imaginary parts of the complex poles are the same as in Fig.\[fig:DeltaBWMhmax\]. In contrast, the total widths $\Gamma^{\textrm{tot}}_{\textrm{FH}}$ obtained from `FeynHiggs` as the sum of higher-order partial widths are implemented into the Breit–Wigner propagators in the cross section based on the coherent sum of all $h_a$-contributions (black, solid) and the incoherent sum (black, dotted).

The large increase in the widths from $\Gamma_{h_a}^{\textrm{Im}}$ to $\Gamma_{h_a}^{\textrm{tot}}$ has a very significant effect in Fig.\[fig:GammaImTotMhmax\], since the resonant behaviour is smeared out by the larger widths. As a consequence, the separate resonances are less pronounced, and the cross section is suppressed. Here, the incoherent sum without the interference term again overestimates the cross-sections. In addition, it lacks the two-peak structure. This observation emphasizes the importance of including the total width at the highest available precision and to take the interference term into account. One can also see that two resonances might be too close to be resolved if they are smeared by large widths.

Conclusions {#sect:conclusions}
===========

We have shown that the momentum-dependent propagator matrices of systems of unstable particles that mix with each other can be accurately approximated by a combination of Breit–Wigner propagators and wave function normalisation factors, where the latter are evaluated at the complex poles of the propagator matrix. For illustration, we have applied this approach to the example of the neutral Higgs bosons of the MSSM with complex parameters, where the ${\mathcal{CP}}$-violating interactions give rise to a $3\times 3$ mixing between the lowest-order mass eigenstates $h,H,A$ and the loop-corrected mass eigenstates $h_1,h_2,h_3$ (further mixing contributions with unphysical Goldstone bosons and vector bosons can be incorporated separately up to the considered order in perturbation theory). In the special case of ${\mathcal{CP}}$-conservation the mixing between the physical Higgs bosons is reduced to a $2\times 2$ mixing between the states $h$ and $H$.

Analysing the pole structure of propagator matrices, we have shown that for the case of $3\times 3$ mixing each entry of the propagator matrix has three complex poles, while for the $2\times 2$ mixing case each entry of the propagator matrix has two complex poles. Consequently, a single-pole approximation is not sufficient to approximate the momentum dependence of the full propagators. In particular for close-by states with sizeable widths, the different resonance regions may overlap.

We have derived in a process-independent way how the full propagators can be expanded around all of their complex poles. This approximation results in the sum of Breit–Wigner propagators of the corresponding resonances, weighted by wave function normalisation factors which encompass the relation between the lowest-order mass eigenstates and the loop-corrected mass eigenstates. As a key result, we have demonstrated that wave function normalisation factors that have been derived to ensure the correct on-shell properties of external particles in physical processes at their (in general complex) poles are a useful tool also for describing off-shell propagators, i.e. including momentum dependence. It has been shown that the momentum dependence of the full propagators can be accurately approximated by simple Breit–Wigner propagators. The complex wave function normalisation factors properly incorporate the imaginary parts of the self-energies arising from absorptive parts of loop integrals, which in contrast are neglected in effective coupling approaches where the contributing self-energies are evaluated at vanishing external momentum.

In our numerical analysis we have found very good agreement between the approximation presented in this paper and the full propagators. Besides the numerical comparison of the full and approximated propagators around the complex poles or depending on real momenta, we have provided a detailed discussion of the uncertainties involved in the approximation. In particular, the omitted higher-order terms from the expansion of self-energies around the real part of complex momenta, the expansion of the momentum-dependent effective self-energies around the complex poles and numerical effects in the on-shell condition have been quantified. The uncertainty estimates of these sources are in general well reflected in our numerical results. We have also discussed the validity of the approximation far away from the poles.

The formalism of Breit–Wigner propagators and on-shell wave function normalisation factors has several appealing advantages in describing the propagators of unstable particles that mix with each other. It avoids the momentum-dependent evaluation of self-energies and thereby significantly simplifies and accelerates the calculation of higher-order contributions. Besides, it enables the separation of the individual resonant contributions of the mass eigenstates $h_a$ within a given process and the straightforward calculation of their interference term. Moreover, the Breit–Wigner propagator turns into a $\delta$-distribution in the limit of a vanishing width, thus facilitating the separate calculation of the production and decay of an intermediate unstable particle by means of the narrow-width approximation and its generalisation to the case of overlapping and interfering resonances. In addition, the Breit–Wigner formulation allows the implementation of a more precise total width by incorporating important higher-order effects from the partial widths that are not included in the imaginary part of the complex pole with self-energies evaluated at the same perturbative order as for the partial widths. This feature is very useful for phenomenological predictions of processes involving the exchange of unstable particles, benefitting from the use of quantities computed at the highest available order in perturbation theory.

The explicit calculations presented in this paper have been performed in the MSSM, but the introduced method itself can be easily extended to the cases of particles with non-zero spin and mixing among more than three particles, such as, for example, in different models with a non-minimal scalar sector or new vector resonances.

Acknowledgements {#acknowledgements .unnumbered}
================

We thank Alison Fowler for her contributions at an early stage of this work. E.F. thanks the DESY theory group where most of this work was carried out. The work of E.F. was partially supported by the German National Academic Foundation. The work of G.W. is supported in part by the Collaborative Research Centre SFB 676 of the DFG, “Particles, Strings and the Early Universe”, and by the European Commission through the “HiggsTools” Initial Training Network PITN-GA-2012-316704.

[^1]: In order to avoid sign ambiguities, taking the square root of a $Z$-factor, which in general has two solutions in the complex plane, refers here and in the following always to the principal square root, i.e.the solution with a non-negative real part.

[^2]: This index notation differs from the conventions in Refs.[@Frank:2006yh; @Hahn:2006np; @Williams:2007dc; @Fowler:2010eba; @Williams:2011bu].

[^3]: Explicitly, all roots of ${\boldsymbol{\Delta}}_{hHA}{(p^{2})}^{-1}$ are roots of $\det\left[{\boldsymbol{\Delta}}_{hHA}{(p^{2})}\right]^{-1} 
    \equiv \det\left[-{\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}^{-1}\right]^{-1} = 
    -\det\left[{\hat{\boldsymbol{\Gamma}}}_{hHA}{(p^{2})}\right]$.

[^4]: The additional contributions contained in more recent versions are not essential for the numerical comparison carried out here.
