Precision tools and models to narrow in on the 750 GeV diphoton resonance

The hints for a new resonance at 750 GeV from ATLAS and CMS have triggered a significant amount of attention. Since the simplest extensions of the standard model cannot accommodate the observation, many alternatives have been considered to explain the excess. Here we focus on several proposed renormalisable weakly-coupled models and revisit results given in the literature. We point out that physically important subtleties are often missed or neglected. To facilitate the study of the excess we have created a collection of 40 model files, selected from recent literature, for the Mathematica package SARAH. With SARAH one can generate files to perform numerical studies using the tailor-made spectrum generators FlexibleSUSY and SPheno. These have been extended to automatically include crucial higher order corrections to the diphoton and digluon decay rates for both CP-even and CP-odd scalars. Additionally, we have extended the UFO and CalcHep interfaces of SARAH, to pass the precise information about the effective vertices from the spectrum generator to a Monte-Carlo tool. Finally, as an example to demonstrate the power of the entire setup, we present a new supersymmetric model that accommodates the diphoton excess, explicitly demonstrating how a large width can be obtained. We explicitly show several steps in detail to elucidate the use of these public tools in the precision study of this model.

Abstract: The hints for a new resonance at 750 GeV from ATLAS and CMS have triggered a significant amount of attention. Since the simplest extensions of the standard model cannot accommodate the observation, many alternatives have been considered to explain the excess. Here we focus on several proposed renormalisable weaklycoupled models and revisit results given in the literature. We point out that physically important subtleties are often missed or neglected. To facilitate the study of the excess we have created a collection of 40 model files, selected from recent literature, for the Mathematica package SARAH. With SARAH one can generate files to perform numerical studies using the tailor-made spectrum generators FlexibleSUSY and SPheno. These have been extended to automatically include crucial higher order corrections to the diphoton and digluon decay rates for both CP-even and CP-odd scalars. Additionally, we have extended the UFO and CalcHep interfaces of SARAH, to pass the precise information about the effective vertices from the spectrum generator to a Monte-Carlo tool. Finally, as an example to demonstrate the power of the entire setup, we present a new supersymmetric model that accommodates the diphoton excess, explicitly demonstrating how a large width can be obtained. We explicitly show several steps in detail to elucidate the use of these public tools in the precision study of this model.

Introduction
The first data from the 13 TeV run of the Large Hadron Collider (LHC) contained a surprise: ATLAS and CMS reported a resonance at about 750 GeV in the diphoton channel with local significances of 3.9σ and 2.6σ, respectively [1,2]. When including the look-elsewhere-effect, the deviations from standard model (SM) expectations drop to 2.3σ and 1.2σ. [238], HiggsBounds [239,240], HiggsSignals [241] or Vevacious [242] can easily be included in the framework. These powerful packages provide a way to get a thorough understanding of the new models. The main goal of this work is to support the model builders and encourage them to use these tools. We provide several details about the features of the packages in the spirit of making this paper self-contained and bringing the reader unfamiliar with the tools to the level of knowledge necessary to use them. More information can be found in the manuals of each package. We have created a database of diphoton models in SARAH, by implementing 40 among those proposed in recent literature, which is now available to all interested researchers. For each model we have written model files to interface SARAH with SPheno and FlexibleSUSY.
Although in each case we have tried to check very carefully that we implement the model which has been proposed in the literature, it is of course possible that some details have been missed. The original authors of these models are encouraged to check the implementation themselves to satisfy that what we have implemented really does correspond to the model they proposed. In the description of some of the models we state cases where the model has problems or where we find difficulties for the proposed solution. This helps inform potential users about what they may see when running these models through the tools we are discussing here. However especially in these cases we encourage the original authors to check what we have written and let us know if they disagree with any claim we make.
This paper is long but modular, and each section is to a large extent self contained, so the reader can easily jump to the section of greater interest. We have structured it as follows: • In Section 2, we give a list of common mistakes we have found in the literature and emphasise how they are easily avoided by using the tools. This provides the main motivation of this paper and we hope that other model builders will also see the necessity of using these packages.
• In Section 3, we discuss the basic features of SARAH, describing how one can use it to extract all analytical properties of the model and how to generate model files or source code for the numerical tools.
• In Section 4, we give an introduction to SPheno and FlexibleSUSY and explain how these codes can be easily interfaced with MC tools. We also discuss at some length the implementation of the diphoton and digluon effective vertices.
• In Section 5, we give an overview of the models which we have implemented in SARAH and briefly discuss their main features.
• In Section 6, we provide an explicit example of how to quickly work out the details of a model, analytically with SARAH and numerically with the other tools. For this purpose we extended a natural supersymmetric (SUSY) model to accommodate the 750 GeV resonance.
• We conclude in Section 7.
• We provide three appendices containing frequently asked questions (FAQs) aimed at further clarifying the use of the packages.

Motivation
Precision studies in high energy physics have reached a high level of automation. There are publicly available tools to perform Monte-Carlo studies at LO or NLO [243][244][245][246][247][248], many spectrum generators [222,223,[249][250][251][252][253][254][255][256][257][258][259][260] for the calculation of pole masses including important higher order corrections, codes dedicated only to Higgs [261][262][263][264] or sparticle decays [265][266][267], and codes to check flavour [268][269][270][271][272] or other precision observables [273]. This machinery has been used in the past mainly for detailed studies of some promising BSM candidates, like the MSSM, NMSSM or variants of THDMs. There are two main reasons why these tools are usually the preferred method to study these models: (i) it has been shown that there can be large differences between the exact numerical results and the analytic approximations; (ii) writing private routines for specific calculations is not only time consuming but also error prone. On the other hand, the number of tools available to study the new ideas proposed to explain the diphoton excess is still limited. Of these tools, many are not yet widely used largely due to the community's reluctance in adopting new codes. However, we think it is beneficial to adopt this new generation of generic tools like SARAH.
We noticed that several studies done in the context of the 750 GeV excess have overlooked important subtleties in some models, neglected important higher order corrections, or made many simplifying assumptions which are difficult to justify. Using generic software tools in this context can help address these issues: many simplifications will no longer be necessary and important higher order corrections can be taken into account in a consistent manner. In order to illustrate this we comment, in the following subsections, on several issues we became aware of when revisiting some of the results in the literature.

The diphoton and digluon rates beyond leading order
A precise calculation of the diphoton rate is of crucial importance. In the validation process of this work, we identified several results in the literature that deviate, often by an order of magnitude or more, in comparison to our results [65,179,202]. Additionally we observed that in many cases there are important subtleties which we think are highly relevant.
First of all, the choice of the renormalisation scale of the running couplings appearing in the calculation. The majority of recent studies use the electromagnetic coupling at the scale of the decaying particle. However, one should rather use α em (0), i.e. the Thompson limit (see for instance Refs. [274,275]), in order to keep the NLO corrections under control. Taking this into account already amounts to an O(10 %) change of the diphoton rate compared to many studies in the literature. In addition, as we will discuss in Section 2.2.3, an important prediction of a model is the ratio Br(S → gg)/Br(S → γγ), where S is the 750 GeV scalar resonance. It is well known that the digluon channel receives large QCD corrections. If one neglects these corrections the ratio will be severely underestimated.
To demonstrate these effects we show in Fig. 1 the total decay width 1 of the singlet S as a function of the mass M F 1 and coupling Y F 1 for a simple toy model containing only the vector-like fermions Ψ F 1 , presented in Section 5.1. Table 1 contains benchmark points for the partial widths of the digluon and diphoton channels as well as the ratio of these two channels for both CP-even and CP-odd scalar resonances. This table contains the LO calculations performed using SPheno as a comparison to results previously shown in the literature [166]. We also show the partial widths including NLO corrections for the diphoton channel 2 and N 3 LO QCD corrections to the gluon fusion production as implemented in Section 4.4. The discrepancies between the LO calculations arise purely through the choice of the renormalisation scale for the gauge couplings. However, the NLO results clearly emphasise that loop corrections at the considered mass scales are the dominant source of errors. To our knowledge, these uncertainties have thus far not received a sufficiently careful treatment in the literature; we give further discussion of this (and the remaining uncertainty in the SARAH calculation) in Section 4.8.  Figure 1. The approximate total width (sum of the diphoton and digluon channels) of S as a function of the coupling Y F 1 and the mass m F 1 of the vector-like particle F 1 , calculated using SPheno (blue) at LO (dashed) and NLO (solid). The orange contours are the results of the LO calculation from Ref. [166]. Here we assume a single generation of vector-like quarks.

Constraints on a large diphoton width
In order to explain the measured signal, one needs a large diphoton rate of Γ(S → γγ)/M S 10 −6 assuming a narrow width for S, while for a large width one requires Γ(S → γγ)/M S 10 −4 [114]. In weakly-coupled models there are three different possibilities to obtain such a large width: 1. Assuming a large Yukawa-like coupling between the resonance and charged fermions 2. Assuming a large cubic coupling between the resonance and charged scalars 3. Using a large multiplicity and/or a large electric charge for the scalars and/or fermions in the loop  Table 1. Branching fraction ratio, as well as the partial decay widths for the digluon and diphoton channels for the toy model (Section 5.1) containing only the relevant vector-like fermion pair Ψ F i . The above values are for the benchmark points Y F i = 1 and m F i = 1 TeV, where the values are for a CP-even/CP-odd scalar resonances, respectively. The SPheno NLO calculation includes N 3 LO corrections for the digluon channel, while the diphoton decay width is calculated at NLO and LO for a CP-even and odd scalar respectively.
However, all three possibilities are also constrained by very fundamental considerations, which we briefly summarise in the following.

Large couplings to fermions
A common idea to explain the diphoton excess is the presence of vector-like states which enhance the loop-induced coupling of a neutral scalar to two photons or two gluons. This led some authors to consider Yukawa-like couplings of the scalar to the vector-like fermions larger than √ 4π, which is clearly beyond the perturbative regime 3 . Nevertheless, a one-loop calculation is used in these analyses to obtain predictions for the partial widths [209], despite being in a non-perturbative region of parameter space. Moreover, even if the couplings are chosen to be within the perturbative regime at the scale Q = 750 GeV, they can quickly grow at higher energies. This issue of a Landau pole has been already discussed to some extent in the literature [22,36,114,125,205,206], and one should ensure that the model does not break down at unrealistic small scales.

Large couplings to scalars
One possibility to circumvent large Yukawa couplings is to introduce charged scalars, which give large loop contributions to the diphoton/digluon decay. A large cubic coupling between the charged scalar and the 750 GeV one does not lead to a Landau pole for the dimensionless couplings because of dimensional reasons. However, it is known that large cubic couplings can destabilize the scalar potential: if they are too large, the electroweak vacuum could tunnel into a deeper vacuum where U (1) em gauge invariance is spontaneously broken, depending on the considered scenario. The simplest example with such a scenario is the SM, extended by a real scalar S and a complex scalar X with hypercharge Y . The scalar potential of this example is (2.1) In Fig. 2 the dependence of the diphoton partial width as a function of κ and M X is shown, and the stability of the electroweak potential as well as the life-time of its ground state is checked with Vevacious and CosmoTransitions. For more details about vacuum stability in the presence of large scalar cubic terms, we refer to Ref. [205]. The overall conclusion of [205] is that the maximal possible diphoton width, even when allowing for a meta-stable but sufficiently long-lived electroweak vacuum, is not much larger than in the case of vector-like fermions when requiring that the model is perturbative up to the Planck scale. It is therefore essential to perform these checks when studying a model that predicts large cubic scalar couplings.

Large multiplicities
To circumvent large Yukawa or cubic couplings, other models require a large number of generations of new BSM fields and/or large electric charges. As a consequence the running of the U (1) Y gauge coupling, g 1 , gets strongly enhanced well below the Planck scale. Moreover, even before reaching the Landau pole, the model develops large (eventually non-perturbative) gauge couplings. This implies an enhancement of Drell-Yan processes at the LHC, with current data already setting stringent constraints and potentially excluding some of the models proposed to explain the diphoton excess [126,276]. For general studies on running effects in the context of the diphoton excess see [22,36,114,125,206]. We briefly discuss dramatic examples of this class of models proposed in Refs. [158] and [194], which feature approximately ∼ 100 and 6000-9000 generations of doubly-charged scalar fields respectively. In the model of Ref. [194] the SM particle content is enlarged by a vector-like doubly-charged fermion E, a Majorana fermion N R , a singlet scalar S, a singly-charged scalar h + and N k generations of the doubly-charged scalar field k ++ . At the one-loop level the running of g 1 is governed by the renormalisation group equation (RGE) where t = log µ, µ being the renormalisation scale, and is the one-loop β function. It is clear from Eq. (2.3) that a large value of N k necessarily leads to a very steep increase of g 1 with the renormalisation scale, soon reaching a Landau pole. This is shown in Fig. 3, obtained by setting the masses of all the charged BSM states to µ NP = 2.5 TeV, which is already the largest mass considered in this analysis. The running up to µ NP is governed by the SM RGEs, and the result for g 1 is displayed with a black dashed line. For scales above µ NP = 2.5 TeV, the contributions from BSM fields become effective. Fig. 3 shows that a Landau pole can be reached at N k µ Landau 10 2 · 10 13 GeV 100 1.2 · 10 5 GeV 1000 3.8 TeV 6000 2.7 TeV 9000 2.6 TeV Table 2. Energy scale at which a Landau pole in g 1 appears as a function of N k in the model of Ref. [194].
relatively low energies once we allow for such large values of N k . In fact, for N k = 9000, we find that a Landau pole appears already at µ 2.6 TeV. In this specific example the appearance of a Landau pole below 10 16 GeV is unavoidable as soon as N k > 10, as shown in Table 2.  Figure 3. Running of the U (1) Y gauge coupling, g 1 , in the model presented in [194] for N k = 1000 (red), N k = 6000 (green) and N k = 9000 (blue). The black dashed line corresponds to the SM running below µ N P = 2.5 TeV.

How do the tools help?
The tools which we describe in more detail in the following sections can help to address all the above issues: 1. FlexibleSUSY and SPheno can calculate the diphoton and digluon rate including important higher order corrections.
2. Using the effective vertices calculated by FlexibleSUSY/SPheno and the interface to CalcHep or MadGraph, the gluon-fusion production cross-section of the 750 GeV mediator can be calculated numerically and one does not have to rely on analytical (and sometimes erroneous 4 ) approximations.
3. SARAH calculates the RGEs for a model which can be used to check for the presence of Landau pole. 4. Vevacious can be used to check the stability of the scalar potential.
2.2 Properties of the 750 GeV scalar

Mixing with the SM Higgs
It is often assumed that S, although it is a CP-even scalar, does not mix with the SM-like Higgs h. However, if this is done in a very ad hoc way and not motivated by any symmetry, this assumption will not hold when radiative corrections are taken into account. To see this, one can consider, for example, the scalar potential where H is the SM Higgs SU (2) L doublet, which contains the SM Higgs h. This potential in principle has all ingredients to get a large diphoton decay of S via a loop involving the charged scalar X. Note, however, that the potentially dangerous term κ H S|H| 2 has been omitted. One can see immediately that this term would be generated radiatively by the diagram below.
It is straightforward to see that the analytical estimate of the production cross section in Eq. (10) of Ref. [202] is wrong by orders of magnitude: consider the production of a SM-like scalar H with m H = 750 GeV via top-loops. Then, the factor h 2 F m 2 t /m 2 F drops out and one obtains σ = 1440 pb which is too large by roughly three orders of magnitude [277]. The authors of Ref. [179] (which originally made use of this analytic estimate) have revised their results in an updated version of their paper.
Note that it is also not possible to circumvent this decay by forbidding the λ HX term: since H and X are charged under SU (2) L × U (1) Y , also the λ HX |H| 2 |X| 2 term would be generated radiatively via diagrams like X X * H H * Such a mixing has important consequences since it opens the decay channels S → hh and S → ZZ, S → W + W − at tree-level, which are tightly constrained.
Another possibility is that all terms allowed by symmetries are taken into account, but very special relations among them are imposed like in Ref. [186]. When these relations hold, the above-mentioned tree-level decays in SM particles would cancel. However, as long as there is no symmetry behind these relations, they will not be invariant under RGE running. Therefore, immediately the question arises how large the tuning among the parameters must be to have a point that fulfils all constraints. To illustrate this issue, we make small variations in the couplings λ H3 and λ 36 , which cause non-vanishing tree-level couplings between S and the massive vector bosons, and check for which size of the deviations the condition Br(S → W + W − )/Br(S → γγ) < 20 holds. The result is shown in Fig. 4. Here, the diphoton rate was maximized by setting the masses of the vector-like fields to 375 GeV and using a Yukawa coupling of O(1). In principle, one could try to check what this means for the scale dependence of these ratios by calculating the RGEs. However, this cannot really be done for this setup since one obtains the following condition from the relations which have been imposed: Y is needed to maximize the diphoton branching ratio. Thus, f Y of O(1) immediately leads to a huge quartic coupling.
Thus, in general, it is very difficult to justify the assumption that the 750 GeV particles do not mix with the SM-like Higgs if there is no fundamental symmetry to forbid this mixing. However, this can already be forbidden using the CP symmetry: the mentioned problems can be circumvented in models where the diphoton excess stems from a CP-odd particle. While, in the case of a CP-even particle, it is crucial to include the mixing effects and to check at least how large the tuning in parameters must be.  Figure 4. The impact on small deviations from the parameter relations assumed in Ref. [186]. Specifically, the y-axis axis represents the deviation in the coupling λ 36 , while the x-axis represents deviations in the coupling λ H3 . The contour lines show the ratio Br(S → W + W − )/Br(S → γγ). The red line indicates where the model would already be in conflict with current collider limits.

To VEV or not to VEV?
The possibility that the new scalar receives a vacuum expectation value (VEV) is also often neglected. However, as we have just discussed, it often occurs that a H-S mixing will be induced, at least radiatively, in many models. Such radiative effects would immediately lead to a non-zero VEV for the new scalar. Even in cases where there is a symmetry which prevents a mixing with the SM Higgs, the 750 GeV particle will still receive a VEV if it is a CP-even scalar. This arises due to the introduced couplings to charged particles which are necessary to allow diphoton and digluon decays. More specifically, these introduced couplings will generate one-loop tadpole diagrams for S as shown in Fig. 5. Thus, the tadpole equation reads at the one-loop level where T (T ) is the tree-level tadpole, given by Figure 5. Tadpole terms for S generated at one-loop level.
Here, we have parametrised the tree-level expression so that the general form has the solution v S = 0. One finds in general that the one-loop corrections are Taking M Ψ , κ, M X of the order 1 TeV, results in a VEV which is naturally of order 1 TeV 3 /(16π 2 c 1 ). As a result, the simplifying assumption that v S vanishes is in general hard to justify. Therefore, it is important to check how the conclusions made about the model depend on this assumption. Here, the tools discussed in the following sections can really help, as including the non-vanishing v S is no more difficult than assuming the VEV vanishes.

Additional decay channels
Many analyses concentrate only on the decay S → γγ and completely neglect other potential decay channels. However, there are stringent constraints on the branching ratios of S into other SM final states, which are summarised in Table 3. Thus, any e + e − + µ + µ − τ + τ − Zγ ZZ Zh hh W + W − tt bb jj inv. 0.6 6 6 6 10 20 20 300 500 1300 400 Table 3. Limits on Γ(S → X)/Γ(S → γγ) assuming a production of S via gluon fusion or heavy quarks. Values are taken from Ref. [114].
model which tries to explain the excess via additional coloured states in the loop must necessarily worry about limits from dijet searches [278]. Therefore, an accurate calculation of the digluon decay rate is a necessity. As an example that illustrates why both additional channels and the diphoton/digluon width calculation are important we consider the model presented in Refs. [48,99] and considered in more detail here in Sec. 5.2.1.2.
This model extends the SM with a singlet and a scalar SU (2)-doublet colour octet. As an approximation the ratio of the singlet decays to gluons and to photons is (2.8) In [48] this is quoted as 750; before any NLO corrections are applied, we find 700. However, once we include all of the N 3 LO corrections this is enhanced to 1150, near the bound for constraints on dijet production at 8 TeV and significantly squeezing the parameter space of the model.
Additionally in many works we observed that potential decay channels of the resonance were missed. For instance in Ref. [108], the authors, who considered the Georgi-Machacek model [279], missed the decay of the scalar into W ± H ∓ , which can be the dominant mode when kinematically allowed.

How do the tools help?
Many of the assumptions which we criticised were made to keep the study simple. However, when using the public tools presented in the next two section, there is no need for these simplifying assumptions: 1. SARAH automatically calculates all expressions for the masses and vertices, no matter how complicated they are.

2.
FlexibleSUSY and SPheno give numerical predictions for the mass spectrum and the mixing among all states including higher order corrections.
Thus, for the user the study becomes no more difficult when he/she drops all simplifying assumptions but considers the model in full generality. Moreover, there is no chance to miss important effects in the decays of the new scalar: There are several studies which extend an already existing model by vector-like states and then assume that this part of the model is decoupled from the rest. When this assumption is made it is clear that the results from toy models, with the minimal particle content will be reproduced. However, it is often not clear if this decoupling can be done without invoking specific structures in the choice of parameters, and if these assumptions hold at the loop level.
On the other hand, if model-specific features are used to explain the diphoton excess, it is likely that there will be important constraints on the model coming from other sectors. For instance, there might be bounds from flavour observables, dark matter, Higgs searches, neutrino mixing, electroweak precision observables, searches for BSM particles at colliders, and so on. All of that has to be checked to be sure that any benchmark point presented is indeed a valid explanation for all observations. Such a wide range of constraints is much easier to address by making exhaustive use of tools which provide a high level of automation.

Theoretical uncertainties of other predictions
Even if the attempts are made to include the effects of the new states on other sectors of the model, it is important to be aware that there are large uncertainties involved in certain calculations. If the level of uncertainty is underestimated, this can have an impact on what is inferred from the calculation. The large uncertainty in a LO calculation of the diphoton and digluon rate has already been addressed in Section 2.1.1. However, there are also other important loop corrections especially in SUSY models: the accurate calculation of the Higgs mass is a long lasting endeavour where for the simplest SUSY models even the dominant three-loop corrections are partially tackled [280]. The current ball-park of the remaining uncertainty is estimated to be 3 GeV.
However, it is clear that the MSSM cannot explain the excess, hence it must be extended. A common choice is to add additional pairs of vector-like superfields together with a gauge singlet, see Section 6. These new fields can also be used to increase the SMlike Higgs mass. However, this will in general also increase the theoretical uncertainty in the Higgs mass prediction, because these new corrections are not calculated with the same precision as the MSSM corrections. For instance, Ref. [105] has taken into account the effect of the new states on the SM-like Higgs. There, they use a one-loop effective potential approach considering the new Yukawa couplings to be O(1) or below, while also including the dominant two-loop corrections from the stop quark. They assumed that including these corrections is sufficient in order to achieve an uncertainty of 2 GeV in the Higgs mass prediction. One can compare their results encoded in Fig. 7 of Ref. [105] with a calculation including, in addition to the corrections taken into account in the paper, momentum dependence and electroweak corrections at the one-loop level, as well as the additional two-loop corrections arising from all the newly introduced states. These corrections can be important, as was shown for instance in Ref. [281]. The result of the comparison is shown in Fig. 6. We find a similar behaviour,  Figure 6. Comparison of the two-loop Higgs mass calculation of Ref. [105] with the results obtained by SPheno. The parameters are those of Fig. 7 in Ref. [105] and we fixed X κ 10 = 0. The lines are the results from SPheno, while the green and red shades areas are the ranges of κ 10 which predict m h = [123,127] GeV according to Ref. [105]. Red takes X t = 4 and green X t = 2, where X t is defined in Eq. (9) of the reference.
but still there are several GeV difference between both calculations. For κ 10 = 0.8 and X t = 4, the point would be within the interesting range for m h = [123,127] GeV, while the more sophisticated calculation predicts a mass below 120 GeV. Thus, the assumed uncertainty of 2 GeV in Ref. [105], which would even be optimistic in the MSSM, is completely unrealistic without including all the aforementioned higher order corrections.

How do the tools help?
The tools help to ensure that one really considers all aspects of a full model: 3 Using SARAH to understand a model

General
One of the reasons that makes high energy particle physics an exciting field is the vast amount of experimental data available. When proposing a model one first has to check its self consistency, checking for instance the particle mass spectrum and vacuum stability requirements. Then it has to be tested against data related to collider searches, flavour observables, dark matter observations and Higgs measurements. A lot of effort has been devoted to developing an arsenal of specific tools to explore these quantities with high precision for specific classes of models, such as the MSSM, the THDM and the NMSSM to some extent. However, none of them, in their simplest versions, can explain the diphoton excess. For the time being there is no specific model which is clearly preferred over others as an explanation of the excess, as reflected by the large variety of models that different authors have proposed, and it would be impractical to repeat the process of developing a code for each one of them. In the absence of a dedicated tool, the alternative is often to resort to approximations or just to leading order expressions, as described in the previous section, in which case the analysis (in particular for more complicated models) is of limited value. Luckily, a dedicated powerful tool already exists. It is the Mathematica package SARAH [216][217][218][219][220][221], which can perform the most advanced quantum field theory computations and apply them generically to any given model. SARAH has been optimised for an easy, fast, and exhaustive study of renormalisable BSM models: not only can it calculate all the relevant quantities within a given model analytically -it also provides routines to export the derived information in order to use it for numerical calculations with dedicated tools. SARAH can be used for SUSY and non-SUSY models to write model files for CalcHep/CompHep [234,235], FeynArts/FormCalc [282,283], WHIZARD/O'Mega [284,285], as well as in the UFO format [286] which can be handled for instance by MadGraph 5 [236], GoSam [287], Herwig++ [288][289][290], and Sherpa [247,291,292]. The modules created by SARAH for SPheno calculate the full one-loop and partially two-loop corrected mass spectrum, branching ratios and decay widths of all states, and many flavour and precision observables thanks to the FlavorKit [293] functionality. Moreover, an easy link to HiggsBounds [239,240] and HiggsSignals [241] exists. One can as well use another tailor-made spectrum generator, FlexibleSUSY [224], which can handle both SUSY and non-SUSY models generated with SARAH. Finally, SARAH can also produce model files for Vevacious [242], which is a tool dedicated to studying vacuum stability.

Models supported by SARAH
SARAH is optimised for handling a wide range of SUSY and non-SUSY models. The basic idea of SARAH is to give the user the possibility to implement models in an easy, compact and straightforward way. The user simply has to declare symmetries, particle content, and the (super)potential. All steps to derive the full Lagrangian for the gauge eigenstates from this information are then fully automatised. But this is not the final aim: usually one is interested in the mass eigenstates after spontaneous gauge symmetry breaking. To perform the necessary rotations to the new eigenstates, the user has to provide additional information: (i) definition of the fields which develop VEVs, and (ii) definitions of the fields which mix. Using this information, all necessary re-definitions and field rotations are performed by SARAH. In addition, gauge fixing terms for the new eigenstates are derived and ghost interactions are added. Plenty of information can be derived by SARAH for all states in the model, as explained in Section 3.5. Before moving into it, we elaborate on the kind of models and features supported by SARAH.

Global symmetries
SARAH can handle an arbitrary number of global symmetries which are either Z N or U (1) symmetries. In a SUSY model one can also impose a continuous R-symmetry, U (1) R , at the level of the superpotential. However, SARAH will then generate automatically the scalar potential including soft trilinear terms that do not respect the U (1) R . The user has two options: either add "AddTterms = False" to the SARAH model file to forbid the automatic generation of the soft trilinear terms, or set such terms to zero manually within the input files for the spectrum generators. SARAH can also handle approximate symmetries: if the user specifies terms in the Lagrangian/superpotential which violate a global symmetry, the terms will not be forbidden but a warning will be printed.

Gauge sector
Gauge groups SARAH is not restricted to the SM gauge sector, but many more gauge groups can be defined. To improve the power in dealing with gauge groups, SARAH has linked routines from the Mathematica package Susyno [294]. SARAH, together with Susyno, takes care of all group-theoretical calculations: it calculates the Dynkin and Casimir invariants, derives the needed representation matrices as well as the generalised Clebsch-Gordan coefficients. For all Abelian groups one can also give a GUT normalisation; these factors usually arise from considerations about embedding a model in a larger symmetry group such as SU (5) or SO (10).
Gauge kinetic mixing With more than one Abelian gauge group, terms of the form are allowed, where F µν are the field strength tensors of two different Abelian groups A, B [295]. κ is in general an n × n matrix if n Abelian groups are present. SARAH fully includes the effects of kinetic mixing for any number of Abelian groups by performing a rotation to bring the field strength to a diagonal form. The kinetic mixing is then absorbed into the covariant derivatives, which take the form The indices x, y run over all U (1) groups, and g xy are the entries of the gauge coupling matrix G which now also contains mixed gauge couplings. Gauge kinetic mixing is not only included in the interactions with vector bosons, but also in the derivation of the D-terms and in the gaugino interactions for SUSY models.

Matter sector
One can define up to 99 matter fields in a single model in SARAH. Each one of them can come with an arbitrary number of generations and can transform as any irreducible representation with respect to the defined gauge groups.
Supersymmetric models The matter interactions in SUSY models are usually fixed by the superpotential and the soft SUSY-breaking terms. SARAH takes as input the renormalisable terms in the superpotential which the user has to write in the model file, and automatically generates the corresponding soft-breaking terms c L , c M , c T are real coefficients, while the linear, bilinear, and trilinear parameters are treated by default in the most general way by taking them as complex tensors of appropriate order and dimension. If identical fields are involved in the same coupling, SARAH also derives the symmetry properties for the parameter.
In recent years models with Dirac gauginos have been largely explored. They feature mass terms of the form mφ i A D λ A ψ i , where λ A is a gaugino and ψ i the fermionic component of a chiral superfieldφ i in the adjoint representation of the gauge group A. In addition, there are new D-term couplings (see e.g. [296]). To generate Dirac mass terms for all the gauginos, these models always come with an extended matter sector, including at least one singlet, one triplet under SU (2), and one octet under SU (3). Furthermore, these models generate new structures in the RGEs [297]. All these are fully supported in SARAH.
Non-Supersymmetric models For non-supersymmetric models, SARAH supports all general, renormalisable Lagrangians of the form for scalars φ i , and Weyl fermions ψ j . The Lagrangian needs to be defined by the user in the model file. Note, that we have omitted a tadpole term, tφ, for a gauge singlet, as it can always be absorbed in a shift of φ.

Checks of implemented models
After the initialisation of a model, SARAH provides functions to check its (self-)consistency: • Check for gauge anomalies, and mixed gauge/gravity anomalies; • Check for Witten anomalies [298]; • Check if all terms in the (super)potential are in agreement with all global and gauge symmetries; • Check if terms allowed by all symmetries are missing in the (super)potential; • Check if additional mass eigenstates can mix in principle; • Check if all mass matrices are irreducible.
In addition, SARAH performs other formal checks. For instance, it checks if the number of Particle Data Group numbers for a given family of particles fits to the number of generations for each particle class, if L A T E X names are defined for all particles and parameters, and if the position in a Les Houches spectrum file is defined for all parameters.

Analytical calculations performed by SARAH
The full power of SARAH can be unleashed on exhaustive numerical analyses via the dedicated interfaces to other tools, discussed in the next sections. However, since SARAH itself is a Mathematica package, it is capable of many analytical computations beyond the derivation of the Lagrangian. Here we list some of them.

Tree-level properties
Tadpole equations During the evaluation of a model, SARAH calculates 'on the fly' all the minimisation conditions of the tree-level potential, the so-called tadpole equations.
Masses SARAH calculates the mass matrices for the states which are rotated to mass eigenstates. In addition, it calculates the masses of states for which no field rotation has taken place.
Vertices SARAH has functions to extract in an efficient way all tree-level vertices from the Lagrangian. These vertices are saved in different Mathematica arrays according to their generic type.

Renormalization group equations
SARAH calculates the full two-loop RGEs in SUSY and non-SUSY models including the full CP and flavour structure. For this purpose, it makes use of the most sophisticated generic calculations available in the literature.

SUSY RGEs
The calculation of the SUSY RGEs is mainly based on Ref. [299], which however did not cover all possible subtleties which can appear in SUSY models. SARAH has implemented also results from more recent literature: • In the case of multiple U (1) gauge groups, gauge-kinetic mixing can arise if the groups are not orthogonal. Substitution rules to translate the results of Ref. [299] to those including gauge kinetic mixing were presented in Ref. [300] and are used by SARAH; Non-SUSY RGEs SARAH uses the expressions of Refs. [303][304][305][306] for the calculation of the RGEs in a general quantum field theory. These results are completed by Ref. [307] to cover gauge kinetic mixing and again by Refs. [301,302] to include the gaugedependence of the running VEVs also in the non-SUSY case.
3.5.3 One-and two-loop corrections to tadpoles and self-energies One-loop corrections SARAH calculates the analytical expressions for the one-loop corrections to the tadpoles and the one-loop self-energies for all the particles. For states which are a mixture of several gauge eigenstates, the self-energy matrices are calculated. The calculations are performed in the DR scheme using 't Hooft gauge for SUSY models. In the case of non-SUSY models SARAH switches to the MS scheme.
Two-loop corrections It is even possible to go beyond one loop with SARAH and to calculate two-loop contributions to the self-energies of real scalars. There are two equivalent approaches implemented in the SPheno interface of SARAH to perform these calculations: an effective potential approach and a diagrammatic approach with vanishing external momenta. More details about these calculations are given in Section 4.1.1.

Spectrum calculation, Monte-Carlo studies, and more
As mentioned in the last section, SARAH can use the analytical information derived about a model and pass it to other tools. We give in the following an overview about the different possibilities.

SPheno
SARAH writes Fortran source code for SPheno [222,223]  One-loop shifts to pole masses The one-loop mass spectrum is calculated from the information about the one-loop self-energies and tadpole equations. The procedure is a generalisation of the one explained in detail for the MSSM in Ref. [308]. The main features are • Any one-loop contribution in a given model to all fermions, scalars and vector bosons is included • The full p 2 dependence in the loop integrals is included • An iterative procedure is applied to find the on-shell masses m(p 2 = m 2 ).
Two-loop shifts to Higgs pole masses SARAH can also generate Fortran code to calculate the two-loop corrections to the masses of the CP-even scalar states with SPheno 6 . The same approximations usually taken for the MSSM are also applied here: (i) all calculations are performed in the gaugeless limit, i.e. the electroweak contributions are dropped, and (ii) the momentum dependence is neglected. Using these routines, the theoretical uncertainty in the Higgs mass prediction for many models has been shrunk to the level of the MSSM. In general, there are two different techniques to calculate the two-loop corrections with SARAH-SPheno: • Effective potential calculation [232]: SARAH makes use of the generic two-loop results for the effective potential given in Ref. [309]. To get the values for the two-loop self-energies and two-loop tadpoles, the derivatives of the potential with respect to the VEVs are taken numerically, as proposed in Ref. [310].
• Diagrammatic calculation [233]: A fully diagrammatic calculation for twoloop contributions to scalar self-energies with SARAH-SPheno became available with Ref. [233]. The advantage of the diagrammatic approach is that no numerical derivation is needed. This is now the default calculation.
The implementation of the two independent approaches provides a good possibility to double check results. It has been shown that these generic calculations can provide important two-loop corrections for the NMSSM which are not included in dedicated spectrum generators for the NMSSM [225,311].

Decay widths and branching ratios
SPheno modules created by SARAH calculate all two-body decays for fermion and scalar states as well as for the additional gauge bosons. In addition, the three-body decays of a fermion into three other fermions, and of a scalar into another scalar plus two fermions are included. In the scalar sector, possible decays into two particles are calculated at tree level. In case of two quarks in the final state, the dominant QCD corrections due to gluons are included [274]. In addition scalar decays into final states with off-shell gauge bosons (ZZ * , W W * ) are included.
The loop-induced decays into two photons and gluons are calculated up to N 3 LO. More details about this are given in Section 4.4.

Flavour observables
The SPheno modules written by SARAH contain out-of-the-box routines to calculate many quark and lepton flavour violating observables: • Lepton flavour violation: -CR(µ − e, N ) (N=Al,Ti,Sr,Sb,Au,Pb), Br(τ → P ) (with P =π, η,η ) • Quark flavour violation: The calculation is based on the FlavorKit functionality [293], which makes use of the chain FeynArts-SARAH-SPheno. This provides a full one-loop calculation in a given model. In addition, this interface can be used to derive Wilson coefficients for new operators, and to calculate new observables with SPheno using the implemented coefficients.

Fine-Tuning
A widely used measure for the electroweak fine-tuning was proposed in Refs. [312,313] Here, α is a set of independent parameters, and ∆ −1 α is a measure of the accuracy to which the parameters α must be tuned to obtain the correct VEV. The user can choose the set of parameters α in SARAH, and SPheno numerically calculates ∆ F T for that choice using the full two-loop RGEs from the GUT or SUSY breaking scale down to the electroweak scale.

Production cross-sections
SPheno provides an estimate for the production cross sections of all neutral scalars within a certain mass range: values for gluon-fusion and vector-boson fusion are obtained for m ∈ [50, 1000] GeV, while associated production with W , Z and t are considered for m ∈ [50,300] GeV. The results are obtained by re-weighting the SM cross sections with the effective coupling of the scalar in the model to SM states normalised to the SM values. For 7 and 8 TeV the SM cross-sections provided by the Higgs cross section working group are used, while for 13, 14 and 100 TeV the cross sections have been calculated with SusHi 1.5.0 [314].

FlexibleSUSY
FlexibleSUSY is a Mathematica package which uses the SARAH-generated expressions for the mass matrices, self-energies, tadpole equations, vertices and RGEs to create a C++ spectrum generator for both SUSY and non-SUSY models. The spectrum generators created with FlexibleSUSY have the following features: • an interface to GM2Calc [273] in MSSM models without sfermion flavour violation FlexibleSUSY aims to generate spectrum generators which are modular such that components can be easily reused or replaced. This means that it is quite easy to re-use the precision calculations in FlexibleSUSY spectrum generators for other purposes or add additional routines.

Mass calculation with FlexibleSUSY
Determination of the gauge couplings, Yukawa couplings and the Standard Model VEV FlexibleSUSY generates routines allowing for a one-loop calculation of the SU (3) c × SU (2) L × U (1) Y gauge couplings, if they exist in the model under consideration. Furthermore, the Yukawa couplings corresponding to the Standard Model fermions can be calculated in the considered model at the full one-loop level from the known fermion masses. To determine the top quark Yukawa coupling, two-loop Standard Model QCD corrections are also added. FlexibleSUSY also performs a complete one-loop calculation of the running Z and W masses at the low-energy scale in the considered model, which can be used to determine the running Standard Model-like vacuum expectation value, v.
One-loop shifts to pole masses FlexibleSUSY by default performs a full oneloop MS/DR calculation to determine the pole masses of all particles in the model, similar to the procedure presented in Ref. [308] for the MSSM. Thereby, it makes use of the one-loop self-energies and tadpole diagrams generated by SARAH, taking the full momentum-dependence into account. To tune the spectrum generator, the user can choose from three different precision levels, which differ in the way two-loop momentum dependent terms are treated. By default, the Higgs masses are calculated with highest precision, where an iteration over the momentum is performed to determine the pole mass M h at p 2 = M 2 h .
Two-loop and three-loop shifts to Higgs pole masses FlexibleSUSY allows the user to add certain predefined two-loop corrections to the Higgs masses in some specific models: In MSSM-like models the two-loop corrections of the order O((α t + α b ) 2 + α t α s + α b α s + α 2 τ ) from [226][227][228][229][230] can be added to the two CP-even and one CP-odd Higgs boson. In NMSSM-like models with three CP-even and two CP-odd Higgs bosons the two-loop corrections O(α t α s + α b α s ) from [231] plus the MSSM-like contributions O((α t + α b ) 2 + α 2 τ ) [227,230] can be added. In SM-like models with one physical Higgs singlet the two-loop corrections O(α t α s + α 2 t ) from Refs. [315] can be added to the self-energy. In the split-MSSM the three-loop gluino contribution O(α t α 2 s ) from Ref. [316] can be added.

Decay widths and branching ratios
Like SPheno, FlexibleSUSY can calculate the loop-induced Higgs decays into two photons and two gluons up to NNNLO, see Section 4.4.

Mass spectrum calculation: SUSY vs. Non-SUSY
We have outlined that FlexibleSUSY and SPheno can include the radiative corrections to all particles up to the two-loop level in the DR scheme. These corrections are included by default for supersymmetric models. It is known that loop corrections, in particular to the Higgs mass, are crucial. Typically the DR and on-shell calculations are in good agreement. Consequently, the remaining difference between both calculations is often a good estimate for the theoretical uncertainty.
The treatment of non-supersymmetric models in FlexibleSUSY and SPheno is very similar to the treatment of supersymmetric models. The main difference is, that in non-supersymmetric models the parameters are defined in the MS scheme, while in supersymmetric ones the parameters are defined in the DR scheme. In this paper we perform only tree-level mass calculations (if not stated otherwise), in which the definition of the renormalisation scheme is irrelevant. Thus, in the mass spectrum calculations performed in the following, one is allowed to use input parameters which are defined in the on-shell scheme. This is for instance the standard approach in the large majority of studies of the THDM: there are in general enough free parameters to perform a full on-shell renormalisation keeping all masses and mixing angles fixed. We find that the one-loop corrections in the MS scheme can give huge corrections to the tree-level masses in nearly all models presented in the following. Therefore large finetuning of the parameters is necessary once the loop corrections are taken into account. A detailed analysis using a full on-shell renormalisation scheme is possible for each model, but is beyond the scope of this work. Of course, for models where it turns out to be unavoidable that shifts in the masses and mixings appear at the loop-level, the user can simply turn on the loop corrections in FlexibleSUSY and SPheno via a flag in the Les Houches input file.

Calculation of the effective diphoton and digluon vertices in SPheno and FlexibleSUSY
For the calculation of the partial width of a neutral scalar Φ decaying into two gluons or two photons we follow closely [274] for the LO and NLO contributions. The partial widths at LO are given by for X = Φ, A. µ NLO is the renormalisation scale used for these NLO corrections, chosen to be µ NLO = m Φ /2 . In the intermediate range of 100m f > m Φ > 2m f , no closed expressions for the NLO correction exist. Our approach in this range was to extract the numerical values of the corrections from HDECAY [261] and to fit them. For the digluon decay rate, the corrections up to N 3 LO are included and parametrised by  For pseudoscalar we include only corrections up to NNLO as the N 3 LO are not known for CP-odd scalars. In order to check the accuracy of our implementation, we compared the results obtained with SARAH-SPheno for the SM Higgs boson decays with the ones given in the CERN yellow pages [323]. In Fig. 7 we show the results for the Higgs branching ratios into two photons and two gluons with and without the inclusion of higher order corrections. One sees that good agreement is generally found when including higher order corrections.
In Fig. 8 we show the ratio Br(h → gg)/Br(h → γγ) and compare it again with the recommended numbers by the Higgs cross section working group [323]. Allowing for a 10% uncertainty, we find that our calculation including higher order corrections agrees with the expectations, while the LO calculation predicts a ratio which is over a wide range much too small. The important range to look at is actually not the one with m h ∼ 750 GeV because this corresponds to a large ratio of the scalar mass compared to  the top mass. Important for most diphoton models is the range where the scalar mass is smaller than twice the quark mass. In this mass range we find that the NLO corrections are crucial and can change the ratio of the diphoton and digluon rate up to a factor of 2. We also note that if we had used α(m h ) instead of α(0) in the LO calculation, the difference would have been even larger, with a diphoton rate overestimated by a factor (α(m h )/α(0)) 2 (137/124) 2 1.22.

Monte-Carlo studies
SARAH writes all necessary files to implement a new model in different MC tools. First, a short summary of the output formats is given . Then, it is described how parameter values between SPheno/FlexibleSUSY and the MC codes can be exchanged.

Output of models files
CalcHep The model files for CalcHep [234,235] can also be used with MicrOmegas [238] to perform both collider and dark matter studies. CalcHep is able to read the numerical values of the masses and mixing matrices from a SLHA spectrum file based on the SLHA+ functionality [324]. SARAH makes use of this functionality to generate model files by default in a way that they automatically expect to find the input values in spectrum files written by a SARAH-generated SPheno version or by FlexibleSUSY. However, other choices are possible: the parameters can be given via the vars file or tree-level expressions can be calculated internally by CalcHep.   [286]. This is particularly useful to implement new models in MG5_aMC@NLO and perform collider studies. The UFO format is also supported by other tools like GoSam [287], Aloha [325], Herwig++ [288,289] and MadAnalysis 5 [326]. Moreover, the spectrum file written by SPheno or by FlexibleSUSY can be directly used as parameter card in MadGraph.
WHIZARD/O'Mega SARAH writes all necessary files to implement a model in WHIZARD and O'Mega [284,285]. Since the SLHA reader of WHIZARD is at the moment restricted to the MSSM and the NMSSM, SPheno versions generated by SARAH can write all information about the spectrum and parameters in an additional file in the WHIZARD specific format. This file can then be read by WHIZARD. Currently, the handling of general Lorentz structures in WHIZARD and the support of the UFO format are under development. This will provide the possibility to use WHIZARD with the calculated diphoton and digluon vertices as explained in the following.

Interplay SARAH-Spectrum-Generator-MC-Tool
The tool chains SARAH-SPheno/FlexibleSUSY-MC-Tools have one very appealing feature: the implementation of a model in the spectrum generator (SPheno or FlexibleSUSY) as well as in a MC tool is based on just one single implementation of the model in SARAH. Thus, the user does not need to worry that the codes might use different conventions to define the model. In addition, SPheno also provides all widths for the particles so that this information can be used by the MC-Tool to save time.

Effective diphoton and digluon vertices for neutral scalars
The effective diphoton and digluon vertices calculated by SPheno or FlexibleSUSY are directly available in the UFO model files and the CalcHep model files: SARAH includes the effective vertices for all neutral scalars to two photons and two gluons, and the numerical values for these vertices are read from the spectrum file generated with SPheno or FlexibleSUSY. For this purpose, a new block EFFHIGGSCOUPLINGS is included in these files, which contains the values for the effective couplings including all corrections outlined in Section 4.4.
It is important to mention that these effective couplings correspond to the decay of the scalar; if we use CalcHep or MadGraph to compute the decay Φ → gg then the value matches (as closely as possible) the NNNLO value, which includes real emission processes such as Φ → ggg. Therefore, the corrections at NLO and beyond for Φ → gg are not the same as pp → Φ via gluon fusion [317]; the full NNNLO production crosssection includes all processes gg → Φ + jet and is therefore described by a different k-factor to the decay. This k-factor can for instance be obtained via where c Φgg is the ratio squared of the effective coupling between Φ and two gluons at LO in the considered model and the SM. These values can for instance be read off by the block HiggsBoundsInputHiggsCouplingsBosons in the SPheno spectrum file. σ SM (pp → H(M Φ )) is the cross section for a SM-like Higgs with mass M Φ . This value can be calculated for instance with Higlu [327] or Sushi [314] for the considered center-of-mass energy. SPheno also provides values for c Φgg · σ SM (pp → H(M Φ )) for the most common energies in the blocks HiggsLHCX (X=7, 8,13,14) and HiggsFCC, see also Section 4.1.5. On the other hand, this approach is not entirely appropriate for more refined collider analyses where the user would like to actually include, for example, a hard jet in the final state (without the full loop corrections to the effective vertex this is not an infra-red safe quantity). In this case, we note that around 750 GeV the effective vertex output by SARAH gives a fairly accurate result -to within 30% -of the total production cross-section, at least in the Standard Model, when we compute σ SM (gg → Φ + jet) using MadGraph and the standard cuts on momenta. This is illustrated in figure 9.  Figure 9. Comparisons of the total Higgs production cross-section via gluon fusion in the Standard Model as a function of the Higgs mass, computed using the SPheno output from SARAH, the Higgs cross-section working group data, and in MadGraph using our effective vertex.

Checking Higgs constraints
HiggsBounds [239,240,328] and HiggsSignals [241] are dedicated tools to study the Higgs properties of a given parameter point in a particular model. While HiggsBounds checks the parameter point against exclusion limits from Higgs searches, HiggsSignals gives a χ 2 value to express how well the point reproduces the Higgs measurements. In general, HiggsBounds and HiggsSignals can handle different inputs: either the cross sections for all necessary processes can be given at the parton or hadron level, or the effective couplings of the Higgs states to all SM particles are taken as input.
In addition, the masses and widths of all neutral as well as charged Higgs states are always needed. SPheno provides all input for the effective coupling approach. The information is distributed in the SLHA spectrum file and in separated files (called MH_GammaTot.dat, MHplus_GammaTot.dat, BR_H_NP.dat, BR_Hplus.dat, BR_t.dat, effC.dat). While SLHA files can be used with HiggsBounds for models with up to five neutral scalars, the separated files can even be used with up to 99 neutral and 99 charged scalars.

Checking the Vacuum stability
Vevacious is a tool to check for the global minimum of the one-loop effective potential allowing for a particular set of non-zero VEVs. For this purpose, Vevacious first finds all tree-level minima using HOM4PS2 [329]. Afterwards, it minimizes the one-loop effective potential starting from these minima using minuit [330]. If the input minimum turns out not to be the global one, life-time of meta-stable vacua can be calculated using CosmoTransitions [331].
Vevacious takes the tadpole equations, the polynomial part of the scalar potential and all mass matrices as input. All of this information has to be expressed including all VEVs which should be tested. That means that, in case of supersymmetric models, in order to check for charge and colour-breaking minima, both the stop and the stau VEVs must be taken into account everywhere in the scalar potential, the tadpole equations and the mass matrices. Moreover, the possible mixing of all states triggered by the new VEVs must be included. To take care of all that, the corresponding input files can be generated by SARAH.

Accuracy of the diphoton calculation
Before concluding this section, we should draw the reader's attention to the question of how accurate the results are from SARAH in combination with SPheno and FlexibleSUSY. While every possible correction has been included, there are still some irreducible sources of uncertainty, as we shall discuss below.

Loop corrections to ZZ, W W, Zγ production
So far in SARAH, loop-level decays are only computed for processes where the tree-level process is absent. This is to avoid the technical issues of infra-red divergences. If the particle that explains the 750 GeV excess is a scalar, then it must mix with the Higgs and acquire tree-level couplings to the Z and W bosons, and these are fully taken into account. However, due to the existence of such terms, the loop corrections to the decays into Zs and W s are more complicated and are therefore not yet available in SARAH.
To estimate the uncertainty incurred by their absence, let us assume that our 750 GeV resonance S couples to the U (1) Y and SU (2) L gauge bosons via the effective operators SB µν B µν and SW µν W µν . If we can neglect the tree-level contributions to the decays and assume that the dominant contribution originates from a set of particles in the loops, which have the hypercharge Y and the dynkin index D 2 (W ) and dimension of the SU (2) representation d 2 , then the decay widths are approximately given by Figure 10. Examples of potentially important NLO corrections.
where we abbreviated t W for tan θ W . Put together, the uncertainty that we find for the decay S → γγ reads The factor in square brackets is therefore largest for fields that only couple to SU (2) L gauge bosons, giving a factor of ∼ 55, and for SU (2) doublets with hypercharge 1/2 it is 13, although the former case yields too many W bosons (the limit from run 1 searches . Thus, provided that Br(S → γγ) 10 −3 , the relative uncertainty is guaranteed to be less than 10%. In such cases, the proportional error in the total width transfers directly into the proportional error in the total cross-section: On the other hand, for models where the dominant decay channel of the singlet is into gluons, it is not possible to have Br(S → γγ) 10 −3 without violating constraints from dijet production, and the reader should be careful about the possible errors incurred. Fortunately, provided that the loop particles have a hypercharge the error is much smaller, for example in the case that D 2 = 0 the coefficient above is less than one, thus giving an error of ∼ 10 −3 for Br(S → γγ) = 10 −3 .

BSM NLO corrections
As discussed above, SARAH includes the leading-order computation of the diphoton and digluon decay amplitudes including the effects of all Standard-Model and Beyond-the-Standard-Model particles in the loops. Furthermore, it also includes the leading-log corrections to the digluon rate at NLO, NNLO and N 3 LO order in α s in the Standard Model, and some NLO corrections due to diagrams with an extra gluon to both the digluon and diphoton rates. However, the NLO corrections are absent for all other particles, which in the case of large Yukawa couplings or hierarchies could be sizeable. Two examples of such a diagrams are given in Fig. 10; in the context of supersymmetric theories, particularly important are diagrams involving the gluino, which (if it is a Majorana particle) would not couple to a singlet at leading order -naively their contribution is although as we shall discuss below this can be (potentially significantly) an underestimate.

Tree vs pole masses in loops
For consistency of the perturbative series and technical expediency, the masses inside loops (to calculate pole masses and loop decay amplitudes) are MS or DR parameters, not the pole masses of observed particles. The difference between calculations performed in this scheme and the on-shell scheme are at two-loop order, and so is generally small. However, in particular when there are large hierarchies or Yukawa couplings in a model, there can be a large difference between the Lagrangian parameters and the pole masses, and therefore a large discrepancy between the loop amplitudes calculated from these. In principle, this should be accounted for by including higher-order corrections such as the right-hand diagram in Fig. 10 -but applying such a correction to each propagator in the loop would actually correspond to a four-loop diagram. The effect of using the pole mass instead is to essentially resum part of these diagrams, which as is well known is relevant in the case of large hierarchies of masses -and so should give a more accurate result in that case. If we define then for one dominant particle p in the loop, we can estimate the uncertainty as where the factor of 1 or 1/2 assumes that the couplings C ppΦ do not depend upon the mass m p (but the prefactor in r Φ p therefore does). For most values of m p the loop functions are slowly changing (only peaking around τ p = 1) so we will have a proportional uncertainty in the result of As an example, in supersymmetric theories the soft masses of coloured scalarsS acquire a significant decrease from gluino loops: If the scalar is a colour triplet with a pole mass of 800 GeV, then for 2 TeV gluinos and the DR mass is ∼ 1100 GeV; but ∼ 1! This corresponds to a shift of a factor of two in the amplitude, and, if the scalar dominates the total amplitude, a factor of four in Γ(S → gg/γγ); in fact in this case SARAH would be potentially underestimating the diphoton rate. This is a relatively mild example regarding this excess: given that the vast majority of models proposed to explain the diphoton signature contain large Yukawa couplings and many new particles, there is a significant potential for large values of δm 2 m 2 , about which the user should be careful. It is worth noting that this is an effect that would not be significant in the (N)MSSM, where the Higgs couplings to photons/gluons are dominated by the top quark (and, for photons, the W bosons) whose masses are protected by chiral symmetry from large shifts: this issue is a novelty for the 750 GeV excess.
For non-supersymmetric models, due to the fact that (almost) every parameter point is essentially fine-tuned, we have not calculated loop corrections to the masses by default, and this issue does not arise in the same way. The user is then free to regard the result as involving the pole masses of particles instead, if they so desirethe issue then becomes one of tuning the potentially large corrections to the other input parameters.

Models
A large variety of models has been proposed to explain the diphoton excess at 750 GeV. We have selected and implemented several possible models in SARAH. Our selection is not exhaustive, but we have tried to implement a sufficient cross-section which are representative of many of the ideas put forward in the context of renormalisable models. These are the ones that SARAH can handle. Their description is organised in the subsections that follow. Before we turn to this discussion we first want to mention other proposals which we do not deal with in this paper.
Many authors [10,16,20,25,29,35,39,57,80,89,112,114,121,150,152,157,175,189,198,204] have studied the excess with effective (non-renormalisable) models, which is sensible given that there are thus far no other striking hints of new physics at the LHC. As more data becomes available and the evidence for new physics becomes more substantial, one might want to UV complete these models, at which point the tools we are advertising become relevant and necessary. Other authors [23,26,30,40,69,73,78,86,139,140,173,176,182,187,190,191,206] considered strongly coupled models, in which the resonance is a composite state. This possibility would be favoured by a large width of the resonance, as first indications seem to suggest. Another possibility is to interpret the signal in the context of extra-dimensional models [3,5,17,18,24,47,74,85,120,132,185], with the resonance being a scalar, a graviton, a dilaton, or a radion, depending on the scenario. However, some of these interpretations are in tension with the non-observation of this resonance in other channels. In Supersymmetry, the scalar partner of the goldstino could provide an explanation to the diphoton signal [56,91,98,201]. Other ideas, slightly more exotic, include: a model with a space-time varying electromagnetic coupling constant [81], Gluinonia [203], Squarkonium/Diquarkonium [177], flavons [148], axions in various incarnations [4,14,31,149,164,202], a natural Coleman-Weinberg theory [13,181], radiative neutrino mass models [158,193,194], and string-inspired models [12,79,110,144,153].
We turn now to the weakly coupled models, and discuss the ones which we implemented in SARAH. All model files are available for download at http://sarah.hepforge.org/Diphoton_Models.tar.gz and an overview of all implemented models is given in Tabs. 4 and 5. In case of questions, comments or bug reports concerning these models, please, send an e-mail to diphoton-tools@cern.ch which includes all authors.

Toy models
The simplest ideas proposed to explain the diphoton excess extend the SM by a scalar singlet and vector-like fermions, which serve the purpose of enhancing the diphoton rate, and -in the case of coloured states -also the production via gluon fusion. An enhancement of gluon fusion seems to be necessary because a production of the resonance purely by photon fusion is in some tension with 8 TeV data: the increase in the cross section from 8 to 13 TeV is just a factor 2, while a factor of 5 would be needed to make the results from LHC run-I and II compatible. To explore the multitude of possibilities we first consider toy models: they do not contain all possible couplings of the vector-like fermions to SM matter, but are rather engineered to allow one to easily explore the effects of different representations of vector-like matter on the diphoton rate and on the relevant partial widths. Then, from section 5.2 on, we consider complete models, containing all the operators consistent with both field content and symmetries.

Name Section
Ref.
Toy models with vector-like fermions CP-even singlet SM+VL/

Toy models with vector-like fermions
• Reference: [109,133,134,166] • Model names: To begin with we categorise the toy models according to the type of the involved scalar singlet. There are three possibilities: (i) the singlet is a real CP-even scalar, (ii) Model Name Section Ref.   real CP-odd, or (iii) a complex scalar. Each case is considered in separate SARAH model files, where we introduce all possible representations of vector-like fermions. These possibilities, following Tables 3 and 4 of Ref. [166], are shown below in Table 6. This allows one to study combinations of fermion representations or individual choices by giving unwanted fermion representations a mass high enough to effectively decouple them from the model 7 . All mixings between the extra fermions and SM fermions are neglected through the assumption of a discrete Z 2 symmetry. Of course in a realistic model the mixings have to be taken into account, as they allow the necessary decays of the coloured vector-like fermions into SM particles. Table 6. Extra particle content of the toy models. S is either the CP-even, CP-odd or complex scalar. The various fermions Ψ F i ≡ Ψ F iL each come with a right-handed partner Ψ F iR with opposite quantum numbers. These models are based on the collection given in Ref. [166], while the last column contains other works where fermions in these specific representations are used. All SM particles have charge '+' under the additional Z 2 symmetry.
We write the scalar potentials for the three different types of scalars as The Yukawa interactions are given by In the Lagrangian above one should substitute the expression for the relevant scalar field Note that imposing CP conservation forces κ HS and κ S to vanish in the CP-odd potential. For both the CP-even and complex singlet models the CP-even component φ S mixes with the neutral Higgs field φ h at tree-level if κ HS = 0. As discussed in Section 2.2.1, even if one sets κ HS = 0 mixing between the CP-even states is induced at the loop level.

Models based on the SM gauge-group
We now turn our attention to complete models that have been proposed to explain the 750 GeV diphoton excess. To begin with we consider models that are based on the SM gauge group, with or without additional global symmetries. We divide the possible models into two main categories: (i) models with a SM-like Higgs sector and (ii) Two-Higgs-doublet type models.

Portal dark matter model
• Reference: [67,133] • Model name: SM+VL/PortalDM This model proposes that the resonance is produced by a 750 GeV real scalar singlet S, with the diphoton rate boosted through the introduction of vector-like quarks coupling to the new scalar singlet. In this model we have three possible options for the representation of the new vector-like matter. These choices are: (i) the addition of a vector-like up-type quark pair t L/R [67], (ii) in addition to the vector-like up-type quark pair, a vector-like quark doublet pair Q L/R is introduced [67] and finally, (iii) the addition of only the vector-like pair X L/R , also triplet under SU (3) C but with exotic hypercharge [133].
The model also introduces a new real scalar singlet S DM and an additional fermionic singlet Ψ DM as DM candidates, with a Z 2 symmetry to stabilise them. The particle content beyond the SM fields is given in Table 7. In order to avoid mixing with the SM quarks, the fields t L/R and Q L/R are also odd under the Z 2 .
The user can choose between the three model types by setting the couplings of the unwanted fields to zero and choosing their masses to be very large (for example, 10 12 GeV) to decouple them. Originally, the fermionic dark matter is absent from the models described in [67]. The exact settings are given below. The scalar potential for these models reads whereas the three model variants lead to three distinct forms for the Yukawa interactions, given by: where H = iσ 2 H * . In the model variants I and II, the vectorlike quarks decay into SM quarks and the scalar dark matter candidate S DM via the couplings Y S DM Q and Y S DM t . In model III, in turn, X is stable at the level of the Lagrangian and could only decay through higher-dimensional operators which are not included here.
The symmetry breaking pattern of the models is that of the SM, where the neutral component of the Higgs field acquires a VEV, plus the VEV of the scalar singlet S In general, φ S mixes with the SM Higgs. As mentioned previously, the user can choose between the three different models through the following parameter choices:

Scalar octet extension
• Reference: [48,99] • Model name: SM-S-Octet A charged scalar colour octet O coupled to a scalar singlet S was proposed in Refs. [48,99]. Here the singlet is the 750 GeV candidate, while the octet enters the loops that contribute to the generation of the couplings of the singlet to the gauge bosons. While Ref. [99] considers a toy model involving only the term S |O| 2 , Ref. [48] takes the singlet extended Manohar-Wise model [333]. For the SARAH implementation we have used the full model. However, since the cubic and quartic terms in O do not play a significant role, they are turned off by default in the SARAH model file. Table 8. Extra scalar field content of the octet extended SM.
The extra particle content with respect to the SM is a real singlet S and a scalar color octet O which is also charged under SU (2) L × U (1) Y , see Table 8. The isospin components of O are where A = 1, . . . , 8 is the adjoint colour index. The full scalar potential reads Electroweak symmetry-breaking (EWSB) is driven by the VEV of the neutral component of the SM Higgs doublet, which can be decomposed as Here φ H ≡ h is the Higgs boson, to be identified with the 125 GeV state discovered at the LHC. Similarly, the singlet S receives a VEV, and the neutral component of the octet is split into its CP-even and CP-odd eigenstates: We will now briefly discuss the parameter space of the model in order to justify our choice of input parameters. First, we consider the tadpole equations, which can be automatically derived by SARAH. Their solution for µ 2 and κ 1 is The tree-level mass matrix for the CP-even neutral scalars in the (φ H , φ S ) basis is given by We note that, in general, there is singlet-doublet mixing. There are two reasons to consider a small singlet-doublet mixing angle, θ. First, the stringent constraints derived from Higgs physics measurements, and second, the required suppressed decay widths into Higgses, W's and Z's in order to fit the diphoton signal -indeed in [48] values of ∼ 10 −2 were found to be required. If we have a small mixing angle, then we can write This implies λ S > 0, but also However, we also have v 2 S ∼ m 2 750 /8λ S , and so We thus require a tachyonic m 2 S for the SM Higgs mass condition: where in the last step we have taken λ S = 1.2, a rather large value. If we want κ 1 ∼ −1 then we require m 2 S ∼ −(600 GeV) 2 . On the other hand, from the second tadpole equation we have so putting these together we find the narrow window Alternative implementation in SARAH The above discussion suggests to use a different choice for the input parameters of the model in our SARAH implementation: ideally we would like the particle masses, the mixing and only dimensionless couplings to be the inputs. We shall take the input parameters to be In terms of these the other parameters are determined to be The exact version of these equations is implemented in SARAH and can be selected using the InputFile→"SPheno_diphoton.m" option in MakeAll or MakeSPheno.

Octet masses
One further input is taken in [48]: the physical mass of the octet scalars. These are given in terms of the Lagrangian parameters as: The values of λ i are taken to be small and equal in order for the octets to have similar masses, but since this is not the general case, we do not impose this choice in SARAH.
The choice in that paper does however hide the possibility of tachyonic m 2 O (and hence possible charge/colour breaking minima) -indeed, if we insist that m 2 O > 0 we have a lower bound on the masses of Clearly this is violated for m O 0,+ = 600 GeV when κ 2 ∼ 1, λ S 1. On the other hand, this does not guarantee a problem.
The desired vacuum has energy If we instead concentrate on the potential terms containing the octets, where only one component develops a VEV, we find Arranging for the additional minimum of the potential to be higher than the colourbreaking one then places a lower bound on the octet self-couplings, but for the phenomenology of the diphoton excess -when we neglect loop corrections to the mass of the octet -they play no other role.

Comments on fitting the excess
In [48] the authors find that the diphoton excess can be easily fit with octets at 600 or 1000 GeV and κ 2 ∼ 1.5 or 4.5, respectively. The scenario involves merely the simplifying assumption λ 1 = λ 2 = λ 3 so that the octets are of approximately equal mass. The ratio between the digluon and diphoton decay rates is then In [48] this is quoted as 715. In SARAH, before any NLO corrections are applied, the running of the Standard Model gauge couplings yields α s (750 GeV) = 0.091 and we use α(0) 137 −1 , giving a ratio of 700, in good agreement. However, when we include corrections up to N 3 LO, this ratio rises to 1150, putting the model near the boundary of exclusion due to dijet production at 8 TeV. These differences are illustrated in plots produced from SARAH/SPheno in figure 11. To produce these plots, all branching ratios/widths are calculated in SPheno, as is the production cross-section of the resonance at 8 TeV. To calculate 13 TeV cross-sections the 8 TeV cross-section was rescaled by the parton luminosity factor for gluons of 4.693.

Vector-like SU (2) triplet quark model
• Reference: [32] • Model name: SM+VL/TripletQuarks The model introduced in [32] considers a singlet scalar S as the candidate for the 750 GeV resonance. In order to produce the singlet via gluon fusion at the LHC, the model is further extended with the introduction of vector-like quarks, triplet under SU (2) L . Moreover, the singlet is charged under a Z 2 parity, although this is softly broken by the vector-like quark mass terms. Table 9. Extra scalar/fermionic degrees of freedom shown in the top/bottom. All SM particles are even under the imposed discrete symmetry.
This model is based on the SM gauge symmetry, extended with a Z 2 parity. The fermionic and scalar particle content is summarized in Table 9. The vector-like SU (2) L triplet quarks can be expressed in 2 × 2 matrix notation as Here we see that the U 1,2 and D 1,2 components have the same electric charges as the SM up-and down-type quarks, respectively. The components X and Y are exotic coloured states with electric charges 5/3 and −4/3, respectively, and thus they do not mix. The Yukawa Lagrangian of the model can be written as whereas the scalar potential is given by We note that the vector-like mass terms for F 1 and F 2 softly break the Z 2 parity. We assume that the SM Higgs doublet obtains a VEV while the introduced singlet does not, hence The Z 2 discrete symmetry, together with S = 0 ensured by the symmetry, implies that no mixing between the SM Higgs boson φ H and the singlet S appears at tree-level.

Single scalar leptoquark model
• Reference: [27] • Model name: Leptoquarks/ScalarLeptoquarks The model introduced in Ref. [27] considers a singlet scalar, S, as candidate for the 750 GeV resonance. It is based on the scalar leptoquark model in Ref. [334], with the addition of a vector-like fermion multiplet χ transforming in an a priori unspecified representation of SU (2) L . The simplest model, which the authors find to work, requires χ to be a SU (2) L triplet. The model is based on the gauge symmetries of the SM, augmented with either a discrete or gauge symmetry necessary for the DM sector. For the sake of simplicity we choose to realise the model using only a discrete Z 2 symmetry. The particle content beyond the SM is shown in Table 10. The new degrees of freedom are: the scalar leptoquarks φ, the gauge singlet scalar S, which is assumed to be real, and finally the vector-like triplet fermions χ. The vector-like SU (2) L triplet can be expressed in 2 × 2 matrix notation as Due to the discrete Z 2 symmetry the neutral component χ 0 is stable and a suitable DM candidate. In Ref. [27] the hypercharge Y χ and the number of generations N χ are left as free parameters. In order to boost the ratio of decay widths Γ S→γγ /Γ S→gg they claim that the minimal choice is N χ = 2 and Y χ = 0. The model files are implemented with these values.
The Yukawa interactions of the model and the vector-like mass term can be expressed as whereas the scalar potential is given by The Lagrangian in Eq. (5.32) is not generic. One could add more renormalizable operators involving the leptoquark, like d c R u R φ for instance. An operator of this kind, together with those in Eq. (5.32), would lead to rapid proton decay, thus it must be forbidden [334]. This can be achieved by imposing a discrete symmetry (for instance another Z 2 under which the leptoquark and the SM quarks are odd and everything else is even). In the model file we simply omit these dangerous operators. The potential in Eq. (5.33), on the other hand, is generic. It could be simplified by introducing yet another Z 2 to avoid terms with odd powers of S, and we would still have the necessary ingredients to fit the diphoton excess. We assume the EWSB proceeds by the usual Higgs VEV in addition to the singlet VEV v S , where the scalar singlet is expressed as (5.34)

Two scalar leptoquark model
• Reference: [65] • Model name: Leptoquarks/TwoScalarLeptoquarks This model includes four additional scalars, two of which, Φ and Ω are leptoquarks. The other two are a singlet S, which is the 750 GeV resonance, and a scalar Θ which carries hypercharge. The additional quantum numbers for these new field are shown in Table 11.
The Yukawa interactions consistent with these new fields is given by Table 11. Additional scalar particles beyond the SM appearing in the two scalar leptoquark model.
where the generation indices have been suppressed. These indices are important so that one obtains non-zero contractions for the first term in the above Yukawa. The scalar potential of the model is given by where the index i = {S, Θ, Φ, Ω} runs over the introduced scalar fields. The potential has an additional U (1) symmetry which can be gauged or global. In the implementation in SARAH we take it as global. This symmetry can be used to forbid terms such as u R d R Ω † , which could potentially lead to proton decay. One can generate neutrino masses at the two-loop level using a combination of the new Yukawa couplings and the VEV of the singlet. Beyond the neutral components of the Higgs doublet, the singlet also obtains a VEV of the form (5.37)

Georgi-Machacek model
• Reference: [70,108] • Model name: Georgi-Machacek The model originally proposed in Ref. [279] extends the SM by one real scalar SU (2) Ltriplet η with Y = 0 and one complex scalar SU (2) L -triplet with χ with Y = 1, which can be written as The CP-even components of the new scalars mix with the neutral SM Higgs. Usually, the lightest state is the 125 GeV Higgs boson, while the second mass eigenstate is identified with the 750 GeV resonance. The most compact form to write the Lagrangian in a SU (2) L × SU (2) R invariant form is to express the Higgs doublet and the two scalar triplets as In this form, the scalar potential reads τ a and t a are the SU (2) generators for the doublet and triplet representations respectively, while U is given for instance in Ref. [335]. Because of the custodial symmetry, the VEVs for the triplets are identical, and there are no tree-level contributions to electroweak precision observables. The compact form of the potential can not be implemented in SARAH, so we have translated into the form of Eq. (5.38). For example: This introduces more couplings, which have to be identical to preserve the custodial symmetry. A substitution rule to apply these relations is included in the model file.
The triplets receive VEV as which fulfil v η = v χ if the custodial symmetry is preserved.

Two-Higgs doublet models
Two-Higgs-doublet models (THDMs) are among the most common candidates to explain the diphoton excess. In this family of models, the 750 GeV resonance is typically a heavy Higgs, whose diphoton rate is enhanced with the addition of new states that contribute to its production cross-section and/or its decay width into pairs of photons. Several realisations of this idea can be found in the recent literature. Here we review a few representative examples.

THDM generalities
We first comment on some general features and conventions used in THDMs. For a complete review we refer to Ref. [336].

Couplings to fermions
We will consider three types of THDM, depending on the way the two Higgs doublets H 1 and H 2 couple to the SM fermions, as shown in this table: We note that, by convention, the right-handed up-type quarks u R always couple to H 2 , and that the MSSM can be seen as a particular example of a type-II THDM. The type-I and type-II THDMs are usually said to be models with natural flavour conservation. This is because, in contrast to the general type-III THDM, flavour changing neutral currents are completely absent at tree-level, thus satisfying flavour constraints more easily. This is achieved by coupling each fermion species to only one Higgs doublet, which can be enforced using a discrete symmetry. For example, the type-I THDM couplings to fermions can be obtained by imposing that H 1 → −H 1 is a symmetry of the Lagrangian. In our SARAH implementations we will not introduce these discrete symmetries explicitly, but simply allow only for those couplings that characterise the type of THDM under consideration. We point out that SARAH includes also other THDM versions, such as the lepton-specific or flipped ones. These versions can be easily combined with the extensions presented here for type I-III to include additional matter in order to explain the diphoton excess. Finally, we have also not implemented explicitly all representations for the vector-like states proposed in the literature so far. For instance, Ref. [160] introduces quarks with a very large hypercharge. However, it is easy to change the considered quantum numbers by changing the model files accordingly.

Scalar potential
We assume the following scalar potential for all THDMs considered below, In principle, the quartic terms are also allowed by the gauge symmetry. However, in our SARAH implementation we neglect these terms. This corresponds to assuming a global symmetry, which would be softly broken by the m 2 12 H † 1 H 2 term.

Symmetry breaking pattern
The symmetry breaking pattern is the same in all THDM variations, namely

Minimal vector-like THDM
• Reference: [41] • Model name: THDM+VL/min-3 THDM+VL/min-8 Ref. [41] considers a type-III THDM extended with two new vector-like colored fermions, S and Q V , aiming at a simultaneous explanation of the diphoton excess and of the CMS hint for the Higgs lepton flavor violating (LFV) decay h → τ µ [337]. The type-III THDM has already been shown in several papers to be able to accommodate this LFV signal, see e.g. [338,339], but [41] is the first work that addresses the diphoton excess at the same time. The representations of the new vector-like fermions under SM gauge group are shown in Table 12. We choose the explicit realizations R c = 3, Q = 2 (THDM+VL/minTHDM-3) and R c = 8, Q = 2 (THDM+VL/minTHDM-8). In both cases, the additional interaction terms beyond those in the standard THDM are where M V Q and M S are masses for the SU (2) L -doublet and singlet vector-like fermions, respectively, and λ D,S i (with i = 1, 2) are new Yukawa interactions.

THDM with exotic vector-like fermions
• References: [156,215] • Model name: THDM+VL/Type-I-VL THDM+VL/Type-II-VL The THDMs in Refs. [156,215] consider a less minimal framework where three generations of new vector-like fermions are added to the standard THDM scenario. Both type-I and type-II THDMs are studied. Furthermore, the vector-like leptons have exotic hypercharge values shown in Table 13. Table 13. Additional fermion field content for the THDM with exotic vector-like fermions. All SM fields are even under the imposed Z 2 discrete symmetry.
The vector-like lepton states can be decomposed (or denoted) as where we explicitly highlight the presence of doubly-charged leptons. The additional interaction terms beyond those in the standard THDMs are Here M V i are the bare vector-like masses and Y V i are Yukawa couplings. For simplicity, we assume that the vector-like states do not mix with the SM fermions. This can be enforced with a discrete Z 2 symmetry under which the vector-like states are odd and the rest of the states even.

THDM with SM-like vector-like fermions
• Reference: [21] • Model name: THDM+VL/Type-I-SM-like-VL THDM+VL/Type-II-SM-like-VL The THDM proposed in [21] introduces vector-like copies of the SM quarks and leptons, the SU (2) L doublets Q V and L V , as well as the SU (2) L singlets d V , u V , ν V and e V with their respective quantum numbers shown in Table 14. In contrast to other THDMs proposed to explain the diphoton excess, where low values of tan β are taken in order to increase the heavy Higgs coupling to the top quark and induce a large gluon fusion cross-section, this paper considers moderate and large values of tan β. In this case, the heavy Higgs production is mainly induced by the new vector-like colored states. Moreover, the advantage of using largish tan β values is that the total decay width of the heavy Higgs is suppressed, thus allowing for a smaller Γ(H → γγ) to explain the excess.
The new interaction terms take the same form as those in Eq. (5.46). Furthermore, as in the previous models, we assume that the new vector-like states do not mix with the SM fermions. For this reason, we introduce a Z 2 parity under which the vector-like states are odd and the rest of the states even. Finally, Ref. [21] considers two variants of this scenario in what concerns the couplings to the SM fermions: a type-I THDM and a type-II THDM.

THDM with a real scalar septuplet
• Reference: [135,188] • Model name: THDM/ScalarSeptuplet Ref. [135] explores an extension of the type-III THDM with new scalar SU (2) L representations, namely an inert complex Higgs triplet and a real scalar septuplet. A more general analysis can be found in Ref. [188], where a generic real scalar SU (2) L multiplet is considered. In both papers, the 750 GeV resonance is identified with the usual heavy Higgs of the THDM, with the additional scalars being introduced to push the diphoton rate to higher values.
For the realization of this idea here we focus on the septuplet case, closely following Ref. [135]. The introduction of such large scalar multiplets has become fairly popular recently in the context of minimal DM scenarios [340], although [135] does not explore any dark matter implications. It only takes advantage of the multicharged states contained in the septuplet which, due to their couplings to the THDM doublets, lead to a large diphoton rate for the heavy Higgs H.
The type-III THDM is extended with the addition of a real scalar which transforms as a 7 under SU (2) L , see Table 15. It proves convenient to use tensor notation for the septuplet. This is the usual choice, see Refs. [341,342], and the one implemented in SARAH. The septuplet T is represented by a symmetric tensor with six indices, T ijklmn , all of which can take values of either 1 or 2. The relation with the vector notation, employed for example in [135], is given by (5.47) The prefactor i and the sign for each component are introduced in order to satisfy T c = T where T c denotes charge conjugation of the field T , and T 0 is a real scalar.
The new potential terms involving the septuplet T are We note that two independent gauge invariant contractions are possible in case of T 4 , which is reflected by the second term of Eq. (5.48). Finally, the authors of [135] assume a minimum of the scalar potential with T = 0. In this case, the components of the septuplet do not mix with the scalar doublets, but only participate in the H → γγ rate due to the interaction terms λ T 3,4 .

U (1) extensions of the SM
In this section we consider a class of models which extend the SM by a new U (1) X gauge group. One typically introduces, beyond the SM Higgs doublet, new scalars charged under U (1) X , which serve two purposes: (i) a linear combination of them results in the 750 GeV particle, and (ii) via spontaneous symmetry breaking they give mass to the new gauge boson, the Z boson. Typically, one also introduces new fermions charged under the U (1) X symmetry that can either be singlets under the SM gauge group, hence forming a dark or hidden sector, or vector-like under the SM. An advantage of these models is that, through a suitable choice of charge assignments under U (1) X , one can avoid flavour constraints present when allowing the additional quarks to decay. Finally, the presence of a massive Z boson can lead to new collider signals, well studied in the literature, which can serve as a smoking gun for these types of models. We note that the mixing matrix for the neutral gauge bosons can be parametrised by two angles, Θ and Θ , with Θ highly constrained by LEP data. We distinguish two cases in the following: models in which the SM fermions are charged under the new Abelian group, and models in which they are singlet. This model is based on a gauged U (1) X extension of the SM with a dark sector that includes three generations of dark SU (2) L -singlet fermions and a DM scalar candidate.
The properties of the new particles introduced in this model are described in Table 16. We have added three right-handed neutrinos ν i R , neutral under U (1) X , to the original model, since N R (their dark partners) were considered here. Φ is the scalar field responsible for the spontaneous symmetry breaking of U (1) X , while X is the DM candidate. The U (1) X charges a, b are left arbitrary, with their assignment chosen such that they fulfil the anomaly cancellation conditions.  The Yukawa interactions and the scalar potential including new fields in the dark sector are described by where H denotes the SM Higgs doublet. For a = b = 1, an extra term Φ † X 2 would be allowed in the potential, which breaks U (1) X down to a Z 2 subgroup after Φ develops a non-zero VEV. Likewise, for 3a = (a + b), there appears an extra term Φ † X 3 , which breaks U (1) X down to a Z 3 subgroup after Φ develops a nonzero VEV. These possibilities are not considered here. The gauge symmetry is broken after H and Φ get non-zero VEVs, v and v Φ respectively, while X does not receive a VEV 8 . The scalar fields after EWSB can be expressed where φ i and σ i are the CP-even and odd components respectively. The 750 GeV candidate is the mixture of the SM Higgs φ H and φ Φ , which leads to constraints on λ HΦ . It is produced in gluon-gluon fusion via loops of the dark quarks, and it decays into diphotons via loops comprised only of charged dark fermions, assuming that the mixing with the SM-like Higgs is negligibly small.

Hidden U (1)
• Reference: [83] • Model name: U1Extensions/hiddenU1 This is a particularly simple realisation of a gauged U (1) X extension of the SM. As previously, a hidden U (1) X is added to the SM gauge group, under which all SM particles are singlets. A scalar S 1 is added for its spontaneous symmetry breaking, and a further scalar S 2 is added which is a singlet under the entire gauge symmetry of the model. Here we assume that both S 1 and S 2 can develop a VEV, in principle.
The 750 GeV candidate is considered to be predominantly composed by S 2 . To explain the diphoton signal, a vector-like quark is also included, carrying the same charge under U (1) X as the S 1 field. Hence, it couples only to S 2 due to the choice of charge assignments.  In Table 17 the hypercharge of the new vector-like quark Y X is left arbitrary. The implemented case is the most favourable one in terms of the diphoton excess, Y X = 2/3, which allows a mixing with the up-type quarks. We did not considered the case of adding also a vector-like lepton to the spectrum, which would lead to even larger rates. The U (1) X charge of the vector-like quark fields a does not affect the diphoton rate, but it has to be the same as for S 1 , which is relevant for the Z boson mass.
The scalar potential, with the usual doublet Higgs field H is given by This potential leads to a mixing between the three physical neutral scalars. The structure of the new scalar fields is given by where once again φ i and σ i are the CP-even and odd components respectively. As discussed in Section 2.2.2 we have allowed all scalar fields to obtain VEVs. The Yukawa interactions and fermionic mass terms of the hidden sector read The mixing of the vector-like quark X with SM quarks via the interaction with S 1 is kept small, and purely serves the purpose of letting the new quark decay.

Simple U (1)
• Reference: [62] • Model name: U1Extensions/simpleU1 This U (1) X extension of the SM augments its particle content by a pair of exotic vector-like quarks, χ 1 and χ 2 , doublets under SU (2) L and with hypercharge 7/6, and by one complex scalar Σ responsible for the spontaneous breaking of U (1) X . The U (1) Xbreaking Higgs boson is the 750 GeV candidate. The particle content beyond the SM is summarized in Table 18.
The scalar potential of the model is given by while the Yukawa and fermionic mass terms read Note that the original model proposed in [62] contains an effective operator As SARAH cannot handle effective operators this term is dropped from the model and subsequent constraints pertaining to stable charged particles are ignored. Expanding the new scalar field yields where once again φ i and σ i are the CP-even and odd components respectively.

Scotogenic U (1) Model
• Reference: [211] • Model name: U1Extensions/scotoU1 The matter particle content of the Scotogenic U (1) Model is summarized in Table 19, where, in addition to the fields charged under U (1) D , we introduce three copies of right-handed neutrinos ν R which are singlets under the full gauge group. Note, that the U (1) D charge of the H field has been changed to −1 compared to Ref. [211] in order to make the Yukawa interaction terms gauge invariant. The SM fields are not charged under the new U (1) gauge group. The discrete Z 2 symmetry is introduced to stabilize the dark matter candidates N . The scalar potential reads where H = (H + , H 0 ) is the SM Higgs field and q and the SM left-handed quark and lepton doublets. The U (1) D symmetry is eventually broken by the VEV of the scalar field Φ which one can decompose as whereas the CP-even component S is considered by the authors as the candidate for the 750 GeV resonance. The spontaneous symmetry breaking leaves the Z 2 parity intact, and H will therefore not develop any VEV. • Model name: U1Extensions/BL-VL This model is based on Refs. [343,344]. It considers a gauged U (1) B−L extension of the SM with an unconventional B −L charge assignment for the right-handed neutrinos ν R , and further requires 3 extra Dirac neutrinos N . It was originally proposed to explain the smallness of neutrino masses if neutrinos were Dirac particles. Finally, it also contains a scalar DM candidate, stabilised by the residual Z 3 symmetry from the breaking of the B − L symmetry. To fit the diphoton excess, new vector-like quarks are added so that the 750 GeV scalar, a linear combination of the SM Higgs and the χ fields, can be produced in gluon-gluon fusion. The particle content and the quantum numbers for this model are summarised in Table 20. In addition to the SM particle content, the model features three right-handed neutrinos ν i R , three pairs of SU (2) L singlet heavy fermions N i L,R , and two pairs of exotic quarks X L,R , Y L,R which carry color and electromagnetic charges but are singlets under SU (2) L .
The scalar potential of the model reads The residual global Z 3 symmetry protects the singlet scalar χ 2 from acquiring a VEV, i.e. χ 2 = 0. All leptons carry a charge ω = e 2iπ/3 under Z 3 . The CP-even degree of freedom of χ 2 is the DM candidate. The Yukawa interactions of the new sector read In the set up considered in Ref [186], the 750 GeV scalar is given by the combination (χ 6 −χ 3 ), that couples to gluons and photons via loops of X and Y fermions proportional to the Yukawa couplings f X and f Y . The fit to the diphoton signal was studied in a simplified scenario, where special relations are imposed to the scalar parameters. As already described in section 2.2.1, these relations are not protected by symmetries, and hence lead to a large amounts of fine tuning.

Sample of U (1) models based on different charge assignments
• Reference: [66] • Model name: U1Extensions/VLsample The particle content of this model is shown in Table 21. Similarly to all the previous models, a complex scalar S is added to break the U (1) X symmetry. In this paper, the mixing between the scalars is kept small, hence the 750 GeV candidate is predominantly the CP-even part of the U (1) X -breaking field. New doublet and singlet vector-like quarks, charged under U (1) X , are added to fit the diphoton signal strength. In this model, the SM fermions also carry X-charges, which are fixed according to the anomaly cancellation conditions. Only the Higgs doublet is not charged, while the charge of the new scalar field S is double (in absolute magnitude) that of the vector-like quarks.
There are several possibilities to assign the U (1) X charges in an anomaly-free way, with different physical interpretations [66]: Note that in all interpretations a remains a free parameter.
The Yukawa interactions of the extra particle content are given by Finally, the scalar potential is given by the symmetry breaking pattern being while the expansion of the scalar fields in to CP-even and odd components is analogous to Eq. (5.51).

Model with flavour-nonuniversal quark U (1) charges
• Reference: [179] • Model name: U1Extensions/nonUniversalU1 In this model the first generation of left-handed SM quarks carries a U (1) X charge while the second and third generations do not. In this way it is possible to add exotic quarks which are vector-like under the SM gauge group and achieve anomaly cancellation with less than three generations. The scalar sector then needs to be extended with a second Higgs doublet, with a different U (1) X charge compared to the first ones, in order to have Yukawa interactions for all quark families. A further complex scalar field S, which is a singlet under the SM gauge group but charged under U (1) X , is also added for the spontaneous symmetry breaking of the U (1) X symmetry.
The charge assignments that cancel all anomalies are given in [345] and are summarised in Table 22. 9 The Yukawa Lagrangian of the model (ignoring here flavour indices) is given by 66) and the scalar potential is The pattern of EWSB follows where once again φ i and σ i are the CP-even and CP-odd components respectively. In order to obtain a massive pseudo-scalar state, κ HS needs to be non-zero. Hence, for keeping the S − H i mixing small while making the pseudo-scalar massive, either the condition v S v must hold, or κ HS must be small, in conjunction with v 1 v 2 → 0 or ∞. φ S , the CP-even component of S, is then identified with the 750 GeV candidate.

Leptophobic U (1) model
• Reference: [168] • Model name: U1Extensions/U1Leptophobic This model is inspired by an E 6 Grand Unified Theory (GUT), but the authors only consider the low energy version where the SM gauge group is augmented by an extra gauged U (1) X symmetry. This extra U (1) symmetry has zero charges for both left-and right-handed leptons making it entirely lepotophobic. However it is impossible to arrange for this to happen by taking linear combinations of the U (1) χ and U (1) ψ groups that appear from the breakdown of E 6 . Instead, the charges from this extra U (1) can only be obtained by including gauge kinetic mixing, so that the introduced mixture of U (1) Y charges exactly cancel the non-zero leptonic charges. This can be done with the U (1) η gauge symmetry and the charges used in this model correspond exactly to those given in Table I of Ref. [346]. It is these charges rather than the charges of the U (1) η which are set in this model. Of course one may discard the E 6 motivation and treat it as an ad-hoc choice of U (1) charges.
The model contains two Higgs doublets, a complex scalar SM singlet (Φ), charged under U (1) X , plus right handed neutrinos and other new fermions, charged under both the SM and the U (1) X symmetries. The latter are odd under an imposed Z 2 symmetry, so that the lightest neutral fermion is a DM candidate. The particle content is summarized in Table 23, while the scalar potential reads and the Yukawa interactions are given by We assume the following symmetry breaking pattern where once again φ i and σ i are the CP-even and CP-odd components, respectively, and we define tan β ≡ v 2 v 1 . The 750 GeV candidate is taken to be dominantly composed by φ Φ , that is, the real CP-even degree of freedom of the Φ field after U (1) X symmetry breaking. It couples to photons and gluons via loops of the new fermions. The model cannot explain the diphoton excess with Yukawa couplings in the perturbative range, but the authors use values between 5 and 10. As stressed in Sec. 2.1.2.1, this renders the perturbative calculation, and hence the results, very questionable.

U (1) extension with a Z mimicking a scalar resonance
• Reference: [60] • Model name: U1Extensions/trickingLY The idea of this model is that the extra neutral gauge boson decays into Sγ, whereas the scalar S itself decays into a diphoton final state. Because of the high boost, the two photons from S appear to be a single photon in the detector.
Ref. [60] works in a model realization where the third-generation quarks are charged under U (1) whereas the first and second generations are not. While that can be viewed as a toy model to make a point, an actual realisation would either need additional Higgs representations or flavour-universal U (1) charges in order to reproduce the correct CKM matrix. Consequently, we will work with U (1) charges for all three generations of SM quarks and also have to use three generations of each additional exotic particle for anomaly cancellation.
The particle content of the model is summarized in Table 24. For anomaly cancellation, the condition Q 1 − Q 2 = −3 must hold. Furthermore, there are four different valid choices for the hypercharge assignments Y i [60]: In the model implementation at hand, we choose (Y 1 , Y 2 , Y 3 ) = ( 1 2 , 1, 0). The Yukawa interactions including fields beyond the SM read Table 24. Fermionic and scalar particle content of the trickingLY. Here The scalar potential is given by where the U (1) symmetry is broken spontaneously as soon as S receives a VEV according to S = v S / √ 2. Unfortunately, the most interesting vertex for this model, Z − S − γ, only arises at the loop level. This would require a handling via an effective operator which is currently not supported in the automatized tools advertised. Hence, it is not possible to recast the results of Ref. [60] based on this model implementation only.
Note that, in principle, the decay Z → ZS is already possible at tree level due to Z − Z mixing and dominates over the decay into Sγ. Therefore, in order to achieve the desired effect, the mixed gauge interaction term F µν F µν must be forbidden while allowing for SF µν F µν which is hard to justify in general.

Left-right symmetric models
Left-right (LR) symmetric models can potentially provide an interesting explanation of the diphoton excess through the use of an extended scalar sector that is necessary to spontaneously break the enlarged gauge group G L−R ≡ SU (3) C × SU (2) L × SU (2) R × U (1) B−L to the SM gauge group. However, due to the large number of fields, these models are often difficult to analyse even at tree-level. Therefore as a starting point we provide model files for four different left-right models that have been proposed in the literature to explain the diphoton excess. These four models are based on the above mentioned gauge group G L−R with, in two cases, further Abelian gauged or global symmetries.

Left-right symmetric model without bi-doublets
• Reference: [42,92,93] (see also Ref. [347]) • Model name: LRmodels/LR-VL In Ref. [42] the authors explored the possibility of explaining the observed diphoton excess in the context of the minimal left-right symmetric model. They show that it is not possible and that an alternative model is necessary. Therefore they give up on standard Yukawa couplings, for which bi-doublets are necessary, and consider separate SU (2) L -and SU (2) R -Higgs-doublets. In order to be able to introduce Yukawa interactions, new vectorlike SU (2) i -singlet fermions are introduced. After integrating out the vectorlike fermions, the SM fermion masses are generated through a universal seesaw mechanism [348,349].
The particle content for the model is shown in Table 25. Note that the authors consider the second-lightest CP-even Higgs, which should be predominantly singlet-like, as the particle responsible for the observed resonance. The scalar potential given the particle content and consistent with the symmetries (SU (3 The Yukawa interactions can be written as 10 10 Note that our Yukawa interactions differ from the literature: in Ref. [42], they are defined as, e.g., q L H † L U L which contracts to zero because of the implicit left/right projection operators. Moreover, in Refs. [42,92] the 'conjugate' assignments of the H L/R need to be exchanged in order to obtain a gauge-invariant Lagrangian. where N c L is the charge conjugate of N L . Note that we have included both Majorana and Dirac mass terms for the fermionic singlet N L/R . We assume the following symmetry breaking VEVs  (2)-singlet fermions allows for the generation of the SM fermion masses via a see-saw mechanism. Additionally, a parity symmetry is imposed to ensure a vanishingθ parameter at tree-level in the QCD Lagrangian [350], in order to solve the strong CP-problem without introducing an axion.

Field
The particle content of this model and charge assignments under the gauge symmetries are shown in Table 26. The proposed candidate for the 750 GeV resonance is taken to be one of the CP-even scalars associated with the SU (2)-singlet Higgs scalars σ D , σ U and σ E that are responsible for the breaking U (1) L × U (1) R → U (1) B−L . The decays of this state into digluon and diphoton final states are assumed to proceed via loops involving the SU (2)-singlet fermions. The Yukawa interactions consistent with the imposed parity are given by The parity symmetry is taken to be softly broken, so that the part of the scalar potential considered in Ref. [51] is given by The couplings and masses λ, ξ, η, µ 2 φ L , µ 2 φ R , µ 2 ∆ L , µ 2 ∆ R , ρ L and ρ R are taken to be real, with µ 2 φ L = µ 2 φ R , µ 2 ∆ L = µ 2 ∆ R and ρ L = ρ R . Note that Eq. (5.78) and Eq. (5.79) differ from Eq. (6) and Eq. (7) in Ref. [51], which as given are not gauge invariant. One may also include a large number of additional terms that are allowed by the gauge symmetries, given by The full scalar potential that we consider is then V + V . The SU (2)-singlet Higgs scalars are assumed to acquire VEVs of the form are taken to develop VEVs of the form The non-zero VEV of φ 0 R leads to the breakdown of SU (2) L × SU (2) R × U (1) B−L to SU (2) L × U (1) Y , which is subsequently broken by the VEV of φ 0 L . As a result, the triplet scalars also acquire VEVs of the form

Left-right symmetric model with fermionic and scalar triplets
• Reference: [33] • Model name: LRmodels/tripletLR In this model the diphoton signal is produced through a cascade decay, namely, pp → Z → XY → XX(δ 0 → γγ) where X and Y are unspecified soft states and δ 0 is the neutral component of the SU (2) R scalar triplet. However, in order to sufficiently boost the rate three SU (2) R triplet fermion fields are added to the model, T 1 , T 2 and T 3 . The model is based on a SU (3) C × SU (2) L × SU (2) R × U (1) B−L gauge group which is broken to the SM gauge group through the VEV of the triplet field ∆ R whereby EWSB proceeds through the bi doublet Φ VEVs. The entire particle content of the model is illustrated in Table 27.
The scalar fields of the model can be expressed as leading to a scalar potential of the form whereΦ ≡ −σ 2 Φ * σ 2 . The Yukawa interactions of the model are given by where α runs over the quarks and leptons α = Q, L and Ψ L,R = (ψ u L,R , ψ d L,R ) with ψ u = u, ν and ψ d = d, . C is the charge conjugation operator. The VEVs of the scalar fields in the model take the form

Dark left-right symmetric model
• Reference: [95] • Model name: LRmodels/darkLR The main idea of this model is to add an additional symmetry in order to stabilize the DM candidate, namely right-handed neutrinos, so that they cannot decay via a W channel. This additional symmetry takes the form of a global Abelian symmetry labelled as U (1) S . The spontaneous breaking of SU (2) R × U (1) S is such that the combinationL = S − T 3R , where T 3R is the third component of the right-handed isospin, remains unbroken. HereL is interpreted as a generalised lepton number. The full particle content of the model is shown in Table 28. Note that the scalar sector of the model is enlarged to include both left-and right-handed triplets and doublets as well as the usual bi-doublet. These scalar fields can be expressed as Subsequently, the proposed candidate for the diphoton excess is φ 0 R , the neutral component of the SU (2) R doublet. Running in the loop will be W -bosons, as well as δ + R and δ ++ R Higgses. However, these particles are insufficient to boost the rate to diphotons, therefore additional quarks x L and d R are introduced.
The scalar potential of the model is whereΦ ≡ −σ 2 Φ * σ 2 . Note that there are a number of extra terms present in this potential compared to [95], which are allowed under the symmetries of the model. Finally the Yukawa interactions are given by where C is the charge conjugation operator. The structure of the VEVs of the model are

331 models
Models based on the SU(3) c × SU(3) L × U(1) X gauge symmetry [351][352][353][354][355][356][357], 331 for short, constitute an extension of the SM that could explain the number of generations of matter fields. This is possible as anomaly cancellation forces the number of generations to be equal to the number of quark colours. Regarding the diphoton excess, 331 models automatically include all the required ingredients to explain the hint. First, the usual SU (2) L Higgs doublet must be promoted to a SU There are several variants of SU(3) c × SU(3) L × U(1) X models. These are characterized by their β parameter 11 , which defines the electric charge operator as 12 First, in Sec. 5.4.2.1 we consider the model in Ref. [44]. This 331 variant has β = 1/ √ 3, which fixes the electric charges of all the states contained in the SU (2) L triplets and anti-triplets to the usual 0, ±1 values. In Sec. 5.4.2.2 we consider a 331 model with β = − √ 3, a value leading to exotic electric charges. This 331 variant has been discussed in the context of the diphoton excess in [53,147,359]. Although the mechanism to explain the diphoton excess is exactly the same as in [44], the presence of the exotic states leads to slightly different numerical results.

On the SU (3) generators in SARAH
The most common choice for the SU (3) generators is T a = λa 2 , with λ a (a = 1, . . . , 8) the Gell-Mann matrices. However, this is just one of the possible representations. In fact, SARAH uses a different set of matrices, T SARAH a = Λa 2 , following the conventions of Susyno [294]. The relation between the non-diagonal matrices in the two bases is Concerning the diagonal matrices, the usual λ 3,8 Gell-Mann matrices, are replaced by Λ 7,8 , The electric charge operator can be written, using the conventions in SARAH, as This in turn implies that the charge assignments in the SU (3) multiplets must be adapted as well. For example, one can easily check that the electric charges of the first and third components of a SU (3) triplet t are exchanged when going from the usual Gell-Mann representation to the basis choice employed in SARAH, In the following we will use the standard conventions based on the Gell-Mann matrices in order to keep the discussion as close to the original works as possible. However, we emphasize that the implementation of the 331 models in SARAH requires this dictionary between the bases. It should also be noted that in the current implementation in SARAH of the 331 models described below, vertices involving four vector bosons in the generated model files for CalcHep cannot yet be handled correctly. In order to generate model files that will work with CalcHep, one must therefore exclude these vertices from being written out by SARAH, as described in Appendix A.12.

331 model without exotic charges
• Reference: [44] • Model name: 331/v1 The model is based on the SU(3) c × SU(3) L × U(1) X gauge symmetry, extended with a global U (1) L and an auxiliary Z 2 symmetry to forbid some undesired couplings. The fermionic and scalar particle content of the model is summarized in Table 29. In addition, due to the extended group structure, the model contains 17 gauge bosons: the usual 8 gluons; 8 W i bosons associated to SU (3) L and the B boson associated to U (1) X . The fermionic SU (3) L triplets of the model can be decomposed as (5.100) The notation used for the extra quarks that constitute the third components of the SU (3) L triplets Q 1,2,3 L is motivated by the fact that their electric charges are −1/3 and 2/3 for D/S and T , respectively. The scalar multiplets can be written as 3 and S − 1 are electrically charged scalars, the components φ 1,2,3,X , S 2,3 and X are neutral.
The Yukawa Lagrangian of the model can be split as and We definedd R ≡ (D R , S R ). We note that Eq. (5.104) leads to an inverse seesaw mechanism for neutrino masses [360,361]. Here, y a is anti-symmetric while m s is symmetric, whereas the rest of Yukawa couplings are generic matrices, including those in Eq. (5.103). An additional term y s Xψ L sΦ X could be added to Eq. (5.104), but given that Φ X = 0, it does not contribute to neutrino masses and we will drop it for simplicity. Finally, the scalar potential is given by , is required to break unwanted accidental symmetries in the scalar potential. We will assume the following symmetry breaking pattern (5.106)

331 model with exotic charges
• Reference: [53] (see also [147,359] for similar constructions) • Model name: 331/v2 Now, we will consider a 331 variant with β = − √ 3, as discussed in the context of the diphoton excess in [53]. The fermionic and scalar particle content of the model is summarized in Table 30. In addition, the model contains 17 gauge bosons: the usual 8 gluons; 8 W i bosons associated to SU (3) L and the B boson associated to U (1) X .
The fermionic SU (3) L triplet representations of the model can be decomposed as Due to the choice β = − √ 3, the electric charges for the extra quarks that constitute the third components of the SU (3) L triplets Q 1,2,3 L are 5/3, 5/3 and −4/3, respectively. The scalar triplets can be written as (5.108) Therefore, the particle spectrum of the model contains the exotic quarks in Eq. (5.107), as well as the doubly-charged fermion E −− and the scalars ρ ++ and χ −− . The Yukawa Lagrangian of the model can be split as where we have definedd R ≡ (D R , S R ), and We note that the exotic fermions E, D, S and T only couple to the χ scalar triplet, and thus only via its vacuum expectation value (VEV) they will acquire masses. Finally, the scalar potential is given by V = µ 2 1 |ρ| 2 + λ 1 |ρ| 4 + µ 2 2 |η| 2 + λ 2 |η| 4 + µ 2 3 |χ| 2 + λ 3 |χ| 4 + λ 12 |ρ| 2 |η| 2 + λ 13 |η| 2 |χ| 2 + λ 23 |η| 2 |χ| 2 + +λ 12 We will assume the following symmetry breaking pattern In this case, the non-zero VEV of χ is responsible for the breaking SU The requirement that this occurs at a scale much above the EW scale then imposes a hierarchy amongst the VEVs, namely that v 3 v 1 , v 2 . Consequently, one of the CP-even scalar states is predominantly from the χ triplet and decouples from the SM. This scalar is then identified as the candidate for the 750 GeV resonance in this model. The decays of this state into two photons proceed via loops involving the heavy fermions, as well as those involving the charged scalars and additional charged vector bosons.

Other models
Gauged THD model • Reference: [151] • Model name: GTHDM The GTHDM model [362] comes with an additional gauged SU (2) H symmetry and a U (1) X symmetry, which is either global or gauged as well. Since the minimal, gauged version suffers from the fact that two massless vector bosons are present, U (1) X is treated as global symmetry. The scalar and fermion fields are listed in Table 31.
The Lagrangian of the GTHDM contains the SM Lagrangian, extended by the new terms After EWSB, there are three neutral gauge bosons which mix giving rise to the γ, Z, Z mass eigenstates and two charged ones (W , W ) which do not mix. The neutral component of the SM-singlet Φ, φ 0 , is considered to be the candidate for the 750 GeV resonance while its VEV gives mass to the exotic fermions that are needed to run in the loops. As φ 0 will typically mix with the SU (2) L doublet Higgs, this mixing needs to be suppressed by a specific parameter choice in order to avoid the tight bounds from the dijet, ZZ or dilepton channels.

Supersymmetric models
There are several ideas to explain the diphoton excess within a supersymmetric framework. Some of them make use of SUSY models which already exist in the literature, and for which also SARAH model files exist: the MSSM with trilinear R-parity violation [6,100], the simplest models with Dirac gauginos [55], or the model with gauged U (1) L × U (1) B [111]. We will not make any further comment on these models, but concentrate in the following on the models which are newly implemented.

NMSSM extensions with vector-like multiplets
• Reference: [105,129,208,209] • Model name: NMSSM+VL The scalar component of the gauge-singlet superfieldŜ of the Next to Minimal Supersymmetric Standard Model (NMSSM) can explain the 750 GeV resonance, if one adds vector-like SU (5) multiplets to enhance the diphoton rate. The new multiplets are added in pairs of (5,5) and/or (10,10) in order to preserve gauge coupling unification. Typically one also imposes a Z 2 symmetry to forbid mixing of the new vector-like particles with the MSSM particles. The authors of Ref. [209] mention the possibility to interpret the resonance as two nearly degenerate singlet-like bosons, roughly the scalar and pseudoscalar components of the singletŜ. There are some differences in the singlet and Higgs superpotential interactions included in the different papers: • in Ref. [208] the authors assume λŜĤ uĤd and κ/3Ŝ 3 to be present, but µĤ uĤd to be absent; • in Ref. [105] µĤ u H d , λŜĤ uĤd and M S /2Ŝ 2 are present; • in Ref. [129] only M S /2Ŝ 2 is present. This does not cause a mixing between the singlet and Higgs doublet at the tree level, but such a mixing is unavoidable at the loop level.
The SARAH implementations use the most general version of the superpotential: all possible interactions are present. The different limits according to the proposals of Refs. [105,129,208] can be obtained by setting the corresponding parameters to zero in numerical studies. In what follows we describe the models with the vector-like multiplets in different representations of SU (5).

NMSSM with vectorlike top
• Model name: NMSSM+VL/VLtop • Reference: [209] This model is an extension of the NMSSM by a vector-like top. There is a global Z 2 R-parity and a Z 3 symmetry, under which all particles transform as X → exp(i 2π 3 ) X. In this way, only terms with three superfields are allowed in the superpotential. The particle content is given in Table 32. Since this model does not introduce complete multiplets (5 or 10) of SU (5), gauge coupling unification is not achieved. The superpotential is given by Beyond the neutral scalar components of the two Higgs doublets, after EWSB the complex singlet gets a VEV and can be decomposed as The fermionic components ofT ,T mix with the up-type quarks, while the scalar components mix with the up-like squarks.

Pairs of 5 of SU (5)
• Model name: NMSSM+VL/5plets The superfields beyond the MSSM are shown in Table 33. In the current implementation we only have one copy of (5,5) fields, but having at least three copies of them should give a better fit to the diphoton resonance. According to [105] the fit is even better with four copies, however in that case one might hit a Landau pole.
The superpotential is given by Beyond the neutral scalar components of the two Higgs doublets, also the singlet gets a VEV after EWSB and can be decomposed as in Eq. (5.117).

Pairs of 10 of SU (5)
• Model name: NMSSM+VL/10plets The superfields beyond the MSSM are shown in Table 34. The superpotential is given by The symmetry breaking pattern is the same as for the model with 5-plets.

Pairs of 5 and 10 of SU (5)
• Model name: NMSSM+VL/5+10plets This model combines the previous two setups, adding pairs of vectorlike 5 and 10 representations of SU (5). The superfields beyond the MSSM are shown in Table 35. The superpotential is given by The symmetry breaking pattern is the same as for the model with 5-plets.

Pairs of 5 of SU (5) and R-parity violation
• Model name: NMSSM+VL/5plets+RpV One can relax the assumption of the Z 2 symmetry that forbids mixing between vector-like fields and MSSM fields, in which case terms like κ 5ŜLĤu are added to model [105]. Furthermore, this also breaks R-parity.
The superfields beyond the MSSM are shown in Table 33 and the superpotential is given by The inclusion of R-parity violating terms triggers VEVs for the neutral components ofL andL , 122) and causes mixing between the vector-like states and the Higgs components.

Broken MRSSM
• Reference: [58] • Model name: brokenMRSSM In the minimal R-supersymmetric model (MRSSM) the scalar R u (see Table 36) is proposed as an explanation to the 750 GeV resonance. In order to explain the diphoton excess, it is necessary to add explicitly an R-symmetry breaking term to the Lagrangian where T u is a dimensionful trilinear coupling. This source of R-symmetry breaking has several consequences, not discussed in Ref. [58], which however are taken into account in the model implementation:  The superfields beyond the MSSM are listed in Table 36. The superpotential, assumed to conserve R-symmetry, is given by As explained above, because the R-symmetry is broken in the soft sector of the model, all possible tri-and bilinear soft-breaking terms corresponding to the superpotential terms will be generated.
The following VEVs appear after EWSB, beyond those of the neutral scalar components of the two Higgs doublets (5.125) The authors favor to have large stop mixing for a not too large R-symmetry breaking term T u by considering the limit v d ∼ v u and mt L ∼ mt R . However, in this limit the mass of the SM-like Higgs is tiny and often tachyonic: in the MSSM, the Higgs tree-level mass vanishes for tan β → 1, and this model has additional negative contributions to the mass from the new D-terms present in models with Dirac gauginos. It is also questionable if the case with very large T u is a viable scenario because these values are highly restricted by charge and colour breaking minima [363,364], which demands careful checks. This is similar to the vacuum stability issues discussed in Section 2.

U (1) -extended MSSM
• Reference: [155,332] • Model name: MSSM+U1prime-VL  In this model all MSSM fields carry a non-zero U (1) -charge so that anomaly cancellation requires additional superfields (see Table 37), which are also responsible for the spontaneous U (1) breaking. Furthermore, colour-charged and colour-uncharged matter superfields which are vector-like with respect to the MSSM gauge group are introduced. A combination of scalar singlets S and S i is supposed to give the 750 GeV resonance.
The complete superfield content with all gauge quantum numbers is given in Tables 37 and 38. In addition to the usual matter parity, we impose a Z 2 symmetry under which all exotic matter superfields are odd and all other superfields are even in order to reduce the number of superpotential terms and hence reduce the complexity of the model.
The superpotential is given by In addition to the neutral components of the two Higgs doublets, the MSSM singlets get VEVs according to • Model name: SUSYmodels/E6SSMalt E 6 -inspired SUSY models predict extra SM-gauge singlets and extra exotic fermions, so they immediately have the ingredients that many authors have tried to use to fit the diphoton excess. These models are often motivated as a solution to the µ-problem of the MSSM, because the extra U (1) gauge symmetry forbids the µ-term, while when one of the singlet fields develops a VEV at the TeV scale this breaks the extra U (1) giving rise to a massive Z vector boson and at the same time generates an effective µ term through the singlet interaction with the up-and down-type Higgs fields, λŜĤ uĤd . The matter content of the model at low energies fills three generations of complete 27-plet representations of E 6 , which ensures that anomalies automatically cancel. A number of models of this nature have been proposed as explanations of the diphoton excess [64,165,365]. The example we implement here [64] is a variant of the E 6 SSM [366,367]. In this version two singlet states develop VEVs and the idea is that the 750 GeV excess is explained by one of these singlet states with a loop-induced decay through the exotic states.
In E 6 models the extra U (1) which extends the SM gauge group is given as a linear combination of U (1) ψ and U (1) χ which appear from the breakdown of the E 6 symmetry as E 6 → SO(10) × U (1) ψ followed by SO(10) into SU (5), SO(10) → SU (5) × U (1) χ . In the E 6 SSM and the variant implemented here the specific combination is, To allow one-step gauge coupling unification however some incomplete multiplets must be included in the low energy matter content. So in addition to the matter filling complete 27 representations of E 6 there are also two SU (2) multipletsĤ and H , which are the only components from additional 27 and 27 that survive to low energies. All gauge anomalies cancel between these two states, so they do not introduce any gauge anomalies. Furthermore, the low energy matter content of the model beyond the MSSM one includes three generations of exotic diquarks 13 ,D i ,D i , three generations of SM singlet superfieldsŜ i and extra Higgs-like states H u 1,2 and H d 1,2 that do not get VEVs.
The full set of superfields are given in Table 39 along with their representations under SU (3) and SU (2) and the charges of the two U(1) gauge groups and the discrete symmetries, which we will now discuss.
The Z L 2 symmetry plays a role similar to R-parity in the MSSM, being imposed to avoid rapid proton decay in the model. However with this imposed there are still terms in the superpotential that can lead to dangerous flavour changing neutral currents (FCNCs). To forbid these, an approximate Z H 2 symmetry is imposed. In the original E 6 SSM model onlyŜ 3 ,Ĥ d andĤ u were even under the Z H 2 symmetry, however in this variant S 2 is also even under this approximate symmetry.
The full superpotential before imposing any discrete symmetries is given by  However, with the discrete symmetries imposed and integrating out the heavy righthanded neutrinos, the superpotential in this specific variant reduces to 14 , One should remember that the Z H 2 can only be an approximate symmetry as otherwise the exotic quarks could not decay. In this variant the exotic quarks decay therough the Z H 2 violating interactions of W 1 .
In the paper it is assumed that the singlet mixing can be negligible and the numerical calculation was performed under this assumption, neglecting any mixing between the singlet state which decays to γγ via the exotic states and the other CP-even Higgs states from the standard SU (2) doublets. However it is clear that there must be some mixing from the D-terms, and therefore if that is included one important check would be to test whether other decays are sufficiently suppressed. Moreover, the parameters needed to simultaneously get a 125 GeV SM-like Higgs state and a 750 GeV singletdominated state are not given. In this respect we note that the singlet VEVs appear both in the diagonal entries of the mass matrix and in the off-diagonal entries that mix the singlet states with the doublet states.
We finally note that other similar E 6 models have also been proposed in the context of the diphoton excess. These include a model by two authors from the original paper [165], a model with a different U (1) group at low energies [368], and a model that is still E 6 -inspired, but where no extra U (1) survives down to low energies [161].

Study of a natural SUSY explanation for the diphoton excess
We show in this section how one can use the described setup to perform easily a detailed study of a new model that aims at explaining the diphoton anomaly.

The model
We are now going to study a SUSY model which enhances the tree-level Higgs mass due to non-decoupling D-terms. The model is based on that proposed in Ref. [369] as a natural SUSY model which allows for light stops compatible with the measured Higgs boson mass, extended by three generations of pairs of vector-like quarks and leptons. We want to achieve a tree-level enhancement of the SM-like Higgs mass and an explanation of the diphoton excess via the loop-induced decay of a CP-odd scalar. In addition, we will also check whether one can get a broad diphoton resonance in this model. The matter field content is shown in Table 40 and the considered superpotential reads: We will not make the simplifying assumption that mixings between the MSSM states and the new vector-like fields can be neglected. Of course, such mixing could have Table 40. Scalars and fermions in the U (1) X -extended MSSM been forbidden by choosing different U (1) X charges for the new particles. However, in such case there would be a conserved Z 2 symmetry associated to the vector-like states (under which all vector-like superfields are odd and the rest are even) that would make the lightest of them absolutely stable. This would be a problem unless that state is neutral and colourless, and thus this scenario can only be viable if we also consider additional singlet vector-like states, such as vector-like partners for the right-handed neutrinos, and make them lighter than the other vector-like states. Thus, this setup would predict two stable particles to make the dark matter. Such a scenario could also be studied with the tools presented here. However, we decided not to consider this option in the following. The other main ingredients of the model are the general soft-SUSY breaking terms, which read Note that we have included the gaugino mass term M 1X arising from gauge kinetic mixing. All the terms shown in Eq. (6.2) are automatically added by SARAH based on the information provided by the user about the particle content and the superpotential. Several scalar fields acquire VEVs. We decompose them as In addition, the sneutrinos are decomposed with respect to their CP eigenstates, which in general have different masses due to the Majorana mass-term Y X η in the superpotential. Since H 0 d and H 0 u carry charges under both U (1) gauge groups, there will be non-zero Z-Z mixing even in the limit of vanishing gauge kinetic mixing. The list of particle mixings, which go beyond the usual MSSM mixings reads 14) The model files which implement this model in SARAH are discussed in Appendix B.2.
In addition, we provide all files to reproduce the computations that follow at http://sarah.hepforge.org/U1xMSSM_example.tar.gz

Analytical results with Mathematica
Before we perform a numerically precise study of the model, we show how already with just SARAH and Mathematica one can gain a lot of information about a new model.

Consistency Checks
The model is initialised after loading it in SARAH via One can see that SARAH tests all different combinations of gauge anomalies and, given that no warning is printed on the screen, confirms that all of them cancel. Similarly, it also checks that all terms in the superpotential are in agreement with all global and local symmetries. More detailed checks can be carried out by running CheckModel[] when the initialisation is finished. After a few seconds, a message is printed telling that the model is loaded.
All Done . U1xMSSM is ready !

Gauge sector
Before we discuss the matter sector or the scalar potential, we have a brief look at the gauge bosons. We make use of the mass matrices calculated by SARAH during the initialisation of the model. We find a handy expression for the mass matrix of the neutral gauge bosons in the limit of vanishing gauge kinetic mixing (g X1 = g 1X = 0) via As expected, the first two eigenvalues are just the ones of the SM gauge bosons, while the mass of the new gauge boson is given by We will use this relation in the following to replace x by M Z in all equations.

Scalar Sector
Solving the tadpole equations We turn now to the scalar sector of the model. First, we make a list with a few simplifying assumptions which we are going to use in the following Here we assume all parameters to be real, remove any complex conjugation (conj) and use the Landau gauge (RXi[_]->0), then we turn off again gauge kinetic mixing and take the VEVs of η andη to be equal. In the fourth line, we parametrise v d and v u as usual in terms of v and tan β. Finally, we set the parameters κ, T κ , λ, T λ and L ξ to zero. We can now solve the tadpole equations, stored by SARAH in TadpoleEquations We have saved the solution in the variable sol for further usage.
Obtaining a 750 GeV pseudo-scalar We use the solution and our assumptions to get simpler expressions for the mass matrix of the CP-even (called hh) and CP-odd (called Ah) scalars: We omit here the analytical expressions for m 2,X A and m 2,X H because of their length and since they are not needed for the following brief discussion. The mass matrix for the CP-odd states is block-diagonal since the MSSM part is unchanged, while we have mixing in the CP-even sector among all five components. 15 The additional D-Terms can be found in the MSSM block, m 2,MSSM H . This also explains our choice of a pseudoscalar as the resonance behind the diphoton excess: the tree-level mixing between the scalar singlet and the doublets would cause tree-level decays of a 750 GeV scalar into all kinds of SM particles. In particular, those into W W and ZZ are constrained and could easily spoil our setup as an explanation of the excess in this model. Of course, we have to check whether it is possible to obtain a pseudo-scalar of the correct mass, and get the corresponding scalar sufficiently heavy so as to escape detection. For that purpose, we calculate the eigenvalues of the lower 3 × 3 block of the pseudo-scalar mass matrix, and fix B S by demanding to have a pseudo-scalar of the correct mass: we see that the CP-odd state is, as expected, mainly a singlet while the CP-even one is mainly a X-Higgs (composed by φ η and φη). That looks already very promising.
Higgs mass enhancement via non-decoupling D-terms Now, we want to confirm that one gets non-decoupling D-terms in this model which cause an enhancement of the tree-level mass of the SM-like scalar. For this purpose, we define a simple function which calculates the lightest CP-even mass for input values of mη and M Z , and create a contour plot using this function. The result is depicted in Fig. 12, where one sees that for mη M Z it is indeed possible to find a tree-level mass well above 100 GeV, while for mη M Z the tree-level mass approaches M Z . we find a tree-level mass of 38 GeV for the lightest CP-even scalar, which is mainly a mixture of η andη. It will be interesting to see if this scenario is still in agreement with all experimental constraints and how important the loop corrections are.
How to obtain a broad width? So far, we have not considered the total decay width of the 750 GeV scalar. The experimental data shows a slight preference for a rather large width of about 40 GeV, which is not easy to accommodate in weakly coupled models, typically requiring a large branching ratio into invisible states. Therefore, it would be interesting to see if this can be realised in this model. There are three possibilities for invisible decays: (i) neutralinos, (ii) (heavy) neutrinos, (iii) sneutrinos.
We are going to consider the third option here. For this purpose, we have to check two ingredients: can the mass of the sneutrinos be sufficiently light and how can the coupling to the 750 GeV scalar be maximised? To get a feeling for that, we first consider the mass matrix of the CP-even and CP-odd sneutrinos. We assume that flavour and left-right mixing effects are negligible. In that case, it is sufficient to take a look only at the (4,4) entry of the mass matrices:

MassMatrix [ SvIm ][[4 , 4]] MassMatrix [ SvRe ][[4 , 4]]
After some simplifications, we get the following expressions from SARAH: We see that the terms in the second line, ∝ T x M Z and ∝ v S λ C Y x , induce a mass splitting between the CP-even and CP-odd states. Thus, in order to have the decay A →ν IνR kinematically allowed, these terms must be individually small or cancel each other. In addition, one has to compensate the large terms ∼ M Z in order to get sufficiently light sneutrinos. This could be done by assuming a negative m 2 11 . Of course, we must check whether this leads to spontaneous R-parity breaking via sneutrino VEVs, and for this purpose one can use Vevacious, see below.
We can now check the vertex Aν IνR using the same assumptions: and we obtain after some simplification If the pseudo-scalar is a pure singlet, only the term ∝ Z A 35 contributes. This term is independent of v S and T x , i.e. we can reduce the mass splitting between the CP-even and CP-odd sneutrinos by adjusting these parameters without having a negative impact on the coupling strength to the 750 GeV scalar.

Vector-like sector
Before we finish the analytical discussion of the masses, we briefly discuss the extended matter sector. The mass matrices responsible for the mixing between the SM fermions and the vector-like fermions can be obtained from SARAH by calling

MassMatrix [ Fe ] MassMatrix [ Fu ]
which return and obtain by setting all parameters to be diagonal . There is a potentially dangerous term λ e ξ which rapidly increases for increasing λ e . To keep all scalar masses positive, it is necessary to choose a rather large B E as well. Therefore, we are going to choose always in our numerical study to circumvent tachyonic scalars.

RGEs and gauge kinetic mixing
We have so far made the simplifying assumption that gauge kinetic mixing vanishes. However, if the two Abelian gauge groups are not orthogonal, kinetic mixing would be generated via RGE running even if it vanishes at some energy scale. Thus, one of the first checks on the RGEs of the model we can make is whether the two U (1) gauge groups are orthogonal. For this purpose, we first calculate the one-loop RGEs with SARAH via

CalcRGEs [ TwoLoop -> False ];
We have chosen one-loop RGEs only to save time. Without the TwoLoop->False flag, the full two-loop RGEs would have been calculated automatically. We can now check the entries in BetaGauge and find The standard normalisation factor 5/3 for the hypercharge has been included. One can see that the β-functions for g Y X and g XY are non-zero even in the limit g XY , g Y X → 0, i.e. these couplings will be induced radiatively. Thus, in general one has not only two couplings g 1 and g X in this model, but a gauge coupling matrix G defined as In the limit of vanishing kinetic mixing, g Y X = g XZ = 0, the relations g Y Y = g 1 and g XX = g X hold. Even if gauge kinetic mixing is present, one has the freedom to perform a change in basis to bring G into a particular form. The most commonly considered cases are the symmetric basis with g XY = g Y X and the triangle basis with g Y X = 0.
The triangle basis has the advantage that the new scalars do not contribute to the electroweak VEV, and the entire impact of gauge kinetic mixing is encoded in one new couplingg. The relation between g ij (i, j = X, Y ) and g 1 , g X ,g are [370] It is interesting to see how largeg is naturally. With 'naturally' we mean under the assumption that the off-diagonal g Y X and g XY couplings vanish at some high scale Λ and are generated by RGE running down to the SUSY scale. Thus, in this setup, the size of gauge kinetic mixing is a function of Λ and g X at this scale. We can write a simple Mathematica function to get a feeling for the off-diagonal gauge couplings: In the first line, we load the file written by SARAH which provides the RGEs in a form which Mathematica can solve. This file also contains the function RunRGEs that can be used to solve the RGEs numerically. As boundary condition, we used g 1 = 0.45 at the scale 1 TeV. After the running we rotate the couplings to the basis where g XY vanishes. We can make a contour plot via and get the result shown in Fig. 13. Thus, we find that at the SUSY scale the gauge  Figure 13. Gauge kinetic mixingg at the SUSY scale as a function of the high energy scale Λ, where it is assumed to vanish, and of the coupling g X (Λ).
mixing couplingg is negative and not much smaller than an ordinary gauge coupling unless Λ is assumed to be very small.

Boundary conditions and free parameters
For the subsequent numerical analysis we are going to assume some simplified boundary conditions applied at the SUSY scale: In addition, we can set Y ν = 0 since this parameter is highly constrained to be small by the small neutrino masses. In addition, we set the mixing parameters m 2 eE , m 2 uU , M E ,M U ,B E ,B U , M 1X to zero and also assume vanishing λ, κ, and T κ . However, we stress that this is just done to keep the following discussion short and simple. All effects of these parameters can be included without any additional efforts. Thus, the free parameters mainly considered in the following are m SU SY , M λ , µ, B µ , A 0 , tan β, The tadpole equations are solved for m 2 H d , m 2 Hu , m 2 η , m 2 S and ξ, while B E and B U are fixed via Eq. (6.28).

Analysis of the important loop corrections to the Higgs mass
We now turn to the numerical analysis of this model. In the first step, we have written a SPheno.m file for the boundary conditions, see Section 6.2.6, and generated the SPheno code with the SARAH command

MakeSPheno [];
We copy the generated Fortran code to a new sub-directory of SPheno-3.3.8 and compile it via $ cd $PATH / SPheno -3.3.8 $ mkdir U1xMSSM $ cp $PATH / SARAH / Output / U1xMSSM / EWSB / SPheno /* U1xMSSM / $ make Model = U1xMSSM We now have an executable SPhenoU1xMSSM which expects the input parameters from a file called LesHouches.in.U1xMSSM. The SPheno code provides many important calculations which would be very time-consuming to be performed 'by hand' for this model, but could be expected to be relevant. A central point is the calculation of the pole mass spectrum at the full one-loop (and partially two-loop) level. In particular, the loop corrections from the vector-like states are known to be very important. However, the focus in the literature has usually been only on the impact on the SM-like Higgs. We can automatically go beyond that and consider the corrections to the 750 GeV state as well. Moreover, SPheno calculates all additional two-loop corrections in the gaugeless limit including all new matter interactions. Thus, we can check the impact of the vector-like states even at two-loop level. These effects have not been studied in any of the SUSY models proposed so far to explain the diphoton excess. Of course, SPheno also makes a very precise prediction for the diphoton and digluon decay rate of all neutral scalars as described in Section 4.4, and it checks for any potential decay mode. Thus, it is impossible to miss any important decay as sometimes has happened in the literature when discussing the diphoton excess. Finally, there are also other important constraints for this model like those from flavour observables or Higgs coupling measurements. As will be shown in the next sections, all of this can be checked automatically with SPheno and tools interacting with it.
If not mentioned otherwise, we make the following choice for the input parameters In general, the user does not need to worry about these details because SARAH/SPheno take care of them automatically. However, it might be interesting to have an intuitive feeling about the size of the different effects. Since it demands some 'hacking' of the code to disentangle the calculation in that way, we are not making this analysis here, but we briefly summarise the main results of Ref. [281] in Fig. 14. We see that all these effects can alter the Higgs mass by several GeV. Thus, an estimated uncertainty of about 2-3 GeV when using only the one-loop effective potential approximation is usually over optimistic. Furthermore, in models with non-decoupling D-terms the new loop corrections are usually neglected in the literature. Therefore, we are going to check whether this is a good approximation or not. For this purpose we show the SM-like Higgs pole mass at tree and one-loop level as a function of g X for two different values of M Z . Since SPheno performs the two-loop corrections in the gaugeless limit, additional corrections from the extended gauge sector are not included at two-loop, and we concentrate on the one-loop effects here. For this purpose, we use the different flags in the Les Houches input file from SPheno to turn the corrections at the different loop levels on or off:  for both M Z values this shift is compensated to some extent when one-loop corrections are included. Thus, the inclusion of non-decoupling D-terms only at tree-level would overestimate the positive effect on the SM-like Higgs mass by 20-30%. In addition, we also see that off-diagonal gauge couplings of a realistic size as consequence of gauge kinetic mixing reduce the positive effect from the non-decoupled D-terms on the Higgs mass by a few GeV.

Loop corrections to the 750 GeV scalar
There are also important loop corrections to all other scalars in the model if large Yukawa-like couplings are present. We discuss this briefly for the 750 GeV pseudoscalar: in Fig. 16, the mass at tree, one-and two-loop level for varying λ V ≡ λ e = λ u for two different values of m SU SY , 1.5 and 2.5 TeV, is given. For m SU SY = 1.5 TeV there is only a moderate difference between tree-level, one-and two-loop for λ V → 0, but for λ V of O(1) the one-loop corrections cause a shift by 100 GeV and more, which is compensated to some extent by the two-loop corrections. For larger m SU SY we see already a large positive shift for small λ V , which quickly increases and reaches 300-400 GeV for λ V ∼ 0.8. For even larger values of λ V , the difference between tree-level and the loop corrected mass becomes smaller. Still, the overall shift is more than 100 GeV, and this would be highly overestimated by only including one-loop corrections. As we will see in the next section, one needs λ V ∼ O(1) to explain the diphoton signal. For this value, a naive tree-level analysis gives a mass for the lightest CP-odd state which is far off the correct value. Thus, one has to be much more careful with the choice for B S .

Diphoton and digluon rate
We now discuss the diphoton and digluon decay rate of the pseudo-scalar, and its dependence on the new Yukawa-like couplings. As we have just seen, large couplings induce a non-negligible mass shift. Therefore, it is necessary to adjust B S carefully to get the correct mass, 750 GeV, after including all loop corrections. This can be done by SSP, which can adjust B S for each point to obtain the correct mass within 5 GeV uncertainty. The results for the calculated diphoton and digluon rate at LO and with the higher order corrections discussed in Section 4.2.2 are shown in Fig. 17. In order to see the size of the higher order corrections, one can use the flag 521 in SPheno to turn them on and off One finds the expected behaviour: the partial widths rise quadratically with the coupling. For about λ V 1.0 one has Γ(S → γγ)/M S ∼ 10 −6 , which is necessary to explain the observed excess. In Fig. 17 we also show a comparison between a purely LO calculation and the one including the higher order QCD corrections described in Section 4.2.2. There is no change for the decay into two photons, because its NLO corrections for a pseudo-scalar are non vanishing only for m A > 2M F . Instead, the digluon width is enhanced by a factor of 2 when including NLO and NNLO QCD corrections. This also changes the ratio of the digluon-to-diphoton width from about 10 (LO only) to 20 (including higher orders).

Singlet-doublet mixing
So far, we made some strong assumptions about some parameters in this model. In particular, we set the coupling between the singlet and the two Higgs doublets λ = 0. This raises the question how sensitive the results are to this choice. For this purpose, we can test what happens if we slightly deviate from it. The branching ratios of the CPodd scalar of 750 GeV mass, which is nearly a pure singlet, as a function of λ are shown in Fig. 18. For comparison we also show the branching ratios for the CP-even scalar with a mass around 800 GeV. This particle is mainly a mixture of η andη with a small singlet component. For both particles we depict the branching ratios when calculating only tree-level masses and when including loop-corrections. At the tree level we find that the impact of λ on the branching ratios of A is very small. This does not change much when including the loop corrections to the pseudo-scalar rotation matrix. On the other hand, for vanishing λ we already have a large branching ratio of the CP-even scalar into hh even at tree level. Moreover, the decay modes into two massive vector bosons or tt at tree level increase very quickly with λ and for λ > 0.01 they already dominate. At one-loop level, the large dependence on λ is no longer visible, because for very small λ the branching ratios into massive SM vector bosons and fermions are already large. This can be seen in Fig. 19 where we compare the doublet fraction of the two states at tree level and one loop. In general, the behaviour shows that a CP-odd scalar might be a much less fine-tuned candidate for the observed excess.

Constraints from Higgs coupling measurements
We have seen in the Mathematica session that it is possible to obtain two light scalars at tree-level. One question is: is this also possible when including all loop contributions? In order to address this question we change some input parameters to The pole masses and the doublet fraction of two lightest CP even states as a function of tan β X is shown in Fig. 20. We find the very strong dependence on tan β X , known in many U (1) extensions [370,372,373]. One difference here is that, due to the mixing with the singlet, the light state in the extended sector does not become massless for tan β X = 1, but for small deviations from it. We see in Fig. 20 that the SM-like Higgs gets a positive mass shift after the level crossing, while the mass of the lightest state drops very quickly. Of course, it is important to know if such light Higgs-like particles are compatible with all limits from Higgs searches at LEP, Tevatron and the LHC. For this purpose, we can make use of HiggsBounds, which checks whether the decay rates of a scalar into SM states are compatible with the observations at all experiments performed so far. If any of these rates is above 1 (normalised to the SM expectation), such a parameter point would be ruled out. Similarly, we can use HiggsSignals to obtain a χ 2 estimator for each parameter point, based on how well the measured Higgs properties are reproduced. In order to use HiggsBounds and HiggsSignals, we set the flag

1 # Write HiggsBounds file
in the input file for SPheno. In this way, SPheno writes out all files which are necessary to run HiggsBounds and HiggsSignals via the effective couplings input mode (effC). However, there is one caveat: SPheno does not automatically write the file MHall_uncertainties which gives an estimate for the theoretical uncertainty in the mass prediction of all scalars. The reason is that SPheno cannot do such an estimate automatically. However, if this file is missing, HiggsBounds and HiggsSignals would assume that the uncertainty is zero. Therefore, we add this file by hand and assume a 3 GeV uncertainty for all masses. We can now use this setup to make a scan in the (tan β X ,λ) plane, for instance by using the SSP option to automatically call HiggsBounds and HiggsSignals during a parameter scan. The results are shown in Fig. 21. One can see that both, the discovery potential and the χ 2 value, are very sensitive on small changes in these two parameters. The reason is mainly the large dependence of the masses of the two lightest scalars and their mixing. One sees that the best χ 2 value is found close to the tan β X range where the SM-like particle is the second lightest CP-even state, and the lightest one is about 80 GeV. In addition, for a very small stripe close to λ = 0 also points with very light scalars with masses below 40 GeV pass all constraints, but for slightly larger values of λ the mixing already becomes too large and the points are excluded by e + e − → (h 1 )Z → (bb)Z from LEP searches.

Large decay width and constraints from vacuum stability
We have already considered the possibility to enhance the total decay width of the CPodd scalar via decays in pairs of right-sneutrinos. In our tree-level analysis with SARAH we found that one can reduce the mass splitting between the two mass eigenstates by demanding In addition, as already discussed above, one has to use a negative soft-mass for the sneutrinos, 44) to get the states light enough. This immediately raises two questions: (i) how large can the total width be for large values of Y x ? (ii) Is the electroweak vacuum stable or not? First of all, we notice that a negative m 2 ν does not necessarily imply spontaneous R-parity violation, as shown in Ref. [374], in contrast to some claims in this direction in the previous literature. However, the danger of disastrous vacuum decays increases, of course, with decreasing m 2 ν . Therefore, we use Vevacious to check the stability of the potential. For this purpose, we have written a second SARAH model file where we include the possibility of VEVs for the right sneutrinos. We also added in this new implementation those mixings among states which were forbidden by R-parity conservation. This is actually necessary because Vevacious calculates the one-loop corrections to the effective potential and the full mass matrices are required. The Vevacious model file is generated via

MakeVevacious [];
We can now run a point with SPheno. If we turn on The final result is summarised in Fig. 22. 17 To maximise the effect on the total width, we take all sneutrinos to be degenerate and with the same coupling to the 750 GeV scalar. We see that we can get a large total width of the pseudo-scalar for large diagonal entries in Y x . Up to values of Y x of 0.25, which corresponds to a total width of 15 GeV, the vacuum is absolutely stable. One can even reach Y x ∼ 0.36 (Γ ∼ 30 GeV) before the life-time of the correct vacuum becomes too short. The dependence of the tunnelling time on the value of m 2 ν is shown in the middle of Fig. 22. One might wonder how dangerous this vacumm decay is, since spontaneous R-parity violation is not a problem per se. However, we show also in the right plot in Fig. 22 that the electroweak VEV v changes dramatically in the global minimum. Therefore, these points are clearly ruled out.
Even if we cannot reach a width of 45 GeV with the chosen point, we see that the principle idea to enhance the width is working very well. Thus, with a bit more tuning of the parameters, one might even be able to accommodate this value. However, this is beyond the scope of this example. We emphasise that, since the large coupling responsible for the large width is a dimensionful parameter, it will not generate a Landau pole. Thus, the large width hypothesis does not necessarily point to a strongly coupled sector close to the observed resonance.

Dark matter relic density
We have seen in the last section that light sneutrinos are a good possibility in this model to enhance the width of the 750 GeV particle. Of course, it would be interesting to see if they can also be a dark matter candidate. For this purpose, we can implement the model in MicrOmegas to calculate the relic density and to check current limits from direct and indirect detection experiments. In order to implement the model in MicrOmegas, it is sufficient to generate the model files for CalcHep with SARAH via

MakeCHep []
and copy the generated files into the work/models directory of a new MicrOmegas project. SARAH also writes main files which can be used to run MicrOmegas. For instance, the file CalcOmega.cpp calculates the dark matter relic density and writes the result as well as all important annihilation channels to an external file. This information can then be stored when running a parameter scan. The parameters are easily exchanged between MicrOmegas and a SARAH-based spectrum generator by copying the spectrum file into the main directory of the current MicrOmegas project directory. 18 However, it is important to remember that MicrOmegas cannot handle complex parameters. Therefore, one has to make sure, even in the case without CP violation, that all rotation matrices of Majorana fermions are real. This can be done by using the following flag for SPheno:  impact of this small variation on the total width is marginal, but the relic density is clearly affected. Thus, with some tuning of the parameters one can expect that it is possible to explain the dark matter relic density and the total width by light righthanded sneutrinos. However, also finding such a point is again beyond the scope of the example here. Moreover, there are plenty of other dark matter candidates which mainly correspond to the gauge eigenstatesS,X,η,η beyond the ones from the MSSM. The properties of all of them could be checked with MicrOmegas as well. A detailed discussion of neutralino and sneutrino dark matter in U (1) extensions of the MSSM and different mechanism to obtain the correct abundance was given for instance in Ref. [378].

Flavour constraints
As mentioned above, we decided to include in this model mixing terms between the extra vector-like fermions and the MSSM particles in order to let the new states decay. In this way, we have a safe solution to circumvent any potential cosmological problem. If one assumes the new coupling matrices to have a generic form, i.e. large entries of O(1), including off-diagonal ones as well, they can trigger flavour violation effects. For might take several hours, all following points should take no longer than seconds, if no new channels are needed.
instance, let us assume that Y e has the following form with degenerate diagonal entries X, and flavour violating entries α, β, γ. We can now check how strong the constraints on α, β, γ would be for given X. For this purpose, we use the results from SPheno for Br(µ → 3e), Br(τ → 3µ), Br(τ → 3e), and µ-e conversion in Ti and Au, and compare the results with the current experimental limits, see Fig. 24. We find, for instance for X = 0.1, that α must be smaller than ∼ 10 −9 , while the limits on β, γ, obtained from τ decays, can still be as large as O(10 −6 ). If other vector-like states mixing with the left-handed quarks or the right-handed down-like quarks would be present -as would be the case for instance when assuming 5 or 10-plets of SU (5) -there would also be stringent constraints on their couplings: they would cause tree-level contributions to ∆M Bs . Since these observables are also calculated by SPheno, one can easily check the limits on models featuring those states.

Z mass limits
So far, we have picked a Z mass of at least 2.5 TeV. Of course, we have to check that this is consistent with current exclusion limits. Recent exclusion limits for pp → Z → e + e − have been released by ATLAS using 13 TeV data and 3.2 fb −1 [383]. To compare the prediction for our model with these numbers, we can use the UFO model files generated by SARAH via

MakeUFO [];
and add them to MadGraph. For this purpose, we copy the SARAH generated files to a subdirectory models/U1xMSSM of the MadGraph installation. Afterwards, we generate all necessary files to calculate the cross section for the process under consideration by running in MadGraph import model U1xMSSM -modelname generate p p > Zp > e1 e1bar output pp_Zp_ee Note the option -modelname when loading the model. This ensures that MadGraph is using the names for the particles as defined in our model implementation. Using the default names of MadGraph causes naming conflicts because of the extended Higgs sector. One can give the spectrum files written by SPheno as input (param_card.dat) for MadGraph. One just has to make sure that the blocks written for HiggsBounds and HiggsSignals are turned off because the SLHA parser of MadGraph is not able to handle them. This can be done by setting the following flag in the Les Houches input file:

# ( HiggsBounds blocks )
In principle, one could also change the mass directly in the param_card without rerunning SPheno for each point. However, the advantage of SPheno is that it calculates the width of the Z gauge boson including SUSY and non-SUSY states. This usually has some impact on the obtained limits [373,384,385]. We can now scan over M Z for fixed values of g X and compare the predicted cross section with the exclusion limits. In addition, we can also check the impact of gauge-kinetic mixing: as we have seen, these couplings are negative and can be sizeable. Therefore, we compare the results without gauge kinetic mixing and when setting g 1X = − 1 5 g X at the SUSY scale. The results are summarised in Fig. 25. We see that for g X = 0.5 the limit is about 2.8 TeV without gauge-kinetic mixing. Including kinetic mixing, it gets reduced by about 200 GeV. Thus, one sees that kinetic mixing is not necessarily a small effect. This contradicts some claims that sometimes appear in the literature, where it is often argued that kinetic mixing can be ignored. In particular, we emphasise that this is very relevant when discussing a GUT theory with RGE running over many orders of magnitude in energy scale. For the dashed line, we assumed in addition g 1X = − 1 5 g X , while for the full lines gauge kinetic mixing has been neglected. The red line shows the exclusion limit from ATLAS [383].

Summary
We have given an overview on weakly-coupled renormalisable models proposed to explain the excess observed by ATLAS and CMS around 750 GeV in the diphoton channel. We have pointed out that many of the papers quickly written after the announcement of the excess are based on assumptions and simplifications which are often unjustified and can lead to wrong conclusions. A very common mistake is the lack of inclusion of higher order corrections to the digluon and diphoton decay rates, which results in underestimating the ratio typically by a factor of 2. Several authors assume that the new 750 GeV scalar does not mix with the SM Higgs, which is often not justified. Including such a mixing can give large constraints. These and other problems can be easily avoided by using SARAH and related tools which were created with the purpose of facilitating precision studies of high energy physics models. In particular, the link between SARAH and the spectrum generators FlexibleSUSY and SPheno is a powerful approach to obtain the mass spectrum and all the rotation matrices for any given model without neglecting flavour mixing, complex phases or 1st and 2nd generation Yukawa couplings. Optionally, one can also include all the important radiative corrections up to two loops. In addition, we have improved the functionality of FlexibleSUSY and SPheno to calculate the diphoton and digluon decay widths of neutral scalars, including the higher order QCD corrections up to N 3 LO. One can now pass on this information directly to Monte-Carlo tools, like CalcHep and MadGraph, by using the appropriate model files generated with SARAH.
In order to study as many models in as much detail as possible, we have created a database of SARAH model files for many of the ideas proposed so far in the literature. The database is also meant to provide many examples in the context of the diphoton excess with which the novel user can try out to familiarise with SARAH, in order to build up the level of expertise needed to implement their own models in the future.
Finally, we have introduced an attractive SUSY model which combines the idea of non-decoupling D-terms with the explanation of the diphoton excess. We have used this as a new example to show how to use SARAH to first understand the model analytically at leading order. As a second step, we have performed a numerical analysis of the important loop corrections to the different masses, checked limits from Higgs searches, neutral gauge bosons searches, and from lepton flavour violation. We have demonstrated that this model could explain a large width of the 750 GeV scalar, but in this context limits from spontaneous R-parity violation become important. These limits can be checked by using the interface to Vevacious.

A How to use SARAH and related tools
We briefly summarise here the important commands and steps to use many powerful features of SARAH. For pedagogical introductions see [221,386].

A.1 How do I install SARAH?
The installation of SARAH is very simple: the package can be downloaded from http://sarah.hepforge.org After copying the tar file to the directory $PATH, it can be extracted $ cp [ Download -Directory ]/ SARAH -X . Y . Z . tar . gz $PATH / $ cd $PATH $ tar -xf SARAH -X . Y . Z . tar . gz $ ln -s SARAH -X . Y . Z SARAH X.Y.Z must be replaced by the version which was downloaded. In the last line a symbolic link SARAH to the directory SARAH-X.Y.Z is created. There is no compilation necessary, SARAH can directly be used with any Mathematica version between 7 and 10.

A.2 How do I load a model in SARAH?
To load an existing model (called $MODEL in the following) run in Mathematica • Mass matrices : The mass matrix for a particle (Particle) is returned by

MassMatrix[Particle]
• Tadpole equations : The tadpole equation corresponding to a scalar (Scalar) is printed by using

TadpoleEquation[Scalar]
• Vertices : To calculate the vertices for a list of external states (Particles) use

Vertex[{Particles},Options];
The options define the considered eigenstates (Eigenstates -> EWSB/GaugeES) as well as the treatment of dependences among parameters (UseDependences -> False/True). All vertices for a set of eigenstates are calculated at once by

A.4 How do I get the renormalisation group equations for a model?
The calculation of the RGEs at the one-and two-loop level can be performed after the initialization of a model via

CalcRGEs[Options];
The results are saved in three-dimensional arrays: the first entry is the name of the considered parameter, the second entry is the one-loop β-function (×16π 2 ) and the third one is the two-loop β-function (×(16π 2 ) 2 ). For non-SUSY models, the RGEs of the different parameters are saved in The output for SUSY models is saved in the following arrays: $ ./ createmodel --name = $MODEL $ ./ configure --with -models = $MODEL $ make To run the spectrum generator, an SLHA input file has to be provided. The full path to the SLHA input file must be specified using the --slha-input-file= argument, for example: 20 FlexibleSUSY assumes that SARAH can be loaded from within Mathematica using the Needs["SARAH'"] command. Please refer to the installation instructions within the README file in the FlexibleSUSY/ directory for installation instructions.
A.9 How can I check points with HiggsBounds and HiggsSignals?
In the same directory in which the SPheno spectrum file is located, all other input files for HiggsBounds and HiggsSignals are saved by SPheno. The (relative) path to this directory has to be given as the last argument to HiggsBounds when executing it. Thus, working from the directory $PATH, HiggsBounds is started via: The option -modelname is used to keep the names of the particles as given in the model files. This prevents conflicts with internal MadGraph conventions. The spectrum files written by SPheno and FlexibleSUSY can be given as a parameter card to MadGraph (param_card.dat). One must only make sure that the HiggsBounds blocks are not included 21 , because MadGraph cannot parse them and would consider the file to be corrupted. If WHIZARD has not been installed globally in the home directory of the current user, WHIZARD will not be able to find the binaries. Thus, the WO_CONFIG environment variable is used to point explicitly to the binaries. By default, the configure script would install the compiled model into .whizard in the home directory of the user. If the user wants to have several WHIZARD installations or install WHIZARD locally, it might be better to provide a model just for one installation. For these cases the installation path has been defined via the --prefix option of the configure script. More information on the available options is shown with the command ./ configure --help To link WHIZARD and SPheno, all SPheno modules created by SARAH write the information about the parameters and masses into an additional file. This file is written in the WHIZARD specific format and can be directly read by WHIZARD. One just has to make sure that the corresponding flag is turned on in the Les Houches input for SPheno to get this output: To implement the new model in CalcHep, it is sufficient to use the internal "import model" routine from the GUI menu, and give the above absolute path. The model files for CalcHep are also suitable for MicrOmegas, since the latter uses CalcHep to obtain the cross section and all necessary decay widths to evaluate the dark matter abundance. To implement the model in MicrOmegas, a new project has to be created and the files have to be copied in the working directory of this project:

A.13 How do I implement a new model in Vevacious?
The model files for Vevacious are generated by SARAH via

MakeVevacious [];
As soon as the model file is created, it is convenient to copy them to the model directory of the local Vevacious installation. In addition, one can also generate a new subdirectory which contains the SPheno spectrum files for the $MODEL used as input for Vevacious, as well as the output written by Vevacious When SARAH is done with the output, the .tex file are stored in

$PATH / SARAH / Output / $MODEL / EWSB / TeX
The main file which can be compiled with pdflatex is $MODEL_EWSB.tex. SARAH usually also generates Feynman diagrams using the L A T E X package feynmf [387]. In order to compile all Feynman diagrams and the pdf file at once, a shell script MakePDF.sh is generated by SARAH.
B What is necessary to implement a model in SARAH?
All information about the model is saved in three different files: $MODEL.m, parameters.m and particles.m must be located in the subdirectory $MODEL in the Models directory of SARAH. Only the first file, MODEL.m, is absolutely necessary and contains all physical information about the model: the symmetries, particle content, (super)potential and mixings. In parameters.m properties of all parameters can be defined, e.g. L A T E X name, Les Houches block and number, relations among parameters, real/complex, etc. In particles.m additional information about particles are set: mass, width, electric charge, PDG, L A T E X name, output name, and so on. The optional information in parameters.m and particles.m might be needed for the different outputs of SARAH. We give here two examples to show that also more complicated SUSY and non-SUSY models can be defined in SARAH in a rather short form. Detailed information about the meaning and syntax are given in Refs. [221] B.1 Definition of a non-SUSY model As an example how to define a non-SUSY model in SARAH we picked the model with two scalar leptoquarks discussed in Sec. 5 C How can I define the features of a SPheno or FlexibleSUSY version?
Before we can use FlexibleSUSY or SPheno for a model, it is necessary to provide an additional input file which defines the basic setup. In general, there are two different kinds of input versions the user can create which need a different amount of input: • High-scale version: In a high-scale version of SPheno or FlexibleSUSY a RGE running between the electroweak scale, an intermediate renormalisation scale and a high (or GUT) scale is supported. The user can define appropriate boundary conditions at each of these three scales. Furthermore, threshold effects by including additional scales where heavy particles are integrated out can optionally be included. Finally, the user can specify a condition which defines the high-energy scale. The most common choice is the unification scale of gauge couplings, but other choices such as Yukawa unification are possible. In addition, these high-scale versions also include the possibility to define the entire input at the intermediate renormalisation scale and skip the RGE running to the GUT scale. The high-scale version is usually the appropriate option for SUSY models, models with heavy mass spectra, or models which should be studied at very high scales.
• Low-scale version: In a low-scale version usually no RGE running is included. The FlexibleSUSY or SPheno low-scale versions expect all free parameters to be given at the low-energy or the renormalisation scale. This version is usually used for non-SUSY models or models with light mass spectra, which should not be studied at very high scales.