Gauge Dependences of Higher-Order Corrections to NMSSM Higgs Boson Masses and the Charged Higgs Decay $H^\pm \to W^\pm h_i$

In this paper we compute the electroweak corrections to the charged Higgs boson decay into a $W$ boson and a neutral Higgs boson in the CP-conserving NMSSM. We calculate the process in a general $R_\xi$ gauge and investigate the dependence of the loop-corrected decay width on the gauge parameter $\xi$. The gauge dependence arises from the mixing of different loop orders. Phenomenology requires the inclusion of mass and mixing corrections to the external Higgs bosons in order to match the experimentally measured mass values. As a result, we move away from a strict one-loop calculation and consequently mix orders in perturbation theory. Moreover, determination of the loop-corrected masses in an iterative procedure also results in the mixing of different loop orders. Gauge dependence then arises from the mismatch with tree-level Goldstone boson couplings that are applied in the loop calculation, and from the gauge dependence of the loop-corrected masses themselves. We find that the gauge dependence is significant.


Introduction
The discovery of the Higgs boson by the LHC experiments ATLAS and CMS [1,2] structurally completed the Standard Model (SM). Subsequent measurements revealed a very SM-like behavior of the Higgs boson. Due to open questions that cannot be answered within the SM, however, theories beyond the SM are considered, many of which contain extended Higgs sectors. So far, no direct signs of New Physics have been observed. This moves the Higgs sector itself into the focus of searches for indirect manifestations beyond the SM. Due to the very SM-like nature of the Higgs boson, sophisticated experimental techniques together with precise theoretical predictions are required for these investigations to be successful. In particular, higher-order (HO) corrections to the Higgs boson observables, their production cross sections, decay widths, and branching ratios have to be taken into account.
A clear manifestation of extended Higgs sectors would be the discovery of a charged Higgs boson that is not present in the SM. Charged Higgs bosons appear e.g. in the next-to-minimal supersymmetric extension of the SM (NMSSM) [3][4][5][6][7][8][9][10][11][12][13][14] which is the model that we consider in this work. More specifically, we work in the framework of the scale-invariant CP-conserving NMSSM. The main decay channels of the charged Higgs boson are those into fermionic final states, but also decays into a Higgs and gauge boson final state, or into electroweakinos can become numerically important depending on the specific parameter values. In this paper, we compute the electroweak corrections to the decay of the charged Higgs boson into a W boson and a light CP-even Higgs boson. We restrict ourselves to pure on-shell decays. The aim of this paper is not only to quantify the relative importance of the electroweak corrections, but in particular we also highlight problems with respect to gauge dependences that occur in the computation of the HO corrections. In a gauge theory gauge-breaking effects do not appear when the computation is restricted to a fixed order, here the one-loop level. This changes, however, when different loop orders are mixed, see e.g. also the discussions in [15][16][17][18][19][20][21][22][23][24]. We encounter such a mixing when we include loop corrections to the mass of the involved external Higgs boson. Since the tree-level upper bound of the SM-like Higgs boson is below the observed 125.09 GeV [25], loop corrections have to be included to shift its mass to the measured value. This introduces a mismatch between the loop-corrected mass of the external neutral Higgs boson and its corresponding tree-level mass, which is used in the propagators of the internal lines and in the tree-level Higgs-Goldstone boson couplings that occur in the computation of the one-loop amplitude 1 . While the latter can be cured by an appropriate change of the Higgs-Goldstone boson coupling as we will outline below (see also [18,19,[21][22][23] for a discussion), the former cannot be cured by introducing loop-corrected masses for the internal lines since it will break UV finiteness. Furthermore, we encounter additional gauge dependences due to the gauge dependence of the loop-corrected Higgs boson masses and their loop-corrected mixing matrix. The loop-corrected Higgs boson masses are defined through the complex poles of the propagator matrix which are evaluated by using an iterative method. While this method gives precise values of the complex poles, it mixes the contributions of different orders of perturbation theory and therefore introduces a dependence on the gauge parameter. The loop-corrected mixing matrix is used to resum the large corrections that stem from the mixing between different neutral Higgs bosons so that the loop corrections remain small and the higher-order prediction are reliable. In addition, the loop-corrected mixing matrix ensures the on-shell property of the external Higgs boson. This mixing matrix is, however, gauge dependent by definition. With the loop corrections to the light Higgs boson masses and the mixing matrix being substantial also the gauge dependence turns out to be significant and much more important than in the minimal supersymmetric extension of the SM (MSSM) as discussed in this work. The purpose of the present paper is to quantify this effect and to investigate different approximations with respect to their gauge dependences.
The outline of the paper is as follows. In Sec. 2 we introduce the Higgs sector of the NMSSM at tree level and at higher orders, and set our notation. In Sec. 3 we describe our computation of the electroweak one-loop corrections to the charged Higgs decay into a gauge plus Higgs boson final state. In particular we present the decay width at strict one-loop order. We follow up with a general discussion of gauge dependences encountered in the decay before presenting improved effective decay widths that include higher-order-corrected external Higgs states in different approximations. In the numerical results of Sec. 4 we analyze the gauge dependence of the loop-corrected masses themselves and subsequently the decay amplitudes and decay widths. We analyze the latter two with respect to their gauge dependences by including various approximations in the treatment of the loop-corrected external Higgs states. We also compare the results with the the size of the gauge dependences in the MSSM limit. We conclude with a small discussion of the relative size of the electroweak corrections as a function of the relevant NMSSM parameters. In Sec. 5 we summarize our results.

Higgs sector of the NMSSM
In this paper, we calculate within the NMSSM the one-loop corrections to the decays of the charged Higgs boson into a W ± boson and a neutral CP-even Higgs boson. To that end, we briefly introduce the Higgs sector of the NMSSM and set up the notation required both for the calculation of the charged Higgs decays as well as for the discussion of the higherorder-corrected neutral Higgs boson masses. Since we apply the same approximations and renormalization conditions as in our previous works on higher-order corrections to the NMSSM Higgs boson masses and trilinear self-couplings [20,[26][27][28][29][30], we remain here as brief as possible and refer, where appropriate, to the corresponding literature for more information. We work in the framework of an NMSSM wherein a gauge-singlet chiral superfieldŜ is added to the MSSM field content, and the superpotential couplings of this singlet are constrained by a Z 3 symmetry. In terms of the two Higgs doublet superfieldsĤ u andĤ d and the singlet superfield S the NMSSM superpotential is written as with the totally antisymmetric tensor ij (i, j = 1, 2) and 12 = 12 = 1, where i, j denote the indices of the fundamental SU (2) L representation. Working in the framework of the CPconserving NMSSM, the dimensionless parameters λ and κ are taken to be real. The MSSM superpotential W MSSM is expressed in terms of the quark and lepton superfields and their charge conjugates as denoted by the superscript c, i.e.Q,Û c ,D c ,L andÊ c , as For better readability, the color and generation indices have been suppressed, and µ (i.e. the supersymmetric Higgs mass parameter of the MSSM) is set to 0 due to the applied Z 3 symmetry. We neglect flavor mixing so that the Yukawa couplings y u , y d and y e , which in general are 3 × 3 matrices in flavor space, are diagonal.
The soft supersymmetry (SUSY) breaking NMSSM Lagrangian is given in terms of the scalar component fields H u , H d and S by where the summation over all three quark and lepton generations is implicit. Here, we denote byQ andL the complex scalar components of the corresponding left-handed quark and lepton superfields, so that e.g. for the first generation we haveQ = (ũ L ,d L ) T andL = (ν L ,ẽ L ) T . The M k (k = 1, 2, 3) represent the gaugino mass parameters of the bino, wino and gluino fieldsB, W j (j = 1, 2, 3) andG, the m 2 X are the squared soft SUSY-breaking mass parameters of the scalar fields X = S, H d , H u ,Q,ũ R ,d R ,L,ẽ R and A x (x = λ, κ, d, u, e) are the soft SUSY-breaking trilinear couplings. In general, the soft SUSY-breaking trilinear couplings and the gaugino mass parameters can be complex, but in the limit of CP conservation all parameters are taken to be real.

The Higgs Sector at Tree Level
The Higgs potential at tree level reads with g 1 and g 2 being the U (1) Y and SU (2) L gauge couplings, respectively. The two Higgs doublets and the singlet can be expanded around their vacuum expectation values (VEVs) v u , v d and v s as with the ratio between them being defined as such that v u and v d can be expressed in terms of v and tan β. The potential V H is minimized by the VEVs, which implies that the first derivatives of the potential with respect to the Higgs fields must be zero. This leads to the definition of the tadpole parameters t φ , which have to vanish. Since we are working in CP-conserving NMSSM, the tadpole parameters that we have at tree level are given by In the CP-conserving NMSSM, there is no mixing between CP-even and CP-odd Higgs fields so that the bilinear parts of the Higgs potential read with separate mass matrices M H ± , M h and M a for the charged, CP-even and CP-odd Higgs fields, respectively. The explicit expressions of these tree-level mass matrices can be found in [26]. The charged, neutral CP-even and CP-odd mass eigenstates are obtained from the interaction states through rotations with the unitary matrices R H ± , R h and R a as These rotation matrices diagonalize the mass matrices such that The obtained mass eigenstates are ordered by ascending mass so that we have three CP-even Higgs states h i (i = 1, 2, 3) with masses m h 1 ≤ m h 2 ≤ m h 3 , two CP-odd states a j (j = 1, 2) with masses m a 1 ≤ m a 2 and a charged Higgs pair H ± with mass m H ± as the physical Higgs bosons. The fields G 0 , G ± form the massless charged and neutral Goldstone modes 2 . In general, the analytic forms of the mass eigenvalues are rather intricate, but analytic expressions for expansions in special parameter regions can be found in [31]. On the other hand, the squared mass of the charged Higgs boson is at tree level given by the simple expression Note that analogous to the MSSM there is an upper bound on the squared SM-like Higgs boson mass at tree level. In the NMSSM, it is given by In order to comply with the measured value of m H = 125.09 GeV [25], loop corrections therefore have to be included in the computation of the Higgs boson mass.
The Higgs potential in Eq. 4 is parametrized by the parameter set For a physical interpretation, it is convenient to substitute some of these parameters with more intuitive ones, such as e.g. the masses of gauge bosons which are measurable quantities, or the tadpole parameters 3 . We can use Eqs. (6) and (7) to eliminate v u and v d in favour of v and tan β, and Eqs. (9)- (11) to replace the soft SUSY-breaking parameters m 2 H d , m 2 Hu and m 2 S in V H with t h d , t hu and t hs . Furthermore, A λ can be replaced by m 2 H ± using Eq. (15). Finally, g 1 , g 2 and v are substituted by the squared masses M 2 W and M 2 Z of the W ± and Z bosons and the electric charge e via , and e = g 1 g 2 In summary, our set of free parameters in the Higgs sector is given by Finally, the MSSM limit of the NMSSM Higgs sector can be obtained by setting and keeping all other parameters, including and A κ , fixed. In this limit the mixing between singlet and doublet Higgs fields vanishes.

The Loop-Corrected Higgs Sector
For the determination of the loop-corrected Higgs boson masses, the UV-divergent self-energies have to be calculated. The divergent integrals are regularized by the SUSY-conserving dimensional reduction scheme [32,33]. Evaluating the self-energies in D = 4 − 2 dimensions, the divergences can be parametrized by the regulator , leading to poles 1/ in the limit of → 0, i.e. in physical D = 4 space-time dimensions. Also in the one-loop corrections to the process H ± → W ± h i we encounter UV divergences. The UV divergences are cancelled by the renormalization of the Higgs fields and the parameters relevant for the calculation 4 . In order to do so, 3 Whether the tadpole parameters can be called physical quantities is debatable but certainly their introduction is motivated by physical interpretation. 4 Note that we do not renormalize the rotation matrices R H ± , R h , R a . For more details, cf. [26].
the bare parameters p 0 of the Lagrangian are replaced by the renormalized ones, p, and their corresponding counterterms, δp, Analogously, the bare fields φ 0 in the Lagrangian are expressed via the renormalized fields φ and the wave-function renormalization constants (WFRCs) Z φ as In accordance with our previous works on higher-order corrections to the NMSSM Higgs boson masses [26][27][28][29], we apply a mixed on-shell (OS) and DR renormalization scheme to fix the parameter and WFRCs. The free parameters of Eq. (19) are defined to be either OS or DR parameters as follows The renormalization scheme for the parameters is chosen such that the quantities which can be related to physical observables are defined on-shell, whereas the rest of the parameters are defined as DR parameters 5 . In addition, we introduce the WFRCs for the doublet and singlet fields before rotation into the mass eigenstates as We apply DR conditions for the WFRCs of the Higgs fields. We introduce a WFRC for the W boson field, needed in the computation of the loop corrections to the charged Higgs boson decay, as The WFRC δZ W is defined through the OS condition where Σ T W W denotes the transverse part of the W boson self-energy. The Higgs boson masses and hence the mixing matrices receive large radiative corrections. Therefore it is necessary to include these corrections at the highest order possible to improve the theoretical predictions. Recently, we completed the two-loop order O(α 2 t ) corrections to the neutral Higgs boson masses in the CP-violating NMSSM [29], thus improving our previous results, which were available to two-loop order O(α t α s ) [28]. The Higgs boson masses corrected up to two-loop order O(α t α s + α 2 t ) require also the renormalization of the top/stop sector at one-loop order. The computation of the two-loop corrections together with the renormalization of the parameters in the above defined mixed OS-DR scheme has been described in great detail in [28,29] 6 , hence we do not repeat it here and instead refer to these references for details. The CP-conserving limit of these results given in the CP-violating NMSSM is straightforward, further information can also be found in [26] where the one-loop calculation is presented for the real NMSSM. We review here, however, important points and highlight differences for the purpose of discussing the parameter dependence. In the following, we focus on the CP-even Higgs bosons. Their loop-corrected masses are defined as the real parts of the poles of the propagator matrix. These complex poles are the zeros of the determinant of the renormalized two-point correlation functionΓ(p 2 ), where p 2 denotes the external squared four-momentum. The renormalized two-point correlation function is expressed as 7 where the renormalized self-energyΣ h i h j (p 2 , ξ) of the transition h i → h j (i, j = 1, 2, 3) is given byΣ Here,Σ 1l (p 2 , ξ) denotes the full one-loop self-energy with full momentum-dependent contributions computed in general R ξ gauge, where ξ stands for the gauge parameters ξ W , ξ Z 8 . The last two terms are the two-loop corrections of order O(α t α s ) [28] and O(α 2 t ) [29], respectively, which are evaluated in the approximation of vanishing external momentum. These contributions do not introduce additional gauge-dependent terms in the renormalized self-energies as they are evaluated in the gaugeless limit. We want to point out that the full one-loop renormalized self-energiesΣ 1l (p 2 , ξ) in general R ξ gauge are newly computed by us and implemented in NMSSMCALC. We computed them both in the standard tadpole scheme and in the Fleischer-Jegerlehner scheme 9 and the results are identical. We apply the iterative procedure described and applied in [26] to extract the zeros of the determinant. In the first iterative step for the calculation of the n th CP-even Higgs boson mass, the squared external momentum in the renormalized self-energies is chosen to be at its squared tree-level mass value, i.e. p 2 = m 2 hn . The mass matrixM 2 (p 2 , ξ) is then diagonalized, yielding the n th eigenvalue as a first approximation to the squared mass of the n th CP-even Higgs boson. In the next step of the iteration, this value is taken as the new OS value for p 2 , and the mass matrix is again diagonalized. This procedure is repeated until the n th eigenvalue changes by less than 10 −9 GeV 2 between two consecutive steps of the iteration. The resulting complex pole is denoted byM 2 Hn and the loop-corrected Higgs mass by M Hn = Re(M 2 Hn ). The loop-corrected CP-even Higgs masses are sorted in ascending order M H 1 < M H 2 < M H 3 . Note that we denote the loop-corrected masses and Higgs mass eigenstates by capital letters M and H i , respectively, whereas the corresponding tree-level values and eigenstates are denoted by lowercase letters, i.e. m and h i . The iterative procedure automatically mixes different orders of perturbation theory and therefore introduces intricate gauge dependences 10 . This will be investigated in more detail in the numerical section.
Besides the computation of the loop-corrected masses, the code NMSSMCALC allows us to compute the loop-corrected mixing matrices in several approximations which are discussed also in [57]. First, we introduce the matrix R 0 for the rotation of the mass matrix in the approximation of vanishing external momentum, The corresponding loop-corrected mass eigenvalues are denoted by an index 0, hence M 2 0,H i (i = 1, 2, 3). In this approximation the mixing matrix R 0 is unitary, but does not capture the proper OS properties of the external loop-corrected states as momentum-dependent effects are neglected.
A different approach leads to the rotation matrix R mtree . Here we diagonalize the mass matrix with the elements evaluated at fixed momentum squared which is given by the arithmetic average of the squared masses, We hence have and the corresponding mass values denoted by M 2 mtree,H i are obtained through rotation with the matrix R mtree as By this approach we approximately take into account the momentum dependence of the selfenergies (see [22] for a discussion).
The correct OS properties of the loop-corrected mass eigenstates are obtained by following the procedure described in [15], i.e. by multiplying the tree-level matrix R h with the finite wavefunction normalization factors given by the Z matrix [15] for external OS Higgs bosons at higher orders, 10 The gauge dependence of the electroweakino masses determined by the iterative method has been discussed in [56].
with the indices i, j = H 1 , H 2 , H 3 . The prime denotes the derivative with respect to p 2 . The quantity is evaluated at the complex polesM 2 i . In contrast to the rotation matrices R 0 and R mtree , which are unitary matrices, the Z matrix is not as it contains resummed higher-order contributions. If we want to compute the decay width at exact one-loop level, we therefore have to expand the Z matrix and take only the one-loop terms where (Σ 1l (p 2 )) = ∂ p 2Σ 1l (p 2 ). Note that theΣ 1l are evaluated at the tree-level mass values m 2 h i , since using loop-corrected masses would introduce higher-order effects.
3 Electroweak One-Loop Corrections to H ± → W ± h i

Decay Width at Leading Order
The decay of the charged Higgs boson H ± into a W ± boson and a CP-even neutral Higgs boson h i (i = 1, 2, 3) depends on the coupling The leading-order (LO) decay width can be written as with λ(x, y, z) = x 2 +y 2 +z 2 −2xy−2xz−2yz denoting the usual Källén function and M tree the reduced matrix element, which for the tree-level decay is given by We remind the reader that by m h i we denote the tree-level mass of the final state Higgs boson.

Decay Width at Strict One-Loop Order
The next-to-leading order (NLO) decay width Γ NLO is given by the sum of the LO width Γ LO , the virtual corrections Γ virt and the real corrections Γ real as The virtual corrections contain the counterterm contributions that cancel the UV divergences. The IR divergences in the real corrections cancel those in the virtual corrections. Γ virt is given by Note that in Eq.
we set the external Higgs boson on its tree-level mass shell. We do the same in the LO decay width and the real corrections. In this way, we ensure that the NLO decay width Γ NLO remains at strict one-loop order by avoiding admixtures of higher orders through loop corrections to the mass. M virt consists of the one-particle-irreducible (1PI) diagrams depicted in Fig. 1. They show the external leg corrections M ext,h and M ext,H ± to the neutral Higgs boson and to the charged Higgs boson, respectively, and the genuine vertex corrections M vert . They already include the counterterms and are hence finite. The corrections to the W boson leg vanish due to the OS renormalization of the W boson. We hence have The amplitudes for the external leg contributions to the neutral and charged Higgs bosons, M ext,h and M ext,H ± , respectively, can be factorized into the tree-level amplitude and the propagator corrections to the external legs. For M ext,h we obtain Note that here we apply δZ 1l ij at strict one-loop order, defined as with the Z matrix at strict one-loop order, Z 1l , given in Eq. (40). The charged Higgs WFRC is determined in the DR scheme so that there are finite contributions to the LSZ factor Ẑ H ± at one-loop order, where the prime denotes the derivative with respect to the squared four-momentum.
The mixing between H ± and G ± can be related to the mixing between H ± and W ± by using the Slavnov-Taylor identity for the renormalized self-energies, whereΣ H + W − (m 2 H ± ) denotes the renormalized truncated self-energy of the transition H + → W + . The correction to the H ± propagator thereby results in where the coupling is given by The genuine vertex corrections M vert are given by the diagrams depicted in Fig. 2 plus the corresponding counterterm contributions that are not shown here. The vertex corrections comprise the 1PI diagrams given by the triangle diagrams with scalars, fermions and gauge bosons in the loops, and the diagrams involving four-particle vertices. The counterterm amplitude is given by in terms of the WFRCs δZ H d , δZ Hu and δZ W and the counterterm δg 2 for the SU (2) L gauge coupling constant g 2 , which in terms of the counterterms of our input parameters, cf. Eq. (24), reads The vertex diagrams also contain IR divergences. These arise from the exchange of a soft virtual photon between the external legs (cf. diagrams 21 and 28 of Fig. 2). Also δM 2 W , δZ W andẐ H + H − are IR-divergent. These soft singularities in the virtual corrections are canceled by the IR-divergent contributions from real photon radiation [58,59] in the process H ± → W ± h i γ. The process is independently gauge-invariant and most easily calculated in unitary gauge. The diagrams that contribute in unitary gauge to the process are shown in Fig. 3. They consist of the proper bremsstrahlung contributions, where a photon is radiated from the charged initial and final state particles, and the diagram involving a four-particle vertex with a photon. The decay width for the real emission is given by where E W , E h i and E γ denote the energies of the corresponding particles and the plus sign corresponds to the outgoing momenta p W and p h i while the minus sign belongs to the incoming momentum p H ± . No upper (lower) index j k (i l ) means that the corresponding ±2p γ p j k (±2p γ p i l ) in the numerator (denominator) of the fraction is replaced by 1. The total NLO width Γ NLO is then both UV-and IR-finite. Furthermore, since Γ NLO has been calculated at strict one-loop level with p 2 h = m 2 h i , it is also independent of the gauge parameter as we explicitly checked. We finish with the remark that for the computation of the loop-corrected decay width we used a FeynArts-3.10 [61] model file for the NMSSM generated by SARAH-4.12.3 [62][63][64][65]. The various pieces of the one-loop corrected decay width were obtained with the help of FormCalc-9.6 [66], FeynCalc-8.2.0 [67,68] and LoopTools-2.14 [66]. Both for the computation of the loopcorrected Higgs boson masses and the decay widths two independent calculations have been performed which are in full agreement.

The Issue of Gauge Dependence
Phenomenology requires that the loop corrections to the masses of the external Higgs bosons should be taken into account. This is particularly important when the external Higgs boson is the SM-like scalar, as its upper mass bound at tree level is well below the measured value of 125.09 GeV. Therefore, in the decay of the charged Higgs boson H ± into a final state with neutral Higgs bosons, we should consider the loop-corrected Higgs states H i (i = 1, 2, 3) with the corresponding loop-corrected masses M H i . For the decay H ± → W ± h i , this means that we should set the external momentum to p 2 = M 2 H i . However, this introduces contributions beyond the one-loop order in the one-loop decay width, which has two implications. First, it invalidates the tree-level relation between the couplings of the charged and neutral Higgs bosons with a charged Goldstone boson or a W ± boson, i.e. between the couplings g h i H − G + and g h i H − W + . This relation needs to be satisfied, however, in order to cancel the IR divergences occuring in the decay H ± → W ± h i at one-loop order. Additionally the relation between g h i a j G 0 and g h i a j Z is spoiled, but this does not influence IR divergence. Second, the introduction of loop-corrected neutral Higgs masses leads to mixing of different orders of perturbation theory, and the gauge independence of the matrix element is no longer guaranteed, and it is indeed violated as will be shown in the following. The problem of gauge dependences arising from loop-corrected Higgs masses, and their restoration via the inclusion of partial two-loop terms in the case of neutral Higgs decays in the NMSSM with complex parameters was recently discussed in [23]. In this work, the gauge dependence arose from the mixing of the neutral Higgs bosons with the Z boson. This does not apply to our work as we are working in the CP-conserving NMSSM and we are considering only decays into CP-even neutral Higgs bosons. The gauge dependence in our case originates from other sources and cannot be remedied easily, if at all.
In this paper, we want to investigate the impact of this gauge dependence on the treatment of the wave-function normalization factors and on the parameters of the model. We also investigate the issue of the gauge dependence of the loop-corrected Higgs masses themselves. As long as there is no recipe on how to achieve gauge-independent results 11 , the value of the gauge parameter applied in the computation of a specific Higgs observable needs to be specified in order to consistently relate measured observables with the parameters of the underlying model.

Decay Width at Improved One-Loop Order
In the following, we look more closely into the relation between the gauge dependence of the loop-corrected decay H ± → W ± h i and the treatment of the external Higgs boson, in particular the treatment of the Z matrix. While curing the IR divergence beyond strict one-loop order is fairly straightforward, the intricacies of gauge dependence in our calculation with respect to setting p 2 h = M 2 H i are much more involved. In order to study this in more detail, we proceed in two steps: 1. In the first instance, we modify our result obtained in Sec. 3.2 by changing p 2 h from the treelevel value m 2 h i to the loop-corrected one M 2 H i , and ensure that all IR divergences cancel by enforcing the correct relations between the gauge couplings g h i H − G + and g h i H − W + beyond tree level. However, we retain the use of the one-loop diagrammatic expansion of the Z matrix as applied in Eq. (47). It is clear that in order to get correct OS properties for the external neutral Higgs boson, we need to make use of the resummed Z matrix defined in Eq. (37). However, it is instructive to demonstrate the breaking of the gauge symmetry that occurs simply by using loop-corrected masses M 2 H i , before we discuss the full result obtained by using the resummed Z matrix.
2. For the next step, we set p 2 h = M 2 H i and apply the resummed Z matrix in our calculation, treating it as a part of the LO amplitude. This means that we no longer need to explicitly include external leg corrections M ext,h . As a result of this modification we will be required to include Z factors also for the real corrections, as well as to modify the gauge coupling relation between g h i H − G + and g h i H − W + , in order to obtain an IR-finite result.
Step 1: The first modification of our strict one-loop decay width consists of calculating Γ virt the real decay width Γ real is separately gauge independent even when computed at p 2 h = M 2 H i . The virtual amplitude M virt , however, which is gauge independent when computed at strict one-loop order, acquires a dependence on the gauge parameters due to higher-order mass effects. Moreover, the cancelation of the IR divergences in Γ NLO H ± →W ± h i only takes place when the relation is satisfied. Due to the gauge structure of the Lagrangian, the relation holds, with m 2 h i being the tree-level CP-even Higgs mass calculated from Eq. (14), cf. [69]. We can enforce the relation Eq. (57) beyond tree level by modifying the coupling g h i H − G + where necessary 13 such that it is expressed in terms of the loop-corrected masses, This is equivalent to an effective potential approach [69]. This modification ensures an IRfinite NLO width. The decay width obtained with these modifications will be referred to as "off-shell" and denoted as Γ NLO,off-shell H ± →W ± h i . This nomenclature points towards the fact that the external loop-corrected neutral Higgs boson does not have the correct OS properties yet. We emphasize that while Γ NLO,off-shell H ± →W ± h i is UV-and IR-finite, the modification of the coupling constants in Eq. (59) does not restore gauge independence. The global modification of the Goldstone couplings g G + G − h i , g G + W − h i , g G 0 G 0 h i and g G 0 a j h i is not possible while keeping the result UVfinite, such that gauge independence cannot be restored by a modification of these couplings. Additionally, we have to deal with the gauge dependence of the loop-corrected Higgs masses themselves.
Step 2: The corrections from M ext,h to the "off-shell" Γ NLO H ± →W ± h i can be large. In order to obtain numerically stable NLO corrections 14 and to ensure that the external neutral Higgs bosons have proper OS properties, we need to use the full Z matrix defined in Eq. (37) which resums higher-order contributions to the external leg corrections, and treat it as part of the LO amplitude. The second modification to our strict one-loop computation therefore consists of not only using p 2 h = M 2 H i , but also employing the resummed Z factors. The LO and the virtual amplitude Eq. (46) are now computed as and Including the resummed Z factors in the LO amplitude means that M ext,h does not have to be calculated anymore, and that the virtual NLO amplitude contains contributions that are formally of two-loop order and higher. We refer to these amplitudes as "improved" amplitudes and to their corresponding widths as "improved" widths. They are given by The Z matrix and the loop-corrected masses are obtained from the program NMSSMCALC [26][27][28][29][70][71][72][73] with the new implementation of the full one-loop renormalized self-energies in general R ξ gauge as discussed in subsection 2.2.
The inclusion of resummed higher-order corrections to the external neutral Higgs boson via the full Z factor also needs to be accounted for in the real corrections, so that the IR divergences cancel properly. This means that we have Finally, since we set p 2 h = M 2 H i , we need to use the same modified gauge couplings that we introduced in Eq. (59) in order to cure the breaking of IR finiteness caused by using loopcorrected masses. We refer to this width as "improved" real corrections. The complete width obtained after these modifications will henceforth be referred to as the "improved" NLO width and is given by

Numerical Results
In this section, we will investigate in detail the gauge dependence of our results. We start by studying the gauge dependence of the loop-corrected neutral Higgs boson masses and then investigate the gauge dependence of the virtual corrections to the decay H ± → W ± h i before studying the gauge dependence of the complete NLO width, by applying various treatments of the external Higgs bosons. We do this for parameter points obtained from a scan in the NMSSM parameter space as described in the following.

The NMSSM Parameter Scan
In order to find scenarios that are compatible with the recent experimental constraints for the purpose of our numerical analysis, we perform a scan in the NMSSM parameter space. We apply the same procedure as in [74][75][76], where also further details can be found. The parameters tan β, λ and κ are varied in the ranges  Table 1: Scan ranges for the NMSSM scan. All parameters are varied independently between the given minimum and maximum values, denoted by "min" and "max", respectively.
so that we obey the rough perturbativity constraint The scan ranges of further parameters are listed in Table 1 . We set and the mass parameters of the first and second generation sfermions are chosen to be The soft SUSY-breaking trilinear couplings of the first two generations are set equal to the corresponding values of the third generation. We follow the SUSY Les Houches Accord (SLHA) format [77,78], which means that the soft SUSY-breaking masses and trilinear couplings are understood as DR parameters at the scale The SM input parameters have been chosen to be [79,80] We calculate the spectrum of the Higgs particles including corrections up to two-loop order O(α t α s + α 2 t ) [29] with the recently published NMSSMCALC version 3.00. For the scan, M H ± has been used as input parameter, and not A λ which is also provided as an option by NMSSMCALC. We choose OS renormalization for the top/stop sector (see [28,29] for details). One of the neutral CP-even Higgs bosons is identified with the SM-like Higgs boson and its mass is required to lie in the range 122 GeV ≤ m h ≤ 128 GeV . We demand the total χ 2 computed by HiggsSignals with our given effective coupling factors to be compatible with the total χ 2 of the SM within 2σ. The input required for HiggsSignals is calculated with NMSSMCALC.
Furthermore, the most relevant LHC exclusion bounds on the SUSY masses are taken into account. These constrain the gluino mass and the lightest squark mass of the second generation to lie above 1.8 TeV, see [85]. The stop and sbottom masses in general have to be above 800 GeV [85,86], and the slepton masses above 400 GeV [85].
For the numerical analysis we have chosen two sample parameter points among all parameter sets obtained by our scan. For the first scenario, denoted by P1, the lightest CP-even Higgs boson is singlet-like and the second-lightest CP-even Higgs boson is the SM-like Higgs boson. In this scenario the mixing between the singlet-and the SM-like states is quite significant. The second point is denoted by P2 and features an SM-like Higgs boson that is the lightest CP-even Higgs boson with the mixing between the singlet and the SM-like state being non-significant.

Gauge Dependence of the Neutral Higgs Boson Masses
As it has been discussed in Sec.  The Higgs boson masses and their main compositions in terms of singlet/doublet and scalar/pseudoscalar components at tree level, one-loop level as well as two-loop O(α t α s ) level and two-loop O(α t α s +α 2 t ) level are given in Table 2 for OS and in Table 3 for DR renormalization in the top/stop sector. They have been computed with NMSSMCALC in the 't Hooft-Feynman gauge (ξ = 1). In the Table, lowercase h i refers to the tree-level and capital H i to the loop-corrected mass eigenstates. In our chosen parameter point, the h s -like and the h u -like Higgs boson masses are light and receive significant higher-order corrections. We call the Higgs boson singlet-like in case its dominant contribution to the mass eigenstate stems from the singlet admixture. The second-lightest Higgs boson is dominantly h u -like and has a mass of 125 GeV at O(α t α s + α 2 t ) for OS renormalization in the top/stop sector. It reproduces the LHC production rates which proceed dominantly through gluon fusion for small values of tan β. Since the LO process is dominated by top-quark loops the Higgs coupling to the tops must be substantial, as is the case for a Higgs boson with large h u admixture.
We first vary the gauge parameter ξ of the general R ξ gauge, which we set throughout the section ξ W = ξ Z ≡ ξ, while all other parameters are kept fixed. The masses of the h s -and the h u -like Higgs bosons depend significantly on ξ. All remaining Higgs bosons have masses larger than 600 GeV and show a very small gauge dependence. This is due to the fact that only the light Higgs bosons receive significant loop corrections. In Fig. 4, we show these dependences for the mass M hs of the CP-even singlet-like Higgs boson in the upper left plot and for the mass M hu of the SM-like Higgs boson in the upper right plot including one-loop (black lines), O(α t α s ) two-loop (blue lines) and O(α t α s + α 2 t ) two-loop (red lines) corrections. These corrections are obtained for the OS (full lines) and DR (dashed lines) renormalization schemes of the top/stop sector. The two lower plots display the relative difference between the masses in general R ξ gauge and in the 't Hooft-Feynman gauge ξ = 1,    fixed loop order, calculated in general R ξ gauge (M x (ξ)) and in the 't Hooft-Feynman gauge (M x (ξ = 1)). We remind the reader that only the renormalized one-loop Higgs self-energies are calculated in general R ξ gauge and therefore depend on ξ, while the renormalized two-loop Higgs self-energies at order O(α t α s ) and O(α 2 t ) do not depend on ξ as they are computed in the gaugeless limit. Note that the tree-level masses for the h s -and h u -like Higgs bosons are 9.8 GeV and 91.38 GeV, respectively. The loop corrections to their masses are positive. From the plots, we can infer that the loop-corrected masses decrease with increasing ξ, which we chose to lie between 0 and 100. We can, however, increase ξ to a larger value and find that for ξ = 1089 the  Next, we fix the value of ξ to 10 and vary λ and κ at the same time to very small values (10 −5 ) while we keep the ratio λ/κ constant. In this way we smoothly approach the MSSM limit, where the singlet and doublet Higgs bosons do not mix. In Fig. 5 we show the thus obtained loop-corrected masses of the h s -like (left) and h u -like (right) Higgs boson as function of λ. The line and color codes are the same as in Fig. 4. The lower plots show ∆ M ξ=10 , i.e. the deviation of the h s -and h u -like masses, respectively, calculated for ξ = 10 from the value obtained in the 't Hooft-Feynman gauge.
As expected, when λ and κ are close to zero, the h s -like Higgs boson decouples and the loop corrections to the h s -like Higgs boson mass vanish in this limit. Therefore all lines in the left plots converge at the left endpoint where λ ≈ 0. As we can see from Fig. 5 (right), the ξ dependence of the SM-like loop-corrected Higgs boson masses does not vanish in the MSSM limit. For our chosen parameter point the relative deviation of the masses for ξ = 10 and ξ = 1 even slightly increases. The kink around λ = 0.29 appears at the threshold where the loopcorrected h s mass is twice the tree-level h s mass.
In order to investigate the influence of the NMSSM-specific contributions to the mass corrections at O(α t α s + α 2 t ) and their ξ dependence we simultaneously vary the couplings λ and κ and show in Fig. 6 (upper)   depict the corresponding relative ξ dependence. For h s the two-loop corrected mass value and the ξ dependence decrease with smaller values of the NMSSM-like couplings, as the singlet-like Higgs boson decouples from the spectrum. The h u -like mass shows the expected behavior and decreases with smaller singlet admixtures. 16 The relative ξ dependence increases, however, for the chosen parameter point. The increasing contribution to the mass corrections for larger λ, κ values from the h u − h s admixture mixes with the doublet-like mass corrections and diminishes their gauge dependence.
The gauge dependence strongly depends on the chosen approach to determine the loopcorrected masses as we will show next. Figure 7 displays the mass corrections (upper plots) at O(α t α s +α 2 t ) for OS renormalization in the top/stop sector and their relative ξ dependence (lower plots) determined through the iterative method to extract the zeros of the determinant [26] (red lines) as well as when we apply the zero momentum approximation ('R 0 -method' called in the following, blue lines), cf. Eq. (33), and when the mass matrix is diagonalized at the arithmetic squared mass average ('R mtree -method', black lines), cf. Eqs. (35) and (36). The gauge dependence becomes very small for the latter in contrast to the two former methods. This is because of the fact that the dependence of the renormalized self-energyΣ h i h j evaluated at the arithmetic squared mass average on the gauge parameter is small for i being different from j and vanishes completely for i being identical to j. Their behavior as a function of ξ depends on the difference between the tree-level mass and the squared momentum at which the mass matrix is diagonalized resulting in ∆ M ξ values up to about 1% (4%) for the R 0 -method (iterative method) for the h s -like mass and 10% (7%) for the h u -like mass when ξ is varied up to values of 100.   Fig. 8, we show these curves for the strict one-loop calculation as described in Sec. 3.2. This means that the external leg corrections to h u are accounted for diagrammatically using Eq. (40), as opposed to using the resummed Z matrix. Moreover, we use the tree-level mass m hu for the external momentum such that p 2 h = m 2 hu . Here and in the following plots we use the mixing matrix that diagonalizes the tree-level mass matrices in the computation of the couplings as otherwise the result will not be UV-finite. In this strict one-loop computation, the virtual corrections are gauge independent, as can be checked explicitly numerically, cf. the solid black curve of Fig. 8 (left): while each individual component of the virtual correction is gauge dependent, their sum, resulting in M virt , is gauge independent. Actually, the relative gauge dependences of the external leg contributions to the charged and the neutral Higgs mass, M ext,H ± and M ext,h , and the ones of the vertex corrections, M vert , come with opposite sign (not visible from the plot as we show the absolute values). The reason for the kinks in the red (dotted) and green (dot-dashed) curves is the following. The masses of the Goldstone bosons depend on the gauge parameter ξ, and these kinks occur when a production threshold for the Goldstone boson is reached, i.e. at the points In Fig. 8 (right), we investigate how this gauge independence of the strict one-loop computation Step 1' that we denoted 'off-shell'. The means, we use loop-corrected masses for the external leg corrections to the neutral Higgs boson h u , i.e. we set p 2 h = M 2 hu . Note, however, that the Z matrix is evaluated at pure one-loop order, as defined in Eq. (40). The masses are calculated at O(α t α s + α 2 t ) for OS renormalization in the top/stop sector by NMSSMCALC in general R ξ gauge as described in subsection 2.2. Going from the strict one-loop calculation to the 'off-shell' one, we see that the dependence of the individual components of the virtual corrections on the gauge parameter ξ changes, such that their sum M virt (solid, black curve) is no longer gauge independent. The overall gauge dependence does not cancel any more, as we move away from the strict fixed-order calculation and include partial higher-order effects coming from the loop-corrected Higgs mass M 2 hu , such that with increasing ξ values, the NLO amplitude becomes arbitrarily large. The relative change of the total virtual amplitude is of up to O(10%) for ξ values up to 100.
In Fig. 9 (left), the curves corresponding to the right plot of Fig. 8 are plotted in the MSSM limit of the chosen parameter point. This limit is taken by setting λ, κ → 0 and keeping the ratio λ/κ constant (actually λ, κ = O(10 −8 )). From Fig. 9 we see that in the MSSM limit the resulting gauge dependence of M virt has a numerically small effect. It varies up to 0.7% 17 , although we are using gauge-dependent loop-corrected masses for the external momentum p 2 h = M 2 hu . This is illustrated once again in the right plot of the figure, where we show M virt and its relative gauge dependence alone for the chosen parameter point (red) and after varying λ and κ to half their values (blue) and to λ/100 and κ/100 (black). The relative gauge dependence decreases successively from 10% to 1% at ξ = 100. While the gauge dependence of the masses increases in the MSSM limit for our chosen parameter point the opposite is hence the case for the loop corrections to the decay. This again shows that the singlet admixtures play an important role for the gauge dependence of the parameters and observables and do not follow a simple law.

The Complete Loop-Corrected Decay Width
In the following, we study for the parameter point P1 the gauge dependence of the complete loop-corrected partial decay width of the decay H ± → W ± h u . In the upper plots of Fig. 10, the black (solid) curve shows, as a function of ξ, the improved LO decay width for H ± → W ± h u , i.e. we apply Eq. (62) as denoted with 'Step 2' in Sec. 3.4. This means that we set the external momentum to the loop-corrected Higgs boson mass M hu , p 2 h = M 2 hu , which is calculated at O(α t α s + α 2 t ) with OS renormalization in the top/stop sector. Additionally, we include in the external-leg corrections to the neutral Higgs boson the resummed Z matrix defined in Eq. (37) in order to ensure the correct OS properties. The blue (dot-dashed) curve displays the corresponding improved NLO width, given by Eq. (65). The left plots are for the NMSSM parameter point P1, whereas the right plots are for the MSSM limit of the same benchmark point. The NMSSM widths show a stronger dependence on ξ than the ones in the MSSM limit (note that the scales of the two plots are different). In the lower plots we show the relative gauge dependence of the LO and NLO widths, respectively, as a function of ξ, as defined by  Figure 11: Analogous to Fig. 10, but using the R 0 -method instead of the iterative approach to extract the loop-corrected masses and mixing matrix.
Here Γ ξ denotes the decay width calculated in general R ξ gauge at fixed loop order, i.e. Γ LO, impr or Γ NLO, impr , and Γ ξ=1 the width calculated in the 't Hooft-Feynman gauge. In the NMSSM the relative gauge dependence is larger for the NLO width than for the LO one while in the MSSM, where the singlet admixture vanishes, the opposite is the case. For the complete NLO width of scenario P1 the relative corrections ∆ Γ ξ can become as large as 200% for ξ = 100. Such a strong gauge dependence is clearly unacceptable for making meaningful predictions for the decay widths. For phenomenological investigations, on the other hand, the interesting quantity is the branching ratio. In order to make meaningful predictions, this requires the inclusion of the electroweak corrections to all other charged Higgs boson decays, so that the total width entering the branching ratio is computed at higher order in the electroweak corrections. This is beyond the scope of the present paper and left for future work. Even if one argues not to introduce new mass scales in the process and to remain below ξ values of 100 the ξ dependence is large, in particular it is far beyond the relative size of the loop corrections which is about 11% for ξ = 1. In the MSSM limit, depicted in the right plot of Fig. 10, the relative gauge dependence is smaller with values of up to about 20% for ξ = 100.
In Fig. 11, we show the corresponding curves analogous to Fig. 10, however now using the R 0 method to extract the loop-corrected mass values and mixing matrix 18 . The LO and NLO widths are then calculated by applying Eqs. (60) to (64), but with the Z matrix replaced by the R 0 matrix, defined in Eq. (33). We denote the corresponding widths with the superscript R 0 as Γ R 0 . Figure 12 shows the corresponding results if the masses are extracted at the arithmetic squared mass average such that the Z matrix is replaced by the R mtree matrix, defined in Eqs. (35,36). The corresponding widths are denoted by the superscript R mtree .
The comparison of Figs. 11 and 12 with Fig. 10 shows that for this parameter point the 18 As remarked above, in the couplings we always use the tree-level mixing matrix elements, however.   Figure 12: Same as Fig. 10, but using the R mtree -method instead of the iterative approach to extract the loop-corrected masses and mixing matrix.
gauge dependence is smallest in the R 0 approximation. The relative change of the complete NLO width with ξ compared to its value for ξ = 1, i.e. ∆ Γ ξ , is about -18% for ξ = 100, while in the R mtree approximation it is about +28%, which is still smaller than if the Z matrix is applied. The corresponding values in the MSSM limit are 6.5% (R 0 ) and -6% (R mtree ). The method of extracting the mixing matrix elements has a strong influence on the ξ dependence of the NLO width and also on the sign of this ξ dependence. For the parameter point P1 the h u -like Higgs boson has a strong singlet admixture. From previous analyses [49,51], we know already that the mixing matrix elements are then very sensitive to changes in the approximation of the loop calculation. Since the mixing matrix elements enter the Higgs couplings, the computed observables, in this case the decay width, become very sensitive to the applied approximation. This is confirmed by our results on the ξ dependence but also by the values of the widths themselves for the various approximations. 19 Overall, we found that the gauge dependence of the loop-corrected mass of the external neutral Higgs boson has a much smaller influence on the gauge dependence of the NLO width than the matrix that is used to set the external Higgs boson OS. The strength of this effect sensitively depends on the chosen parameter set, as can be inferred from Fig. 13. The figure displays the partial decay widths of the decay H ± → W ± h u (left plot) and H ± → W ± h s (right plot) both at LO and NLO as a function of tan β. All other parameter values are fixed to those of scenario P1. Shown are the results for the pure LO and strict one-loop width (black lines) and the ones when we calculate the improved LO and NLO widths applying the Z matrix (red) and the R 0 matrix (green). The lower plots show the corresponding K factors, defined as the  Figure 13: Parameter point P1: partial decay width of the decay H ± → W ± h u (left) and H ± → W ± h s (right) at LO (full) and NLO (dashed) at strict LO and one-loop order (black), in the improved approach applying the Z matrix (red) and the R 0 matrix (green), as a function of tan β. Lower panels: corresponding K factors.
ratio of the NLO width and its corresponding LO width As can be inferred from the left plot, the value of the decay width Γ(H ± → W ± h u ) strongly depends on the applied approximation for our chosen parameter point, i.e. for tan β = 3.11, while the K factor is approximately the same for the improved widths, with a value around 1.1. The K factor for the pure one-loop result largely differs from the improved ones as it does not take into account any resummation of higher orders. For values of tan β between about 3.26 and 3.52 the improved NLO results are rather close, but differ otherwise. In the singlet-like final state shown in Fig. 13 (right), i.e. Γ(H ± → W ± h s ), the K factors for the improved widths differ by less than 5% over the whole scanned range.

Analysis for Scenario P2
In order to further investigate the impact of the gauge dependence, we analyze the gauge dependence of the Higgs boson masses and the charged Higgs decay widths for a second parameter point, P2, defined in Eqs. (75) and (76). We summarize the Higgs boson masses obtained in the OS renormalization scheme of the top/stop sector in Table 4 and in the DR scheme in Table 5, at tree level, at one-loop level and at two-loop level including the O(α t α s ) and the O(α t α s + α 2 t ) corrections, respectively. We have deliberately chosen this scenario in which the h u -like Higgs boson is the lightest state with mass around 125 GeV at O(α t α s + α 2 t ) in the OS renormalization scheme of the top/stop sector while the CP-even singlet-like Higgs boson is the second-lightest state with mass around 433 GeV. In this scenario we analyze only the mass of the h u -like Higgs boson since only its mass is affected substantially by the change of the gauge parameter ξ. We     Fig. 14 the h u -like Higgs boson mass as function of ξ. The left plot shows its two-loop mass at O(α t α s + α 2 t ) for OS (full) and DR (dashed) renormalization in the top/stop sector for three different singlet admixtures. This means we start with the λ and κ values of our original scenario P2 (red lines) and compare with the results when we take half (blue lines) and 1/100 their values (black lines), where the latter means that we are close to the MSSM limit. As can be inferred from the plot, the gauge dependence ∆ M ξ shown in the lower plot is only mildly dependent on the renormalization scheme and on the singlet admixture, and amounts up to values of about 7 to 8 % at ξ = 100.
In the right plot of Fig. 14 we show the ξ dependence of M hu when we apply different approximations to determine the loop-corrected Higgs mass eigenstates with OS renormalization of the top/stop sector, namely through the iterative method (red line), by applying the rotation matrix R 0 to the mass matrix in the zero momentum approximation (blue), or finally by applying R mtree to the mass matrix evaluated at its arithmetic squared mass average (black). The ξ dependence of the iterative and the zero momentum procedure is about the same, with ∆ M ξ amounting to 8 and 9%, respectively, at ξ = 100. For the arithmetic squared mass average method, however, we again find that the ξ dependence is very small. Overall, the gauge dependence of the h u -like mass in scenario P1 is larger than in P2.  Right: The h u -like Higgs boson mass of scenario P2 as a function of ξ by applying iterative (red), the R 0 (blue) and the R mtree -method (black). Lower panels: The corresponding ∆ M ξ values in percent, as a function of ξ. Figure 15 depicts the gauge dependence of the partial widths of the decays H ± → W ± h u (left) and H ± → W ± h s (right) at LO (full lines) and NLO (dashed lines) by applying in Eqs. (62) and (65), respectively, the three different approximations for the matrices that diagonalize the corresponding loop-corrected mass matrices, namely the Z matrix (red), R 0 (blue) and R mtree (black). The corresponding decay widths are denoted by the superscripts 'impr', 'R 0 ' and 'R mtree '. The lower panels display the corresponding ∆ Γ ξ values. The inspection of the plots shows that the ξ dependence of the NLO widths not always decreases compared to the LO one. Moreover, there is no pattern for the two decays that allows to decide which of the three approximations induces the smallest gauge dependence in the NLO widths. A closer investigation reveals that the mixing between h u and h d is responsible for the gauge dependence of H ± → W ± h u and the mixing between h s and h d for the one in H ± → W ± h s . Overall, however, the gauge dependence of the partial widths is much smaller than for the parameter point P1, with maximum values of ∆ Γ ξ around −22% for H ± → W ± h u (for Γ NLO,impr ) and 14% for H ± → W ± h s (for Γ LO,R mtree ). The relative NLO corrections at ξ = 1 amount to -20% (for Γ impr ) for the former and to -12% for the latter decay (for Γ R mtree ), however, so that the gauge dependence is of the order of the loop correction.

Conclusions
In this paper we investigated the influence of the gauge parameter both on the higher-order corrections to the NMSSM Higgs boson masses and the partial decay width of H ± → W ± h i , calculated in general R ξ gauge. The gauge dependence enters through the mixing of loop or-  Figure 15: P2: Decay widths for H ± → W ± h u (left uppper plot) and H ± → W ± h s (right uppper plot) as functions of ξ using the iterative (red), the R 0 (blue) and the R mtree method (black). Lower panels: Corresponding ∆ Γ ξ in percent, as function of ξ.
ders: for the masses, this happens due to the application of an iterative method to determine the loop-corrected mass values. This is transferred to the decay width as phenomenology requires the inclusion of the mass corrections to the external Higgs bosons in order to match the experimentally measured values. These are calculated up to two-loop order including higher-order terms through the application of the iterative procedure. Gauge dependence then enters the process through different mechanisms. On the one hand there is a mismatch between the use of tree-level masses in the propagators of the internal lines and in the Higgs-Goldstone boson couplings appearing in the computation of the loops, and the use of the higher-order-corrected Higgs mass for the external Higgs bosons. The latter prevents the cancelation of IR divergences when adding up the virtual and real corrections. While this can be cured by an appropriate adaption of the involved couplings, the second source of the gauge dependences persists: it stems from the resummation of higher orders that enters both through the external Higgs boson mass and the mixing matrix applied to set the external Higgs boson OS. The latter has a particularly large impact as we found by applying different approximations to determine the higher-order masses and mixing matrix elements. The relative gauge dependence can then largely exceed the relative size of the loop correction itself, so that for the interpretation of the results the specification of the used gauge is crucial. By analyzing different parameter sets, we found that the impact of the gauge dependence depends on the chosen parameter point and can vary substantially depending on the applied parameter set.
work is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.01-2017.78.