Towards Higgs masses and decay widths satisfying the symmetries in the (N)MSSM

In models with an extended Higgs sector, such as the (N)MSSM, scalar states mix with one another. Yet, the concept of Higgs mixing is problematic at the radiative level, since it introduces both a scheme and a gauge dependence. In particular, the definition of Higgs masses and decay amplitudes can be impaired by the presence of gauge-violating pieces of higher order. We discuss in depth the origin and magnitude of such effects and suggest two strategies preserving or restoring gauge invariance. In addition, the intuitive concept of mixing and the simplicity of its definition in terms of two-point diagrams can make it tempting to include higher-order corrections on this side of the calculation, irrespectively of the order achieved in vertex diagrams. Using the global $SU(2)_{\mathrm{L}}$-symmetry in the decoupling limit, we show that no improvement can be expected from such an approach at the level of the Higgs decays, but that, on the contrary, the higher-order terms may lead to numerically large spurious effects.


Introduction
Many models of new physics suggest the existence of an extended Higgs sector. In such a context, the Higgs boson discovered at the LHC [1][2][3] and presenting characteristics that are approximately consistent with a Standard Model (SM) interpretation [4][5][6], would be regarded as only one of many scalar states. The absence of conclusive evidence for the additional Higgs bosons admittedly constrains the available parameter space but continues to spare multiple scenarios where, in general, the decoupling of light new-physics states from the SM particles requires a careful handling of the Higgs mixing (see e. g. Ref. [7]). On the other hand, the production of heavy states at colliders is kinematically suppressed, hence leaving little constraints on new Higgs bosons beyond the TeVrange [8]. The prototype of such extensions of the SM is the Two-Higgs-Doublet Model (THDM) [9] but further enlargement through supersymmetric (SUSY) sectors [10,11] or singlet fields are easy to motivate.
The presence of extended Higgs sectors opens up the possibility for mixing between the Higgs fields. This concept appears as relatively intuitive at the tree level, but is in fact ill-defined when considering radiative corrections. Indeed, corresponding definitions depend on the renormalization scheme and on the chosen gauge. The effective-potential approximation [12][13][14][15][16] offers a popular definition of the loop-corrected mixing, preserving the unitarity of the mixing matrix, but the missing momentumdependent corrections make it inappropriate-or insufficient-for a consistent description of external legs in Feynman amplitudes, in particular when considering the decays of such mixed states. On the other hand, one may directly use the LSZ reduction formula in order to define the Higgs mixing in terms of the loop corrections applying on an external Higgs leg in a physical amplitude: this has been described in e. g. Refs. [17][18][19]. Among the advantages of this approach, loop diagrams applying on the external leg should be automatically contained within the mixing matrix, making it formally suitable for a 'pseudo-on-shell' treatment of the external legs. On the other hand, unitarity of the mixing matrix is lost in the inclusion of momentum-dependent pieces.
In this paper, we outline several shortcomings in the use of a mixing matrix as a substitute to the loop expansion derived from the LSZ reduction formula. These problems rest less in the principle of the procedure than in its technical implementation. Indeed, the formalism described in e. g. Ref. [17] is designed in such a way that it should coincide with the LSZ expansion, at least at the order of the calculation. However, the separation of the contributions of one-loop order between mixing effects on one side, and vertex corrections (e. g. in the case of a two-body Higgs decay) on the other, lends itself to a misleading step, namely working with different orders (or numerical parameters) on each side. While the resulting mismatch is formally a higher-order effect, it can be numerically significant due to imperfect cancellations. Reasons for distrusting these partial higher-order contributions appear clearly when they violate a symmetry that is expected to hold, either exactly or approximately in a given regime. The electroweak gauge symmetry or the global SU (2) L -symmetry in the limit of Higgs masses far above the electroweak scale are examples of such handles on the validity of the calculation, providing a guideline for the resolution of the unphysical effects or at least an estimate of the associated uncertainties.
In practice, we work in the context of the Next-to-Minimal Supersymmetric Standard Model [20,21] and aim at improving our previous work on the Higgs decays at one-loop electroweak order in this model [18,22,23]-we refer the reader to e. g. [19,[24][25][26][27][28][29] for similar projects. However, our discussion in this paper should be valid for a large class of extensions of the SM based on a THDM framework (in particular the MSSM). In fact, as the symmetry arguments control the properties of doublet, but not of singlet states, we focus below on scenarios with doublet-dominated, MSSM-like Higgs bosons. In addition, the SUSY context induces some additional complications related to the connection between the gauge and the quartic Higgs couplings. Indeed there are then too few degrees of freedom available to simultaneously renormalize all the Higgs masses on-shell. One advantage is the gain in predictivity, since the mass of the SM-like Higgs boson cannot be set to an arbitrary value, while the mass-splitting between heavy doublet states is determined by electroweak effects. On the other hand, this causes a difficulty in evaluating Feynman amplitudes involving a loop-corrected mass as kinematical input and simultaneously preserving gauge invariance. A priori, such an issue also exists in a THDM, but there it can be easily evaded when working in an on-shell scheme. Consequently, we study in some detail the gauge dependence in the observables and suggest possible means of restoring a manifest gauge invariance. All our one-loop calculations are performed with the assistance of FeynArts [30], FormCalc [31,32] and LoopTools [32]. Results at the two-loop order are derived with the help of TwoCalc [33], TSIL [34] and TLDR [35].
In the following sections, we first analyze the electroweak gauge dependence in Higgs masses and decays. This provides us with a first formal argument to disfavor the naive inclusion and/or resummation of higher-order corrections in the Higgs-propagator matrix. Then, we consider the global SU (2) L -symmetry in the limit of massive doublet-like Higgs states and compare the consequences of this symmetry at the analytical and the numerical level for heavy-Higgs mass-splittings and several decay channels. Finally, we conclude as to a cautious and consistent use of mixing formalisms and the combination of mixing contributions with the vertex corrections in decay amplitudes.

Aspects of gauge invariance on the determination of Higgs masses and decays
The loop-induced mixing in the Higgs sector is defined at the level of the Higgs self-energies. However, the latter are not gauge-invariant objects, in general, highlighting the artificiality of the loop-corrected mixing matrix. Below, we detail how gauge invariance is ensured at the one-loop order in observable quantities such as the Higgs masses and decays, or how one could attempt to remedy its violation by terms of higher order.

Self-energies and gauge dependence
We first consider the gauge-dependent contributions to the self-energies of one-loop order between two external Higgs legs h i and h j (a priori mass-eigenstates at the tree level) that convey an external momentum p. The calculation applies to extensions of the SM of type THDM-details concerning the THDM Higgs sector are provided in appendix A-though additional fields (e. g. scalar singlets) can be mixed to this sector as long as they do not generate additional breaking of the electroweak symmetry. For simplicity, we focus on the terms associated with the electroweak charged current (the generalization to the neutral current is straightforward). The terms that depend on the gauge-fixing parameter ξ read All the g xyz (x, y, z ∈ {h i,j , G ± , H ± , W ± }) represent Higgs couplings, while M W , m H ± , and m h i symbolize the W , charged Higgs, and neutral Higgs masses, respectively. We checked that Eq. (1) agrees with expressions available in the literature, e. g. Eq. (3.1) of Ref. [36] (see also Refs. [37,38]). Here, we have not explicitly written the momentum-independent terms in A 0 (ξ M 2 W ) and just collected them within C A . In particular, we expect additional gauge-dependent terms of this form from the counterterms, so that the expression for C A before renormalization is of limited interest.
In order to exploit Eq. (1), it is useful to consider the various relevant couplings, which are fully determined by the gauge symmetry: Here, m h i or m H ± represent the tree-level Higgs masses (in contrast with loop-corrected masses M h i or M H ± that we consider later on). The symbol g 2 is the SU (2) L -gauge coupling and X R/I iu/d (we also define X ik ≡ X R ik + ı X I ik ) correspond to the components of the external Higgs h i projecting onto the real (R) or imaginary (I) part of the doublet fields H d or H u that take on a vacuum expectation value (v.e.v.) v c β or v s β , respectively, with v = M W /( √ 2 g 2 ). Then, the first three lines of Eq. (1) obviously vanish as long as h i (or h j ) is orthogonal to the neutral Goldstone boson. The other terms amount to The dependence of this expression on the tree-level Higgs masses originates in the Higgs-Goldstone couplings, except for the A 0 term. For this latter term, one systematically expects contributions of the same form from other origins, e. g. tadpoles. In order to clarify, we may consider the example of a SUSY framework with on-shell renormalization conditions for the W , Z and charged-Higgs masses, and vanishing tadpoles (see e. g. Ref. [39] for details in the MSSM): in such a 'physical' setup, all the mass parameters in the doublet sector are determined by observable quantities; then the tadpole and W -, Z-, H ± -mass counterterms contribute terms ∝ A 0 (ξ M 2 W ) that cancel outC A at the level of the renormalized diagonal self-energies. Alternatively, still in the SUSY case, we could employ DR conditions for the W , Z, H ± masses, so that a non-trivialC A would persist at the level of the renormalized self-energies. We call such a scheme 'unphysical' because the renormalized treelevel mass parameters are not directly observable quantities then, but implicitly gauge-dependent parameters. This would also apply in a THDM with MS conditions.

Gauge dependence and masses
Since we are interested in the gauge dependence emerging from the electroweak one-loop order, we restrict ourselves to a determination of the masses (and later decays) at the strict one-loop (1L) level, i. e. all two-loop (2L) contributions (whether included in the calculation or not) are regarded as being objects of higher order. Obviously, it is possible to include fully known two-loop orderse. g. O α t,b α s , α 2 t,b , etc. -in the picture as well, but what we comment about the gauge-dependence only applies as long as the two-loop electroweak order is not considered.
In order to discuss gauge invariance, it is convenient to work in a 'physical' scheme, as defined above, i. e. a renormalization scheme where tree-level parameters are directly related to observable quantities, making them gauge-independent objects. Then, the dependence on the gauge-fixing parameters is fully explicit in the radiative contributions and should explicitly vanish (at the considered order) in any observable quantity. A similar analysis would be possible in an 'unphysical' scheme as well, but only after extracting the implicit gauge dependence contained in the tree-level parameters, i. e. after relating them to observable quantities. Thus, we assume below that we are working in a 'physical' scheme, and that e. g. mass parameters are fixed by on-shell renormalization conditions. From the explicit form of the gauge-dependent terms in Eq. (3), we observe that the ξ-dependent contribution to Σ h i h i vanishes at p 2 = m 2 h i -up to theC A term, which only disappears after renormalization in a 'physical' scheme. Thus, no gauge-dependent term contaminates the one-loop correction to the mass of h i : withΣ denoting the renormalized self-energy. If, on the contrary, the Higgs mass is determined by the implicit condition consisting of the (real) pole mass M h i and the width Γ h i (and derived through e. g. an iterative procedure), then the mismatch between M 2 h i and m 2 h i generates a gauge-dependent piece of two-loop order due to the terms in Eq. (3): this is formally a contribution of higher order-as is exhibited in the expansion of Eq. (5b)-and, though a source of uncertainty, can be dismissed in principle in virtue of the expansion. In addition, the off-shell renormalization of the Higgs self-energy requires the introduction of field renormalization, introducing an explicit (and unphysical) dependence on the latter at the level of the Higgs masses as defined in Eqs. (5); this was already discussed in Ref. [40]. We apply DR-conditions on the Higgs fields, so that the dependence on field renormalization is going to translate into a dependence on the renormalization scale.
It is actually common practice [17,18,39,[41][42][43][44][45][46][47][48][49] to include the full propagator matrix in the determination of the Higgs masses. The latter are then obtained from a complex pole search using the following condition on the characteristic polynomial of the inverse propagator matrix: The impact of off-diagonal one-loop self-energies on the Higgs-mass calculation is in general of twoloop order. Yet, such off-diagonal terms intervene at one-loop order in the case where they mix nearly degenerate tree-level states, justifying a more intricate procedure. Let us first consider the case where there is no approximate degeneracy with the state h i (i. e. |m 2 Then, the condition defining the Higgs mass amounts to the following expansion: From the expressions in Eq. (3), we observe that the off-diagonal terms generate new gauge-dependent contributions of two-loop order. In fact, even though the self-energies would be evaluated at the treelevel mass m 2 h i , these gauge-dependent terms would persist, indicating the need for a full electroweak calculation of two-loop order in order to control gauge invariance at this level. Now, let us consider the case where h i and h j are nearly degenerate. Then, , the gauge-violating piecessee Eq. (3)-originating in the off-diagonal terms are still of subleading order at the level of the Higgs mass.
The presence of gauge-dependent pieces of two-loop order in the one-loop-corrected Higgs masses obtained from the pole search can be exploited in order to estimate (a lower bound on) the uncertainty associated with such a determination. This was already performed in e. g. Ref. [49], and running quark masses pole quark masses ξ ξ Figure 1: The gauge dependence of the MSSM Higgs masses is shown for m H ± = 1 TeV, t β = 10. The horizontal green lines correspond to the expansion of one-loop order as in Eq. (4). The blue curves correspond to the ξdependence originating from the diagonal self-energy (the off-diagonal terms are set to 0) as in Eqs. (5). In the orange curves, the momentum in the diagonal self-energy is evaluated at the tree-level mass, so that no gauge-dependence is introduced by this term, and the ξ-dependence originating from the off-diagonal terms-see Eq. (7)-is studied. The plots on the left-hand side are obtained for the renormalization scale µ dim = m t , while those on the right-hand side use µ dim = m SUSY ≡ 1.5 TeV. The input parameters (t β ) have been accordingly transformed so that the two schemes consider the same point in parameter space. Finally, for the Higgs masses at the one-loop order, it is formally equivalent to employ Yukawa couplings and quark masses defined on-shell (dashed curves) or QCD-corrected DR running masses, with a scale corresponding to the Higgs mass (solid curves): the difference between these mass predictions provides an estimate of the magnitude of leading twoloop order effects.
we only consider one example point for illustration in Fig. 1. 1 There, we consider the MSSM with m H ± = 1 TeV, t β = 10, squark masses of ∼ 1.5 TeV for the third generation (∼ 2 TeV for the other two) and electroweakino masses in the range of a few 100 GeV. The ξ-dependence from the diagonal self-energy is found to dominate the ξ-variation (blue curve), when setting the external momentum to the loop-corrected mass value. The effect is sizable for the SM-like Higgs and much smaller for the heavy-doublet states (though in fact of comparable magnitude at the level of the self-energies). The off-diagonal contribution to the ξ-dependence (orange curve) is heavily suppressed by the clear hierarchy in the CP-even sector. The somewhat more-pronounced ξ-dependence from the off-diagonal terms in M h 2 at ξ ≥ 30 is due to the crossing of thresholds in loop-functions (e. g. originating from the G + -G − loop). In addition, all these definitions of the mass-from Eqs. (5) or Eq. (7)-explicitly depend (at two-loop order) on the renormalization of the Higgs fields, as we commented above. This introduces an explicit dependence on the renormalization scale, as can be observed by comparing the left-and right-hand side of Fig. 1. In contrast, the tiny scale dependence in the masses obtained from Eq. (4) (green curves) is implicit and originates in higher orders in the conversion of t β between the two choices of renormalization scale. From this analysis, we see that there is no gain in precision, at the one-loop order, in including the shifts of two-loop order of Eqs. (5) or Eq. (7) since these introduce a non-physical behavior for the would-be observable masses: namely an explicit dependence on ξ and the field renormalization. The impact of these artifacts on the mass determination competes with the magnitude of the leading two-loop corrections of O(α t,b α s )-estimated in Fig. 1 through a variation of the numerical input for the Yukawa couplings (solid vs. dashed lines). Admittedly, the result obtained for low values of ξ (e. g. equal to 1) with Eqs. (5) looks comparatively close to the gaugeindependent result of Eq. (4). However, the scale variation alone evidences a shift of ∼ 2 GeV at the level of the SM-like Higgs mass, hence a sizable uncertainty due to the explicit dependence on the field renormalization; this value exceeds the estimate of Refs. [50,51], obtained without or through a very narrow scale variation. However, we believe that the SUSY scale should be as legitimate as the electroweak scale for the calculation of the Higgs spectrum. Thus, this uncertainty originating in the prescription for the Higgs-pole determination-as already discussed in Ref. [40]-should be added to the previous assessments. We thus believe that the violation of the gauge symmetry by the inclusion of an incomplete electroweak two-loop order indicates that the reliability of the prediction is not improved, but rather that the uncertainty is inflated.

Gauge dependence and decays
At the level of the Higgs decays, loop corrections on the external Higgs leg appear explicitly unless the Higgs fields are renormalized on-shell-in which case corresponding contributions are entirely shifted to the vertex counterterms. Contrarily to the case of the Higgs-mass determination, where only the diagonal self-energy was formally needed at the one-loop order, both diagonal and off-diagonal selfenergies intervene at this order in the decays. As we discussed above, the off-diagonal self-energies are gauge-dependent at the one-loop order even when we set p 2 = m 2 h i : this indicates that these objects-or the associated Higgs-mixing matrix at the one-loop order-are not physical observables (obviously, they are also scheme-dependent), but simply intermediate steps in the calculation of the decay width.
For definiteness, we can consider the example of the decay h i → tt. The vertex corrections (including the on-shell renormalization of the top-quark fields) produce the following gauge-dependent terms (still restricting ourselves to the charged currents): with the vertex functions Here, h i is identified with the field (mass eigenvector) that is associated with the tree-level mass m h i . Yet we allow its external momentum p to be 'free', with p = p t + pt (p t and pt are the external momenta associated with the on-shell quark lines). In the expression above, m 2 h i originates in the Higgs-Goldstone couplings, while p 2 appears in scalar products of external momenta. Obviously, the ξ-dependence in the three-point vertex functions I 1,2 only disappears if the mass appearing in the Higgs-Goldstone coupling coincides with the kinematical one. Thus, employing a loop-corrected mass under the claim that it provides a better description of the kinematical situation would also generate a gauge-violating contribution of two-loop order from these terms: it is consequently arguable whether this choice brings any actual improvement at the numerical level.
The gauge-dependent B 0 and A 0 terms in Eqs. (8) must be combined with the contributions from loop corrections on the external Higgs leg, as prescribed by the LSZ reduction formula. Then, the terms of Eq. (3) generate Adding Eq. (8a) and Eq. (10a), we observe the cancellation of the B 0 and A 0 terms up to the remainder in A mix 1 of Eq. (10b). The latter (generically) disappears only if its prefactor is zero, i. e. the kinematical mass p 2 and the tree-level mass m 2 h i appearing in the Higgs-Goldstone couplings coincide.
The ξ-dependence from the electroweak neutral current vanishes in a similar way; additional terms from the mixing of the external Higgs with internal Z and G 0 cancel out separately up to contributions proportional to (p 2 − m 2 h i ). This cancellation has already been discussed, e. g. in Refs. [17,22]. So far, we have checked how the LSZ reduction formula ensures gauge invariance in the decay amplitude, up to higher-order terms ∝ (p 2 − m 2 h i ). Now, let us turn to the mixing formalism of Refs. [17,18]. The loop-corrected field is defined as for the associated eigenvalue M 2 H k and satisfying the normalization condition [δ ij +Σ ij (M 2 H k )] Z ki Z kj = 1 [18]. By construction, in the absence of degeneracies, at this order. If degeneracies are present, the same formal expansion as above applies (though ill-converging). Gauge invariance is thus satisfied at the strict one-loop order in all the cases. However, there remain gauge-dependent pieces of higher order. The terms ∝ (p 2 − m 2 h i ) of Eqs. (8) and Eqs. (10) which is an object of one-loop order, hence generates gauge-violating contributions of two-loop order. In addition, there are gauge-violating effects beyond those contained in the LSZ expansion, due to the resummation of mixing effects in the mixing matrix and the inclusion of terms from the product of mixing and vertex corrections. The restoration of gauge invariance will thus prove more difficult in this formalism.
In Fig. 2, we show the ξ-dependence in the decay widths Γ[h i → bb] in the MSSM for the point considered in Fig. 1. The latter is threefold, originating firstly in the explicit ξ-dependence of the decay width from the terms of Eqs. (8) and Eqs. (10) when the Higgs mass is set to a loop-corrected value, secondly in a possible processing of mixing and vertex corrections on different footings, thirdly in the implicit ξ-dependence of the mass when it is derived from e. g. Eq. (7). In order to illustrate these features, we plot several definitions of the decay widths: • the solid black curves correspond to the (inclusive) QCD-corrected tree-level decay widths, where however the kinematical factors employ the loop-corrected masses of Eq. (4); • a first version of the decay widths of full one-loop order is shown in solid green: there, the amplitudes are evaluated at the tree-level Higgs mass, while the kinematics employs the mass determination of Eq. (4), leading to an explicitly gauge-independent result (for the heavy-doublet states, this curve is hardly distinguishable from the blue and orange ones); the difference with the solid black curves provides the magnitude of the (non-QCD) radiative corrections; • for the dot-dashed blue curves, the ξ-independent Higgs mass of Eq. (4) is employed everywhere, leading to explicit ξ-dependent decay widths due to the terms of Eqs. (8) and Eqs. (10); this gauge-dependence is found to be rather mild in the example, showing only at the level of the light Higgs through threshold effects; • the dashed red curves show the impact of keeping a term |A vert | 2 in the decay width: the ξ-dependence is sizable for all the states, competing with the absolute magnitude of the electroweak effects; • finally, the dotted orange curve corresponds to a decay width of one-loop order employing a ξ-dependent loop-corrected Higgs mass, i. e. adding to the explicit gauge dependence of the decay width the implicit one contained in Eqs. (5): the effect is mostly relevant for the light Higgs, as could be expected from the impact at the level of the mass determination.
This comparison in particular shows that while the explicit gauge dependence of the decay width at strict (truncated) one-loop order remains rather mild, the predictivity of the calculation can be wasted when vertex and mixing contributions are not consistently combined so as to neutralize the gauge dependence. In addition, we also vary the renormalization scale between m t (left) and m SUSY = 1.5 TeV (right): the associated effects appear to be mostly driven by the dependence on ξ Figure 2: The ξ-dependence in the Higgs decay widths into bb is shown in the scenario of Fig. 1. The solid black curve represents the QCD-corrected tree-level width, while the solid green line corresponds to a full oneloop calculation where the one-loop functions are evaluated at the tree-level value of the Higgs mass. In both cases, the kinematics employs the loop-corrected mass of Eq. (4). For the dot-dashed blue curve, this same loop-corrected mass is also used in the one-loop functions. In addition, a |A vert | 2 piece is kept in the squared amplitude for the dashed red curve. In the dotted orange curve, no explicit contribution to the decay of two-loop order is included, but the momentum is set to the ξ-dependent Higgs mass of Eqs. (5). The renormalization scale is set to m t on the left-hand side, and to m SUSY on the right-hand side. The decay widths for the pseudoscalar state have been omitted since they are essentially identical to those of the heavy CP-even state (lower row).
the Higgs mass (dotted orange curves), hence essentially affect the decays of the light Higgs (where the mass prediction is most sensitive to the variations in ξ). There, even at small ξ, the fluctuations associated with the scale dependence compete in magnitude with the absolute size of the electroweak corrections. For the heavy-doublet Higgs states, the percent-level shift of the decay widths of strict one-loop order (solid green curve) should be seen as belonging to the higher-order uncertainty and remains much smaller than the magnitude of the one-loop corrections, which are dominated by effects of Sudakov type as explained in Ref. [23]. The shift of the tree-level widths (solid black curves) by 10% is associated with the scheme conversion, i. e. the modified value of t β .

Restoring gauge invariance in the non-degenerate case
From the perspective of a strict order-counting, it can appear superfluous to worry about the gaugeviolating pieces described above, since they correspond to a higher order in the expansion. In fact, gauge dependence can be exploited as a means to estimate part of the theoretical uncertainty, as we showed at the level of the Higgs masses in Fig. 1. On the other hand, the gauge-violating effects can numerically dominate the radiative contributions to the decays, because internal cancellations caused by symmetries are no longer enforced. Thus, the violation of the Ward identity in h i → γγ [22,52] can sizably unsettle the corresponding estimate of the decay width. In h i → W W , infrared (IR) divergences do not cancel between virtual QED corrections and soft photon radiation [22,24,53]. We thus believe that it is meaningful for reliable predictions to employ decay amplitudes that satisfy the symmetry principles. A first obvious method neutralizing explicit gauge dependence would simply consist in expanding the loop functions in terms of the external Higgs masses appearing as argument, in the vicinity of the tree-level value, then truncating the expansion at the order achieved in the calculation. In this way, the ξ-dependent terms in Eqs. (5), Eqs. (8) and Eqs. (10) would explicitly vanish indeed. We stress that only the 'amplitudes' or 'form-factors' need to be expanded and truncated in this fashion: the kinematical factors continue to be written in terms of kinematical (i. e. loop-corrected) masses, as e. g. in the green curves of Fig. 2. The usual prejudice against this procedure arises from the fact that the thus shifted argument of the loop-functions displaces or even obstructs internal effects, e. g. thresholds. This is especially true in the case of a light Higgs state (with mass at or below the electroweak scale), since radiative effects can compete with the tree-level. Nevertheless, as we argued above, it is unlikely that directly inputting a loop-corrected mass in the loop functions actually improves the reliability of the calculation, since it then introduces gauge-violating pieces and explicit dependence on the field renormalization. 2 Thus, if one chooses to inject a loop-corrected mass in the decay amplitudes, gauge dependence should be carefully analyzed, either for an estimate of the associated uncertainties or for an attempt at restoring the symmetry. The latter obviously requires the addition of a two-loop order piece absorbing the ξ-dependent terms in Eqs. (5), Eqs. (8) and Eqs. (10). Below, we continue to focus on a 'physical' scheme, since it makes the analysis of gauge dependence more convenient. Here, we observe that: • the ξ-dependence in the three-point functions only appears at the level of the vertex diagrams; thus, this gauge dependence must be neutralized separately; • the ξ-dependence in the two-point functions appears both in vertex diagrams and self-energies; therefore, both objects must be combined consistently in order to neutralize this form of gauge dependence; • the ξ-dependence in the one-point functions appears in the vertex and self-energy diagrams, as well as counterterms (e. g. tadpoles); • the ξ-dependence in the two-and three-point functions originates in the mismatch between the kinematical mass and the tree-level mass appearing in the Higgs-Goldstone couplings, while additional sources intervene at the level of the A 0 functions.
(2)-by substituting a kinematical mass to the tree-level one. This solution works at the level of the three-point functions and sets IR divergences under control. However, applied to the two-point functions of Eqs. (5) or Eqs. (10), it shifts the amplitude by an ultraviolet (UV)-divergent piece (as already noticed in Ref. [49]). Such UV-divergences can admittedly be regularized in an ad-hoc fashion, but this means that UV-logarithms are thus arbitrarily introduced. In addition, the gauge-dependent A 0 term is not removed by this procedure. In order to simultaneously work with loop-corrected masses, preserve gauge invariance, and keep control over the UV-divergences, a more elaborate and consistent procedure needs to be constructed. The cancellation of the two-and three-point gauge-dependent terms makes it clear that the promotion of the Higgs-Goldstone couplings discussed above is a necessary step if one aims at restoring ξindependence. However, we note that this 'upgrade' is just a subset, restricted to these specific couplings, of a larger transformation of the Higgs potential that would impose the loop-corrected Higgs mass as a tree-level value of the new potential. Such a reshaping of the Higgs potential is straightforward to implement as the THDM parameters (or at least specific linear combinations) can be expressed in terms of the masses and mixing angles [54]-see Ref. [55] for such a reconstruction in the context of a THDM with an additional singlet. In addition, it is not possible in general to preserve the properties of a SUSY tree-level Higgs sector when upgrading the Higgs masses to their loopcorrected value, because this operation spoils the connection of the quartic Higgs parameters with the electroweak gauge couplings. Properties of the tree-level spectrum, such as m 2 H ± = m 2 A + M 2 W , are violated at the radiative level (though still constrained by the electroweak symmetry). Relaxing these relations appears as a necessary sacrifice in order to restore gauge invariance in a controlled fashion. As a consequence, the gauge counterterms appearing in the Higgs self-energies in the SUSY context lose any sort of meaning in the generalized framework, and are insufficient in order to absorb the UV-divergences. Therefore, the calculation of the Higgs self-energies needs to be performed in the new framework (after re-definition of the Higgs couplings), that of a THDM (with singlet) with SUSY matter content.
A detailed procedure allowing to map the MSSM onto a THDM+SUSY framework is provided in appendix B and can be applied to the recursive determination of the Higgs masses of Eqs. (5). Then, the Feynman amplitudes employ the value of the loop-corrected Higgs masses as tree-level input parameters both explicitly-when the Higgs states appear in propagators-and implicitly-in the cubic and quartic Higgs couplings. The generated shift is still formally of two-loop order with respect to the original calculation in the SUSY context and, indeed, allows one to restore gauge invariance. Nevertheless, as explained in appendix B, the extension of the renormalization conditions of the (N)MSSM to the THDM (with singlet) framework for the parameters of the Higgs potential is not unambiguous in general. Several choices may appear as 'natural', such as restoring the logarithms of the SUSY self-energies at tree-level on-shell external momenta, or employing logarithms of the same form as those appearing in the gauge couterterms. Failing to identify a physical principle determining these logarithms, we must concede that, while the UV-divergences are now under control, the added UV-logarithms are still largely arbitrary. This arbitrariness can be exploited in the form of a scale dependence (which we denote as µ map below) as a measure of the uncertainty introduced in the mapping. This form of uncertainty replaces that of the field renormalization in the original SUSY calculation with off-shell external momentum, while the ξ-dependence has been neutralized: the result satisfies the symmetry principle.
On the left-hand side of Fig. 3, we show the mapping-scale dependence in the recursive mass determination for the same scenario as in Fig. 1. A comparison of the range of variation with that of Fig. 1 is not meaningful, since the ξ-dependence in the latter is polynomial, hence problematic as ξ → ∞,   while the scale dependence in Fig. 3 is logarithmic and has been freed of symmetry-violating effects. In terms of predictivity, there is no obvious gain with respect to the mass determination of Eq. (4), except perhaps in the control of the uncertainty associated with electroweak corrections of higher-order.
Let us now assume that the physical Higgs masses are calculated in a consistent fashion (either via the truncation method or through the mapping procedure) and are gauge-independent objects. We turn to the question of the Higgs decays. In order to define ξ-independent transition amplitudes, we can employ the same strategy as for the mass determination, i. e. work in a THDM+SUSY framework where the physical masses are tree-level masses. However, such a framework where tree-level and loop-corrected masses coincide is that of a THDM with on-shell renormalization conditions for the Higgs masses: indeed, the Higgs potential of the THDM possesses enough degrees of freedom to allow for an on-shell definition of the masses. In contrast to the mass calculation, there is no arbitrariness in the renormalization conditions for the definition of the decay amplitude (except in the case of Higgsto-Higgs transitions). The precise implementation of this on-shell THDM+SUSY is discussed in appendix C: it can be viewed as a simple switch of renormalization scheme with respect to the MSSM. The arbitrariness of renormalization that we mentioned at the level of the mass determination is now lifted (though the associated uncertainty is still hidden in the input values for the Higgs masses). Corresponding results are shown on the right-hand side of Fig. 3 in the case of neutral Higgs decays into bottom quarks. The explicit uncertainty associated with the matching procedure (variations of the blue curves) is very small. The implicit dependence on the mapping scale (orange curves) when the input values of the Higgs masses are defined by the recursive condition in the THDM mostly matters for the SM-like state. Finally, we note that the predictions for the decay widths obtained with this matching procedure are very close to that of the truncation method in the MSSM (green curve) and spread far less than the ξ-dependent results of Fig. 2.
Nevertheless, there is another source of uncertainty in this mapping/matching procedure of the MSSM onto a THDM, because the Higgs potential of the THDM is not fully determined by the identification of the (four) Higgs masses, leading to an arbitrariness in the choice of the (seven) λ i parameters (out of which three are complex). The simplest choice consists in applying λ 5,6,7 ! = 0, as in the MSSM, and can be justified formally. However, if the radiative corrections to the Higgs masses involve large effects of λ 5,6,7 -type, such as a large splitting between the two neutral heavy-doublet Higgs bosons, the inadequate mapping of these effects onto λ 1,2,3,4 can lead to possibly large spurious Higgs-to-Higgs corrections: this means that the exact form of the THDM must be carefully chosen in order to be consistent with the Higgs spectrum. As yet, we do not have a systematic recipe to optimize this selection. The latter would require assessing several Higgs-to-Higgs transitions in order to shape the radiative Higgs potential in a more realistic way.

Gauge-dependence in the (near-)degenerate case
The procedures that we discussed above, restoring explicit gauge invariance in the renormalized diagonal self-energies or the transition amplitudes, are well-defined only in scenarios where Higgs states do not receive large mixing effects at the radiative level (i. e. in scenarios where the LSZ expansion applies). Now we consider the near-degenerate case. The main difficulty consists in defining the mixed state in a gauge-invariant way, because of the gauge dependence present in the off-diagonal self-energies. Let us focus on a two-dimensional near-degenerate subspace, generated by the tree-level fields h i and h j (larger degenerate sectors follow the same logic). Then, the pole equation reads The right-hand side-of two-loop order-is not neglected because both factors in the left-hand side are of one-loop order each (since m 2 h i − m 2 h j ∼ Σ h i h j (p 2 ) and the difference between p 2 , m 2 h i and m 2 h j is of one-loop order when p 2 coincides with the pole mass). At this leading (non-trivial) order, one can freeze the momentum in the self-energies, e. g.
Then, the only gauge-dependent pieces in Eq. (11) appear in the off-diagonal self-energies and are due to the term ∝ (p 4 , such terms are formally of three-loop order and could be explicitly set to 0 in the self-energy. This would again imply an arbitrary regularization of the associated UV-divergence, which, as before, may be translated into a scale uncertainty-equivalently, the ξ-dependence could be kept and varied in order to estimate the associated uncertainties. The corresponding object (with possibly neutralized ξ-dependence) is denoted asΣ Another feature in the choice of momenta as presented above is the independence of the thus defined Higgs masses from field renormalization, hence the absence of corresponding uncertainties. We can then consider the effective mass matrix in the degenerate sector, which is symmetric, hence diagonalizable in an orthogonal basis with eigenvaluesm 2 H k , and eigenvectors H k = S ki h i + S kj h j . This definition of the masses and fields at the one-loop order essentially extends the truncation procedure to the degenerate case. We stress that the inclusion of the offdiagonal elements is only meaningful because m 2 Indeed, in the non-degenerate case, the off-diagonal element is a piece contributing at an incomplete higher order, hence it only increases the uncertainty in the determination of the Higgs properties. Setting its argument to the average squared mass, or the explicit neutralization of its ξ-dependence also become transformations of this object by one-loop effects, depriving it of any quantitative meaning. Of course, assessing the exact point of transition between the near-and non-degenerate regimes remains largely arbitrary.
Further momentum-dependent corrections (though formally subleading) may be included as 'diagonal' effects according to is again re-defined by a shift of two-loop order absorbing such dependence. This can be achieved in the context of a THDM+SUSY with tree-level fieldsH k absorbing both the tree-level rotation by U n and the loop-level rotation by S of the MSSM. Still, the imaginary parts in S induce further complications for the mapping, so that the complex rotation S is conveniently replaced by a real rotation S minimizing the size of the off-diagonal term in S · M 2 eff · S T . Again, the form of the Higgs potential should be carefully chosen so that mass corrections are appropriately mapped.
Having defined the external states in an explicitly gauge-invariant way, we may now consider the decay amplitudes. In the 'truncation' approach, one can define the gauge-independent object as if both h l and h i belong to the same degenerate subsystem, and Σ eff h i h l = 0 otherwise. Alternatively, it is again possible to embed the Higgs masses in an on-shell model with SUSY matter content. The only difference with respect to the non-degenerate case is that in the on-shell model, the S ( ) -rotation must be included in the definition of the mixing matrix U n → S ( ) · U n , so that the tree-level fields differ in the SUSY model and its on-shell counterpart. TeV, ϕ At ∈ 0, π 2 (upper row), then ϕ At = 3 π 20 ≈ 0.47. Upper left: the mass-splitting vs. mixing is shown: in green, the mass-splitting between the CP-even and CPodd heavy scalars at the tree level; in blue, the mass-splitting between the two diagonal entries of the effective mass matrix, as defined in Eq. (12); in orange, the mixing term of Eq. (12) (in absolute value). Upper right: the mass-splitting vs. mixing is shown after rotation by an appropriate real orthogonal matrix: in green, the mass-splitting between the diagonal elements of the rotated effective mass matrix; in orange, the magnitude of the subsisting off-diagonal element-the effective mass matrix being complex, it cannot be fully diagonalized by a real orthogonal matrix, but the off-diagonal element can be minimized.  12). The dashed lines are obtained with self-energies that are calculated in a THDM where the Higgs potential takes tree-level masses that coincide with the loop-corrected mass value that are injected as external momenta and tree-level fields including the S rotation. The parameters of this THDM are determined from the conditions of Eq. (46) together with the (arbitrary) constraints λ 1 ≡ λ 2 , λ r,i 6 ≡ 0, λ r 7 ≡ 0. The dependence on µ map originates in the freedom of scheme for the fixing of the renormalization conditions of these Higgs-potential parameters and can be seen as a symmetry-conserving uncertainty on the mass-determination.
We illustrate this discussion for the degenerate case with Fig. 4. We consider a CP-violating scenario with t β = 30, m H ± = 0.5 TeV, m SUSY = 1 TeV, A t = 3 TeV, ϕ At ∈ 0, π 2 -likely to exhibit several phenomenological shortcomings (e. g. direct searches for heavy Higgs bosons [8], B-meson decays [56], electric dipole moments [57]), but solely presented here in order to exemplify the impact of gauge dependence in degenerate systems. The phase of the trilinear stop couplings induces a mixing between the (tree-level) CP-even and CP-odd neutral heavy-doublet components, which are close in mass due to the SU (2) L -symmetry. This CP-violating mixing strictly appears at the radiative level (the tree-level MSSM Higgs sector is still CP-conserving) and does not involve the gauge sector: the mixing term is thus automatically ξ-independent, without need of further manipulations. We then consider the effective mass matrix M 2 eff of Eq. (12) and compare the mixing with the masssplitting between the diagonal entries: this is depicted in the plot in the upper left-hand quadrant. The degeneracy at the tree level (green curve) is partially lifted by the diagonal loop corrections (blue curve), but the mixing (orange curve) still competes in the vicinity of ϕ At ∼ 0.5. We then introduce the rotation matrix S (we omit ) that minimizes the size of the off-diagonal entry of M 2 eff . It is clear that a complex matrix S could fully diagonalize M 2 eff . However, in order to recast the system onto an effective THDM, it is more convenient to introduce a real orthogonal S, so that a subsidiary offdiagonal piece remains as a tribute to the imaginary parts in M 2 eff . As is shown on the upper righthand side of Fig. 4, this term is subleading and can be neglected in the mass determination. Next, we focus on the point ϕ At = 3 π 20 , where S shows a mixing angle of about π 4 . The plots in the lower row of Fig. 4 illustrate the ξ-dependence in the mass determination. The dotted lines correspond to the treelevel masses. The solid lines correspond to the eigenvalues of the (2×2) mass system while the dashed lines are obtained from the diagonal entries after rotation by the (gauge-independent) matrix S: the good agreement between both approaches proves the reliability of the second one. Then, the horizontal lines are obtained from the effective mass matrix M 2 eff and are, by construction, ξindependent. The curves showing a ξ-variation are obtained by retaining full momentum dependence in the MSSM self-energies: as before, this approach is both gauge dependent and renormalizationscheme dependent-the latter is made obvious by the impact of the renormalization scale µ dim . If the amplitude of the variations in ξ and µ dim is interpreted as the uncertainty on the electroweak corrections, we see that the latter competes with the mass shift between tree-level and one-loop result, which is problematic. At low ξ, the scale variation of the masses still reaches ∼ 50% of the size of the loop-corrected mass-splitting.
In Fig. 5, we extend the calculation to a THDM+SUSY framework as explained in appendix B, so that the self-energies can receive loop-corrected mass values as external momenta without violating the gauge symmetry. Due to the degeneracy in the MSSM, the S rotation is included in the definition of the tree-level Higgs fields of the THDM. The parameters of the tree-level Higgs potential are defined from the loop-corrected mass values according to Eq. (46), together with the (arbitrary) constraints λ 1 ≡ λ 2 , λ r,i 6 ≡ 0 and λ r 7 ≡ 0. The scale µ map encodes the arbitrariness in the associated renormalization conditions and can be seen as a measure of the uncertainty associated with the mapping procedure. The resulting self-energies in the THDM are ξ-independent and differ from the MSSM ones by a shift of two-loop order. The associated masses are obtained at the strict one-loop order from the THDM self-energies projecting onto the tree-level Higgs-field directions of the THDM (the off-diagonal term is subleading by construction): they are displayed in dashed orange/brown and show a variation of ∼ 1 GeV for µ map ∈ [M W , m SUSY ]. This uncertainty, from which electroweakviolating effects have been cleansed, is much smaller than the one observed in Fig. 4 which is triggered by ξ-dependent contributions. The iteration procedure is stable (the shift is hardly noticeable after one iteration). Finally we note that the mass values are consistent with those of the truncation method (solid blue lines).
Another mixing scenario in the NMSSM is presented in appendix D.
As a concluding remark on the gauge dependence, we have seen that the higher-order terms introduced in the determination of observable quantities, via e. g. not expanding and truncating the momenta in the radiative corrections, are liable to spoil the quality of these calculations. We note, however, that sizable deviations from gauge-independent definitions only appear for large values of the gauge-fixing parameter. In particular, the choice ξ = 1 always seems to produce results that are 'accidentally' close to the gauge-independent definitions, so that the impact of the gauge dependence should remain small at the numerical level for calculations performed in the 't Hooft-Feynman gauge. Nevertheless, the explicit dependence on the field renormalization-associated with the UV-regularization of self-energies away from their mass shell-obviously makes noticeable contributions to the theoretical uncertainty.

Logarithms in loop-order Higgs mixing
The analysis of the gauge dependence in Higgs decays has shown that splitting the loop corrections between mixing and vertex contributions to the decay amplitudes is a purely artificial procedurewhich was already clear from the fact that it is scheme dependent. While this separation is not a fundamental problem per se, it could lead to misleading results through the inclusion of incomplete higher orders. We aim to illustrate this below by comparing decay amplitudes in the limit of heavy-doublet states, at a scale where the electroweak-violating effects are subleading, hence the global SU (2) L -symmetry should be approximately satisfied (up to breaking terms scaling with v and suppressed by the heavy Higgs mass).

Fermion-loop contributions to the Higgs mixing
For simplicity, we consider the MSSM with decoupled SUSY particles and m H ± much above the electroweak scale. This is akin to a THDM framework of type II in the decoupling limit. In particular, the mixing in the tree-level Higgs sector can be approximated by a β-angle rotation: where h 0 u,d , a 0 u,d and H ± u,d are the CP-even, CP-odd and charged Higgs fields (respectively) in the gauge-eigenstate basis. In addition, m H ± ≈ m H ≈ m A m h ≈ M Z . We will exploit this approximation below for the purpose of deriving simple analytical formulae capturing the main features of the calculation, although our numerical results still consist in a full calculation of one-loop order.
We further target corrections associated with quark Yukawa couplings of the third generation Y t,b , i. e. top and bottom loops. For a heavy SUSY sector, the sfermion loops contributing at the same order can be regarded as constant with respect to variations of the Higgs external momentum p, and we do not document them further (though they may involve large logarithms of the form ln m 2 SUSY /M 2 EW ). As before, we work in a scheme where the charged-Higgs mass, as well as the W -, Z-, and fermion masses are renormalized on-shell, while other parameters receive a DR-renormalization with the ultraviolet regulator set to m t .
The leading contributions to the renormalized Higgs self-energies can be summarized as follows: Here, terms ∝ v 2 are neglected, unless p 2 is set to a value of electroweak magnitude (like inΣ hh (p 2 ) for the determination of M h ). The scale M 2 EW ∼ M 2 W ∼ M 2 Z ∼ m 2 t is characteristic of the scheme. From Eqs. (15), we can extract the leading contributions to the heavy Higgs masses and mixing. In the considered scheme (with the charged Higgs renormalized on-shell), these one-loop corrections to the masses are of subleading order O Y 2 t,b M 2 EW /(16 π 2 ) in the limit M 2 EW m 2 H ± . On the other hand, the contributions to the wave-functions and mixing of the heavy-doublet states can be sizable and matter at the level of the Higgs decays: It is common practice [17,22] to include the wave-function corrections and the mixing in the CP-even sector within a loop-corrected mixing matrix, while the mixing with the Goldstone bosons is kept separately at the diagrammatic level. One could regret that the SU (2) L -symmetry is thus explicitly broken by the formalism, but this should not be cause for any inconsistency, in principle. The mo-tivation behind absorbing the Higgs mixing within the definition of the external field rests with a resummation of mixing effects in the case where tree-level states are almost degenerate.
Similarly, the fields of the CP-odd and charged Higgs receive loop corrections according to A → Z A A and H ± → Z H ± H ± with Z A Z H Z ± H , while the mixing with the Goldstone bosons is kept apart. Below, we continue to use the notations h ∼ h 1 or H ∼ h 2 indifferently.

Two-loop O Y 4 q corrections to the Higgs masses
The choice of processing the mixing in the CP-even sector differently than in the CP-odd or charged sectors intervenes first at the level of the mass calculation, when corrections of two-loop order are considered. Indeed, the off-diagonal corrections of one-loop order contribute only from this order on, but the two-loop effects of O Y 4 q belong to those that are commonly included in the mass determination. From the perspective of an expansion up to the two-loop order one has We expect the mass-splitting between heavy-doublet states to be protected by the SU (2) L -symmetry, i. e. contributions of order m 2 H ± ln 2 m 2 H ± /M 2 EW , ln m 2 H ± /M 2 EW , 1 should disappear in the mass difference, leaving only terms of the form M 2 EW ln 2 m 2 H ± /M 2 EW , ln m 2 H ± /M 2 EW , 1 . 3 As the charged Higgs mass is renormalized on-shell in our scheme, the radiative corrections should satisfy the above property directly at the level of the renormalized self-energies-i. e. the counterterms should automatically remove the SU (2) L -symmetry-violating terms.
Let us first examine the genuine diagonal two-loop pieceΣ 2L ii . 4 As the two-loop electroweak calculation is incomplete in the MSSM, 5 we work in the gaugeless limit in this subsection. The full two-loop self-energyΣ 2L ii can be decomposed into several contributions, starting with genuine oneparticle irreducible (1PI) two-loop self-energy diagrams Σ 2L,1PI ii , 1PI one-loop self-energy diagrams with counterterm insertion Σ 1L×CT ii , two-loop counterterm diagrams involving a pair of one-loop counterterms Σ CT×CT ii and the genuine two-loop counterterm contribution Σ 2L,CT ii . Each piece can be expanded in the heavy-mass limit, providing the following leading terms of O Y 4 q : 6 O(Y t 4 ): terms ∝ p 2 log 2 (p 2 ) + terms ∝ p 2 log(p 2 ) where δ OS,1 mt = 1 = δ OS,1 m b are related to the DR/on-shell conversion of fermion masses. The convergence of this expansion is illustrated on the left-hand side of Fig. 6 for both Σ 2L,1PI ii and Σ 1L×CT ii , considering only the terms of order Y 4 t , at t β = 3. TwoCalc [33] and TLDR [35] have been utilized for analytical manipulations of the two-loop diagrams, and numerical results for the two-loop integrals are obtained via TSIL [34]. These terms are exactly mirrored by the corresponding contributions to the charged Higgs self-energy which, in our renormalization scheme, appears in the two-loop counterterms to the neutral self-energies. The resulting cancellation is shown on the right-hand side of Fig. 6, resizing the diagonal contributions to the renormalized Higgs self-energies to electroweak magnitude (as we expected). Now, let us turn to the one-loop squared terms of Eq. (18). The diagonal one,Σ 1L ii dΣ 1L ii /dp 2 , obviously satisfies the property that we stated before: when replacing the renormalized one-loop self-energies by their approximate expressions of Eqs. (15), we find only terms of electroweak size.
The off-diagonal term (last term of Eq. (18)) behaves differently, however, and generates a leading contribution of order m 2 H ± ln 2 m 2 H ± /M 2 EW : This applies to the CP-even sector (h i ≡ H, h j =i ≡ h), but also to the CP-odd (h i ≡ A, h j =i ≡ G 0 ) and charged (h i ≡ H ± , h j =i ≡ G ± ) ones, so that the SU (2) L -symmetry is still not (strongly) violated.
In the renormalization scheme under consideration, with an on-shell charged Higgs, this also means that the two-loop counterterm contribution Σ 2L,CT ii should contain the term −Σ 1L G + H − /m 2 H ± , and balance the one-loop squared term directly at the level of the renormalized self-energy. To be explicit, the renormalization condition e Σ H + H − m 2 H ± = 0 should be applied to the full expression, including off-diagonal contributions, at each order of the loop expansion. At the two-loop order, this fixes the on-shell counterterm to The renormalization scheme at the one-loop order is not repeated here (see e. g. Ref. [39]). We note that-although e Σ 1L is requested to vanish-the term of the momentum expansion, , still contributes a finite shift due to the imaginary parts (see e. g. Ref. [48] for details). In addition, m 2 G ± = 0 in the gaugeless limit, while, in the presence of gauge contributions, the combination of gauge and Goldstone self-energies also restores the denominator 1/m 2 H ± , as explained in Refs. [17,22] for the neutral case. 7 A problem arises when the off-diagonal contribution to the charged-Higgs mass counterterm is overlooked, as seems to have been the case in some earlier calculations-see e. g. Eq. (11) in [48]-, leading to an explicit violation of the SU (2) L -symmetry and an inconsistent order O Y 4 q . While this issue is made obvious for heavy-doublet states by the analysis of logarithms, it actually points at a conceptual shortcoming in the mass calculation for all Higgs states. Even in a calculation of one-loop order, the partial inclusion of the CP-even off-diagonal contribution, but not of the corresponding Higgs-Goldstone mixing in the charged and/or CP-odd sector, in fact worsens the quality of the mass prediction for the heavy states with respect to a simple truncation at strict oneloop order, since it introduces unwarranted SU (2) L -violating effects. Nevertheless, the charged Higgs counterterm has a limited impact on the mass prediction for the SM-like Higgs in general, due to the decoupling of this latter state from the heavy doublet. Therefore, only in scenarios where m H ± ∼ v and a large h-H mixing develops could this issue have any numerical consequences for M h , but we note that there is no longer any logarithmic enhancement in such scenarios either. In addition, the partial result obtained in the gaugeless limit is not really suited to explore such configurations.
In Fig. 7, we first consider the squared mass-splitting between the heavy-doublet states of the scenario of Fig. 6. The plot on the left-hand side of the first row shows this squared mass-splitting for the heavy CP-even and charged Higgs. The result of FeynHiggs-2.16.0 with the setting FHSetFlags[4,4,0,1,0,0,0,3] is displayed in orange and shows a comparatively large and growing squared mass-splitting. Our result is the solid blue curve, corresponding to considerably smaller values. We are able to roughly reproduce the result of FeynHiggs when omitting the off-diagonal  self-energy-squared term in the mass counterterm of Eq. (21) (dotted blue line; also the CP-even mixing is then included in the pole search). This mismatch is thus explained by the large SU (2) Lviolating terms of order O Y 4 q scaling like m 2 H ± ln 2 m 2 H ± /M 2 EW , ln m 2 H ± /M 2 EW , 1 that have been erroneously introduced in the calculation. It is remarkable that these 'wrong' two-loop effects then dominate the one-loop corrections, which, in our prediction, explain most of the squared masssplitting. The plot on the right-hand side examines our result more closely. The bulk of the variation both in M 2 H − M 2 H ± and M 2 A − M 2 H ± originates in one-loop effects of gauge type, dominated by a term scaling linearly with the mass of the heavy doublet, ∼ α M Z m H ± (plotted in green). The following plots further analyze the remainder in the CP-even (second row) and the CP-odd cases (third row). On the left, we compare the magnitude of the remaining one-loop contributions (i. e. beyond the gauge term scaling linearly; in green), the two-loop corrections of order O Y 4 q in the effectivepotential approximation (black curves) and in the gaugeless limit (blue curves). It thus appears that the genuine effects of O Y 4 q are considerably smaller than what the off-diagonal contribution in the CP-even sector hinted. In addition, we note that the two-loop terms obtained in the effective potential approximation are ∼ 30% away from those obtained in the a priori more complete gaugeless limit. 8 Finally, the plots on the right-hand side of Fig. 7, second and third  As a closing remark, we note that, while the order Y 4 q is known, it may not be the most relevant type of two-loop contribution for heavy-doublet states. Indeed, already at the one-loop order we mentioned that the radiative corrections to the squared masses are dominated by a linear term ∼ α M Z m H ± due to electroweak gauge effects. It is thus likely that two-loop corrections to the masses of the heavy doublets are dominated by contributions of the same type. This does evidently not spoil the relevance of our previous remarks as to a consistent treatment of the one-loop squared term, since the 'erroneous' terms of order Y 4 q m 2 H ± ln 2 m 2 H ± /M 2 EW , ln m 2 H ± /M 2 EW , 1 first need to be put under control before the linear contribution becomes apparent.

Application to the decays into quarks
Now, we focus on the decays H, A → tt, H + → t RbL , processes that are mediated through Y t at the tree level. The vertex corrections at the same order as considered in Eqs. (15) read Following the LSZ reduction, we then add the contributions from self-energy diagrams on the Higgs leg, provided by Eqs. (16), and obtain at the one-loop order: 8 We stress that these two-loop contributions contain large logarithms of the type ln(m 2 q /M 2 EW ) in our scenario with a heavy SUSY spectrum; the threshold corrections for the heavy SUSY sector are explained for a SM-EFT matching in Refs. [63][64][65][66][67][68][69], and for a THDM-EFT matching in Refs. [70][71][72][73][74][75], but this feature goes beyond the scope of this article.
The propagator in 1/p 2 for the Goldstone bosons is justified by the combination with the weak gauge-boson effects (see e. g. section 4.3 in Ref. [17]). We observe that (at the considered order) the radiative corrections preserve the SU (2) L -symmetry for the heavy-doublet states: their decay widths are identical at the tree level (up to corrections of O M 2 EW /m 2 H ± ) and continue to be so at the one-loop order. This is indeed what we expect for M 2 EW m 2 H ± : massive states are hardly sensitive to the electroweak symmetry-breaking effects.
Obviously, we would formally obtain the same expansion as that of Eqs. (23) at this order when employing the mixing formalism instead of the LSZ reduction. Therefore, any deviation from comparable values in the decay widths in the CP-even, CP-odd or charged sector would have to originate in higher-order terms. Nevertheless, on the left-hand side of Fig. 8, we observe a sizable mismatch between the decay widths of the CP-even state (blue curves), on the one hand, and the CP-odd (green) and charged Higgs (orange), on the other hand, when employing the mixing formalism. This discrepancy increases at large masses, where we would expect the SU (2) L -violating effects to become smaller. The results of FeynHiggs (dotted curves) show a similar behaviour but differ from ours (solid lines) because of several higher order pieces (|A vert | 2 in FeynHiggs, resummation of Sudakov double logarithms on our side, etc.). At this level, we are forced to consider the difference between the blue and green curves of Fig. 8 as setting the magnitude of the higher-order uncertainty in the calculation, an uncertainty close to 100% at sufficiently high mass! In fact, the origin of this issue entirely rests with the different procedures that are employed in the CP-even, CP-odd and charged sectors in order to include the loop corrections on the external Higgs leg. As expected, these variants formally deviate by two-loop effects of order O Y 2 q α s , which however prove to be numerically significant. Indeed, the Higgs mixing in the CP-even sector has been defined at the same level as the Higgs mass, i. e. the Yukawa couplings are given at the electroweak scale, as prescribed by the renormalization scheme. On the contrary, in the CP-odd and the charged sectors, the Higgs mixing is incorporated at the same time as the vertex corrections, i. e. explicitly for the decay. At this level, the QCD analysis (see e. g. Refs. [76,77]) makes it clear that QCD logarithms can be resummed in the calculation of the inclusive width (i. e. including gluon radiation) through the incorporation of the QCD running in the Yukawa couplings up to the scale of the decaying Higgs state: QCD logarithms are solely of ultraviolet type at the level of the inclusive width. This justifies the use of running Yukawa couplings defined at the high scale in the decay width, and in particular in the CP-odd and charged mixing contributions. As announced above, the latter thus differ from the mixing contribution implemented in the CP-even sector by terms of order O Y 2 q α s (and higher orders). Consequently, the resulting SU (2) L -breaking effect is purely artificial and of higher order. However, the QCD analysis indicates that the recipe employed in the CP-odd or charged Higgs decays is more reliable in this specific case than the mixing procedure in the CP-even sector.
On the right-hand side of Fig. 8, the same decay widths are shown, but with a consistent combination of mixing and vertex contributions for all the states. Then the decay widths for the CPeven and CP-odd states are roughly identical. A small shift persists between neutral and charged states. The latter has a physical meaning: it already appears in the coefficients of the Sudakov double-logarithms [22], which indeed differ at the level of the exclusive widths (i. e. discarding Wand Z-radiation, though inclusive with respect to QCD and QED radiation), as a tribute to the SU (2) L × U (1) Y → U (1) em breaking.

Application to the Higgs decays into weak gauge bosons
The coupling of the heavy-doublet states to weak gauge bosons vanishes in the decoupling limit. It is indeed difficult to build an electroweakly invariant operator coupling exactly one doublet scalar to two triplet or singlet vectors: this requires a breaking of the SU (2) L -symmetry, which, as we argued before, should appear as a subleading effect for states with a mass substantially larger than the electroweak scale. Consequently, the decay widths for H, A → W W, ZZ and H ± → W ± Z exactly or approximately (in the CP-even case) vanish at tree-level. On the other hand the vertex corrections involve fermion loops that lead to unsuppressed logarithms: The absence of a logarithmic contribution in the pseudoscalar case points at an apparent deviation from the SU (2) L -correspondence between the heavy-doublet states. In addition, the existence of logarithmic terms for the CP-even and charged states contradicts our previous comment that, in the SU (2) L -conserving limit, no operator mediating these decays can be written. However, the result is as yet incomplete at the considered one-loop order and should also include the mixing contribution. Following Eqs. (15), we have As expected, the contributions to the vertex in Eqs. (24) and to the mixing in Eqs. (25) exactly cancel each other. On the other hand, if mixing and vertex are not processed at the same level, e. g. due to the use of different parameters at the loop order, then the cancellation is imperfect and large spurious effects, though formally of higher order, develop. Such misleading effects are exacerbated for heavy states as the fine cancellation between mixing and vertex encompasses orders of magnitude. Unluckily, beyond the higher orders that are inherent to the mixing formalism, explicit effects of higher order, e. g. O Y 2 t,b α s , that are considered in the mass calculation are also routinely included within the mixing in the CP-even sector. Such terms are going to cause an imbalance in the cancellation between vertex and mixing as long as the vertex corrections are not known to the same order. We illustrate this fact below.
In the upper row of Fig. 9, we present the decay width of the heavy CP-even Higgs into W + W − for t β = 10, considering only corrections of O Y 2 q . The plot on the left-hand side shows the magnitude of the individual contributions to the decay width in the strict LSZ expansion. The tree-level prediction (dotted black) is smaller by orders of magnitude than the mixing (orange) and vertex (dashed green) contributions of one-loop order, which however largely compensate one-another. As the sum of oneloop contributions (dot-dashed red) is negative and larger in absolute value than the tree level, the decay width at truncated one-loop order would be negative at large Higgs masses. In this context, it is legitimate to include the one-loop squared piece (dotted blue), which is expected to supersede the remaining contributions of two-loop order. In the plot on the right-hand side, we display the tree-level decay width (dotted black), and the width of one-loop order including the 1L 2 term (solid red). In dashed orange, the two-loop O Y 2 t,b α s corrections to the mixing have been included: the associated prediction is in excess by a factor 10 with respect to the one-loop result (red): however, this enhancement is most likely not a genuine effect, but an artifact caused by the non-cancellation of mixing and vertex contributions of order O Y 2 t,b α s -since the corresponding order has not been included in the vertex corrections. The purely radiative Γ[A → W + W − ] (in green) is shown as a reference. This comparison indicates that the inclusion of a partial order O Y 2 t,b α s in the mixing in fact worsens the quality of the prediction for the decay width by introducing large symmetryviolating effects (though formally of higher order). It thus appears as misleading to process vertex and mixing contributions in a decoupled way and introduce partial higher orders, which are liable to violate the symmetries.
However, the mixing formalism introduces further higher-order terms due to the factorization and resummation of mixing effects in a loop-corrected mixing matrix. These contributions are liable to blur further the fine cancellation resulting from the symmetry requirements. In the lower row of Fig. 9, we no longer restrict to the order O Y 2 q , but include the full electroweak and SUSY corrections-the SUSY spectrum is still at a scale of ∼ 100 TeV. We consider both the decay widths for H → W + W − (left) and H → ZZ (right). The one-loop corrections (green curves) can reach the magnitude of the Born-level amplitudes, and the one-loop squared term may dominate the widths (red curves). However, these decay widths obtained from the straightforward application of the LSZ reduction all remain comparatively suppressed, as a tribute to the electroweak symmetry. On the other hand, if vertex and mixing contributions are not linearly added, but rather combined via a  loop-corrected mixing matrix, then the imperfect cancellation of these contributions to the decay amplitudes generate pieces of higher-order that break the electroweak symmetry strongly and come to dominate the width at high masses (dashed orange curves). The dramatic enhancement of these spurious effects is due to the kinematic prefactor M 3 H /M 2 V , that the decay amplitude needs to balance by a careful scaling A ∝ M 2 V /M 2 H , which is spoiled by the separate processing/factorization of mixing contributions. Again, this argues against an indiscriminate use of the mixing formalism in physical transitions. Finally, we stress that, even though the consistent combination of mixing and vertex contributions leads to decay widths that are compatible with symmetries, the latter still come with a sizable uncertainty for such rare processes: for, instance, exchanging pole quark masses by QCDrunning masses at the scale of the decaying Higgs-a legitimate shift at the order controlled in the calculation-in the scenario of Fig. 9 typically leads to a reduction of the widths by a factor O (10). This points at the necessity to include two-loop contributions for a reliable assessment of these symmetry-violating channels.
The analysis of the one-loop radiative corrections in the decoupling limit, in a limit where the SU (2) L -symmetry should hold, thus indicates that misleading large effects purely associated with partial higher-order contributions may develop as a consequence of considering vertex and mixing diagrams on a different footing. In particular, the deliberate inclusion of higher-order effects in the mixing without the corresponding terms in the vertex appears as an unfruitful effort (as already suggested by the analysis of the gauge dependence).

Conclusions
In this paper, we have analyzed how terms of higher order introduced in the calculation of the radiative corrections to observables in extended SUSY Higgs sectors could lead to spurious effects in view of the symmetries. As noticeable numerical variations accompany these artifacts, it appears necessary to consider them seriously in the uncertainty estimates. On the other hand, the associated behaviour is unphysical and needlessly burdens the error budget; therefore, we regard it as meaningful to attempt and avoid such symmetry-violating pieces of higher order.
We first discussed gauge dependence in Higgs-mass determinations and decays and explained how setting the arguments of the loop functions away from the tree-level mass values generates pieces that depend both on the gauge-fixing parameter and the field renormalization. We then presented two possible strategies avoiding such undesirable effects: the first one consists in systematically expanding and truncating the amplitudes at the relevant order controlled in the calculation; the alternative one extends the SUSY model by a more flexible structure where the Higgs potential can adjust to the values of the loop-corrected masses. While the former method is more conventional and straightforward to implement, the latter one can be interesting in order to assess uncertainties in a gauge-conserving context, although redundancies in the definition of the Higgs potential limit its efficiency. We also discussed how to define a mixing matrix in a gauge-invariant way in the case of near-degenerate states, a recipe that fails in the non-degenerate case, where gauge invariance is most efficiently enforced by a strict application of the LSZ reduction. In particular, we saw that gauge invariance required a careful combination of vertex and mixing contributions to the decay amplitudes, so that a separation of both, e. g. through the definition of a mixing matrix at radiative order, can be source of inconsistent behaviours.
Then, we focussed on the decoupling limit of the MSSM, where the SU (2) L -symmetry still controls the dominant properties of the heavy-doublet states. We illustrated this analytically, with expressions for the radiative corrections of Yukawa type, as well as numerically, with a full calculation of oneloop order. These arguments allowed us to spot several issues, first in the calculation of the twoloop corrections to the Higgs masses of O(Y 4 q ), then in the implementation of Higgs decays when mixing and vertex corrections are included at different orders or via a mixing matrix defined at the radiative level. The corresponding inconsistencies become large, admittedly because of the choice of a heavy-doublet spectrum, but still point at shortcomings in the general implementation of these observables in all regimes. Observables measuring SU (2) L -breaking effects, such as mass-splittings among SU (2) L -partners or heavy-Higgs decays into electroweak gauge bosons or lighter Higgs states, are particularly sensitive to the introduction of spurious SU (2) L -symmetry-violating terms of higher order, so that a proper control on the symmetries appears as imperative for a meaningful study of such channels.
As we took the restoration of symmetries as our guiding principle, it proved more convenient to work with electroweakly-charged states. However, what we learnt with doublet states can be extended straightforwardly to more exotic Higgs spectra, and even to other fields that are not necessarily renormalized on-shell (such as sfermions or electroweakinos; see e. g. Refs. [78,79]). Furthermore, several additional issues that are not constrained by symmetries can appear in connection with a careless use of mixing formalisms at the radiative order. Double-counting of mixing corrections can thus emerge from e. g. hybridizing the mixing formalism with that of effective couplings for integratedout SUSY sectors (for the latter formalism, see e. g. Ref. [80] and references therein): indeed, SUSY corrections on the external Higgs line are then potentially double-counted.
As a concluding word, we believe that, while it is of course still possible to employ a mixing formalism in calculations of radiative corrections to the Higgs sector, corresponding results should be critically analyzed in order to verify whether the effects that they produce are genuine or just the outcome of symmetry-violating artefacts.
For the off-diagonal block matrix one finds After application of the minimization conditions, the Goldstone boson is obtained as the linear combination G 0 = c β a 0 d − s β a 0 u . The other (three) tree-level mass-eigenstates may be writ- where U n is a unitary matrix. In the absence of CP-violation, the pseudoscalar sector can be diagonalized separately with a rotation of angle β 0 , which is found to coincide with β after application of the minimization conditions. The pseudoscalar state A 0 ≡ s β 0 a 0 d + c β 0 a 0 u takes on the mass The Higgs couplings appearing in loop amplitudes can be summarized as follows: • the neutral symmetric triple-Higgs couplings • the hermitian couplings of one neutral and two charged fields • weak gauge-boson masses: the W -and Z-boson masses are renormalized on-shell, so that the associated counterterms cancel out with the transverse W and Z self-energies evaluated at • the angle β: t β is renormalized in the DR-scheme, The associated contributions are purely of Yukawa type, see Ref. [81] for details.
• the THDM shifts δ( i v 2 ), for i = 1, . . . , 7, are as yet undetermined. In the MSSM, they are set equal to 0, since they are not needed to achieve renormalizability (and could spoil the SUSY relations for non-shifted momenta).

A.3. Counterterms in the neutral Higgs sector
The renormalized self-energies (ignoring field renormalization for now) for the neutral Higgs bosons can be expressed in the following way (in the gauge eigenbase): The renormalized diagonal self-energies of the neutral sector in the mass basis are obtained after rotation by U n :Σ From the explicit calculation of the self-energies and tadpoles, it is possible to extract the ξdependence of these expressions: are UV-finite determine the UV-divergences of the δ( i v 2 ) counterterms as well.
of Eqs. (46) are fundamental in ensuring that the massesm h 0 i (identified with the loop-corrected masses in the MSSM) are tree-level masses in the THDM.

B.2. Definition of the THDM counterterms
Now our problem rests with the determination of the countertermsδ( i v 2 ) in the THDM, or at least of those linear combinations that enter the Higgs mass matrix. The extension should be smooth in the sense that eachδ( i v 2 ) remains an object of two-loop order from the perspective of the counting in the MSSM. The renormalization scheme should also ensure gauge invariance at the level of the renormalized Higgs self-energies evaluated at the corresponding tree-level Higgs mass in the THDM, so that the definition of the Higgs masses is gauge invariant.
It is obvious that all the counterterms δC that are already fixed (masses, tadpoles and t β ) differ by terms of two-loop order between the two models: in the MSSM, the renormalization constants are defined in terms of one-loop one-or two-point functions Q MSSM i that are evaluated with parameters λ MSSM i and m 2 h i . In the THDM, the Q THDM i are computed from λ i andm 2 h i that differ from the MSSM parameters by terms of one-loop order. Formally: Similarly, the Higgs self-energies Σ THDM h i h j (m 2 h i ) and Σ MSSM h i h j (m 2 h i ) differ by a shift of two-loop order. If we wish to preserve the ξ-independence of the Higgs self-energies evaluated in the THDM at on-shell values of the external momentum, then Eqs. (45) explicitly constrain the ξ-dependence of theδ( i v 2 )-s. Similarly, the UV-divergences can be worked out explicitly from the condition that the renormalized self-energies (without field renormalization) evaluated at on-shell external momentum are UV-finite in both the MSSM and the THDM. The differences between the self-energies in the two models originate from different sources: shifted external momentum; shifted Higgs couplings; shifted Higgs masses in internal lines. It is convenient to split these contributions between bosonic diagrams where the modified Higgs sector directly intervenes, and fermionic contributions where only the shifted external momentum matters: The matrices H ξ and C were defined in Eqs. (45) while the (symmetric) others are determined by the following entries in the basis (h 0 d , h 0 u , a 0 ): (H 6 ) dd = r 6 (6 1 + 3 3 + 4 4 + 5 r 5 ) s 2β − r 2 6 (5 + 7 c 2 β ) − 2 r 6 r 7 s 2 β + 5 i (H 6 ) ud = r 6 6 1 c 2 β + 3 3 + 2 4 (1 + c 2 β ) − (7 r 6 + 5 r 7 ) s β c β + ( r • UV-divergences from fermionic contributions (only terms of the third generation are displayed): Matching the MSSM with an on-shell THDM+SUSY In the previous section, we have discussed how it was possible to extend the Higgs self-energies of the MSSM by a shift of two-loop order embedded within a THDM with SUSY matter content. The corresponding procedure was meant to restore gauge invariance in the determination of the loop-corrected Higgs masses. Now, we assume that these masses have been determined-either in the MSSM via the truncation method or in a THDM context via the method of appendix B-and are gauge-invariant quantities, and we wish to consider particle scattering and decays involving the Higgs bosons of the MSSM in a fully on-shell context. This is not possible in the MSSM stricto sensu since the quartic scalar couplings are associated to gauge couplings. Yet, since this structure of the couplings is not preserved by the radiative corrections, we can instead choose to work in an effective field theory (EFT) with identical field content and symmetries, but with Higgs fields that are renormalized on-shell. This EFT is obviously a THDM+SUSY with λ i -parameters satisfying the conditions of Eqs. (46). However, contrarily to the case discussed in the previous section, theδ i -s of this EFT are no longer shifts of two-loop order, since they absorb the full radiative corrections to the Higgs masses (on-shell condition). In addition, this EFT is not suitable for the determination of the Higgs masses in the MSSM since, in the on-shell THDM+SUSY, these masses are free input. Yet, the matching conditions between the MSSM and the on-shell THDM+SUSY are akin to the usual scheme-matching conditions between two renormalization schemes applied to the same model.

C.2. Matching conditions
Now we turn to the matching conditions. On the side of the MSSM, the parameters of the Higgs sector are M 2 W , M 2 Z , m 2 H ± , t β and G F (encoding the electroweak v.e.v. v). For the THDM, we have m 2 H ± , λ 1,··· ,7 , ϕ 5,6,7 ,t β , G F . The mixing angles (in the CP-even sector and/or CP-violating) are an output of the potential and thus a priori different between the two models. The matching conditions at the tree level are trivial and provide λ THDM, (0) i = λ MSSM i , with λ MSSM i given in Eqs. (27). At the one-loop order, the requirement that the physical Higgs masses coincide in both models provide the six matching conditions of Eqs. (46). Furthermore, we may consider the transition H ± → t LbR as a condition determiningt β : and similarly for the MSSM witht β → t β . Then, we observe that A THDM 1L and A MSSM 1L are identical up to a shift of two-loop order: indeed, the only difference between the two amplitudes comes from the modified Higgs potential, but in a quantity of one-loop order, we can employ λ THDM, (0) i without spoiling the expansion. Therefore, we can chooset β = t β as a matching condition of one-loop order. The same analysis in e. g. h 0 i → bb shows that we may also identify the mixing angles provided the off-diagonal mass-counterterms are defined as in Eq. (61)-otherwise, the tree-level contributions would disagree by an effect of one-loop order, forcing a different choice of mixing angles.
At this point, we still have four unconstrained degrees of freedom in the THDM Higgs sector. Fixing them would require considering Higgs-to-Higgs transitions, e. g. h 0 2 → h 0 1 h 0 1 or H + h 0 i → H + γ. However, the same argument as for the mixing angles shows that this choice is arbitrary. Indeed, since λ THDM i − λ MSSM i is formally of one-loop order, the distribution of corresponding finite effects between tree level and counterterms only amounts to a formal shift of higher order. For instance, if one chooses λ THDM 5,6,7 ≡ 0, it is always possible to defineδλ 5,6,7 so that the Higgs-to-Higgs transitions chosen as matching conditions are satisfied. Therefore, we are left with the same arbitrariness as in the previous section with respect to the choice of potential in the EFT. On the other hand, the renormalization conditions are well-defined in this scheme-conversion approach, so that the Higgs transitions that are studied in this framework are no longer explicitly subject to the renormalizationscale dependence (this dependence or the gauge one is still implicitly present within the Higgs masses that are used as input for the matching).
As a final remark, we stress that the transition amplitudes evaluated in this on-shell THDM+SUSY framework are automatically gauge-invariant and still a variation of two-loop order with respect to the corresponding MSSM transition amplitudes.
Our renormalization scheme in the doublet sector is still determined by 'physical' conditions, connecting the Higgs sector to the masses of the electroweak gauge bosons. On the other hand, the singlet and singlet-doublet-mixing parameters are renormalized DR. In these circumstances, it is not completely trivial whether gauge invariance in the Higgs masses can be discussed without relating the DR parameters to physical quantities. As it turns out, this subtlety does not matter for parameters associated with a singlet state since the latter does not couple to the electroweak gauge sector. Technically, looking at Eq. (3), DR-renormalization conditions could only spoil the cancellation ofC A at the level of the renormalized self-energies. However, this cancellation is fully ensured by the counterterms in the doublet sector, so that, while our scheme is not 'physical' in the strict sense, it still displays the same properties in view of the gauge dependence. The variation of |A κ | modulates the diagonal mass input for the singlet component, hence the mass-splitting with the SM-like component. The radiative mixing also varies with the spectrum. In Fig. 10, we first study the magnitude of the mixing effect between the two light CP-even states (upper left-hand corner). To this end, we consider the effective mass-matrix of Eq. (12). Its diagonal entries are gauge invariant as explained in section 2.2. The off-diagonal self-energy contains a gauge dependence of three-loop order, which we will not attempt to neutralize here (its effect is numerically negligible). Then we plot |mixing/diagonal splitting| 2 : we observe that this quantity is 'large' for 100 GeV ≤ |A κ | ≤ 260 GeV, hence that the mixing generated at radiative order is important in this range. It can be safely neglected in the mass determination outside-even though the mass-splitting is of electroweak order in the full scenario. In the upper right-hand corner of Fig. 10, we show the magnitude of the singlet composition of the Higgs fields at the tree level (dashed curves) and after diagonalization of M 2 eff (using a real orthogonal matrix S; solid curves): we see that the maximal singlet-doublet mixing is displaced from |A κ | ∼ 150 GeV to ∼ 250 GeV by the radiative effects.
In the lower left-hand quadrant of Fig. 10, we plot various (almost) ξ-independent definitions of the Higgs masses. The dotted lines correspond to the tree-level values: the mass of the mostly-singlet state decreases with increasing |A κ | and crosses the mass of the SM-like state at |A κ | ∼ 150 GeV, leading to substantial mixing. The (square root of the) diagonal entries of the effective mass matrixexactly gauge independent-are shown with dashed curves: radiative corrections displace the mass values. Remarkable points are |A κ | ∼ 150 GeV, where the tree-level mixing is 'un-mixed' (leading to an inversion of the hierarchy between M 2 eff 11 and M 2 eff 22 ), and |A κ | ∼ 250 GeV, where the actual crossing of singlet and doublet masses takes place at the one-loop order. As argued before, this definition of the masses at one-loop order is sufficient for |A κ | ≤ 100 GeV and |A κ | ≥ 260 GeV, but not in the intermediate regime where the mixing at radiative order competes with the diagonal mass-splitting. The masses in the mixing formalism are shown with solid curves and 'meet' at |A κ | ∼ 250 GeV. As explained above, a small gauge dependence remains present due to the offdiagonal self-energy. The impact of the latter on the mass-determination is of order ∼ 2-10 MeV when varying the gauge-fixing parameter between 0.1 and 100, hence completely negligible in view of higher-order corrections. The 'spikes' of the orange curve in the upper plot can now be understood: they correspond to the points of maximal tree-level mixing-where a large radiative mixing is needed to counteract the 'fake' tree-level mixing-and maximal mixing at the radiative order (with a very narrow mass-splitting).
The plot on the lower right-hand corner of Fig. 10 displays the dependence of the mass determination on the gauge-fixing parameter ξ for the point |A κ | 243 GeV with near-maximal mixing. The dashed lines are obtained with Eq. (12) and show negligible variation. The solid lines employ an iterative diagonalization procedure with variable external momentum-set to one eigenvalue of the mass matrix at each iteration in an attempt to solve Eq. (6): this procedure is highly dependent on ξ and the difference at low ξ with respect to the eigenvalues of Eq. (12) cannot be viewed as a genuine improvement in accuracy.