Generic calculation of two-body partial decay widths at the full one-loop level

We describe a fully generic implementation of two-body partial decay widths at the full one-loop level in the SARAH and SPheno framework compatible with most supported models. It incorporates fermionic decays to a fermion and a scalar or a gauge boson as well as scalar decays into two fermions, two gauge bosons, two scalars or a scalar and a gauge boson. We present the relevant generic expressions for virtual and real corrections. Whereas wave-function corrections are determined from on-shell conditions, the parameters of the underlying model are by default renormalised in a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\text {DR}}$$\end{document}DR¯ (or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\text {MS}}$$\end{document}MS¯) scheme. However, the user can also define model-specific counter-terms. As an example we discuss the renormalisation of the electric charge in the Thomson limit for top-quark decays in the standard model. One-loop-induced decays are also supported. The framework additionally allows the addition of mass and mixing corrections induced at higher orders for the involved external states. We explain our procedure to cancel infrared divergences for such cases, which is achieved through an infrared counter-term taking into account corrected Goldstone boson vertices. We compare our results for sfermion, gluino and Higgs decays in the minimal supersymmetric standard model (MSSM) against the public codes SFOLD, FVSFOLD and HFOLD and explain observed differences. Radiatively induced gluino and neutralino decays are compared against the original implementation in SPheno in the MSSM. We exactly reproduce the results of the code CNNDecays for decays of neutralinos and charginos in R-parity violating models. The new version SARAH 4.11.0 by default includes the calculation of two-body decay widths at the full one-loop level. Current limitations for certain model classes are described.


Introduction
While the large hadron collider (LHC) has completed the standard model (SM) of particle physics with the discovery of a scalar which has all expected properties of the long searched for Higgs boson [1][2][3], there is no indication for new physics up to now. This has lead to impressive exclusion limits for particles predicted by either supersymmetry (SUSY) or other extensions of the SM which were proposed to resolve the open questions of the SM. However, these exclusion limits for beyond the standard model (BSM) particles depend strongly on the decay properties of these particles. For instance, it is well known that the often cited limits for SUSY squarks and gluinos of 1.8 TeV and more hold only in vanilla models where these states decay to 100 % into a given final state [4][5][6][7]. Once realistic decay patterns for the particles are used, the limits become much weaker [8][9][10][11]. Thus, a precise knowledge of the branching ratios of BSM states is necessary to be able to draw firm conclusions from the null results. On the other hand, once a new particle is discovered, precise calculations become especially important to extract the underlying parameters and compare against the predictions of many different models.
There has been a lot of effort to improve the predictions of the decay widths for new Higgs-like scalars not only in the minimal supersymmetric standard model (MSSM)  and the next-to-minimal supersymmetric standard (NMSSM) [40], but also in several singlet and doublet extensions of the SM [41][42][43][44][45][46][47]. These results are implemented in public tools such as HDECAY [48,49], FeynHiggs [50][51][52] or NMSSMCALC [53]. However, for the plethora of other states, tree-level results are often used. Exceptions are the MSSM, where one-loop corrections to all sfermions and gauginos were discussed in Refs. ; and neutralino and chargino decays in the NMSSM [75,76]. For other SUSY models with R-parity violation and CP violation, only a few selected decay modes were discussed so far in Refs. [77,78]. The available codes to study decays at the oneloop level in the MSSM are SDECAY [79], SUSY_HIT [49] and SFOLD [80] for sfermion decays, FVSFOLD for flavour violating squark as well as gluino decays, and SloopS [76] and CNNDecays [75,77] for neutralino and chargino decays without and with R-parity violation.
This limited number of codes and supported models has to be seen in contrast to the increasing number of models which are currently studied. With the increasing limits on the SUSY masses within the MSSM, other ideas for new physics are seeing more and more attention. In order to be able to also give more accurate predictions for the decays in non-minimal SUSY models or also in non-supersymmetric extensions of SM, a high-level of automatisation is needed. A very powerful ansatz to obtain robust results for BSM models has been established with the Mathematica package SARAH [81][82][83][84][85][86]: SARAH derives from a short model file all analytical properties of a given model. This information together with generic expressions for various observables is then used to generate Fortran code for SPheno [87,88] which can be used to obtain numerical results. Up to now, one-and twoloop masses [89][90][91], one-loop flavour and precision observables [92], as well as two-and three-body tree-level decays could be obtained via this setup. We have now enhanced the decay calculation to the next level by a generic ansatz to calculate two-body decay widths at the full one-loop level. These extensions are now available with SARAH 4.11.0. In this paper we give all necessary details about the calculation, including the renormalisation scheme; the generic expressions for virtual and real corrections; and the handling of ultraviolet and infrared divergences.
While wave-function corrections are determined from onshell conditions, the default settings use a DR (or MS) renormalisation for the parameters of the underlying model. However, the user can also define model-specific counter-terms in SARAH to be used in the numerical evaluation in SPheno. For now the self-energies of all particles of the underlying model are available for this purpose. Since many particle species receive significant higher-order corrections to their masses and mixing beyond tree level, we also allow the inclusion of mass and mixing corrections for the involved external states. This needs a careful treatment of the infrared divergences, for which we add an infrared counter-term making use of modified Goldstone boson vertices. The setup also supports loop-induced decays. An extension to models with CP violation or additional charged and massive, coloured vector particles is left for future work. Higgs boson decays, which are very sensitive to corrections of the external states, will be more thoroughly addressed in the future.
The paper is organised as follows: In Sect. 2 we discuss the technical details of the implementation employing external tree-level masses. The incorporation of higher-order corrections for the external states is lined out in Sect. 3. In Sect. 4 we explain how the new features of SARAH and SPheno can be used. In Sect. 5 we present some results obtained with the new machinery: we first show the implementation of counter-terms for two SM examples and then compare our implementation in SARAH with other public codes as SFOLD, HFOLD, CNNDecays. We conclude in Sect. 6. The appendix contains all relevant generic expressions for virtual and real corrections as well as a derivation of the employed Goldstone boson vertices.

Calculation of decay widths at the full one-loop level
In this section we discuss the technical details of the calculation of two-body decay widths at next-to-leading order for decays that are mediated through a tree-level diagram X → Y 1 Y 2 . Our implementation can handle the decays S → SS, S → SV , S → V V , S → F F, F → F S and F → F V , where S denotes a scalar, F a fermion and V a heavy gauge boson. For loop-induced processes V can also be a photon or gluon. At next-to-leading order we include full QCD and electroweak corrections. Thus, apart from ultraviolet divergences, which we need to address through the renormalisation of the parameters of the underlying model, infrared divergences due to massless photons and gluons have to be taken care of. For loop-induced decays the subsequent discussion simplifies substantially, since neither ultraviolet nor infrared divergences have to be tamed, i.e. also the detailed renormalisation of parameters is not of relevance. We continue as follows: we describe the generic form of unpolarised squared matrix elements for two-body decays in the subsequent subsection and thereafter present the various ingredients in terms of tree-level and one-loop amplitudes. This includes vertex and wave-function corrections as well as a discussion of counter-terms. Then in Sect. 2.6 we discuss the relevant real corrections being 1 → 3 processes, before we combine the results in Sect. 2.7. Finally, we list the limitations of our implementation in Sect. 2.8.

Generic unpolarised squared matrix elements
For any two-body decay we write its partial width in the form where m X , m Y 1 and m Y 2 are the masses of the mother and daughter particles in the initial and final state, respectively.
We denote their momenta with p 0 , p 1 and p 2 , respectively. The sum runs over all helicities (h) and polarisations ( p) in the initial and final state. A symmetry factor C S and colour factor C C have to be employed. The symmetry factor is C S = 1 by default. For X = F we have C S = 1 2 , if Y 1 = Y 1 and Y 2 = Y 2 . For X = S it is C S = 1 2 , if Y 1 = Y 1 and Y 2 = Y 2 and Y 1 = Y 2 . Therein Y denotes the antiparticle of Y . The colour factor C C for a decaying colour singlet is equal to the dimension of the final states under SU(3) C . For example, for a colour octet decaying into triplets, it yields C C = 1 2 , while for more complicated colour configurations C C can easily be extracted from the colour-dependent part of the vertex triggering the decay: the colour of the initial state is fixed and a sum over all possible colour combinations in the final state is performed. The Källén function λ is given by For decay modes with fermions and gauge bosons in the initial and/or final state the matrix elements are a sum over Lorentz structures; we label these with a lower index as M i and therefore split the total squared amplitude in sums of contributions M i M * j , which are multiplied with different kinematic dependences obtained from helicity and polarisation sums. The structures and their sums are given by For S → V V we split the squared amplitude as follows: Lastly the squared amplitude for F → F V is given by (2.8) We implemented special cases for final states with vanishing masses, which are not given here. 1 1 In the calculation of one-loop decays we introduced a minimal allowed mass for fermions and scalars of 10 −15 GeV. Smaller masses 2.2 Tree-level amplitudes For the two-body decays at tree level the contributions to the matrix elements M T i can be directly identified with the (leftand right-handed) couplings as follows: The conventions for the parametrisation of the vertices are summarised in Appendix A. For F → F V M T 3 and M T 4 and for S → V V M T 2 vanish at tree level, but contributions are generated at the one-loop level.

One-loop amplitudes
Before discussing the detailed form of vertex and wavefunction corrections, we show their combination with the previously presented results. Once the amplitude due to vertex corrections M V and due to wave-function corrections M W are split into M V i and M W i , which encode the various contributions to different combinations of helicities and polarisations, they can be added to the tree-level amplitudes as follows: For an exact next-to-leading order calculation the complexconjugated part of this amplitude M * i is inserted in the complex-conjugated matrix elements of Eqs. (2.3)-(2.8), whereas M T i is used for the non-conjugated ones. The total partial width is obtained from the real part of the full expressions in Eqs. ( Our results are expressed in terms of Passarino-Veltman integrals obtained through dimensional reduction (DR). 2 Thus, the ultraviolet divergences can be split off in terms of = 1 −γ E +ln(4π), where regularises the divergence and equals the difference with four dimensions d = 4−2 and γ E is the Euler-Mascheroni constant. Terms denoted by r need to be set to zero in dimensional reduction and correspond to the difference with respect to dimensional regularisation, i.e. it yields r = 1 for calculations performed in the minimal subtraction scheme (MS). By default SARAH sets r = 0 in SUSY models and r = 1 in non-SUSY models. In order to match mass dimensions correctly in less than 4 dimensions, dimensional reduction also introduces a new scale Q, the renormalisation scale. The generic particle U denotes a Faddeev-Popov ghost. In Feynman-'t Hooft gauge their masses are identical masses to the gauge bosons masses, i.e. also ghosts obtain the subsequently discussed regulator mass. Infrared divergences due to massless photons and gluons are regularised through a finite regulator mass. There are no diagrams that contain both photons and gluons. The cancellation of infrared divergences will be addressed in Sect. 2.6, whereas the cancellation of ultraviolet divergences is obtained by adding the subsequently discussed corrections M W .
The combinatoric part to populate the generic diagrams with all possible field insertions in a given model is done by SARAH. SARAH also checks for possible symmetry factors which appear if in the topologies 2-4 in Fig. 1 two real and identical particles are in the loop. In addition, it calculates relevant colour factors to be multiplied with the interference terms M T (M V ) * .

Wave-function corrections
The amplitude M W contains the corrections due to wavefunction normalisation as well as the counter-term for the tree-level coupling. They cancel the ultraviolet divergences of the vertex corrections M V and are mostly determined through renormalisation prescriptions, in contrast to M V . Omitting the complication of fermions and gauge bosons for a moment the amplitude M W i jk for a vertex of the form with the counter-term δc of the tree-level coupling c and the wave-function corrections δ Z for the three particles involved. In the last three terms a sum over l has to be performed. In the following we will first describe the derivation of the wavefunction corrections and then comment on the counter-term for the tree-level coupling. For the wave-function corrections we employ an on-shell scheme for the three fields S, V and F. For the fermions we distinguish left-and right-handed components F L and F R . In all cases we allow for mixing among particles induced through loop effects, such that the wave-function corrections are generally matrices. We have 14) In order to determine the wave-function corrections δ Z from on-shell conditions, we need the self-energies for our three particle species. Their notation can be read off from the inverse propagators at the one-loop level, which we write as follows: The renormalised self-energies are indicated throughˆ and compared to the unrenormalised ones and , which are of relevance for the subsequent discussion. V V and SS are the self-energies of the gauge bosons and scalars, respectively. For the gauge bosons we are only interested in the transverse part V V . The only mixing induced between the gauge bosons of the SM is among the photon and the Z boson. P L and P R are the left-and right-handed projection operators, which split the self-energies of the fermions in L , R , SL and S R . The topologies which can contribute are shown in Fig. 2. Moreover, all possible generic diagrams contributing to the fermion, scalar and vector bosons self-energies are shown in Appendix A in Figs. 12, 13 and 14. We give also in Appendix A the expressions for the generic amplitudes for the self-energies and their derivatives. Note that the above structure for the gauge bosons implies the usage of Feynman-'t Hooft gauge. The derivatives of the wavefunction corrections, denoted with˙ and˙ , are defined as follows: (2.21) Demanding on-shell conditions for the external states fixes the wave-function corrections. Their derivation can for example be found in Refs. [75,96] and results in similar expressions for scalars and gauge bosons: For the fermions we need to distinguish four cases: (2.24) By Re we indicate that and entering δ Z include only the real parts of the loop functions, whereas couplings enter with real and imaginary components. In case of CP violation the definition of wave-function corrections usually distinguishes between in-and outgoing particles in order to correctly multiply absorptive parts of self-energies with complex couplings, see the appendix of Ref. [66]. This is beyond our implementation. We note that despite the fact that we employ on-shell conditions to determine the wave-function corrections our external particles are not necessarily on-shell particles, see the discussion at the end of this section.
With this setup at hand we can also define counter-terms to be used for tree-level rotation matrices, which at lowest order transform gauge into mass eigenstates. Those counterterms enter the counter-term of the tree-level coupling. For any particle species in mass eigenstates, which is obtained from gauge eigenstates through i = R i j j , the counterterm is given by For fermions left-and right-handed states are rotated with two matrices, such that two counter-terms employing leftand right-handed wave-function corrections also need to be defined. For Majorana fermions we refer to Ref. [75]. It is well known that the definition of such counter-terms for mixing matrices based on the wave-function corrections needs a proper treatment of Goldstone boson tadpole contributions in order to achieve gauge invariance; see Ref. [75] for a more detailed discussion. Since we work in Feynman-'t Hooft gauge we can completely omit these Goldstone boson tadpole contributions, since they ultimately cancel between the wave-function corrections and the counter-term of the mixing matrices. As for the vertex corrections, SARAH inserts all combination of particle species in the generated code, and includes colour as well as symmetry factors. The non-trivial and mostly non-automatisable part of the calculation of two-body partial decay width is the renormalisation prescription used for the bare parameters of the underlying theory, which enter the tree-level coupling counter-term of the two-body decay under consideration. The counterterms are usually chosen depending on the model and process. However, a simple DR (or MS) prescription for the renormalisation of the parameters of the underlying theory is always easily applicable: from the β functions and anomalous dimensions used for the renormalisation group equations implemented in SARAH [97][98][99][100][101][102][103][104] we can define all counterterms of the parameters of the underlying theory to be just proportional to the pure ultraviolet divergence only. We will refer to this scheme as DR (or MS) scheme in the following. It is well known that this scheme will not perform well in various cases. Therefore, the user of SARAH can define their own counter-terms; see Sect. 4.2 for a more detailed discussion. We also add an example of a proper renormalisation of the electric charge in Sect. 5.
A consequence of the application of the DR (or MS) scheme is that our partial decay widths are left with a dependence on the renormalisation scale Q introduced through the regularisation of ultraviolet divergences. This is most prominent in the running of the parameters that enter the tree-level coupling obtained from the renormalisation group equations, which is not cancelled at the one-loop level. In the generated code the scale Q is by default set to the average stop mass √ mt 1 mt 2 in supersymmetric models and the top-quark mass m t in non-supersymmetric models. However, the user can control the scale Q in the input file, either throughout SPheno or only for the calculation of the decays at one-loop level; see Sect. 4.2. A common choice for the renormalisation scale is also Q ∼ m X close to the mass of the decaying particle X . We refrain from making it the default option, since Q ∼ m X slows down the numerical evaluation substantially.
In this case loop contributions need to be evaluated multiple times. We recommend to vary the scale to check the stability of the partial decay width calculation, as we demonstrate in Sect. 5 for the decay of the SM Higgs boson into bottom quarks. If the scale is changed throughout SPheno keep in mind that also masses and thus kinematics can change. For a full on-shell calculation the scale dependence also completely vanishes. We demonstrate this for the decay of the top-quark in Sect. 5. In order to achieve a renormalisationscale-independent result, external states have to have fixed masses and mixing, which for gauge bosons and fermions can be achieved through the settings explained in Sect. 4.2. Until now we ignored the fact that particles receive higherorder mass corrections. 3 By construction we have to employ the mass values at lowest order throughout the calculation. We will discuss in Sect. 3 how for the external states mass corrections and mixing beyond tree level can be incorporated into our calculation. If we allow for mass corrections we limit the discussion to the inclusion of DR (MS) corrections to the masses, whereas full on-shell prescriptions for BSM particles (as it would be appropriate) are left for future work. With the outlined procedure in the previous subsections, we obtain a gauge-independent and ultraviolet finite result for the partial width X → Y 1 Y 2 , which in the most general case, however, is scale dependent. As mentioned, the cancellation of infrared divergences is addressed in the next section.

Real corrections
In the previous calculation of vertex and wave-function corrections we regularised infrared divergences through the introduction of a finite, but small regulator mass for the photon and/or gluon. The artificial dependence of the cross section calculation on that mass is cancelled by adding the real emission of a photon and/or gluon to the two-body decay, i.e. by adding three-body decays. For a soft photon and/or gluon a divergence is induced, which can again be regularised through a mass and cancels the mass dependence from the vertex and wave-function corrections. With the help of FeynArts [93] and FormCalc [94] we generated generic results for the emission of one additional photon γ or gluon g for the previously discussed two-body processes S → SS, 26) where k denotes the momentum of the photon or gluon and momenta with upper index 0 equal the zeroth component of the corresponding four vector. External momenta are set to We then have C S = C S for S → V V and S → SV , otherwise C S = 1 2 C S . The charge and colour structure is encoded in matrices C i j explained in Appendix B, which is why Eq. (2.26) only contains an additional c for colour to be summed over. It is clear that the real corrections due to photon emission and gluon emission can be calculated individually and summed up afterwards. By rewriting denominators in terms of eikonal factors, the above integrals can be mapped onto where p i, j ∈ {p 0 , p 1 , p 2 } and the minus signs refer to cases where p i, j equals the momentum p 0 of the initial particle X . The notation follows Ref. [96], where also results for the relevant integrals are shown. Only integrals with double lower indices are infrared divergent and thus dependent on the regulator mass in addition. We present our results in Appendix B. Through our procedure we calculate the full soft-and hard emission of such photons and gluons and thus for the threebody decay S → SV + γ /g also include the four-point interaction, which does not diverge as the regulator mass approaches zero. The correct charge and colour factor assignments in the real corrections are done by SARAH as explained in Appendix B. Where possible we compared to the analytic results for real corrections implemented in SFOLD [80] and HFOLD [107] as well as the result presented in Ref. [75]. Apart from finite contributions in S → SV we found complete agreement. We avoid additional collinear divergences by keeping finite masses for all three particles in the initial and final state of our two-body decay calculation, if they interact with photons or gluons. Thus, this problem does not arise for e.g. final-state neutrinos. For fully massless charge-and colour neutral particles in the final state we implemented dedicated routines for F → F Sγ and S → F Fγ , where one final-state fermion F can be massless. Keep in mind that if final-state charged or coloured particles are very light, large collinear logarithms can induce a bad numerical behaviour of our routines. This should not cause problems in practical applications unless charged or coloured states with very small masses ( keV) are present.
Lastly note that since the real correction decay widths are gauge independent, we performed the calculation in unitary gauge for simplicity. This ensures that the results depend only on the gauge couplings and the original tree-level vertex, and we are not obliged to include would-be Goldstone boson vertices as we do for the corresponding loop corrections. The exception is the decay S → SV + γ /g, where gauge invariance fixes the form of the four-point vertex in terms of the three-point one, and we implicitly assume this relation.

Combination of results
The partial width at next-to-leading order is thus obtained as follows:

Current limitations
In the approach described so far, we made some assumptions which make the results not applicable to all models which are currently supported by SARAH.
• While complex parameters in all calculation can be handled in principle, the setup is not yet supposed to be used for CP violation. The reason is that, for decays of real particles into complex final states, only the decay mode Y 1 Y 2 is calculated, while Y 1 Y 2 is assumed to have the same partial width. Also note that in the case of CP violation a common approach is to define extended wave-function corrections as discussed in the appendix of Ref. [66]. • The calculation of the real divergences has neglected the possibility of massive, coloured vector bosons as for instance in Pati-Salam, deconstructed, or trinification models. • When using loop-corrected external masses as described in the next section, we need to cure all infrared divergences through a proper treatment of Goldstone boson vertices. Currently we assume that the W boson is the only massive, charged vector boson, such that models with a W cannot be used with loop-corrected external masses. • Gauge boson decays are not implemented yet. This is partially due to the previously two mentioned limitations.
On the other hand for decays of the gauge bosons of the SM our framework can easily be extended, which we leave for future work.

Higher-order corrections to the external states
Our previous discussion was based on the usage of tree-level masses for internal as well as external particles. However, the masses of various particle species receive significant higherorder contributions. One way to address this is to adopt an on-shell scheme throughout the calculation, but pure on-shell schemes are not always the best choice for such calculations, as is well known from the Higgs sector of the MSSM. Also even if the calculation is performed in terms of on-shell states in particular in supersymmetric models the limited number of renormalisable parameters in the Lagrangian does not allow for a renormalisation procedure where all on-shell masses correspond to their tree-level values. Instead, if we do not want to change our renormalisation prescription, we should use the LSZ reduction formula to connect S-matrix elements of on-shell states with Feynman diagrams in our scheme. This results in external normalisation factors and external loop-corrected masses (which need to be distinguished from the previously discussed tree-level masses and mixing matrices). A schematic picture is shown in Fig. 3. For Higgs bosons the approach is discussed in detail in Ref. [108], where such wave-function normalisation factors are denoted Z -factors. In particular, as noted there, since we are working at only one loop, there are different truncations of the perturbative series that we can make and different approximations can also be made for expediency. Here we outline the choice(s) that we have made.
Firstly, in the gauge boson sector we recommend to use onshell values for gauge bosons, see the discussion in Sect. 4. 3 Schematic picture of our method to include higher-order mass and mixing corrections. The calculation presented in Sect. 2 is combined with external normalisation factors and external loop-corrected masses For scalars and fermions, however, we introduce matrices that we denote U . Let us introduce our notation for the example of n scalars S i , which are mass eigenstates obtained from gauge eigenstates S i = R S i j S j at tree level and mix at higher orders: following Eq. (2.18) the mass matrices beyond tree level take the form with the tree-level masses m i . Let us suppose that for the calculation of the external masses a MS or DR scheme is preferred, i.e. the self-energiesˆ andˆ are renormalised such that only the corresponding ultraviolet divergent part is omitted. In a first approximation we set p 2 = 0 and diagonalise the obtained mass matrix M(0) through a unitary (n × n) matrix U 0 . The tree-level mass eigenstates S i are thus rotated into statesS i = U 0 i j S j with masses m i . This matrix U 0 incorporates the additional mixing induced at higher orders and in principle corresponds to the Z -factors in the p 2 = 0 approximation of Ref. [108]. 4 It is used to rotate the tree-level, vertex and counter-term corrections uniformly by applying it at the amplitude level. For the decayS i →S jSk we for example shift the amplitudes by for M = M T , M V and M W . Also we define rotated treelevel couplings c i jk in the same manner to be used in the calculation of tree-level ampltiude and real corrections as discussed subsequently. This concept can be very similarly employed for fermions, where again left-and right-handed mixing matrices U have to be introduced. Before we discuss the cancellation of ultraviolet and infrared divergences let us note that we also implemented two more methods to obtain the mixing matrix U : instead of setting p 2 = 0 an alternative choice is to use p 2 = m 2 i . This results in the mass eigenvalue m i , which is used to repeat the procedure iteratively with p 2 = m i until the mass determination stabilises. Our default choice is that the relative change between the masses of two iterations should be below 10 −6 . This procedure needs to be performed for each mass eigenstateS i separately and the matrix U p is determined row by row and is thus not unitary any more, as is also well known from the general form of Zfactors. Lastly a possible choice is p 2 = m 2 1 , i.e. the external momenta is chosen to be equal to the lightest mass eigenstate. In this case U m 1 is again a unitary matrix. We note that the outlined procedure to determine m and U can be performed beyond one-loop level, i.e. for supersymmetric Higgs boson masses corrections at the two-loop level can be incorporated.

Ultraviolet and infrared divergences
The application of external masses m different from the treelevel values m and mixing matrices U in addition to tree-level mixing induces a problem with the cancellation of ultraviolet and infrared divergences. The first problem can be solved easily. We employ tree-level masses m for all propagators of loop functions as well as external momenta entering loop functions. This applies to vertex and counter-term corrections and guarantees the cancellation of ultraviolet divergences.
The infrared problem is more demanding. In order to achieve the cancellation of infrared divergences we define infrared counter-terms. These counter-terms encode the mismatch between the masses and mixings of internal and external states and are formally of higher order. These counterterms are used to shift the wave-function and vertex corrections: The aim is to cancel the infrared divergences stemming from where IR(M) takes only the infrared divergent part of the amplitude M. The definition of such counter-terms in supersymmetric models is a common strategy: Ref. [66] introduces an infrared counter-term for the decayst 1 →b i W + , see Eq. (191) in Ref. [66], which exactly encodes the difference between on-shell, i.e. loop-corrected, and tree-level masses due to the limited number of renormalisable parameters in the stop and sbottom sectors. Reference [109] discusses the introduction of such counter-terms for the heavy Higgs boson decay H → W + W − . Reference [110] discusses these aspects in detail in the context of a generalised narrow-width approximation, where also infrared divergent parts in the loop contributions are sorted out and evaluated at a common mass scale. It is clear that the subtraction and re-addition of infrared divergent logarithms as discussed before induces a spurious dependence on other masses, namely on the masses being the counterpart of the regulator mass(es) in the logarithms. This is unavoidable, however, numerically of minor relevance.
In practice, the following procedure is applied: 1. We calculate the virtual corrections using tree-level masses. 2. We extract the infrared divergences of all two-and threepoint functions using the results given in Appendix C. These result are used to obtain IR M V,W . 3. We use loop-corrected masses m throughout all infrared divergent diagrams for the external legs and the particles in the loop. We take again the infrared divergent parts of these amplitudes to obtain IR M V,W . 4. The calculation of the kinematics as well as of the helicity and polarisations sums for both, the virtual and real corrections, is done with loop-corrected masses. 5. Lastly, the usage of an additional external mixing applied through the mixing matrices U , named U -factors, works as follows: we rotate the amplitudes of the treelevel, wave-function and virtual corrections according to Eq. (3.2). Instead for the contribution of the infrared counter-term we use rotated tree-level couplings c rather than rotated amplitudes M. Those rotated couplings also enter the calculation of the real corrections. In this context we note that by construction the infrared counter-term always contains exactly one occurrence of the coupling c of the tree-level two-body decay.
These steps give for most cases infrared finite results. However, there is one complication: if the infrared counter-term contains loops with massive gauge bosons, then there will necessarily also be related diagrams with charged Goldstone bosons, and the gauge symmetries require several relationships between the couplings-and masses-of the internal and external particles in order for the infrared divergences to cancel. If we were to apply the above procedure then the infrared counter-term would not be gauge invariant; for these diagrams we therefore use loop-corrected masses and couplings. Note that if we used unitary gauge we could avoid a discussion of corrected couplings in Goldstone boson vertices. Denoting a would-be Goldstone boson by G, massive gauge bosons by V G with masses m G V and massless ones by γ a , real scalars as S i with masses m i and Dirac fermions as F I with masses m I , the relevant couplings are The couplings c a i j , c a GG , c a,L/R I J are just generators of the unbroken gauge group in the appropriate representation multiplied by the unbroken gauge coupling. We find that we must enforce the following relations, which we derive in Appendix D: (3.10) The implementation in SARAH currently assumes that the gauge sector is that of the SM; so there are no infrared divergent diagrams with neutral Goldstone bosons, and we do not shift their couplings to loop-corrected masses. In practice, a new set of Goldstone vertices is derived by the following relations which is then used in the calculation of the IR shifts: Note that in Eq. (3.11) we explicitly assume no CP violation.
Employing the outlined procedure we obtain partial decay widths at next-to-leading order with full cancellation of ultraviolet and infrared divergences. Though the application of loop-corrected masses in the infrared counter-term can induce a spurious higher-order gauge dependence, for phenomenological purposes this is, however, small, see e.g. Refs. [75,77]. Note that for external heavy gauge bosons of the SM we give the option to put the heavy gauge bosons onshell, such that the cancellation of a gauge dependence in the real corrections among internal gauge bosons and Goldstone bosons is always guaranteed.

Mixing of species
Particular attention is needed in the calculation of processes where self-energy diagrams allow for the mixing between different particle species beyond tree level. As an example (CP-odd) Higgs bosons including the neutral Goldstone boson can mix with the Z boson and even the photon. Then wave-function corrections to the two-body decay come with an internal propagator with a state different from the external state. Such diagrams potentially need to sum up correctly to ensure a gauge-independent partial width. For this purpose, in order to avoid unphysical poles the momenta flowing through the propagators have to match. As an example, Ref. [36] keeps tree-level masses in diagrams mixing Higgs bosons and the Z boson/Goldstone boson in the calculation of Higgs decays to Higgs bosons. The Slavnov-Taylor identities then ensure that the sum of the Z and Goldstone contributions give zero; see also Refs. [46,111]. Generally such diagrams are beyond our generic implementation described here and require a process-dependent treatment, i.e. they are not included. Still, in our setup they can easily be added.

Loop-induced decays
We finish with a remark about loop-induced decays like F i → F j γ . Since infrared divergences do not appear for these processes at the one-loop level, there are fewer restrictions on which masses should be used to calculate the vertex one-loop diagrams. As a default setting we therefore use loop-corrected masses everywhere. The reason is that these decays are of particular importance in regions of kinematical thresholds. Thus, the mass splitting between the two massive states should be taken properly into account in the one-loop calculation.

SARAH-SPheno interface
The possibility to calculate one-loop decay widths is available from SARAH 4.11.0. This is a new feature of the SARAH interface to SPheno which was established with SARAH 3.0.0: SARAH generates Fortran code which can be compiled together with the standard SPheno package to obtain a spectrum generator for a given model. The main features of a spectrum generator obtained in that way are a precise mass spectrum calculation including two-loop corrections to real scalars [89][90][91], a prediction for many precision and flavour observables [92] and up to now the calculation of two-and three-body decays mainly at tree level.
The general procedure to obtain the SPheno code for a given model starts with the download of the most recent SARAH version from HepForge: This creates a new binary bin/SPheno$MODEL which reads all input parameters from an external file. SARAH writes a template for this input file which can be used after filling it with numbers by typing: The output is written to SPheno . spc .$MODEL and contains all running parameters at the renormalisation scale, the loop-corrected mass spectrum, as well as all other observables calculated for the given model and parameter point.
The time for generating the Fortran code for the oneloop two-body decays as well as the compilation time of SPheno are extended by these new routines. Therefore, in the case that the user is not interested in the loop-corrected two-body decays, they can be turned off via: They can be permanently turned off for a given model by adding 1

SA ' A d d O n e L o o p D e c a y = False ;
to SPheno.m. Usually, the calculation of the one-loop decays triggers also the calculation of the RGEs even when using the option "OnlyLowEnergySPheno = True;" to generate the SPheno code. The reason is that the βfunctions are used to check the cancellation of ultraviolet divergences. However, for non-supersymmetric models, in particular in the presence of many quartic couplings in the scalar potential, the RGE calculation can be very timeconsuming. In this case the option 1

SA ' N o R G E s f o r D e c a y s = True ;
skips the RGE calculation. Of course, the verification of the cancellation of ultraviolet divergences will not be performed with this setting.

Definition of counter-terms
We included the possibility to define counter-terms to be used in the calculation of the one-loop decays. This is done in SPheno.m via the new array RenConditionsDecays. For instance, the standard renormalisation conditions for the electroweak gauge couplings are set via: We give an example for the application of the above electroweak counter-terms and their derivation in Sect. 5.1. If RenConditionsDecays is not defined, a pure MS/DR renormalisation for the bare parameters of the underlying model is performed. The counter-terms can also be turned on/off in the numerical session via new flags in the SPheno input file as explained in the next subsection. The conventions are: • The names for the counter-terms are the names of the corresponding parameter starting with d. • For a rotation angle X, no counter-term for the angle itself is introduced, but for the trigonometric functions involving that angle. Those are called dCosX, dSinX and dTanX.
The following objects can be used to define the counterterms: • Parameters of the model: the internal SARAH names must be used. • Masses of particles in the model: those are called MX where X is the name of the particles in SARAH. • Self-energies for scalars and vector bosons: those are called PiX where X is the name of the particles in SARAH. • The derivatives of self-energies of scalars and vector bosons: those are called derPiX where X is the name of the particles in SARAH. When SARAH is finished generating the SPheno output, a list of all self-energies and their derivatives which are available in SPheno is stored in SA'SelfEnergieNames, and the names for all counter-terms are saved in SA'ListCounterTerms. One needs to be careful when using self-energies or their derivatives for particles which come with several generations. In this case, the objects defined above are arrays with three indices. The last two indices give the involved generations, the first one the external momentum, e.g.
When defining the counter-terms, commands for matrix or tensor operators should already have been evaluated in Mathematica. Although we offer the possibility to the user to define counter-terms in that way, we want to stress that it has not been tested in practice beyond the examples given in this paper. Thus, this option should be used carefully and the results should be tested throughout, e.g. the ultraviolet finiteness of the partial decay widths is a first test to be performed. Again we emphasise that such counter-terms are for now only applied in the calculation of decay widths. Thus, on-shell prescriptions for the calculation of masses as e.g. known from the neutralino and chargino sector, see Refs. [65,75,105,106], cannot yet be incorporated.

Options for the evaluation with SPheno
There are several options to steer the performed oneloop calculations which can be controlled via the block DECAYOPTIONS in the Les Houches input file for SPheno.
In practice the most important options are:

... # Use loop -corrected masses for loop -induced decays
The following settings are possible: • DECAYOPTIONS[10XY]: the one-loop decays for each particle can individually be turned on (1) or off (0) via these flags. The particle to which a given flag corresponds to is written as comment by SARAH. The default value is 1. The default value is 1.
In addition, the following options exist which are mainly supposed to be used for testing and validation of the virtual and real corrections: The following settings are possible: • Although this block DECAY1L is not officially supported by the Les Houches conventions, there are the following reasons not to overwrite the results of the 'old' calculation: • The sizes of the one-loop corrections are immediately apparent. • The results given in DECAY are not only pure tree-level decay widths, but include in particular for the Higgs decays crucial higher-order corrections adapted from literature. Those are beyond the one-loop corrections which we can provide in the new automatised framework at the moment. • The 'old' calculations also include tree-level three-body decays. We leave the choice of how to combine them with the two-body decay widths obtained at the one-loop level to the user.

Numerical results
We start this section with two examples for the calculation of two-body decay width in the SM, where we demonstrate the relevance of model-and process-dependent counter-terms. Our default implementation makes use of an MS or DR renormalisation of all parameters of the underlying theory. However, for many processes different schemes are actually better suited. This is particularly true for the calculation of electroweak corrections. For this purpose the user of SARAH can define their own counter-terms, as outlined in Sect. 4.2. We show two simple examples in the SM, namely the calculation of the partial decay width t → W b and H → bb.
In the first example we discuss different schemes for the renormalisation of the electric charge, in the second example we show that our MS renormalisation for the bottom-quark Yukawa coupling is actually sufficient. After these examples we continue with a detailed comparison of our implementation with existing codes, among them SFOLD, HFOLD and CNNDecays. Whereas SFOLD and HFOLD are also based on a DR renormalisation of the parameters of the MSSM, the code CNNDecays calculates neutralino and chargino decays in the MSSM, NMSSM and in models with R-parity violation again renormalising the electric charge in the Thomson limit. Thereafter, we compare loop-induced decays with the original implementation of SPheno and lastly show the effect of U -factors in the calculation of two-body decay widths. A more thorough comparison for Higgs boson decays is left for future work.

Renormalisation of α and the top-quark width
First we perform a calculation of the top-quark partial width in the decay t → W b including electroweak and QCD corrections using a SM version of SPheno. Since this process is mediated through the gauge coupling g 2 of SU(2) L at tree level, we will discuss the renormalisation of g 2 in this context. We choose the following input parameters: We neglect quark mixing (i.e. the CKM matrix is approximated by the identity matrix). Note that in a more general approach the renormalisation prescription introduced in Eq. ( 19 12 Next, we provide the decay width for the renormalisation of the electric charge in the Thomson limit of the f f γvertex, i.e. at zero momentum transfer [96]. The counterterms for the electric charge are given by (see Ref. [112] for an overview) where we distinguish two schemes: At NLO we can make use of the very precise value of α(0) together with the corresponding counter-term δ Z e (0) or we employ α(m Z ) and compensate for the shift through the additional terms in δ Z e (m Z ). The relevant self-energies include only contributions from light fermions. Ultimately we also need the renormalisation of the weak mixing angle, which is given by such that in schemes (2) and (3) we obtain In this section we keep α and α s fixed as a function of the renormalisation scale Q and therefore schemes (2) and (3) lead to a scale-independent partial decay width, whereas scheme (1) develops a scale dependence. Also, in scheme (1) it is a priori not clear whether α(0) or α(m Z ) is the better choice, so we provide both values. The results are shown in Table 1. We include the LO partial width as well as the NLO partial width only including electroweak and including electroweak and QCD corrections. We also provide (absolute and) relative corrections in brackets. For NLO,EW+QCD t→W b they are with respect to the partial width NLO,EW t→W b . For scheme (1) evaluated with α(m Z ) we obtain NLO,EW+QCD t→W b = 1.434 GeV at Q = 90 GeV, which demonstrates that the renormalisation scheme dependence is not very pronounced. The absolute QCD correction remains constant (for fixed values of α and α s ) and yields ∼ −9% as expected [113,114]. It is apparent that schemes (2) and (3) yield very comparable results at NLO despite the different input values for α. This is due to the compensation through the shift in the counter-term, which guarantees that the electric charge is renormalised in the Thomson limit. In contrast scheme (1) shows a significant dependence on the input value, where not surprisingly the choice α(0) comes closer to the results in schemes (2) and (3).
We conclude that for precision predictions the proper renormalisation of certain parameters is rather important. The relevant counter-terms for scheme (3) can be defined by the user in the SARAH framework as discussed in Sect. 4.2 through RenConditionsDecays. Note that such counter-terms will only apply at the moment to the calculation of decay widths, not to the calculation of masses.

Renormalisation of Yukawa couplings and fermionic
Higgs decays We also briefly discuss the calculation of H → bb in the SM, which is mediated through the bottom-quark Yukawa coupling Y b , and it turns out that the MS renormalisation of Y b is the preferred choice. A priori, we would expect that the calculation of NLO electroweak corrections [115][116][117] would again be optimally performed using an on-shell renormalisation of all parameters involved. The counter-term of the bottom-quark Yukawa coupling in the on-shell case is given by such that renormalisation prescriptions for δm b and δv are needed. Whereas δm b can be obtained from the selfenergies of the down-type quarks, the on-shell renormalisation of the vacuum expectation value depends on other parameters: one requires that renormalised tadpoles vanish as well as the on-shell renormalisation of the Higgs mass and the Higgs self-coupling, see e.g. Ref. [117]. Also, such counter-terms can be implemented in principle through RenConditionsDecays. However, it turns out that electroweak corrections are small (∼ 1%) for a Higgs mass of m H = 125 GeV. In contrast QCD corrections are much larger and for them the renormalisation of the Yukawa coupling in the MS scheme is more convenient, since it resums large logarithms [118,119]. We demonstrate this effect in Table 2 We see from Table 2 that the electroweak corrections are indeed small. The QCD corrections coincide with the term found in the literature, being 5.667α s (m H )/π ∼ 20.4% [120]. The depicted scale dependence can be used to estimate the remaining uncertainties, which can be significantly reduced by including higher-order QCD corrections beyond one-loop level.

Comparison with other codes
In order to further validate our calculations and implementations, we compared the obtained results for the MSSM and in R-parity violating models against other public tools.
Our comparison is twofold: we compared neutralino and chargino decays into neutralinos and charginos and heavy gauge bosons with CNNDecays, where we employ a full onshell scheme for the gauge couplings, but work with tree-level neutralino and chargino masses. We use the counter-terms for the electroweak sector as outlined in Sect. 4.2. Given that we adjust all input parameters to be identical, we can therefore exactly reproduce the results of CNNDecays.
Secondly, we made use of the three codes SFOLD, HFOLD and FVSFOLD which also use a DR scheme for the renormalisation of the parameters of the MSSM to calculate one-loop corrections to two-body decays. Since these codes do not make use of external U -factors, we have turned them off in our evaluation with SPheno. In addition, we forced all codes to use tree-level (DR) masses in all loops and for the kinematics. Thus, the partial widths and the size of the loop corrections presented in the following might be of limited physical interest since the inclusion of U -factors and external loopcorrected masses can change the results substantially, see also Sect. 5.4. Thus in this section our aim is to only demonstrate the agreement (and disagreement) between the codes. For the comparison with SFOLD, HFOLD and FVSFOLD we have chosen a parametrisation for the general MSSM which depends only on one dimensionful parameter m as follows: All other soft-terms are set to zero. This parametrisation has no physical motivation but was chosen in a way to open many different decay channels to be compared among the codes. In addition, we fixed tan β = 10. We show results for the predicted partial widths when varying m from 300 to 2500 GeV. We also performed a comparison for the loop-induced neutralino and gluino decays which were already implemented in SPheno. The details of this comparison and the outcome are summarised in Sect. 5.3.5.

Neutralino and chargino decays in the MSSM and in bilinear R-parity violation: CNNDecays vs. SARAH
We compared the decay modesχ ± i →χ 0 j W ± andχ 0 i → χ ∓ j W ± as well asχ 0 i →χ 0 j Z andχ ± i →χ ± j Z in the Rparity conserving MSSM and adding bilinear R-parity violation. Bilinear R-parity violation allows an explanation of neutrino masses, but makes the lightest supersymmetric particle (LSP) unstable. Since its decay modes are related to the R-parity violating parameters, which are small in order to explain the size of neutrino masses, the decay width of the LSP is also small. We remain with real parameters, both in the MSSM as well as for the R-parity breaking parameters. We adjust the tree-level masses and mixing as well as the gauge couplings g 1 and g 2 to be exactly identical in both codes and also ensure to choose the same renormalisation scale (through DECAYOPTIONS[1205]), namely Q = m Z . Since we employ the renormalisation of the electric charge in the Thomson limit as outlined in Sect. 4.2, the partial decay widths are in principle all renormalisationscale independent, however, we also want to compare wavefunction and vertex corrections individually. We find full agreement between both codes, i.e. numerically identical results beyond 8 digits in the MSSM. In particular this is also true for the vertex and wave-function corrections individually as well as the individual pieces to the counter-terms. Also for the R-parity violating decaysχ 0 4 → l ∓ W ± in bilinear R-parity violation we find agreement at the per mille level; the smallness of couplings and masses makes those decay modes more sensitive to numerical errors (factors too small or large for the precision of the code). Decay modes into light neutrinos and a gauge boson or a scalar like e.g. χ 0 4 → ν Z or ν S, which are of relevance for R-parity violating scenarios, suffer from bad numerical errors. Therefore, neutrino masses, which are analytically zero at tree level, are set to zero in the calculation of one-loop decays. In contrast, the mixing matrix of neutrinos (and neutralinos) remains exact, such that the associated error is small. As we already explained a detailed check of CP-violating scenarios as discussed in Refs. [67,70,78] is left for future work for the reasons explained in Sect. 2.8.

Sfermion decays in the MSSM: SARAH vs. SFOLD
Next we turn to the comparison of the two-body decays of sfermions in the MSSM. For this purpose, we compared our results against the public code SFOLD 1.2. We have applied several modifications to the code SFOLD: • The variable to use loop masses either in the loops or in the kinematics are set to 0 in SFOLD.F: This is done to ensure that the two codes use exactly the same masses everywhere. • We find a disagreement for the bremsstrahlung routines for S → SV decays. Therefore, we add to line 440 in Bremsstrahlung.F of SFOLD the terms: 1 -2* g1 **2* gt * gtC * Ii -2* g1 **2* gt * gtC * I0up1 • We find for S → SS and S → SV decays huge numerical loop-corrections that could even cause a negative width. We could trace back the problem to diagrams with two massive vector bosons in the loop. 5 The problem is avoided by setting 1 Xipart1 = 0 d0 in Decay.F of SFOLD. With this choice, it is no longer possible to change the gauge in SFOLD, but the results are only valid in Feynman-'t Hooft gauge, sufficient for our comparison. • We find that SFOLD uses a different renormalisation prescription for the rotation matrices: it includes only the divergent parts for the counter-terms, while SARAH calculates the counter-terms from the wave-function renormalisation constants using Eq. (2.25). In particular for S → SV decays this can induce large differences in the one-loop corrections. If the finite parts for the counterterms of the rotation matrices are included, a cancellation between the wave-function corrections and the counterterm correction appears which in sum gives much smaller one-loop corrections. Therefore, we added at the end of the file CalcRenConst99.F of SFOLD the following re-definitions of the counter-terms: 5 In a private discussion with one author of SFOLD the origin of the problem could not be identified. It might be a numerical problem with the high-rank loop integrals which appear when performing the calculation in R ξ gauge. The results for some representative decays for the light and heavy stop are shown in Fig. 4. Here and in the following we give the partial widths at LO and NLO as well as the relative size of the one-loop corrections defined as (5.10) With our described adjustments we find an excellent agreement for the heavy stop decays into a Higgs or a gauge boson and a stop or sbottom. While the corrections for the decays into gauge bosons are comparably small and only of order of a few per-cent, the situation changes if the finite parts for the counter-terms described above are not included. In that case, i.e. when using SFOLD out of the box, the corrections for the decays with a Z or W boson in the final state can be a factor of 10 larger. For the decays into a pair of fermions we also find very good agreement with only very small differences for small values of m. Similarly we show the results for the light and heavy sbottom decays in Fig. 5. Here, the results are very similar to those of the stop decays. We do not add figures for stau or τ -sneutrino decays, or the decays of first and second generation sfermions; they would look very similar to the ones for stop and sbottoms, only the overall size of the loop corrections being smaller. Thus, in total we found a very good agreement between SARAH and SFOLD for all kinds of two-body decays of sfermions.

Gluino decays in the MSSM: SARAH vs. FVSFOLD
In this section we compare the decays of gluinos in the MSSM obtained with SARAH and SPheno against the results generated with the code FVSFOLD. We also performed similar adjustments in FVSFOLD as done for SFOLD for our comparison. However, FVSFOLD already includes the finite parts of the counter-terms of the squark rotation matrices, i.e. it was not necessary to add those. Therefore without any larger adjustments, we find a very good agreement between SARAH and FVSFOLD as shown in Fig. 6 Thus, SARAH reproduces also the result of Ref. [121], namely that the one-loop corrections to gluino decays reduce the decay width by about 10%.

Heavy Higgs decays in the MSSM: SARAH vs. HFOLD
SARAH also makes predictions for the one-loop corrections of Higgs boson decays. However, it must be clearly stated that those predictions have to be interpreted with some care: the automatised calculations are not yet optimised for the calculation of Higgs boson decays, in particular for the SMlike Higgs boson. For such decays, we leave an appropriate definition of counter-terms, following our explanations in Sect. 4.2, to future work. One reason is that for consistency it will be necessary to use the counter-terms in the calculation of the mass spectrum as well. This is, however, not yet possible. We want to stress that SARAH already calculates the light Higgs into SM particle decays by adapting higher-order corrections (even beyond NLO) for the SM and MSSM from literature. Thus, the 'old' results obtained with SARAH are expected to be more accurate.
On the other hand, for the decays of heavy Higgs bosons, whose mass corrections are usually much smaller, and/or for decays into BSM states the applied NLO corrections are expected to work well, and the obtained results supersede the pure tree-level calculations often done for these decay modes. In order to validate these results, we compared them against the code HFOLD which also makes predictions for the one-loop corrections of Higgs decays in the MSSM. Here, we made the same adjustments as for FVSFOLD: on-shell masses in loops and kinematics have been turned off. In addition, we needed to turn off all improvements for the 'old' calculation in SPheno to obtain equivalent LO results to HFOLD. The results are summarised in Fig. 7 where we compare our  Fig. 7 Comparison between SARAH and HFOLD for the heavy CP-even and -odd Higgs. On the left, the loop-corrected partial widths are shown. On the right, the relative size of the loop correction is given. Blue lines are obtained with SARAH, red lines with HFOLD results for the decays of the neutral heavy Higgs states H and A 0 into SUSY particles and SM-like Higgs bosons. Without further modifications we find that the predictions of the size of the one-loop corrections of both codes agree rather well in particular in the dominant decay modes. However, we find that for decay channels with small partial widths also sizeable differences can be present. A detailed investigation of the remaining differences and also a comparison with other Higgs boson decay widths calculations is left for a dedicated work. Such a future investigation should also focus on the detailed derivation and incorporation of the U -factors, which admix the Higgs bosons beyond tree level. This is particularly crucial when comparing to codes such as NMSSMCalc or FeynHiggs, where the definition of their Z -factors is different [40,108]. We briefly discuss the relevance of Ufactors for Higgs boson decays in Sect. 5.4.

Radiative neutralino and gluino decays: SARAH vs.
SPheno SARAH does not only calculate the one-loop corrections to tree-level two-body decays, but also calculates the LO result for loop-induced decay widths. 6 The most important applica- For our comparison, we choose parameter points with a light mass splitting between: (i) a bino and wino LSP and NLSP, respectively; (ii) the gluino and all three kinds of neutralinos. For the case of the neutralino decay, the result for the obtained width as a function of the wino mass parameter M 2 as well as the relative difference between SPheno and SARAH as a function of the mass splitting are shown in Fig. 8. We show the SARAH results for three different choices of the U -factors: (i) without U -factors, (ii) using the rotation matrices obtained with the momentum being the mass of the lightest neutralino in all vertices, (iii) using p 2 -dependent U -factors. The second option corresponds to the procedure applied in SPheno and thus we find a reasonable agreement within 10%. The results without U -factors are very similar and only very close to the level crossing visible differences occur. However, when using the p 2 -dependent U -factors, the obtained width is significantly smaller. This is due to a cancellation between the vertex and wave-function corrections, which is most efficient when including the p 2 dependence in the U -factors. For the decays of the gluino into a neutralino and gluon, we find very good agreement between SPheno and SARAH for all three kinds of neutralinos, see again Fig. 8. Note that throughout the calculation of loop-induced decays loop-corrected masses are inserted.

Impact of external U -factors
Before we conclude, we want to give some impression of the numerical impact induced by the inclusion of U -factors. For this purpose, we show in Figs. 9, 10 and 11 the size of the oneloop corrections for selected stop, sbottom and Higgs decays, respectively, using the three available options to calculate the U -factors. We apply loop-corrected DR masses for all cases and particles, i.e. DECAYOPTIONS[1115]=1. We focus on two effects: First we apply U -factors both at LO and NLO equally and, in the left figures, show the relative correction induced by the NLO corrections. Second, in the right figures, we show the effect of including the U -factors in the NLO calculation compared to the NLO decay width calculation without external U -factors. More precisely, the value U shown is defined as Here, 0 are the decay widths without applying U -factors. U encodes the difference in the relative correction factor from LO to NLO when applying U -factors in contrast to not applying U -factors. It thus encodes the effect of U -factors in the one-loop correction, factoring out their effect already present at tree level. Depending on the particle species the effect at tree level can already be pronounced, and thus was already included in previous SARAH and SPheno versions. We therefore focus on the effect of the U -factors in the relative NLO correction.
For the sbottom decays into gauginos depicted in Fig. 9 the changes due to the inclusion of U -factors are moderate. From the left figures it is apparent that the size of NLO corrections is mostly independent of the inclusion of U -factors. From the right figures we deduce that the effect of U -factors on the relative NLO correction remains below 10% for all choices. The reason is that the left-right mixing in the sbottom sector is in general small and nearly identical at tree and loop level. Thus the U matrices are almost diagonal. This is different for the decays of stops shown in Fig. 10 where the left-right mixing is more pronounced. This mixing receives also a sizeable radiative correction which is encoded in the U -factors. Consequently, there is also a larger sensitivity on how this matrix is calculated and incorporated as shown in the right two figures. We find that the results without momentumdependence can differ from the other two options by 30% for the considered decays. For the heavy stop, this effect is even more pronounced. However, for the decay widtht 2 →χ 0 4 t, where the relative NLO corrections encoded in U differ by more than 100%, the absolute NLO correction almost vanishes, as can be seen from the left figure. Thus, in all examples for stop and sbottom decays, the inclusion of U -factors gives only a moderate change in the relative NLO corrections once (momentum-dependent) U -factors are taken into account compared to the calculation without U -factors. This is slightly different for the heavy Higgs decays shown in Fig. 11. As shown in the left figures all three options for the U -factors can alter the size of the relative NLO corrections significantly. In the right figures it is apparent that even for the relative NLO correction differences of 50% and more compared to the calculations without U -factors are easily possible, a fact which is well known for Higgs bosons. This shows the need to properly include these factors for Higgs boson decays even if the radiative corrections to the masses are moderate and the particles are clearly separated in their masses. Further studies for Higgs boson decays and a comparison of the U -factors to Z -factors as discussed in Refs. [40,108] are in order in future work.

Conclusions
In this paper we described a fully generic implementation of the calculation of two-body decay widths at the full oneloop level in the SARAH and SPheno framework, which can be used in a wide class of supported models. We presented the necessary generic expressions for virtual and real corrections. Wave-function corrections are determined from on-shell conditions. On the other hand, the parameters of the underlying model are by default renormalised in a DR (or MS) scheme. We described how higher-order corrections for the external states can be taken into account. We also explained how we restore gauge invariance as well as ultraviolet and infrared finiteness when setting the external masses to their loop-corrected values. We commented on the drawbacks compared to a full on-shell approach which is model and process dependent.
We have shown how the new features of SARAH and SPheno can be used and how the user can implement own counter-terms to be used for the calculation of two-body decay widths. We studied the impact and relevance of such counter-terms for two examples in the SM, namely the decay to the top-quark and the SM Higgs boson decay into bottom quarks. In addition, we compared our implementation for sfermion and gluino decays within the MSSM against other available codes, namely SFOLD, HFOLD and FVSFOLD, which also employ a DR renormalisation for the MSSM parameters. After a few described adjustments in those codes we found an overall excellent agreement. For the MSSM and R-parity violating models we also compared chargino and neutralino decays against CNNDecays, which uses a full on-shell scheme for masses and couplings and found numerically identical results.
The new extension is included in SARAH 4.11.0 and makes it possible to study radiative corrections to two-body decay modes in many different supersymmetric and nonsupersymmetric models. However, models with CP violation and/or (additional) massive gauge bosons charged under U(1) em × SU(3) c are not yet supported. This is left for future work. Other future extensions aim at necessary improvements to better handle Higgs boson decays, in particular for the decays of the SM-like Higgs boson to SM particles and the inclusion of external higher-order mass and mixing corrections. Lastly the inclusion of decays of gauge bosons is in order, but it is left for future work.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: Conventions and expressions for loop contributions
In this section we present our conventions for vertices and the generic expressions for the loop contributions to the wavefunction corrections and the vertex corrections (Figs. 12,13,14,15,16,17,18,19,20). We factor out the Lorentzdependent part of the vertices and work with the following conventions.
(a) (b) Fig. 12 Generic one-loop diagrams for the fermion self-energy For the loop corrections SARAH inserts the various particle species of the model under consideration also taking into account additional symmetry and colour factors, which are not depicted here. The various contributions are then summed up using All results in the following are expressed in terms of Passarino-Veltman integrals. The scalar loop functions A 0 , B 0 and C 0 are calculated numerically in SPheno according to the standard recipe of Ref. [124]. Tensor integrals are related to the scalar functions according to the famous techniques developed in Ref. [125]. Explicit expressions for derivatives of two-point functions,Ḃ 0 ,Ḃ 1 , which are used by SPheno are given in the appendix of Ref. [75].
For Majorana fermions, an additional overall factor 1 2 is present. 19 Generic diagrams contributing to S → SV decays

A.2.5: Scalar to scalar and gauge boson decays
A.2.6: Scalar to two gauge boson decays Following the notation introduced in Sect. 2.6, we present here the formulae employed for the calculation of real corrections (massless gauge boson emission) to the six generic decays under consideration. Common to all processes is the calculation of the group-theory factors, which we describe in the first subsection.

B.1: Group-theory factors
For real three-body decays, we start with a two-body matrix element M 0 where a particle having four-momentum p 0 decays to two others with four-momenta p 1 , p 2 , and we can attach a photon in four ways: either to the incoming or outgoing legs, or to the vertex itself. This latter case is only possible when the tree-level vertex involves two scalars and a massive gauge boson; in SARAH such vertices are stored in the form If we attach a massless gauge boson A b μ to the above term then we can have the coupling We should now understand that the indices for V, A correspond to the unbroken gauge symmetry under which A b μ transforms. By insisting on gauge invariance up to linear order we find three sets of conditions: Summation is only implied over the primed indices. The first two equations arise from requiring simple gauge invariance of the individual vertices, while the third involves the derivatives of the gauge transformation cancelling against the A a μ → A a μ + ∂ μ a + · · · part. Importantly, Eq. (B.5) completely determines c i jab .
If the gauge group is that of the SM, then the quartic coupling can only be relevant for a W boson and a photon, for which case the first two equations simply become charge conservation and Eq. (B.5) implies The above logic can also be used to compute the grouptheory factors for the Bremsstrahlung decay cross-sections, where a particle with momentum p 0 decays to final states with momenta p 1 , p 2 and a photon/gluon with momentum k.
We first compute all of the relevant processes using only the primitive vertices stripped of group-theory factors, and split up the squared amplitude depending on which leg the photon or gluon propagator has been attached. We first make a general definition of the coupling c i jk in the same was as we did in Eqs. (B.1) and (B.2) and strip off any Lorentz indices, so for fields (scalars, fermions or vectors) i in representations of the gauge group R 0 → R 1 , R 2 we have the coupling Then the matrix element is where the functions f μ i ( p i ) may have spinor or other Lorentz indices, and include the wave-function factors, as appropriate to the final states; and the index i labels the leg to which the photon/gluon is attached: 0 for incoming, 1, 2 for the outgoing legs and 3 to denote the vertex in the case of scalar to scalar plus gauge boson decays. We have explicitly included the intermediate propagators to show the potential infrared divergences-and also because they allow easy identification of the appropriate diagram.
Then defining x 0 ≡ −2 p 0 · k, x 1 ≡ 2 p 1 · k, x 2 ≡ 2 p 2 · k, x 3 ≡ 1, a generic squared matrix element can be expressed as Now all of the group-theory factors are encoded inC i j . However, whileC ji =C * i j , it is also clear from the hermiticity of the generators that theC i j are real-for example using d(R 0 ) as the dimension of representation R 0 and therefore we can write Using the gauge invariance of the Lagrangian terms, we obtain for i < 3 where for S → SV decays we define the heavy gauge boson as particle 2. These are almost enough to completely determine the colour factor of the amplitude in terms of the twobody decay factor C, because we can writẽ (B.14) These conditions are sufficient to determine all of the group factors in terms of C, the quadratic Casimir C 2 (R i ), and one remaining colour factor. However, we can also use the gauge invariance of the Lagrangian term to second order to obtain 0 =C 00 +C 11 +C 22 − 2C 01 − 2C 02 Hence all of the group-theory factors are proportional to C, and we can therefore define giving and we can write for a generic process In the case of a U(1) gauge boson (the photon), we have C i j = Q i Q j and we define Q 3 ≡ Q 0 +Q 1 . This also follows from the above expressions, for example The real corrections for the fermionic decays F → F S are given by The coupling c g is the electromagnetic coupling e or the strong coupling g s for photon and gluon emission, respectively. c L and c R are the left-and right-handed coupling of the tree-level vertex of F → F S. We identify F in = X , F out = Y 1 and S = Y 2 with the indices 0, 1 and 2, respectively. The individual contributions AB i j are given by We also implemented real corrections for which Y 1 is taken to be massless, if Y 1 is charge and colour neutral.

B.3: S → F F
The real corrections for the scalar decays S → F F are given by (B.33) We follow the notation introduced in the previous subsection. However, we now identify S = X , F = Y 1 and F = Y 2 with the indices 0, 1 and 2, respectively. The individual contributions AB i j are given by We also implemented real corrections where one of the finalstate fermions can be massless. Again this fermion has to be charge and colour neutral.

B.4: F → F V
The real corrections for the fermionic decays F → F V are given by (B.46) The coupling c g is the electromagnetic coupling e or the strong coupling g s for photon and gluon emission, respec-tively. c L and c R are the left-and right-handed coupling of the tree-level vertex of F → F V . The final-state gauge boson Y 2 is fixed to one of the two heavy gauge bosons, Z or W . Clearly only the photon couples to the W boson, such that we subsequently present results independently for F → F W + γ on the one side and F → F Z + γ /g and F → F W + g on the other side. For gluon emission C 00 = C 11 = C 01 = C 2 (R F ), C i2 = 0; for photon emission C i j = Q i Q j , where we again identify F in = X and F out = Y 1 with indices 0 and 1, respectively. For F → F W + γ the individual contributions are given by The coupling c g is the electromagnetic coupling e or the strong coupling g s for photon and gluon emission, respectively. In contrast to previous decays, S → SV +γ /g also has a contribution from a four-point interaction, which accordingly does not contribute to the infrared divergent part of the real corrections, but yields a finite contribution. The sum over the gauge factors thus includes four indices 0, . . . , 3.
Since we assume a non-extended gauge sector in this work, the four-point vertex is only present for photon emission, so we can simplify the calculation of the group-theory factors to C i j = Q i Q j for i, j ∈ {0, 1, 2}, where the charges of the incoming and outgoing scalars are Q 0 and Q 1 and the charge of the outgoing gauge boson is Q 2 ; the factors involving the four-point interaction are then For the emission of a gluon we can use Q 2 = 0 and Q 0 = Q 1 = 1. However, we stress that the results below are true in general for any gauge groups; for extended gauge sectors we would just need to employ Eq. (B.17). The individual contributions i j are given by we are working with real scalars, the t G i j are antisymmetric and real (to translate to complex fields we require roughly t G i j → i T G i j ). Then we find that the Goldstones G are defined by 3) The vectors should be chosen to be orthogonal with normal- Looking at the vector-mass term we find for the scalars S 0 i having vacuum expectation values v i and we see that where m G V is the mass of the vector, and we define the SSV coupling when we diagonalise the scalars to mass eigenstates through S 0 i ≡ R ik S k (where R G j = N G α G j ): For the SGV coupling and SV V couplings, we read off This is the relationship that we enforce between the on-shell vertices to ensure that infrared divergences are cancelled. We can also read off the GV V coupling For any unbroken gauge groups, the Goldstones must transform as and we therefore identify the α G i t a i j α G j factor in Eq. (D.8) with T a GG to obtain For the photon, this just becomes the familiar vertex L ⊃ −em W G + γ a μ W − μ + h.c.. (D.11) Now the V V γ coupling will be given just in terms of the gauge coupling of the unbroken gauge group, so the electromagnetic coupling e here; indeed from decomposing the kinetic terms of the gauge bosons we will find Of these, the first two terms vanish in the limit where the massless gauge boson is soft, while we identify the last term as the conventional gauge current and c aGG = − gt a GG .
(D. 13) This leads to the relation between the γ GV and γ V V vertices of simply a factor of m G V .