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

TL;DR

This paper presents a comprehensive framework for calculating two-body partial decay widths at the full one-loop level across various models, incorporating virtual, real, and higher-order corrections with IR divergence cancellation.

Contribution

It introduces a fully generic implementation in SARAH and SPheno for one-loop decay width calculations, supporting diverse decay channels and model-specific renormalization schemes.

Findings

Validated against existing codes for MSSM decay processes.

Reproduced results for R-parity violating neutralino decays.

Compared and contrasted with other public decay codes.

Abstract

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 wavefunction corrections are determined from on-shell conditions, the parameters of the underlying model are by default renormalised in a DR (or 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 infra-red divergences for such cases, which is achieved through an infra-red 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. Radiative 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.

Tables

Table 1. Table 1: Partial decay width t → W ​ b → 𝑡 𝑊 𝑏 t\to Wb in different schemes, see text for details.

Scheme	${\Gamma }_{t\to Wb}^{\text{LO}}$ [GeV]	${\Gamma }_{t\to Wb}^{\text{NLO}\text{,EW}}$ [GeV]	${\Gamma }_{t\to Wb}^{\text{NLO}\text{,EW+}\text{QCD}}$ [GeV]
$\left(1\right)$, $\alpha \left(0\right)$, $Q=173$ GeV	$1.443$	$1.487$ $\left[+3.0\%\right]$	$1.352$ $\left[-0.135\right]\left[-9.1\%\right]$
$\left(1\right)$, $\alpha \left(m_Z\right)$, $Q=173$ GeV	$1.546$	$1.596$ $\left[+3.2\%\right]$	$1.452$ $\left[-0.144\right]\left[-9.0\%\right]$
$\left(2\right)$, $\alpha \left(0\right)$	$1.443$	$1.519$ $\left[+5.3\%\right]$	$1.384$ $\left[-0.135\right]\left[-8.9\%\right]$
$\left(3\right)$, $\alpha \left(m_Z\right)$	$1.546$	$1.522$ $\left[-1.6\%\right]$	$1.378$ $\left[-0.144\right]\left[-9.5\%\right]$

Table 2. Table 2: Partial decay width H → b ​ b ¯ → 𝐻 𝑏 ¯ 𝑏 H\to b\bar{b} for different values of Q 𝑄 Q , see text for details.

Scale	${\Gamma }_{H\to b\overline{b}}^{\text{LO}}$ [MeV]	${\Gamma }_{H\to b\overline{b}}^{\text{NLO}\text{,EW}}$ [MeV]	${\Gamma }_{H\to b\overline{b}}^{\text{NLO}\text{,EW+}\text{QCD}}$ [MeV]
$Q=125$ GeV	$1.959$	$1.972$ $\left[+0.6\%\right]$	$2.376$ $\left[+20.5\%\right]$
$Q=62.5$ GeV	$2.188$	$2.215$ $\left[+1.5\%\right]$	$2.473$ $\left[+11.6\%\right]$
$Q=250$ GeV	$1.778$	$1.783$ $\left[+0.3\%\right]$	$2.280$ $\left[+27.8\%\right]$

Equations

$${\Gamma }_{X\to Y_1{Y}_2}=\frac{1}{16\pi m_X^3}\lambda (m_X^2,m_{Y_1}^2,m_{Y_2}^2)C_S{C}_C\sum _{h,p}|\mathcal{M}{|}^2,$$

$$\lambda (p_0^2,p_1^2,p_2^2)=\sqrt{p_0^4+p_1^4+p_2^4-2p_0^2{p}_1^2-2p_1^2{p}_2^2-2p_0^2{p}_2^2}.$$

$$F\to FS:$$

$$\mathcal{M}$$

$$\sum _{h,p}|\mathcal{M}{|}^2$$

$$S\to FF:$$

$$\mathcal{M}$$

$$\sum _{h,p}|\mathcal{M}{|}^2$$

$$S\to SV:$$

$$\mathcal{M}$$

$$\sum _{h,p}|\mathcal{M}{|}^2$$

$$S\to SS:$$

$$\mathcal{M}$$

$$\sum _{h,p}|\mathcal{M}{|}^2$$

$$\mathcal{M}\equiv$$

$$\sum _{h,p}|\mathcal{M}{|}^2=$$

$$+\frac{1}{8m_{Y_2}^2{m}_{Y_3}^2}{\left[m_X^4+(m_{Y_1}^2-m_{Y_2}^2{)}^2-2m_X^2(m_{Y_2}^2+m_{Y_3}^2)\right]}^2{M}_2{M}_2^{\ast }$$

$$+\frac{1}{4m_{Y_2}^2{m}_{Y_3}^2}[m_X^6-3m_X^4(m_{Y_1}^2+m_{Y_2}^2)-(m_{Y_1}^2-m_{Y_2}^2{)}^2(m_{Y_1}^2+m_{Y_2}^2)$$

$$+m_{Y_1}^2(3m_{Y_1}^4+2m_{Y_1}^2{m}_{Y_2}^2+3m_{Y_2}^4)](M_1{M}_2^{\ast }+M_2{M}_1^{\ast })$$

$$\displaystyle\mathcal{M}\equiv\epsilon_{\mu}^{*}(p_{2})\bigg{(}M_{1}\bar{v}(p_{0})\gamma^{\mu}P_{L}v(p_{1})+M_{2}\bar{v}(p_{0})\gamma^{\mu}P_{R}v(p_{1})$$

$$\displaystyle\qquad\qquad+M_{3}p_{0}^{\mu}\bar{v}(p_{0})P_{L}v(p_{1})+M_{4}p_{0}^{\mu}\bar{v}(p_{0})P_{R}v(p_{1})\bigg{)}$$

$$\sum _{h,p}|\mathcal{M}{|}^2=\frac{1}{2m_{Y_2}^2}\left[m_X^4+m_{Y_1}^4+m_{Y_1}^2{m}_{Y_2}^2-2m_{Y_2}^4+m_X^2(-2m_{Y_1}^2+m_{Y_2}^2)\right](M_1{M}_1^{\ast }+M_2{M}_2^{\ast })$$

$$+\frac{1}{8m_{Y_2}^2}\left[(m_X^2+m_{Y_1}^2-m_{Y_2}^2)(m_X^4+(m_{Y_1}^2-m_{Y_2}^2{)}^2-2m_X^2(m_{Y_1}^2+m_{Y_2}^2))\right](M_3{M}_3^{\ast }+M_4{M}_4^{\ast })$$

$$-3m_X{m}_{Y_1}(M_1{M}_2^{\ast }+M_2{M}_1^{\ast })$$

$$-\frac{1}{4m_{Y_2}^2}\left[m_{Y_1}(m_X^4+(m_{Y_1}^2-m_{Y_2}^2{)}^2-2m_X^2(m_{Y_1}^2+m_{Y_2}^2))\right](M_1{M}_3^{\ast }+M_3{M}_1^{\ast }+M_2{M}_4^{\ast }+M_4{M}_2^{\ast })$$

$$-\frac{1}{4m_{Y_2}^2}\left[m_X(m_X^4+(m_{Y_1}^2-m_{Y_2}^2{)}^2-2m_X^2(m_{Y_1}^2+m_{Y_2}^2))\right](M_1{M}_4^{\ast }+M_4{M}_1^{\ast }+M_2{M}_3^{\ast }+M_3{M}_2^{\ast })$$

$$+\frac{1}{4m_{Y_2}^2}\left[m_X{m}_{Y_1}(m_X^4+(m_{Y_1}^2-m_{Y_2}^2{)}^2-2m_X^2(m_{Y_1}^2+m_{Y_2}^2))\right](M_3{M}_4^{\ast }+M_4{M}_3^{\ast })$$

$$F\to FV:M_1^T=i{c}_R,M_2^T=i{c}_L,$$

$$S\to FF:M_1^T=-i{c}_R,M_2^T=-i{c}_L,$$

$$S\to SV:M^T=-2ic,$$

$$M_i=M_i^T+2M_i^V+2M_i^W.$$

$$M_{ijk}^W=i\left(\delta c_{ijk}+\frac{1}{2}{c}_{ljk}\delta Z_{X_l{X}_i}+\frac{1}{2}{c}_{ilk}\delta Z_{Y_{1l}{Y}_{1j}}+\frac{1}{2}{c}_{ijl}\delta Z_{Y_{2l}{Y}_{2k}}\right)$$

$$V_i^{\mu ,0}$$

$$S_i^0$$

$$F_i^{L0}$$

$$F_i^{R0}$$

$${\Gamma }_{S_i{S}_j}(p^2)$$

$${\Gamma }_{V_i{V}_j}^{\mu \nu }(p^2)$$

$${\Gamma }_{F_i{F}_j}(p)$$

$$\dot{\Pi }(k^2)={\frac{\partial }{\partial p^2}\Pi (p^2)|}_{p^2=k^2}\text{and}\dot{\Sigma }(k^2)={\frac{\partial }{\partial p^2}\Sigma (p^2)|}_{p^2=k^2}$$

$$\delta Z_{S_i{S}_i}=$$

$$\delta Z_{S_i{S}_j}=$$

$$\delta Z_{F_i{F}_i}^L=$$

$$\delta Z_{F_i{F}_i}^R=$$

$$\delta Z_{F_i{F}_j}^L=$$

$$\delta Z_{F_i{F}_j}^R=$$

$$\delta R_{ij}^{\Phi }=\frac{1}{4}\sum _k\left(\delta Z_{{\Phi }_i{\Phi }_k}-(\delta Z_{{\Phi }_k{\Phi }_i}{)}^{\ast }\right)R_{kj}^{\Phi }.$$

Full text

1. 1 Introduction

2. 2 Calculation of decay widths at the full one-loop level

3. 2.1 Generic unpolarised squared matrix elements

4. 2.2 Tree-level amplitudes

5. 2.3 One-loop amplitudes

6. 2.4 Vertex corrections

7. 2.5 Wavefunction corrections

8. 2.6 Real corrections

9. 2.7 Combination of results

10. 2.8 Current limitations

11. 3 Higher-order corrections to the external states

12. 3.1 Ultraviolet and infra-red divergences

13. 3.2 Mixing of species

14. 3.3 Loop-induced decays

15. 4 Implementation in SARAH

16. 4.1 SARAH–SPheno interface

17. 4.2 Definition of counter-terms

18. 4.3 Options for the evaluation with SPheno

19. 4.4 Output of SPheno

20. 5 Numerical results

21. 5.1 Renormalisation of $\alpha$ and the top-quark width

22. 5.2 Renormalisation of Yukawa couplings and fermionic Higgs decays

23. 5.3 Comparison with other codes

24. 5.4 Impact of external $U$-factors

25. 6 Conclusions

26. A Conventions and expressions for loop contributions

27. A.1 Wavefunction corrections

28. A.2 Vertex corrections

29. B Expressions for real corrections

30. B.1 Group theory factors

31. B.2 $F\to FS$

32. B.3 $S\to FF$

33. B.4 $F\to FV$

34. B.5 $S\to VV$

35. B.6 $S\to SS$

36. B.7 $S\to SV$

37. C Infra-red divergent parts of Passarino-Veltman integrals

38. D Goldstone boson vertices

1 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) [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39] 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. [54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74]; and neutralino and chargino decays in the NMSSM [75, 76]. For other SUSY models with $R$-parity violation and $\mathcal{CP}$ violation, only a few selected decay modes were discussed so far in Refs. [77, 78]. The available codes to study decays at the one-loop 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 two-loop 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 infra-red divergences.

While wavefunction corrections are determined from on-shell conditions, the default settings use a $\overline{\text{DR}}$ (or $\overline{\text{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 infra-red divergences, for which we add an infra-red counter-term making use of modified Goldstone boson vertices. The setup also supports loop-induced decays. An extension to models with $\mathcal{CP}$ violation or additional charged and massive, coloured vector particles is left for future work. A more thorough discussion of Higgs boson decays, which are very sensitive to corrections of the external states, will be addressed in the future.

The paper is organised as follows: In Section 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 Section 3. In Section 4 we explain how the new features of SARAH and SPheno can be used. In Section 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 Section 6. The appendix contains all relevant generic expressions for virtual and real corrections as well as a derivation of the employed Goldstone boson vertices.

2 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\to Y_{1}Y_{2}$. Our implementation can handle the decays $S\to SS$, $S\to SV$, $S\to VV$, $S\to FF$, $F\to FS$ and $F\to FV$, 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, infra-red 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 infra-red 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 wavefunction corrections as well as a discussion of counterterms. Then in Section 2.6 we discuss the relevant real corrections being $1\to 3$ processes, before we combine the results in Section 2.7. Finally, we list the limitations of our implementation in Section 2.8.

2.1 Generic unpolarised squared matrix elements

For any two-body decay we write its partial width in the form

[TABLE]

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}=\frac{1}{2}$, if $\overline{Y}_{1}=Y_{1}$ and $\overline{Y}_{2}=Y_{2}$. For $X=S$ it is $C_{S}=\frac{1}{2}$, if $\overline{Y}_{1}=Y_{1}$ and $\overline{Y}_{2}=Y_{2}$ and $Y_{1}=Y_{2}$. Therein $\overline{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}=\frac{1}{2}$, while for more complicated colour configurations $C_{C}$ can be easily 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 $\lambda$ is given by

[TABLE]

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:

[TABLE]

For $S\to VV$ we split the squared amplitude as follows:

[TABLE]

Lastly the squared amplitude for $F\to FV$ is given by:

[TABLE]

We implemented special cases for final states with vanishing masses, which are not given here.111In the calculation of one-loop decays we introduced a minimal allowed mass for fermions and scalars of $10^{-15}$ GeV. Smaller masses are set to zero to stabilise numerics, see Section 5.3.1 for a discussion in the context of $R$-parity violation.

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 (left- and right-handed) couplings as follows:

[TABLE]

The conventions for the parametrisation of the vertices are summarised in Appendix A. For $F\to FV$ $M^{T}_{3}$ and $M^{T}_{4}$ and for $S\to VV$ $M^{T}_{2}$ vanish at tree-level, but contributions are generated at the one-loop level.

2.3 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 wavefunction 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

[TABLE]

For an exact next-to-leading order calculation the complex-conjugated part of this amplitude $M_{i}^{*}$ is inserted in the complex-conjugated matrix elements of Eq. (2.1) to Eq. (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 Eq. (2.1) to Eq. (2.8). When squaring the amplitude one needs to be careful if external coloured particles are involved. For these cases, SARAH calculates individual colour factors for tree- and loop-level contributions and sums them up congruently. For loop-induced decays the sum of the amplitudes of the vertex corrections $M^{V}_{i}$ and wavefunction corrections $M^{W}_{i}$ is inserted into all occurrences of matrix elements in Eq. (2.1) to Eq. (2.8).

2.4 Vertex corrections

In general, there are four different topologies contributing to the vertex corrections $M^{V}$ which are shown in Fig. 1. For decays involving fermions only the first topology is of relevance. Depending on the considered decay, different generic diagrams are associated with these topologies. They are depicted for the different decays in Appendix A in Figs. 15–20. The results are a function of internal masses $m_{1},m_{2}$ and $m_{3}$ entering the loop diagrams, but also of the external momenta $p_{0},p_{1}$ and $p_{2}$. In this Section 2 their squared values correspond to the squared external masses of the particles, i.e. $m_{X}^{2}$, $m_{Y_{1}}^{2}$ and $m_{Y_{2}}^{2}$, respectively. The calculation of the generic amplitudes $M^{(k)}_{i}$ for each diagram is straightforward and all results are given in Appendix A. They are obtained with FeynArts [93] and FormCalc [94] in Feynman-’t Hooft gauge, i.e. charged and neutral Goldstone bosons are included in the calculation and cancel the unphysical contributions from heavy gauge bosons.

Our results are expressed in terms of Passarino-Veltman integrals obtained through dimensional reduction ($\overline{\text{DR}}$).222FormCalc works with constrained differential renormalisation [95], which equals dimensional reduction at the one-loop level [94]. Thus, the ultraviolet divergences can be split off in terms of $\Delta=\tfrac{1}{\epsilon}-\gamma_{E}+\ln(4\pi)$, where $\epsilon$ regularises the divergence and equals the difference to four dimensions $d=4-2\epsilon$ and $\gamma_{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 ($\overline{\text{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. Infra-red 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 infra-red divergences will be addressed in Section 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})^{*}$.

2.5 Wavefunction 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}_{ijk}$ for a vertex of the form $c_{ijk}X_{i}Y_{1j}Y_{2k}$ for the process $X_{i}\to Y_{1j}Y_{2k}$ yields

[TABLE]

with the counter-term $\delta c$ of the tree-level coupling $c$ and the wavefunction corrections $\delta 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 wavefunction 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 wavefunction corrections are generally matrices.

[TABLE]

In order to determine the wavefunction corrections $\delta 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:

[TABLE]

The renormalised self-energies are indicated through $\hat{\Pi}$ and $\hat{\Sigma}$ compared to the unrenormalised ones $\Pi$ and $\Sigma$, which are of relevance for the subsequent discussion. $\Pi_{VV}$ and $\Pi_{SS}$ are the self-energies of the gauge bosons and scalars, respectively. For the gauge bosons we are only interested in the transverse part $\Pi_{VV}$. 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 $\Sigma^{L}$, $\Sigma^{R}$, $\Sigma^{SL}$ and $\Sigma^{SR}$.

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–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 $\dot{\Pi}$ and $\dot{\Sigma}$, are defined as follows:

[TABLE]

Demanding on-shell conditions for the external states fixes the wavefunction corrections. Their derivation can for example be found in Refs. [96, 75] and results in similar expressions for scalars and gauge bosons:

[TABLE]

For the fermions we need to distinguish four cases:

[TABLE]

By $\widetilde{{\rm Re}}$ we indicate that $\Pi$ and $\Sigma$ entering $\delta Z$ include only the real parts of the loop functions, whereas couplings enter with real and imaginary components. In case of $\mathcal{CP}$ violation the definition of wavefunction corrections usually distinguishes between in- and outgoing particles in order to correctly multiply absorbative 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 wavefunction 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 counter-terms enter the counter-term of the tree-level coupling. For any particle species $\Phi$ in mass eigenstates, which is obtained from gauge eigenstates $\Phi^{\prime}$ through $\Phi_{i}=R^{\Phi}_{ij}\Phi^{\prime}_{j}$, the counter-term is given by

[TABLE]

For fermions left- and right-handed states are rotated with two matrices, such that two counter-terms employing left- and right-handed wavefunction 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 wavefunction 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 wavefunction 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, that enter the tree-level coupling counter-term of the two-body decay under consideration. The counter-terms are usually chosen depending on the model and process. However, a simple $\overline{\text{DR}}$ (or $\overline{\text{MS}}$) prescription for the renormalisation of the parameters of the underlying theory is always easily applicable: From the $\beta$ 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 counter-terms of the parameters of the underlying theory to be just proportional to the pure ultraviolet divergence only. We will refer to this scheme as $\overline{\text{DR}}$ (or $\overline{\text{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 Section 4.2 for a more detailed discussion. We also add an example of a proper renormalisation of the electric charge in Section 5.

A consequence of the application of the $\overline{\text{DR}}$ (or $\overline{\text{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 $\sqrt{m_{\tilde{t}_{1}}m_{\tilde{t}_{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 Section 4.2. A common choice for the renormalisation scale is also $Q\sim m_{X}$ close to the mass of the decaying particle $X$. We refrain from making it the default option, since $Q\sim 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 Section 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 Section 5. In order to achieve a renormalisation-scale 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 Section 4.2.

Until now we ignored the fact that particles receive higher-order mass corrections.333Also on-shell schemes can lead to finite shifts at higher orders, see e.g. Refs. [105, 106, 65] for examples in the neutralino and chargino sector. By construction we have to employ the mass values at lowest order throughout the calculation. We will discuss in Section 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 $\overline{\text{DR}}$ ($\overline{\text{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\to Y_{1}Y_{2}$, which in the most general case however is scale dependent. As mentioned, the cancellation of infra-red divergences is addressed in the next section.

2.6 Real corrections

In the previous calculation of vertex and wavefunction corrections we regularised infra-red 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 wavefunction corrections. With the help of FeynArts [93] and FormCalc [94] we generated generic results for the emission of one additional photon $\gamma$ or gluon $g$ for the previously discussed two-body processes $S\to SS$, $S\to SV$, $S\to VV$, $S\to FF$, $F\to FS$ and $F\to FV$. We denote the real corrections for $X\to Y_{1}Y_{2}+\gamma/g$ with

[TABLE]

where $k$ denotes the momentum of the photon or gluon and momenta with upper index [math] equal the zeroth component of the corresponding four vector. External momenta are set to $p_{0}^{2}=m_{X}^{2}$, $p_{1}^{2}=m_{Y_{1}}^{2}$, $p_{2}^{2}=m_{Y_{2}}^{2}$ and $k^{2}=0$. We then have $C^{\prime}_{S}=C_{S}$ for $S\to VV$ and $S\to SV$, otherwise $C^{\prime}_{S}=\tfrac{1}{2}C_{S}$. The charge and colour structure is encoded in matrices $C_{ij}$ explained in Appendix B, which is why Eq. (2.30) 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

[TABLE]

where $p_{i,j}\in\{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 infra-red 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 three-body decay $S\to SV+\gamma/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\to 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\to F^{\prime}S\gamma$ and $S\to F^{\prime}F\gamma$, where one final state fermion $F^{\prime}$ 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 shouldn’t cause problems in practical applications unless charged or coloured states with very smalles masses ($\ll$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\to SV+\gamma/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.

2.7 Combination of results

The partial width at next-to-leading order is thus obtained as follows

[TABLE]

where $\Gamma_{X\to Y_{1}Y_{2}}$ is obtained from Eq. (2.1) with the squared amplitudes given in Eq. (2.1) to Eq. (2.8). The individual parts $M_{i}$ are taken from Eq. (2.12) and Eq. (2.9) for the complex conjugated and non-complex conjugated squared amplitudes, respectively. $\Gamma_{X\to Y_{1}Y_{2}+\gamma/g}$ are the real corrections from Eq. (2.30), which are calculated individually for photons and gluons and summed up. For loop-induced decays the virtual contributions in $M^{V}$ are by definition ultraviolet finite, still we also include the described wavefunction corrections. We note that ultraviolet finiteness can be checked through a variation of $\Delta$ defined in Section 2.4 and infra-red finiteness through a variation of the regulator mass for the photon and/or gluon, see Section 4 for a description how to access them.

2.8 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 $\mathcal{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 $\overline{Y}_{1}\overline{Y}_{2}$ is assumed to have the same partial width. Also note that in case of $\mathcal{CP}$ violation a common approach is to define extended wavefunction 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 infra-red 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^{\prime}$ 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 be easily extended, which we leave for future work.

3 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 higher-order 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 the gauge boson sector we recommend to use on-shell values for gauge bosons, see the discussion in Section 4. 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}_{ij}S^{\prime}_{j}$ at tree-level and mix at higher orders: following Eq. (2.18) the mass matrices beyond tree-level take the form

[TABLE]

with the tree-level masses $m_{i}$. Let’s suppose that for the calculation of the external masses a $\overline{\text{MS}}$ or $\overline{\text{DR}}$ scheme is preferred, i.e. the self-energies $\hat{\Pi}$ and $\hat{\Sigma}$ are renormalised such that only the corresponding ultraviolet divergent part is omitted. In a first approximation we set $p^{2}=0$ and diagonalize the obtained mass matrix $\mathcal{M}(0)$ through a unitary $(n\times n)$ matrix $U^{0}$. The tree-level mass eigenstates $S_{i}$ are thus rotated into states $\tilde{S}_{i}=U^{0}_{ij}S_{j}$ with masses $\widetilde{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].444It differs by the prescription how the self-energies $\Pi$ are renormalised. It is used to rotate the tree-level, vertex and counter-term corrections uniformly by applying it at the amplitude level. For the decay $\tilde{S}_{i}\to\tilde{S}_{j}\tilde{S}_{k}$ we for example shift the amplitudes by

[TABLE]

for $M=M^{T},M^{V}$ and $M^{W}$. Also we define rotated tree-level couplings $\widetilde{c}_{ijk}$ 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 infra-red 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_{i}^{2}$. This results in the mass eigenvalue $\widetilde{m}_{i}$, which is used to repeat the procedure iteratively with $p^{2}=\widetilde{m}_{i}$ until the mass determination stabilizes. 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 eigenstate $\tilde{S}_{i}$ separately and the matrix $U^{p}$ is determined row by row and is thus not unitary any more, as it is also well-known from the general form of $Z$-factors. Lastly a possible choice is $p^{2}=m_{1}^{2}$, 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 $\widetilde{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.

3.1 Ultraviolet and infra-red divergences

The application of external masses $\widetilde{m}$ different from the tree-level values $m$ and mixing matrices $U$ in addition to tree-level mixing induces a problem with the cancellation of ultraviolet and infra-red 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 infra-red problem is more demanding. In order to achieve the cancellation of infra-red divergences we define infra-red 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 counter-terms are used to shift the the wavefunction and vertex corrections:

[TABLE]

The aim is to cancel the infra-red divergences stemming from $2M^{T}(M^{V}+\delta M^{V}+M^{W}+\delta M^{W})^{*}$ against the ones from the real emission calculated with loop-corrected masses. The counter-terms $\delta M^{V}$, $\delta M^{W}$ are defined to be the difference of the infra-red divergences of our default scheme and the one with loop-corrected masses $\widetilde{m}$

[TABLE]

where $\text{IR}(M)$ takes only the infra-red divergent part of the amplitude $M$. The definition of such counter-terms in supersymmetric models is a common strategy: Ref. [66] introduces an infra-red counter-term for the decays $\tilde{t}_{1}\rightarrow\tilde{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. Ref. [109] discusses the introduction of such counter-terms for the heavy Higgs boson decay $H\to W^{+}W^{-}$. Ref. [110] discusses these aspects in detail in the context of a generalised narrow-width approximation, where also infra-red 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 infra-red 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:

We calculate the virtual corrections using tree-level masses. 2. 2.

We extract the infra-red divergences of all two and three-point function using the results given in Appendix C. These result are used to obtain $\text{IR}\left(M^{V,W}\right)$. 3. 3.

We use loop-corrected masses $\widetilde{m}$ throughout all infra-red divergent diagrams for the external legs and the particles in the loop. We take again the infra-red divergent parts of these amplitudes to obtain $\text{IR}\left(\widetilde{M}^{V,W}\right)$. 4. 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. 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 tree-level, wavefunction and virtual corrections according to Eq. (3.2). Instead for the contribution of the infra-red counter-term we use rotated tree-level couplings $\widetilde{c}$ rather than rotated amplitudes $\widetilde{M}$. Those rotated couplings also enter the calculation of the real corrections. In this context we note that by construction the infra-red counter-term always contains exactly one occurrence of the coupling $c$ of the tree-level two-body decay.

These steps give for most cases infra-red finite results. However, there is one complication: if the infra-red 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 infra-red divergences to cancel. If we were to apply the above procedure then the infra-red 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_{V}^{G}$ and massless ones by $\gamma^{a}$, real scalars as $S_{i}$ with masses $m_{i}$ and Dirac fermions as $F_{I}$ with masses $m_{I}$, the relevant couplings are

[TABLE]

The couplings $c^{a}_{ij},c_{GG^{\prime}}^{a},c_{IJ}^{a,L/R}$ 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:

[TABLE]

The implementation in SARAH currently assumes that the gauge sector is that of the SM; so there are no infra-red 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.

[TABLE]

Note, that in Eq. (3.11) we explicitly assume no $\mathcal{CP}$ violation.

Employing the outlined procedure we obtain partial decay widths at next-to-leading order with full cancellation of ultraviolet and infra-red divergences. Though the application of loop-corrected masses in the infra-red 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 on-shell, such that the cancellation of a gauge dependence in the real corrections among internal gauge bosons and Goldstone bosons is always guaranteed.

3.2 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 ($\mathcal{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. Slavnov-Taylor identities then ensure that the sum of the $Z$ and Goldstone contributions give zero, see also Refs. [111, 46]. 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 be easily added.

3.3 Loop-induced decays

We finish with a remark about loop-induced decays like $F_{i}\to F_{j}\gamma$. Since infra-red 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.

4 Implementation in SARAH

4.1 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:

⬇

http://sarah.hepforge.org/

Then the user should copy the tar-file into a directory called $PATH in the following and extract it through:

⬇

tar -xf SARAH-4.11.0.tar.gz

Afterwards, start Mathematica, load SARAH, run a model $MODEL and generate a SPheno version through

⬇

<< $PATH/SARAH-4.11.0/SARAH.m;

Start["$MODEL"];

MakeSPheno[];

The last command initialises all necessary calculations and writes all Fortran files into the output directory of the considered model. These files can be compiled together with SPheno version 3.3.8 or later. SPheno is also available at HepForge:

⬇

http://spheno.hepforge.org/

The necessary steps to compile the new files are:

⬇

tar -xf SPheno-4.0.2.tar.gz

cd SPheno-4.0.2

mkdir $MODEL

cp -r $PATH/SARAH-4.11.0/Output/MODEL/EWSB/SPheno/* MODEL

make Model=$MODEL

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:

⬇ ./bin/SPheno$MODEL$MODEL/LesHouches.in.$MODEL

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 one-loop two-body decays as well as the compilation time of SPheno are extended by these new routines. Therefore, in case that the user in not interested in the loop-corrected two-body decays, they can be turned off via:

⬇

MakeSPheno[IncludeLoopDecays->False];

They can be permanently turned off for a given model by adding

⬇

1SA‘AddOneLoopDecay = 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 $\beta$-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 time-consuming. In this case the option

⬇

1SA‘NoRGEsforDecays=True;

skips the RGE calculation. Of course, the verification of the cancellation of ultraviolet divergences will not be performed with this setting.

4.2 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:

⬇

1RenConditionsDecays={

2{dCosTW, 1/2*Cos[ThetaW] * (PiVWm/(MVWM^2) - PiVZ/(mVZ^2))},

3{dSinTW, -dCosTW/Tan[ThetaW]},

4{dg2, 1/2g2(derPiVPheavy0 + PiVPlightMZ/MVZ^2

5 - 2dSinTW/Sin[ThetaW] + (2PiVZVP*Tan[ThetaW])/MVZ^2)},

6{dg1, dg2Tan[ThetaW]+g2dSinTW/Cos[ThetaW]

7 - dCosTWg2Tan[ThetaW]/Cos[ThetaW]}

8};

We give an example for the application of the above electroweak counter-terms and their derivation in Section 5.1. If RenConditionsDecays is not defined, a pure $\overline{\text{MS}}$/$\overline{\text{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 counter-terms:

• •

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.

• •

Self-energies and their derivatives mixing vector bosons: those are called PiXY respectively derPiXY where X and Y are the names of the particles in SARAH and $p^{2}=m_{Y}^{2}$.

• •

Special self-energies for vector bosons containing only light/heavy states:

• –

PiVPlight0/derPiVPlight0: only light degrees of freedom are included in the loops; external momentum $p^{2}=0$.

• –

PiVPlightMZ/derPiVPlightMZ: only light degrees of freedom are included in the loops; external momentum $p^{2}=m_{Z}^{2}$.

• –

PiVPheavy0/derPiVPheavy0: only heavy degrees of freedom are included in the loops; external momentum $p^{2}=0$.

• –

PiVPheavyMZ/ derPiVPheavyMZ: only heavy degrees of freedom are included in the loops; external momentum $p^{2}=m_{Z}^{2}$.

• •

The different parts of the fermion self-energies and their derivatives: those are called SigmaLX, SigmaRX, SigmaSLX, SigmaSR respectively DerSigmaLX, DerSigmaRX, DerSigmaSLX, DerSigmaSR, 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.

[TABLE]

When defining the counter-terms, commands for matrix or tensor operators should already be 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 emphasize 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. [105, 106, 65, 75], can not yet be incorporated.

4.3 Options for the evaluation with SPheno

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

⬇

1Block DECAYOPTIONS #

2...

31001 ... # One-loop decays of particle X

41002 ... # One-loop decays of particle Y

5...

61114 ... # U-factors (0: off, 1:p2_i=m2_i, 2:p2=0, p3:p2_i=m2_1)

71115 ... # Use loop-corrected masses for external states

81116 ... # OS values for W,Z and fermions

9 (0: off, 1:g1,g2,v 2:g1,g2,v,Y_i)

101117 ... # Use defined counter-terms

111118 ... # 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.

• •

DECAYOPTIONS[1114]: this defines the choice for the external $U$-factors:

• –

0: no $U$-factors are applied.

• –

1: the $U$-factors including the full $p^{2}$ dependence are used ($U^{p}$).

• –

2: the $U$-factors calculated for $p^{2}=0$ are used ($U^{0}$).

• –

3: the $U$-factors are calculated from the loop-corrected rotation matrix for the lightest mass eigenstate ($U^{m_{1}}$).

The default value is 1.

• •

DECAYOPTIONS[1115]:

• –

0: the kinematics is done with tree-level masses.

• –

1: the kinematics is done with loop-corrected masses.

The default value is 1.

• •

DECAYOPTIONS[1116]:

• –

0: $\overline{\text{MS}}$/$\overline{\text{DR}}$ parameters are used for gauge couplings, $v$ and Yukawa couplings.

• –

1: $g_{1}$, $g_{2}$ and $v$ are set to reproduce the measured values of $M_{Z}$, $\alpha_{ew}(M_{Z})$ and $\sin\Theta_{W}$.

• –

2: same as 1, but in addition the Yukawa couplings are set to reproduce the measured values of SM fermions.

The default value is 0.

• •

DECAYOPTIONS[1117]:

• –

0: the counter-terms defined in RenConditionsDecays are not used.

• –

1: the counter-terms defined in RenConditionsDecays are used.

The default value is 0.

• •

DECAYOPTIONS[1118]:

• –

0: for loop-induced decays tree-level masses are used.

• –

1: for loop-induced decays loop-corrected masses are used.

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:

⬇

1Block DECAYOPTIONS #

2...

31101 ... # Only ultraviolet divergent parts of integrals

41102 ... # Ultraviolet divergence

51103 ... # Debug information

61104 ... # Only tree-level decay widths

7...

81201 ... # Photon/Gluon regulator mass

91205 ... # Renormalisation scale

The following settings are possible:

• •

DECAYOPTIONS[1101]: this option can be used to check the cancellation of ultraviolet divergences.

• –

0: One-loop functions employed in the calculation of one-loop decay widths return the finite part and the ultraviolet divergence defined in DECAYOPTIONS[1102].

• –

1: Only the ultraviolet divergence defined in DECAYOPTIONS[1102] is returned.

The default value is 0 .

• •

DECAYOPTIONS[1102]: this option can be used to check the cancellation of ultraviolet divergences. X sets the value used for the ultraviolet divergence $\Delta$ defined in Section 2.4. The default value is 0.

• •

DECAYOPTIONS[1103]:

• –

0: No debug information is shown.

• –

1: Additional information is shown on the screen. This includes the individual contributions from vertex, wavefunction and real corrections, which are useful to check the cancellation of ultraviolet and infra-red divergences.

The default value is 0.

• •

DECAYOPTIONS[1104]: this option can be used to check the consistency between the tree-level and one-loop calculation of decay widths.

• –

0: The one-loop routines return decay widths at NLO.

• –

1: The one-loop routines only return the tree-level decay widths, which can be compared to the tree-level results contained in Block DECAY.

The default value is 0.

• •

DECAYOPTIONS[1201]: this option can be used to check the cancellation of infra-red divergences. X defines the value (in GeV) used for the photon/gluon mass. The default value is 1.0E-5. Note that this option does not work with loop-corrected masses. The user should ensure that the regulator mass dependence of vertex and wavefunction corrections cancels against the one of the real corrections and yields a regulator-mass independent decay width. In order to show the individual contributions DECAYOPTIONS[1103] should be set to 1.

• •

DECAYOPTIONS[1205]: this option can be used to check the renormalisation scale dependence. If defined, X sets the value (in GeV) used for the renormalisation scale $Q$ in all one-loop functions employed in the calculation of decay widths. The default option is to use the same renormalisation scale as used in the calculation of masses, see Section 2.5.

4.4 Output of SPheno

The results of the one-loop calculation of decay widths are written in the SPheno output. For this purpose, we introduced the keyword DECAY1L beside the standard Block DECAY which lists the results of the ‘old’, i.e. leading order, calculation. Thus, for an arbitrary MSSM point, the output file contains:

⬇

1DECAY 1000001 5.03001929E+01 # Sd_1

2*# BR NDA ID1 ID2*

3 2.91393772E-01 2 6 -1000024 # BR(Sd_1 -> Fu_3 Cha_1)

4 1.70527978E-01 2 6 -1000037 # BR(Sd_1 -> Fu_3 Cha_2)

5 2.29399038E-04 2 3 1000023 # BR(Sd_1 -> Fd_2 Chi_2)

6 7.62216405E-03 2 5 1000022 # BR(Sd_1 -> Fd_3 Chi_1)

7 1.47930395E-01 2 5 1000023 # BR(Sd_1 -> Fd_3 Chi_2)

8 1.87674294E-03 2 5 1000025 # BR(Sd_1 -> Fd_3 Chi_3)

9 2.08779423E-03 2 5 1000035 # BR(Sd_1 -> Fd_3 Chi_4)

10 3.78314020E-01 2 1000002 -24 # BR(Sd_1 -> Su_1 VWm )

11...

12DECAY1L 1000001 5.07518318E+01 # Sd_1

13*# BR NDA ID1 ID2*

14 2.86487000E-01 2 6 -1000024 # BR(Sd_1 -> Fu_3 Cha_1)

15 1.63886304E-01 2 6 -1000037 # BR(Sd_1 -> Fu_3 Cha_2)

16 2.35164668E-04 2 3 1000023 # BR(Sd_1 -> Fd_2 Chi_2)

17 6.82861774E-03 2 5 1000022 # BR(Sd_1 -> Fd_3 Chi_1)

18 1.50793873E-01 2 5 1000023 # BR(Sd_1 -> Fd_3 Chi_2)

19 1.91193810E-03 2 5 1000025 # BR(Sd_1 -> Fd_3 Chi_3)

20 2.14376064E-03 2 5 1000035 # BR(Sd_1 -> Fd_3 Chi_4)

21 3.87696880E-01 2 1000002 -24 # BR(Sd_1 -> Su_1 VWm )

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 to the user how to combine them with the two-body decay widths obtained at the one-loop level.

5 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 $\overline{\text{MS}}$ or $\overline{\text{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 Section 4.2. We show two simple examples in the SM, namely the calculation of the partial decay width $t\to Wb$ and $H\to b\bar{b}$. In the first example we discuss different schemes for the renormalisation of the electric charge, in the second example we show that our $\overline{\text{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 $\overline{\text{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.

5.1 Renormalisation of $\alpha$ and the top-quark width

First we perform a calculation of the top-quark partial width in the decay $t\to Wb$ 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

[TABLE]

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. (2.29) can be applied to quark mixing. Subsequently we employ external tree-level masses without running, i.e. we effectively calculate with on-shell masses for all three involved particles (setting flag SPHENOINPUT[61]$=$0 to disable the RGE running for the parameters and and flag DECAYOPTIONS[1116]$=$2 to use on-shell mass values). This also fixes $g_{1},g_{2}$ and $v,m_{W}$ from $G_{F},m_{Z}$ and $\alpha(m_{Z})$. Our simple $\overline{\text{MS}}$ scheme for the renormalisation of $g_{2}$ (named scheme $(1)$) yields

[TABLE]

Next, we provide the decay width for the renormalisation of the electric charge in the Thomson limit of the $ff\gamma$-vertex, i.e. at zero momentum transfer [96]. The counter-terms for the electric charge are given by (see Ref. [112] for an overview)

[TABLE]

where we distinguish two schemes: At NLO we can make use of the very precise value of $\alpha(0)$ together with the corresponding counter-term $\delta Z_{e}(0)$ or we employ $\alpha(m_{Z})$ and compensate for the shift through the additional terms in $\delta 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

[TABLE]

such that in schemes $(2)$ and $(3)$ we obtain

[TABLE]

In this section we keep $\alpha$ and $\alpha_{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 $\alpha(0)$ or $\alpha(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 $\Gamma^{\text{{NLO},EW+{QCD}}}_{t\to Wb}$ they are with respect to the partial width $\Gamma^{\text{{NLO},EW}}_{t\to Wb}$.

For scheme $(1)$ evaluated with $\alpha(m_{Z})$ we obtain $\Gamma_{t\to Wb}^{\text{{NLO}{},EW+{QCD}{}}}=1.434~{}{\rm 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 $\alpha$ and $\alpha_{s}$) and yields $\sim-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 $\alpha$. 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 $\alpha(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 Section 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.

5.2 Renormalisation of Yukawa couplings and fermionic Higgs decays

We also shortly discuss the calculation of $H\to b\bar{b}$ in the SM, which is mediated through the bottom-quark Yukawa coupling $Y_{b}$, and it turns out that the $\overline{\text{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

[TABLE]

such that renormalisation prescriptions for $\delta m_{b}$ and $\delta v$ are needed. Whereas $\delta m_{b}$ can be obtained from the self-energies of the down-type quarks, the on-shell renormalisation of the vacuum expectation value depends on other parameters: one requests 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 ($\sim 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 $\overline{\text{MS}}$ scheme is more convenient, since it resums large logarithms [118, 119]. We demonstrate this effect in Table 2, which is obtained with the SM version of SPheno setting flag DECAYOPTIONS[1116]$=$1 and flag SPHENOINPUT[61]$=$1. Through these settings $Y_{b}$ as well as the gauge couplings, in particular $\alpha_{s}$, are evaluated at the renormalisation scale $Q$ and for $Y_{b}$ the $\overline{\text{MS}}$ scheme is employed. We again depict the LO as well as the NLO partial width with only electroweak as well as electroweak and QCD corrections including relative corrections. The most relevant parameters are $m_{H}=125$ GeV, $m_{b}(m_{b})=4.18$ GeV and $\alpha_{s}(m_{Z})=0.1187$. The running to $Q=m_{H}$ yields $Y_{b}\propto m_{b}(m_{H})=2.781$ GeV and $\alpha_{s}(m_{Z})=0.1133$.

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\alpha_{s}(m_{H})/\pi\sim 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.

5.3 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 two-fold: we compared neutralino and chargino decays into neutralinos and charginos and heavy gauge bosons with CNNDecays, where we employ a full on-shell 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 Section 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 $\overline{\text{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 ($\overline{\text{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 loop-corrected masses can change the results substantially, see also Section 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:

[TABLE]

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\beta=10$. We show results for the predicted partial widths when varying $m$ from $300$ GeV 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 Section 5.3.5.

5.3.1 Neutralino and Chargino decays in the MSSM and in bilinear $R$-parity violation: CNNDecays vs. SARAH

We compared the decay modes $\tilde{\chi}_{i}^{\pm}\to\tilde{\chi}_{j}^{0}W^{\pm}$ and $\tilde{\chi}_{i}^{0}\to\tilde{\chi}_{j}^{\mp}W^{\pm}$ as well as $\tilde{\chi}_{i}^{0}\to\tilde{\chi}_{j}^{0}Z$ and $\tilde{\chi}_{i}^{\pm}\to\tilde{\chi}_{j}^{\pm}Z$ in the $R$-parity 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 Section 4.2, the partial decay widths are in principle all renormalisation-scale 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 wavefunction corrections individually as well as the individual pieces to the counter-terms. Also for the $R$-parity violating decays $\tilde{\chi}_{4}^{0}\to l^{\mp}W^{\pm}$ 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. $\tilde{\chi}_{4}^{0}\to\nu Z$ or $\nu 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 $\mathcal{CP}$-violating scenarios as discussed in Refs. [67, 70, 78] is left for future work for the reasons explained in Section 2.8.

5.3.2 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 [math] in SFOLD.F:

⬇

1 osextmassesOn = 0

2 osloopmassesOn = 0

This is done to ensure that both codes use exactly the same masses everywhere.

• •

We find a disagreement for the bremsstrahlung routines for $S\to SV$ decays. Therefore, we add to line 440 in Bremsstrahlung.F of SFOLD the terms:

⬇

1 -2g1**2gtgtCIi - 2g1**2gtgtCI0up1

• •

We find for $S\to SS$ and $S\to 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 555In 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_{\xi}$ gauge.. The problem is avoided by setting

⬇

1 Xipart1 = 0d0

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 wavefunction renormalisation constants using Eq. (2.29). In particular for $S\to SV$ decays this can induce large differences in the one-loop corrections. If the finite parts for the counter-terms of the rotation matrices are included, a cancellation between the wavefunction corrections and the counter-term 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:

⬇

1 Do i1=1,3

2 Do i2=1,3

3 dUSf1(:,:,i1,i2) = 0.25*MatMul(dZSf1(:,:,i1,i2) &

4 & - Conjg(Transpose(dZSf1(:,:,i1,i2))),USf(:,:,i1,i2))

5 dUSfC1(:,:,i1,i2) = Conjg(0.25*MatMul(dZSf1(:,:,i1,i2) &

6 & - Conjg(Transpose(dZSf1(:,:,i1,i2))),USf(:,:,i1,i2)))

7 End Do

8 End Do

9 dVCha1 = MatMul(dZChaL1 - Conjg(Transpose(dZChaL1)),VCha)/4

10 dUCha1 = MatMul(dZChaR1 - Conjg(Transpose(dZChaR1)),UCha)/4

11 dZNeuRM1 = MatMul(dZNeuL1 - Conjg(Transpose(dZNeuL1)),ZNeu)/4

12 dVChaC1 = Conjg(dVCha1)

13 dUChaC1 = Conjg(dUCha1)

14 dZNeuRMC1 = Conjg(dZNeuRM1)

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

[TABLE]

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 $\tau$-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.

5.3.3 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$%.

5.3.4 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 SM-like Higgs boson. For such decays, we leave an appropriate definition of counter-terms, following our explanations in Section 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 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, that 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 [108, 40]. We shortly discuss the relevance of $U$-factors for Higgs boson decays in Section 5.4.

5.3.5 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 widths666Even if SARAH will use the new routines to obtain also loop-induced decay width for Higgs states, one should still use the old results for the diphoton and digluon rate. The latter also include the full model dependence at LO but in addition also higher-order QCD corrections, see Ref. [122]. The most important applications for these routines are radiatively induced decays of BSM particles. The main candidates for such decays in the MSSM are $\tilde{\chi}^{0}_{2}\to\tilde{\chi}^{0}_{1}\gamma$ and $\tilde{g}\to\tilde{\chi}^{0}_{1}g$. Those decays were already implemented in SPheno based on the results of Ref. [123]. 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 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 wavefunction 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 kind of neutralinos, see again Fig. 8. Note that throughout the calculation of loop-induced decays loop-corrected masses are inserted.

5.4 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–11 the size of the one-loop corrections for selected stop, sbottom and Higgs decays, respectively, using the three available options to calculate the $U$-factors. We apply loop-corrected $\overline{\text{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 shown value $\Delta_{U}$ is defined as

[TABLE]

Here, $\Gamma_{0}$ are the decay widths without applying $U$-factors. $\Delta_{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 momentum-dependence 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 width $\tilde{t}_{2}\to\tilde{\chi}_{4}^{0}t$, where the relative NLO corrections encoded in $\Delta 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. [108, 40] are in order in future work.

6 Conclusions

In this paper we described a fully generic implementation of the calculation of two-body decay widths at the full one-loop 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. Wavefunction corrections are determined from on-shell conditions. On the other hand, the parameters of the underlying model are by default renormalised in a $\overline{\text{DR}}$ (or $\overline{\text{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 infra-red 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 $\overline{\text{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 non-supersymmetric models. However, models with $\mathcal{CP}$ violation and/or (additional) massive gauge bosons charged under U$(1)_{\rm em}\times\text{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 left for future work.

Acknowledgements

We are grateful to Dominik Stöckinger for discussions at the early stage of this project. We thank Helmut Eberl for help when comparing to the codes SFOLD and HFOLD, and Thomas Hahn for his support concerning FeynArts/FormCalc. SL acknowledges support by Deutsche Forschungsgemeinschaft through the SFB 676 “Particles, Strings and the Early Universe”. MDG acknowledges support from French state funds managed by the Agence Nationale de la Recherche (ANR), in the context of the LABEX ILP (ANR-11-IDEX-0004-02, ANR-10-LABX-63), and the young researcher grant “HiggsAutomator” (ANR-15-CE31-0002). FS is supported by ERC Recognition Award ERC-RA-0008 of the Helmholtz Association.

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. We factor out the Lorentz dependent part of the vertices and work with the following conventions

• •

$FFS$ vertex ($\bar{F}_{1}F_{2}S$)

[TABLE]

• •

$FFV$ vertex ($\bar{F}_{1}F_{2}V^{\mu}$)

[TABLE]

• •

$SSV$ vertex ($S^{*}_{1}S_{2}V^{\mu}$)

[TABLE]

• •

$SVV$ vertex ($S_{1}V_{1}^{\nu}V_{2}^{\mu}$)

[TABLE]

• •

$VVV$ vertex ($V^{1}_{\mu}V^{2}_{\nu}V^{3}_{\sigma}$)

[TABLE]

• •

$VVVV$ vertex ($V^{1}_{\mu}V^{2}_{\nu}V^{3}_{\sigma}V^{4}_{\rho}$)

[TABLE]

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

[TABLE]

and

[TABLE]

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, $\dot{B}_{0}$, $\dot{B}_{1}$, which are used by SPheno are given in the appendix of Ref. [75].

A.1 Wavefunction corrections

A.1.1 Fermion

(a) $FS$ diagram

[TABLE]

[TABLE]

(b) $FV$ diagram

[TABLE]

[TABLE]

For Majorana fermions, an additional overall factor $\frac{1}{2}$ is present.

A.1.2 Scalar

(a) $FF$ diagram

[TABLE]

with

[TABLE]

(b) $SS$ diagram

[TABLE]

(c) $UU$ diagram

[TABLE]

(d) $VV$ diagram

[TABLE]

(e) $SV$ diagram

[TABLE]

with

[TABLE]

(f) $S$ diagram

[TABLE]

(g) $V$ diagram

[TABLE]

A.1.3 Gauge boson

(a) $FF$ diagram

[TABLE]

with

[TABLE]

(b) $SS$ diagram

[TABLE]

(c) $UU$ diagram

[TABLE]

(d) $VV$ diagram

[TABLE]

with

[TABLE]

(e) $SV$ diagram

[TABLE]

(f) $S$ diagram

[TABLE]

(g) $V$ diagram

[TABLE]

A.2 Vertex corrections

A.2.1 Fermion to fermion and scalar decays

(a) $SFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{2}^{2},p_{0}^{2},p_{1}^{2},m_{3}^{2},m_{2}^{2},m_{1}^{2})$.

(b) $FSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(c) $VFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{2}^{2},p_{0}^{2},p_{1}^{2},m_{3}^{2},m_{2}^{2},m_{1}^{2})$.

(d) $FSV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(e) $FVS$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{3}^{2},m_{2}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(f) $FVV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

A.2.2 Fermion to fermion and gauge boson decays

(a) $SFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{2}^{2},p_{1}^{2},p_{0}^{2},m_{2}^{2},m_{3}^{2},m_{1}^{2})$.

(b) $FSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(c) $VFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{2}^{2},p_{1}^{2},p_{0}^{2},m_{2}^{2},m_{3}^{2},m_{1}^{2})$.

(d) $FSV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(e) $FVS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(f) $FVV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

A.2.3 Scalar to two fermion decays

(a) $FFS$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{0}^{2},p_{2}^{2},p_{1}^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2})$.

(b) $SSF$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{0}^{2},p_{2}^{2},m_{3}^{2},m_{1}^{2},m_{2}^{2})$.

(c) $FFV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{0}^{2},p_{2}^{2},p_{1}^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2})$.

(d) $SVF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{3}^{2},m_{2}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{0}^{2},p_{2}^{2},m_{3}^{2},m_{1}^{2},m_{2}^{2})$.

(e) $VSF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{3}^{2},m_{2}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{0}^{2},p_{2}^{2},m_{3}^{2},m_{1}^{2},m_{2}^{2})$.

(f) $VVF$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{0}^{2},p_{2}^{2},m_{3}^{2},m_{1}^{2},m_{2}^{2})$.

A.2.4 Scalar to two scalar decays

(a) $FFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(b) $SSS$ diagram

[TABLE]

(c) $UUU$ diagram

[TABLE]

(d) $VVV$ diagram

[TABLE]

(e) $SSV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(f) $SVS$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{3}^{2},m_{2}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(g) $VSS$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(h) $SVV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(i) $VSV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(j) $VVS$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{3}^{2},m_{2}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(k) $SS$ diagram

[TABLE]

(l) $VV$ diagram

[TABLE]

(m) $SS$ diagram

[TABLE]

(n) $VV$ diagram

[TABLE]

(o) $SS$ diagram

[TABLE]

(p) $VV$ diagram

[TABLE]

A.2.5 Scalar to scalar and gauge boson decays

(a) $FFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(b) $SSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(c) $UUU$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(d) $VVV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(e) $SSV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(f) $SVS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(g) $VSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(h) $SVV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(i) $VSV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(j) $VVS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(k) $SV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{0}^{2},m_{1}^{2},m_{2}^{2})$.

(l) $SV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{1}^{2},m_{1}^{2},m_{2}^{2})$.

(m) $SS$ diagram

[TABLE]

(n) $VV$ diagram

[TABLE]

A.2.6 Scalar to two gauge boson decays

(a) $FFF$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(b) $SSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(c) $UUU$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(d) $VVV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(e) $SSV$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(f) $SVS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(g) $VSS$ diagram

[TABLE]

with $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(h) $SVV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(i) $VSV$ diagram

[TABLE]

with $B_{i}=B_{i}(p_{2}^{2},m_{2}^{2},m_{3}^{2})$ and $C_{i}=C_{i}(p_{1}^{2},p_{2}^{2},p_{0}^{2},m_{1}^{2},m_{3}^{2},m_{2}^{2})$.

(j) $VVS$ diagram

[TABLE]

(k) $SS$ diagram

[TABLE]

(l) $VV$ diagram

[TABLE]

(m) $SV$ diagram

[TABLE]

(n) $SV$ diagram

[TABLE]

Appendix B Expressions for real corrections

Following the notation introduced in Section 2.6, we present here the formulas 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 $\mathcal{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

[TABLE]

If we attach a massless gauge boson $A_{\mu}^{b}$ to the above term then we can have the coupling

[TABLE]

We should now understand that the indices for $V,A$ correspond to the unbroken gauge symmetry under which $A_{\mu}^{b}$ transforms. By insisting on gauge invariance up to linear order we find three sets of conditions:

[TABLE]

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_{\mu}^{a}\rightarrow A_{\mu}^{a}+\partial_{\mu}\Lambda^{a}+\ldots$ part. Importantly, Eq. (B.5) completely determines $c_{ijab}$.

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

[TABLE]

The above logic can also be used to compute the group theory 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_{ijk}$ in the same was as we did in Eq. (B.1) and Eq. (B.2) and strip off any Lorentz indices, so for fields (scalars, fermions or vectors) $\Phi_{i}$ in representations of the gauge group $R_{0}\rightarrow R_{1},R_{2}$ we have the coupling

[TABLE]

Then the matrix element is

[TABLE]

where the functions $f_{i}^{\mu}(p_{i})$ may have spinor or other Lorentz indices, and include the wavefunction factors, as appropriate to the final states; and the index $i$ labels the leg to which the photon/gluon is attached: [math] 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 infra-red divergences – and also because they allow easy identification of the appropriate diagram.

Then, defining $x_{0}\equiv-2p_{0}\cdot k,x_{1}\equiv 2p_{1}\cdot k,x_{2}\equiv 2p_{2}\cdot k,x_{3}\equiv 1$, a generic squared matrix element can be expressed as

[TABLE]

Now all of the group-theory factors are encoded in $\tilde{C}_{ij}$. However, while $\tilde{C}_{ji}=\tilde{C}_{ij}^{*}$, it is also clear from the hermeticity of the generators that the $\tilde{C}_{ij}$ are real – for example using $d(R_{0})$ as the dimension of representation $R_{0}$

[TABLE]

and therefore we can write

[TABLE]

Using the gauge invariance of the Lagrangian terms, we obtain for $i<3$:

[TABLE]

where for $S\rightarrow 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 two-body decay factor $C$

[TABLE]

because we can write

[TABLE]

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

[TABLE]

Hence all of the group theory factors are proportional to $C$, and we can therefore define

[TABLE]

giving

[TABLE]

and we can write for a generic process

[TABLE]

In the case of a U$(1)$ gauge boson (the photon), we have $C_{ij}=Q_{i}Q_{j}$ and we define $Q_{3}\equiv Q_{0}+Q_{1}.$ This also follows from the above expressions, for example

[TABLE]

B.2 $F\to FS$

The real corrections for the fermionic decays $F\to FS$ are given by

[TABLE]

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\to FS$. We identify $F_{\text{in}}=X$, $F_{\text{out}}=Y_{1}$ and $S=Y_{2}$ with the indices $0,1$ and $2$, respectively. The individual contributions $\Omega_{ij}^{AB}$ are given by:

[TABLE]

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\to FF$

The real corrections for the scalar decays $S\to FF$ are given by

[TABLE]

The notation follows 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 $\Omega_{ij}^{AB}$ are given by:

[TABLE]

We also implemented real corrections where one of the final state fermions can be massless. Again this fermion has to be charge and colour neutral.

B.4 $F\to FV$

The real corrections for the fermionic decays $F\to FV$ are given by

[TABLE]

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\to FV$. 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\to FW+\gamma$ on the one side and $F\to FZ+\gamma/g$ and $F\to FW+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_{ij}=Q_{i}Q_{j}$, where we again identify $F_{\text{in}}=X$ and $F_{\text{out}}=Y_{1}$ with indices [math] and $1$, respectively.

For $F\to FW+\gamma$ the individual contributions are given by:

[TABLE]

For $F\to FZ+\gamma/g$ and $F\to FW+g$ the individual contributions are given by:

[TABLE]

B.5 $S\to VV$

The real corrections for the decays $S\to VV$ are given by

[TABLE]

The coupling $c$ is the tree-level coupling of the vertex of $S\to VV$. Given that both final state particles $V$ are heavy gauge bosons only photon emission is of relevance for this process. We identify $S=X$, $V=Y_{1}$ and $V=Y_{2}$ with the indices $0,1$ and $2$, respectively. The individual contributions $\Omega_{ij}$ are given by:

[TABLE]

B.6 $S\to SS$

The real corrections for the decays $S\to SS$ are given by

[TABLE]

The coupling $c_{g}$ is the electromagnetic coupling $e$ or the strong coupling $g_{s}$ for photon and gluon emission, respectively. $c$ is the coupling of the tree-level vertex of $S\to SS$. Therein we identify $S=X$, $S=Y_{1}$ and $S=Y_{2}$ with the indices $0,1$ and $2$, respectively. The individual contributions $\Omega_{ij}$ are given by:

[TABLE]

B.7 $S\to SV$

The real corrections for the decays $S\to SV$ are given by

[TABLE]

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\to SV+\gamma/g$ also has a contribution from a four-point interaction, which accordingly does not contribute to the infra-red divergent part of the real corrections, but yields a finite contribution. The sum over the gauge factors thus includes four indices $0,\ldots,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_{ij}=Q_{i}Q_{j}$ for $i,j\in\{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

[TABLE]

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 $\Omega_{ij}$ are given by:

[TABLE]

Appendix C Infra-red divergent parts of Passarino-Veltman integrals

Subsequently we present the infra-red divergent parts of the Passarino-Veltman integrals relevant for our purposes regularized through a photon or gluon mass $m_{\Lambda}$:

[TABLE]

Therein we use the Källén function given in Eq. (2.2). The result for $C_{0}$ is consistent with the infra-red divergent part of Eq. (B.5) of Ref. [126]. For the cases $p_{1}^{2},p_{2}^{2}\ll p_{0}^{2}$ and $p_{2}^{2}\ll p_{1}^{2},p_{0}^{2}$ which is of particular relevance for scalar decays into light fermions, we also implemented formulas

[TABLE]

These are equivalent to Eqs. (B.8) and (B.9) of Ref. [126].

Appendix D Goldstone boson vertices

For decays into scalars and fermions, the cancellation of infra-red divergences is straightforward: they correspond to summing the real emission of a massless gauge boson with the virtual process of the same massless gauge boson in the loop. Because the current of the unbroken gauge symmetries are necessarily flavour diagonal, the only vertices involved are the original two-body decay vertex and the gauge couplings of the external states. Therefore, when we want to use loop-corrected external masses to cancel the infra-red divergences for this case it is straightforward to either put all masses of internal and external states to the loop-corrected values, or (as we do here) to subtract off the infra-red divergent parts separately from the Bremsstrahlung and virtual corrections before multiplying by a kinematic factor employing loop-corrected masses.

However, for decays with a massive gauge boson in the final state, the cancellation between the real and virtual infra-red divergences is a little subtle. The reason is that a would-be Goldstone boson can propagate as an internal state. While we perform the Bremsstrahlung calculation in the unitary gauge (so there are no Goldstone bosons), for the virtual corrections for practical purposes the default choice is Feynman-’t Hooft gauge, and we must therefore sum the diagrams with an internal massive gauge boson with those having a massive Goldstone boson. We show the relevant virtual corrections for the processes $F\rightarrow FV$, $S\rightarrow SV$ and $S\rightarrow VV$ (for decays to massive vectors) in Figs. 21, 22 and 23 respectively. Denoting the heavy gauge boson as “$V$”, the massless one as “$\gamma$,” and the Goldstone boson as “$G$,” we see that we have the gauge coupling of the unbroken group appearing in $FF\gamma,SS\gamma$ and $VV\gamma$ vertices, but we must maintain a relationship between the $FFV$ and $FFG$ vertices, between the $SSV$ and $SSG$ vertices, the $VV\gamma$ and $VG\gamma$ vertices, in order for these cancellations to occur. For the $S\rightarrow VV$ process, we require a relationship between the three sets of vertices $SVV$, $SGV$ and $SGG$.

The required relations follow from Slavnov-Taylor identities, or alternatively we could examine the infra-red divergent part of the loop amplitudes. However, here we will more simply derive the conditions imposed by symmetry when inspecting the lagrangian at tree level. We first identify the would-be Goldstone boson by writing the symmetry transformations of real scalars $S_{i}^{0}$ (before spontaneous symmetry breaking, and not necessarily in a mass-diagonal basis) as

[TABLE]

where $\epsilon^{G}$ is the gauge transformation parameter for each broken direction $G$ and $a_{ij}^{G}$ are numbers for a linearly realised gauge symmetry (as always assumed in SARAH). Then

[TABLE]

where $t_{ij}^{G}$ are the generators of the gauge symmetry that will be broken by appropriate vacuum expectation values; since we are working with real scalars, the $t_{ij}^{G}$ are antisymmetric and real (to translate to complex fields we require roughly $t_{ij}^{G}\rightarrow iT^{G}_{ij}$). Then we find that the Goldstones $G$ are defined by

[TABLE]

The vectors should be chosen to be orthogonal with normalisation $1$, so $N_{G}=1/\sqrt{\sum_{i}(\alpha_{i}^{G})^{2}}$.

Looking at the vector mass term we find for the scalars $S_{i}^{0}$ having vaccum expectation values $v_{i}$

[TABLE]

we see that

[TABLE]

where $m_{V}^{G}$ is the mass of the vector, and we define the $SSV$ coupling when we diagonalise the scalars to mass eigenstates through $S_{i}^{0}\equiv R_{ik}S_{k}$ (where $R_{Gj}=N_{G}\alpha_{j}^{G}$):

[TABLE]

For the $SGV$ coupling and $SVV$ couplings, we read off

[TABLE]

This is the relationship that we enforce between the on-shell vertices to ensure that infra-red divergences are cancelled.

We can also read off the $GVV$ coupling

[TABLE]

For any unbroken gauge groups, the Goldstones must transform as

[TABLE]

and we therefore identify the $\alpha_{i}^{G^{\prime}}t_{ij}^{a}\alpha_{j}^{G}$ factor in Eq. (D.8) with $T_{GG^{\prime}}^{a}$ to obtain

[TABLE]

For the photon, this just becomes the familiar vertex

[TABLE]

Now the $VV\gamma$ 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

[TABLE]

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

[TABLE]

This leads to the relation between the $\gamma GV$ and $\gamma VV$ vertices of simply a factor of $m_{V}^{G}$.

To find the relationship between the Goldstone coupling to scalars and the gauge boson coupling, we can use Eq. (2.32) of Ref. [127] to find the derivative of the (effective) potential $V$ with respect to the scalars:

[TABLE]

When we diagonalise the masses this gives

[TABLE]

Again we enforce this relationship between the $GSS$ coupling and the $SSV$ coupling when we use loop-corrected masses in order to cancel infra-red divergences.

The $GGS$ coupling is a special case of the above, with the understanding that we must use zero for the Goldstone boson mass to obtain

[TABLE]

We also find the a similar relationship for fermions, which can in that case also be derived from the fact that the masses are Yukawa couplings; writing in terms of Weyl fermions we have:

[TABLE]

The gauge invariance of the Yukawa coupling gives

[TABLE]

and so

[TABLE]

Now we write down the gauge coupling

[TABLE]

and we diagonalise into $\psi_{I}^{0}=R_{IJ}\psi_{J}$ so we obtain

[TABLE]

If we now split the fermions into left and right-handed states with separate rotation matrices $L$ and $R$, then we the gauge couplings are given in Dirac notation by

[TABLE]

while the Goldstone couplings are

[TABLE]

If we work in a basis where the fermion masses are complex and the rotation matrices $R,L$ are real, then we find the relationship between the Goldstone boson couplings and the gauge couplings

[TABLE]