Supersymmetric and non-supersymmetric models without catastrophic Goldstone bosons

The calculation of the Higgs mass in general renormalisable field theories has been plagued by the so-called “Goldstone Boson Catastrophe,” where light (would-be) Goldstone bosons give infra-red divergent loop integrals. In supersymmetric models, previous approaches included a workaround that ameliorated the problem for most, but not all, parameter space regions; while giving divergent results everywhere for non-supersymmetric models! We present an implementation of a general solution to the problem in the public code SARAH, along with new calculations of some necessary loop integrals and generic expressions. We discuss the validation of our code in the Standard Model, where we find remarkable agreement with the known results. We then show new applications in Split SUSY, the NMSSM, the Two-Higgs-Doublet Model, and the Georgi–Machacek model. In particular, we take some first steps to exploring where the habit of using tree-level mass relations in non-supersymmetric models breaks down, and show that the loop corrections usually become very large well before naive perturbativity bounds are reached.


Introduction
The Large Hadron Collider has opened a new era of precision physics. Following the discovery of the Higgs, the measurement of its properties -in particular its mass -have now been performed with an astonishing precision. This is interesting because a precise determination of the Higgs mass is of crucial importance in understanding the fate of the Standard Model (it is used to calculate the Higgs quartic coupling, required to determine whether the electroweak vacuum is metastable) and is especially sensitive to new physics beyond a e-mail: braathen@lpthe.jussieu.fr the SM (BSM). This is particularly important in supersymmetric models, where there is a prediction for the Higgs quartic coupling at tree level in terms of other fundamental parameters of the theory (notably the gauge couplings). There is therefore a long tradition of calculating higher-order corrections to the Higgs mass which was founded at the beginning of the 1990s when the dominant one-loop corrections in the minimal supersymmetric standard model (MSSM) were calculated [1][2][3]. Nowadays, the dominant two-and even threeloop corrections are available for the MSSM in the gaugeless limit, with vanishing external momenta  and the dominant momentum-dependent two-loop corrections were given in [25][26][27].
However, all of these higher-order corrections bypass an intrinsic technical problem of divergences associated with would-be Goldstone bosons of the broken electroweak symmetry. Calculations beyond the Standard Model at two loops and higher have only been performed in Landau gauge in order to decouple ghosts and thus simplify the calculations; however, in this gauge the would-be Goldstone bosons are massless and lead to infra-red divergences. In the MSSM, the gaugeless limit avoids this by turning off the Goldstone boson 1 couplings to the Higgs; and the other (momentumdependent) calculations that have been performed beyond this limit only consider the sector of the theory without the Goldstones.
However, as soon as one considers non-minimal supersymmetric models in which trilinear interactions of the Higgs superfields occur in the superpotential, the gaugeless limit no longer offers much protection against the problem, since the quartic coupling is not determined by the gauge couplings; and this is a generic feature of non-supersymmetric models (such as the Standard Model!). The so-called "Goldstone Boson Catastrophe" was noticed in the first attempt to go beyond the gaugeless limit in the MSSM at more than one loop [15], and leads to divergent values for the Higgs mass at two loops and beyond -it can in fact be a complete obstacle to a precise calculation.
Recently a solution was proposed in the context of the Standard Model [28,29] (see also [30][31][32] for recent related work) and then extended to the MSSM [33] which involved resumming (a subset of) the Goldstone boson propagators. An alternative for the Standard Model based on the 2PI effective action was proposed in [34][35][36], where essentially all particle propagators are resummed. However, both of these approaches are difficult to generalise. Instead, in [37] a general procedure was developed to cure this problem in twoloop Higgs mass calculations, based on setting the Goldstone boson propagators on-shell, which provided a complete set of modified loop functions for the tadpoles and self-energies that were finite. Thus, combining the results of [37] with those of [38][39][40][41] which provide fully generic expressions for the two-loop corrections to real scalar masses in supersymmetric and non-supersymmetric models, all ingredients are present to calculate Higgs masses in any renormalisable model.
The generic expressions of [38][39][40][41] are already used by the Mathematica package SARAH [42][43][44][45][46][47] to calculate in combination with SPheno [48,49] the Higgs masses in supersymmetric models at the two-loop level [41,50,51]. Up to now, the workaround for the Goldstone boson catastrophe in this setup was to introduce finite masses for the electroweak Goldstones by dropping the D-terms in the mass matrices. However, there were many regions of parameter space where the divergences reappeared (see e.g. [52][53][54]) and this does not work at all for non-supersymmetric models, which have no D-term potential! Therefore, to perform this work we have implemented the results of [37], in addition to filling some additional technical gaps which we describe here in Sect. 2 and the appendices; in particular, we complete the basis of required loop functions. The new version of SARAH 4.12.0 therefore now offers the possibility to calculate two-loop masses for neutral scalars in non-supersymmetric models, as well as substantially improving the calculation in supersymmetric ones. As the only non-supersymmetric model for which comparable results exist is the Standard Model, in Sect. 3 we compare our new calculation against the public code SMH [55] and the results of [56], finding excellent agreement (even if our results do not include all of the contributions included in those references). We then illustrate our new routines by computing some new results in Split SUSY in Sect. 5. On the other hand, in Sect. 4 we show how our new approach improves our previous calculation for supersymmetric models through the example of the NMSSM, for which our results should now be considered state of the art.
Momentum-independent renormalisation schemes are the most convenient choices for applying to a large variety of models, and so all mass calculations in SARAH are performed in the MS or DR scheme. In contrast, on-shell schemes might offer some model dependent advantages. This is for instance the case in supersymmetric models with Dirac gauginos and a large mass splitting between the stops and the gluino. It has been shown that in this case an on-shell scheme leads to an improved convergence of the perturbative series [57]. It is also very useful often if a DR and on-shell calculation exists for the same supersymmetric model: the difference between the results can be used as estimate of the missing higher-order corrections; this can now be done for the MSSM and certain classes of NMSSM and Dirac gaugino contributions. On the other hand, there has been hardly any discussion in the literature about radiative corrections to Higgs masses in non-supersymmetric BSM models. One reason for this, besides the technical hurdles, is that the additional freedom in non-supersymmetric models introduces a large number of free parameters, i.e. in some cases it might be possible to absorb any finite correction in the scalar sector into the counter-terms of these parameters. Thus, it is often implicitly assumed that the masses, but also the mixing angles, in the extended Higgs sector in BSM could be kept at their tree-level values. However, this is fraught with danger: (1) not all non-supersymmetric models really have a sufficiently large number of free parameters to absorb all radiative corrections. This is for instance the case in the Georgi-Machacek model. (2) If a low-energy model is combined with an explicit UV completion (such as a GUT theory), the freedom to adjust the couplings is usually lost. (3) Using masses instead of couplings as input hides the presence of huge or even non-perturbative quartic couplings. (4) Even if parameters are checked with respect to simple limits such as λ < 4π or tree-level unitarity bounds, this does not guarantee that the considered parameter point is perturbative or that strongly coupled effects do not appear at lower energies than can be explored at the LHC. Partly motivated by the growing interest in exploring quantum corrections to nonsupersymmetric models, here in Sects. 6 and 7 we explore the corrections to the Two-Higgs-Doublet Model (2HDM) and Georgi-Machacek model (GM), drawing attention to the fact that the corrections pass out of control well before the naive perturbativity or unitarity bounds.
Finally, an MS calculation has the advantage that it can give an impression of the size of the theoretical uncertainty by varying the renormalisation scale. Moreover, to obtain more reliable results for the vacuum stability by considering the renormalisation group equation (RGE) improved effective potential, a translation of masses into MS parameters is necessary. We show in this work how these aspects can be analysed in non-supersymmetric models with the new calculation available now in SARAH.

The Goldstone Boson Catastrophe and its solutions
To calculate the Higgs boson masses in general field theories we require the tadpole diagrams and self-energies. Expressions for the former were given in [41], which were derived from the general expression for the effective potential in the Landau gauge, given in [58]. Hence we must also use the selfenergies in the Landau gauge; these were given in [38] up to order g 2 in the gauge couplings, and so we restrict ourselves to the "gaugeless limit" where we ignore the contributions of broken gauge groups. This has a number of advantages, chiefly simplicity and speed of the calculation; but also the fact that we can compute the one-loop corrections in any gauge desired. Once we have dropped the electroweak contributions, it is also tempting to disregard the momentum dependence of the loop functions, which is typically estimated to contribute at the same order (and indeed is so for the MSSM [26,27]) -hence the popularity of calculations in the effective potential approach.
However, calculations in the Landau gauge/gaugeless limit suffer from the "Goldstone Boson Catastrophe", where the Goldstone bosons lead to ill-defined or divergent loop functions. Let us define the scalar potential in terms of real scalar fields ϕ 0 i and their fluctuations around expectation values v i such that ϕ 0 i ≡ v i +φ 0 i (not necessarily mass diagonal): where t i are tadpoles. Since we define the VEVs to be exact, we must have By defining the potential in terms of fluctuations, we have the MS/DR masses squared m 2 i j for all the scalars in the theory, and these are the values that enter the loop functions. However, the tadpoles are functions of the masses: and so, since these need to be adjusted loop order by order, we must choose some parameters to vary -and the standard choice is the mass-squared parameters, because in this way the couplings are unaffected. So then we define wherem 2 0,i j is the value without loop corrections (so t i = 0) and m 2 i j satisfies the full tadpole equations: where ΔV consists of the loop corrections to the effective potential, and we have written explicitly its dependence on the parameters {m 2 i j }. The Goldstone Boson Catastrophe appears because the mass-squared parameter(s) of the Goldstone boson(s) in the Lagrangian is (are) zero at tree level, but non-zero once we take into account the loop corrections to the potential. Then at two loops and higher we must calculate loop corrections with a small and/or negative mass-squared parameter, which leads to large logarithms and/or phases.
In the context of the Standard Model [28,29] (see also [30][31][32]) and MSSM [33] it was suggested that resumming (a subset of) the Goldstone boson self-energies would cure divergences in the tadpole diagrams; and including external momenta in the self-energies would also be required to cure divergences there. Alternatively, Refs. [34][35][36] proposed using the (symmetry-improved) two-particleirreducible potential to cure the problem in the Standard Model, which provides a consistent theoretical underpinning but unfortunately is particularly difficult to generalise. In the following we shall describe our previous approaches to the problem and the new results and implementation in SARAH.

Previous approaches in SARAH
Up until now, in SARAH the catastrophe appeared in an even more acute form because all of the one-and two-loop tadpoles and self-energies are computed using the tree-level masses in the loops, so without a solution to the problem, the Goldstone bosons are massless and cause several loop functions to diverge. However, for supersymmetric models the original workaround implemented in [41,50] and explored in more detail in [51] relies on the fact that the electroweak gauge couplings appear in the D-term potential. 2 We therefore used the tree-level parameters that are solutions of the full tree-level tadpole equations including the electroweak couplings to calculate the tree-level masses (but set the electroweak gauge couplings to zero in the mass matrices) used in the two-loop routines' loop functions. In other words, the masses in the loop functions are not at the minimum of the potential, and are typically tachyonic, 3 with a size of order the electroweak scale. Since we are neglecting two-loop corrections proportional to these couplings, this error is acceptable. On the other hand, for models beyond the MSSM (in particular, the NMSSM) there are typically regions of the parameter space where the Higgs sector masses still pass near to zero and cause the loop functions to diverge; for example such problems were observed in [52][53][54].
A more recent approach was to introduce regulator masses. All scalar masses in the two-loop routines which are below a certain threshold are set in terms of the renormalisation scale Q and a constant R: This approach was introduced in SARAH to stabilise cases in which the D-term approach fails. This could either be, as demonstrated in an example in Sect. 4, if other scalars artificially become very light, or if the supersymmetric scale is much higher than the electroweak scale. However, in contrast to the D-term solution, this approach violates the symmetries of the theory and can lead to non-zero masses for Goldstone bosons. Furthermore, there is no a priori indication for the optimal size of R; too large and the Goldstone/Higgs contributions are suppressed (because logarithmic contributions including them are artificially reduced), too small and the results become numerically unstable, and the user must use trial and error. Finally, it implicitly assumes that the corrections coming from the Higgs/Goldstone bosons to the Higgs mass are small (so that modifying them is benign). This is not a good approximation in many non-SUSY models, and for this reason the newly implemented solution described in the next section allows non-SUSY models to be studied accurately for the first time.

On-shell Goldstone bosons and consistent tadpole solutions
In [37] a genuine solution was presented for generic field theories: we should treat the Goldstone boson mass as an onshell parameter. A set of modified expressions for tadpoles and self-energies were given -indeed, it was shown that there were a class of loop diagrams that were not made finite purely by including external momenta. In addition, expressions for the "consistent solution" of the tadpole equations were given. 3 Since the mass was tachyonic and generally not small, we then neglected the imaginary part of the self-energies/tadpoles. These two results are closely related, as we shall elaborate a little here. If we take the Goldstone boson mass on-shell, as proposed in [37], then we have two possible ways of calculating the resulting tadpoles and self-energies, which differ in terms of how we solve the tadpole equations. The choice arises because the mass parameters m 2 i j appear on both the leftand the right-hand sides of Eq. (2.5), so we can: 1. Numerically solve Eq. (2.5) to find the m 2 i j exactly. 2. Perturbatively expand the m 2 i j so that and solve for a given loop order.
Since the effective potential ΔV will only be computed to a given loop order, the two approaches are formally equivalent. For the first approach, in practice, this means that we must iteratively solve the tadpole equations; at each iteration we put m 2 i = R ki R li m 2 kl for the tree-level mass parameters, computing a new R each time and therefore modifying the couplings, and then set the Goldstone boson mass to zero in the loop functions and compute the tadpole equations from the expressions in [37]. We find in this case that the couplings are no longer guaranteed to satisfy certain relationships imposed by the broken symmetries; only the full onshell amplitudes will satisfy the appropriate Slavnov-Taylor identities. This is only a problem for the coupling λ GG G between three Goldstone bosons, which is zero at tree level and on-shell; because the parameter in the Lagrangian will in general obtain a small non-zero value (in theories with CP-violation) and yet leads to divergent Goldstone boson self-energies we must impose that this is also on-shell (i.e. zero). Since this coupling does not appear at one loop in the calculation of the Higgs boson mass, taking this coupling to vanish causes no shift at two loops.
On the other hand, if we want to calculate the Goldstone boson self-energy at two loops then we do find a set of shifts when we take this coupling "on-shell": we would need to include the vertex corrections and define a set of shifted loop functions for those contributions (which, of course, only affect the self-energies). We shall return to this in future work.
Instead, in our implementation of the results of [37] in SARAH we take the second approach in the list above: we expand m 2 i j as a series in the couplings, and solve explicitly up to two-loop order in one step without recursion. This was already proposed in [33] for the MSSM, and in [37] explicit formulae for the corrections to the tadpoles and masses with a so-called consistent tadpole solution were given for the general case. Then we can calculate all of our loop functions using the masses (m 2 0,i j ) and couplings in the tree-level Lagrangian, and shifts Δ i j .
However, here we shall also generalise a little the expressions given in [37]: we shall allow Δ i j to be an implicit function of the tadpole shifts, rather than explicitly assuming ; indeed, this equation fails for pseudoscalars, for example. Instead we solve the tadpole equations for some variables {x i } with For example, in the Goldstone model of a single complex scalar Φ having potential when Φ obtains a VEV it decomposes into a real scalar h and a Goldstone boson G as Φ = 1 We solve the tadpole equations for the parameter μ 2 so that However, both the mass of the Goldstone boson and the Higgs are controlled by the μ 2 parameter; (2.11) So in our notation, x h → μ 2 , c 0,h → −λv 2 , c hh → − 1 v and so Substituting the Δ i j into the one-loop tadpole and self-energy expressions then gives a set of two-loop shifts; for scalars these were given in Eqs. (5.2) and (5.3) of [37]. We have implemented this under the assumption that the variables {x i } are dimensionful and there is no explicit dependence of the trilinear/quartic couplings on them (only implicitly through the mixing matrices R); and also we assume that the fermion mass matrices do not depend on these parameters. These assumptions are fulfilled e.g. for {m 2 H u , m 2 H d } in the MSSM, but not {μ, B μ } chosen as parameters to solve the tadpole equations. On the other hand, we give expressions for the shifts to the tadpoles and self-energies when fermion masses depend on the {x i } in Appendix B, and we plan to implement these in the future. Now, since the Goldstone boson is massless at tree level, this then means that we automatically have the Goldstone boson "on-shell." This means that the "on-shell" and "consistent solution" approaches are more closely related than first appears: since the Goldstone boson mass must be zero on-shell and we can identify the Goldstone boson eigenstates using a matrix R kG derived just from the broken symmetries (see e.g. [37]) then the on-shell condition becomes and sincem 2 0,GG = 0 we have (2.14) i.e. the approach of adjusting the loop functions (as we do when setting the Goldstone boson on-shell) or defining a set of shifts to the tadpoles and self-energies involving Δ i j should give the same result when we just consider the shifts to the Goldstone boson masses, even though the expressions look very different.

A complete basis of loop functions and the implementation in SARAH
For the evaluation of tadpoles and self-energies [37] proposed a "generalised effective potential limit," where the selfenergies are expanded in s = −p 2 (= m 2 on-shell) and all terms of order O(s) are neglected (but crucially retaining terms that diverge at s = 0). We therefore require the following basis of loop functions, where {x, y, z, u, v} = 0 are masses squared: Momentum independent: J (x), P SS (x, y), P SS (0, y), I (x, y, z), I (0, y, z), I (0, 0, z), (2.15) All of these functions are implicitly dependent on the renormalisation scale Q, typically containing factors of log x ≡ log(x/Q 2 ). Expressions for all of these functions expanded up to O(1) in the external momenta (or the reference for them) were given in [37]. However, the functionsṼ (x, y, z) were given in terms of the regularised function V (u, 0, y, z) defined in [38]; 4 in Appendix A we derive explicit compact expressions for this function -first with full momentum dependence, and then expanded up to O(1) in the external momenta.
In our practical implementation in SARAH we have extended the available routines for calculating two-loop integrals with the missing ingredients to address the Goldstone boson catastrophe. Moreover, there are three loop functions involving fermions and gauge bosons which needed modification for the MS scheme as used for non-supersymmetric models, as compared to the DR for supersymmetric models; the tadpole and self-energies contain whereV (2) is the two-loop contribution to the effective potential, d(I ), C(I ) are the dimension and quadratic Casimirs of representation I of the gauge group having coupling g, and the loop functions are 4 Note that ourṼ uses one fewer variable than the definition of V in [38] because in our case it always appears in the form V (u, 0, y, z).
Here δ MS is one for MS masses and zero for DR . In addition, routines to calculate the consistent tadpole solution are generated during the output of SPheno code. This is fully automatised beginning with SARAH version 4.12.0 and the user can obtain a SPheno version for nonsupersymmetric models as before -with the difference that two-loop mass corrections are now included. We refer for more detailed explanations of how to use the code to the standard references such as [47]. The only requirements are recent versions of SARAH and SPheno which are available at http://www.hepforge.org.
The new features can now be adjusted in the Block SPHENOINPUT in the Les Houches input file: Note that the solution to the Goldstone boson catastrophe exists only for the diagrammatic calculation (flag 8 → 3), but not for the effective potential calculations using numerical derivatives to obtain the tadpoles and self-energies (flag 8 → 1, 2). By default, the new calculation is used now, but could be turned off if demanded (flag 151 → 0). In this case, it is usually necessary to include a non-zero regulator mass via flag 410 for non-supersymmetric models. In principle, there should not be any reason to revert to the old calculation with regulator masses except for double-checking the result.
The consistent tadpole solution (described in the previous section) is turned off by default but can be turned on by setting flag 150 → 1. This is because, while strictly it is more accurate to include it, there is also the possibility of numerical instability if the shift in the tree-level mass parameters is large; for example, if the expectation values of some scalars are small (such as e.g. the neutral scalar of an electroweak triplet which must have a small expectation value from electroweak precision constraints) then the shift in the mass parameter can be much larger than the tree-level value and the perturbative solution fails. In such cases, it would be better to use a recursive approach which is currently not possible for the reasons given in Sect. 2.2.

A first comparison of our results with existing calculations
Now that two-loop corrections to scalar masses are available in SARAH, free of the Goldstone boson catastrophe, it is important to compare the results we obtain with other computations available in the literature, as a verification of our results and as a way to estimate the impact of missing corrections. We consider in this section the Higgs mass calculations in the Standard Model, and we will compare the results obtained with SPheno with the computations performed at complete two-loop calculation in [56], and the full two-loop (plus leading three-loop) Higgs mass calculation implemented in the public code SMH [55]. These works take into account two-loop electroweak corrections, which are not available for generic theories and are not included in our code, hence we will quantify the size of these effects, together with effects from momentum, and investigate the discrepancy in masses coming from the different determination of the top Yukawa coupling. It is interesting to examine the way that the two calculations avoid the Goldstone Boson Catastrophe. The calculation of Buttazzo et al. [56] was performed in Feynman gauge and using certain parameters on-shell, whereas the results implemented in SMH are in a pure MS scheme and Landau gauge, which is closer to our approach. In the latter paper, some resummation is performed by hand to eliminate the divergence in the mass calculation; it is perhaps surprising that the absence of the function V (0, 0, y, z) from the basis in TSIL was not problematic, but there the calculation was performed by computing the set of integrals explicitly using TARCER [60] rather than starting from a set of generic expressions, so the result was found directly in terms of the other basis functions. In principle this should agree with Eq. (A.10).
For clarity, we recall that we define the (tree-level) Higgs potential as (3.1) A first approach for the comparison between SPheno and [56] is to compute the Higgs mass with the quartic coupling λ ranging in the interval [0.125, 0.130], and only setting the SM inputs to the same values as in [56], which we recall here They furthermore took the experimentally determined central value of the Higgs mass to be 125. 15 GeV, which we shall take as a reference value rather than an input. The use of consistent solutions to the tadpole equations -as derived in [37] -has also been implemented in the SPheno code and this comparison in the context of the SM is a good occasion to study the effect of this additional shift to the tadpoles and mass diagrams, thus we compute the Higgs mass in this first method both with and without using the consistent tadpole solutions. A second approach to computing m h with SPheno, which could potentially improve the comparison, is to use as well the same values for the top-Yukawa y t and electroweak gauge couplings g 1 , g 2 as those given for each order in Table 3 of [56]. We obtain another result for m h with SMH [55], and although this code is made to perform Higgs mass calculations in the Standard Model to partial three-loop order, we use it here with the three-loop corrections always switched off, for the purpose of our comparison with SPheno. We use the routine calc_Mh that gives for a given loop order the value of m h from the inputs of the renormalisation scale Q, the quartic coupling λ, the top-Yukawa y t , the Higgs VEV v, and the gauge couplings g 3 , g, g , all given at scale Q. In order to improve the comparison, we take the same values for the inputs as used at each order in SPheno. We give in Table 1 the values we find for the Higgs mass when taking the same values of λ as found in [56], with the two methods described above for SPheno and with SMH.
At tree level, all the values we find with SPheno and SMH obviously match as the tree-level Higgs mass only depends on λ and v which have almost the same values here, and the divergence from the value of 125.15 GeV is solely explained by the Higgs VEV which is not the same as in [56] since they take it as an on-shell parameter, while we use the MS value as described in [61]. More importantly, the loop-corrected values in the different methods also agree quite well, thanks to the improved determination of the top Yukawa coupling y t (including leading two-loop effects) recently implemented Table 1 Values of the Higgs mass at scale Q = m t for the values of the quartic couplings λ found in [56] at tree level, one loop and two loops, in the two approaches we used for SPheno and with SMH. The first approach was to change only the SM parameter inputs while letting SPheno determine the top-Yukawa and electroweak gauge couplings, and the Higgs mass is computed both with and without the consistent tadpole solutions. The second method was to take the same values of y t , g 1 , g 2 in SPheno as in [56] (and switch off the consistent tadpole routines). For SMH, the values of the input parameters -the top-Yukawa, the electroweak gauge couplings, the Higgs VEV and the strong gauge coupling -were taken from the outputs of the SPheno scans. Computations are made with SARAH-4.12.0, SPheno-4.0.3 and SMH-1.0 [55] Loop order Value of λ found in [56]  The different value of the Higgs VEV is also quite certainly the main reason for the discrepancies between the values we obtain using SMH and those from [56]. I.e. it is because we use the VEV computed in SPheno in SMH, which does not correspond to the same accuracy of parameter extraction as used in [56], which would be required for a fair comparison directly between the two prior approaches: here our aim was to compare our result separately with [56] and SMH.
A further way to compare our results to those of [55] and [56] is to find for each order what value of the quartic Higgs coupling we need to obtain m h = 125.15 GeV, and our results are given in Table 2. We observe that the change of λ between each order of the perturbation expansion is approximately the same in all four methods. Moreover, the value we extract at two loops with SPheno is very close to the value found in [56], only differing by 0.1%.

A detailed comparative study of SPheno and SMH results
After this first comparison, we may now investigate in more depth the effects of the three sources of differences on the Higgs masses listed above, using SPheno and SMH.
To begin with, we should consider the Higgs VEV and its calculation: in SMH, calculations are performed in the Landau gauge, while SPheno is by default set to use the Feynman gauge, and while the Higgs mass should in principle be gauge independent, its vacuum expectation value is not, hence there is an inconsistency coming from the use of a Feynman gauge VEV in SMH. The easiest way to correct this is to switch the SPheno calculation to the Landau gauge -we set in the code the gauge parameter ξ to a very small finite value to approach the limit of the Landau gauge (the current implementation gives a numerical divergence when ξ = 0) -and then to use the new value of the Landau gauge VEV in SMH.
The values we find for m 2 h with the two codes for the two different choices of gauge parameter and fixed values of Q and λ are given in Table 3. The first observation that can be made from these results is that the Higgs mass shows residual dependence on the gaugem 2 h varies by about 50 MeV between ξ = 1 and ξ = 0.01. This is explained mainly 5 by the difference in the calculation of the MS value for the electroweak VEV in SPheno between the Feynman gauge and other gauges: in the case of Feynman gauge one-loop corrections from δ V B as well as two-loop corrections from δ r are included which are not available in general ξ (see Appendix A of [61] for details of the matching in Feynman gauge). On the other hand, in ξ gauge, the VEV is calculated from What is more interesting is that the agreement between the two codes improves greatly once we use the Landau gauge in SPheno; indeed the difference in the Higgs mass results is reduced from approximately 0.4 GeV to less than 0.05 GeV.
A second point we can study is the effect of the two-loop momentum dependence and two-loop electroweak corrections. Let us introduce the notation for calculating the pole mass via where div f (s) denotes all terms in f (s) that diverge as s → 0. Our SPheno code computes the one-loop corrections in any R ξ gauge with full momentum dependence, but the two-loop corrections are performed in a generalised effective potential approachi.e. we keep only the divergent part of the momentum dependence (see Sect. 4 of [37] for more details). The momentum in the two-loop routines is fixed (for speed of calculation) whereas that in the one-loop routines is adjusted to solve the on-shell condition: hh,gaugeless (s) This begs the question of how to compare our result with SMH: ideally, we would like to extract a result from SMH which is comparable to ours. However, this is confounded by several factors: 1. It is impossible to remove or extract the electroweak contributions in SMH, because the individual contributions in the computation diverge as the electroweak gauge couplings become zero; and the total result is not finite at vanishing external momentum. 2. To avoid the Goldstone boson catastrophe and ensure cancellation between Goldstone boson and longitudinal gauge boson diagrams, in the two-loop corrections in SMH the external momentum s has been replaced by 2λv 2 wherever it appears in a pre-factor (but not in the arguments of the loop functions).

The term
which is part of the one-loop correction coming from Goldstone bosons and longitudinal gauge bosons, is moved into the two-loop corrections, with the justification that on-shell s = 2λv 2 + Δ (1) M 2 h so will give a contribution at two-loop order when solving for the onshell mass.
If it were not for point (2) above, it would perhaps have been possible to extract the result for the generalised effective potential approximation for the electroweak corrections. Instead, we will simply compare the results as we vary the momentum in SMH; by modifying slightly the source code, we obtain a version of SMH without the momentum dependence at two loops (but retaining the dependence at one loop). Interestingly, the result of SMH is finite even when s = 0 meaning that the divergence as s → 0 has been removed. It turns out that this is because of the term (3.6), which has the effect of cancelling the divergences as s → 0 (even though this cancellation is fictitious). If we write δ (2) (s) for the missing momentum dependence in SMH from setting the coefficients of loop functions equal to 2λv 2 , then we have (2) hh,gaugeless (s) + Π (2) hh,electroweak (s) + δ (2) The cancellations of the divergences imply that div Π (2) hh,gaugeless (s) + Π (2) hh,electroweak (s) + δ (2) (s) On the other hand, by evaluating the diagrams for the Standard Model in the gaugeless limit retaining only the top Yukawa coupling and the Higgs quartic λ we find div Π (2) hh,gaugeless (s) (3.9) so we can see there are several remaining pieces that must be cancelled by div[Π (2) hh,electroweak (s) + δ (2) (s)]. But if we set s 2 fixed = −Q 2 in our routines we should cancel the divergent part exactly, and leave us only with Π (2) hh,gaugeless (0). We can then determine We find that this residual difference is tiny; at Q = m t = 173.34 GeV with λ = 0.12604, y t = 0.9345, v = 247.07 GeV and the gauge couplings This corresponds to a tiny value of the electroweak corrections; a similar observation was made in [56]. Finally, we compare the more physically meaningful differences between the codes when we take s = m 2 h | tree in our routines. The values of the Higgs mass computed with SPheno after turning off the light SM fermion contributions and with the modified version of SMH is given in Table 3, and strikingly they only differ by 40 MeV when we include the momentum dependence in SMH -in other words, for Table 3 Comparison of two-loop Higgs masses calculated with the codes SPheno and SMH, for different choices of gauge in SPheno and switching on and off the two-loop momentum dependence in SMH. The renormalisation scale is fixed to Q = 173.34 GeV, and the Higgs quartic coupling is λ = 0.12604 and is not varied (the idea being to illustrate the importance of the choice of scale, rather than the stability of the result). All other inputs for SMH are taken to the same values as in SPheno. In SPheno the only two-loop momentum dependence is from pseudo-scalar diagrams and only a generalised effective potential approach (see main text) with s = m 2 h | tree , while in SMH the full two-loop dependence is implemented and is used to find m h iteratively   Fig. 1 shows the difference of the two-loop masses between the two codes -more precisely (m 2 h ) SPheno − (m 2 h ) SMHwith and without momentum, as a function of the renormalisation scale Q (where the MS parameters are extracted by SPheno at each value while keeping λ fixed rather than evolving the parameters: the idea is to show the importance of the choice of scale rather than the stability of the computation). While for large scales the two-loop momentum effects may become large (1 GeV or more), the electroweak corrections represent at most 0.2 GeV and even vanish for a scale close to the MS top mass.

Momentum dependence
Implementing the solution to the Goldstone boson catastrophe in SARAH has required the insertion of external momentum in infra-red divergent loop integrals, and thus we should also investigate the impact of the momentum s = −p 2 on the Higgs mass calculation in SPheno. In practise, we have set for the majority of scans the external momentum for the twoloop calculations to be equal to m 2 h | tree but we will now vary the momentum to study its impact on m h . Table 4 shows the shift to the two-loop Higgs mass -with respect to the value computed with s = (125 GeV) 2 -for external momentum in loops equal to s = α × (125 GeV) 2 , where α ranges from 10 −6 to 10 6 and for λ = 0.126 and λ = 0.130. For all values of the external momentum considered here, the variation of the Higgs mass remains small: at most they become of order ∼ 0.13 × 1 GeV for α = 10 −6 (i.e. √ s = 0.125 GeV), and while this effect is noticeable, it is far from the divergences that could have been feared when approaching the limit of s → 0. All in all, although pole masses -as we compute here -are in principle found as the zero of the inverse propagator, that has to be found iteratively as the self-energy contains momentum dependence, we see from the minute effects of momentum in the range α ∈ [1/2, 100], relevant for scalar masses, that we will not require an iterative solution and that simply taking s = (125 GeV) 2 in the loop diagrams with pseudo-scalars will be a satisfactory approximation. In particular, changing s between m 2 h | tree and 125 GeV causes a difference in m 2 h of less than an MeV. We emphasise, however, that the effect of momentum on Goldstone boson mass diagrams discussed here is only a subset of the general momentum dependence of the two-loop masses, which should in principle be taken into account, as seen in the previous sections.

The NMSSM
As a second check of our new solution, and demonstration of its importance, we shall compare the results for the three different options to solve the Goldstone Boson Catastrophe in the example of the next-to-minimal supersymmetric standard model (NMSSM) -see [62] and the references therein for a detailed description of the model. Indeed, the NMSSM is the first supersymmetric model for which the problems at certain points in the parameter space were found in earlier versions of SARAH. Here we shall show that this is avoided, and have a preliminary look at the impact of the "consistent tadpole solutions".
We start with a test point defined by the following input parameters: 6 λ = 0.7, κ = 0.25, A λ = 1350 GeV, and all other soft-masses set to 2 TeV. The Higgs masses for the following calculations are given in Table 5: 1. D-terms turned off in mass matrices but retained in tadpole solutions (as in previous versions of SARAH), labelled "D" in the table. 2. Regulator masses with R = 10 −5 -10 −1 .    While the new "on-shell" solution of the Goldstone boson catastrophe is optimal, between introducing a regulator R and the previous approach with neglected D-terms in the scalar mass matrix, the latter is preferred because one does not need to check for a suitable choice of R to stabilise the results. However, we can now consider parameter points where the old method fails. This is shown for the point defined by and all scalar soft-masses set to 2 TeV. The lightest scalar tree-level mass with and without the D-terms as a function of λ is shown in Fig. 2. One can see that for λ 0.5, 0.8, the lightest scalar becomes massless in the limit of vanishing D-terms. Thus, for these values, divergences in the two-loop corrections can be expected which are this time not associated with the Goldstone but with the lightest CP-even state. We show the lightest Higgs mass in Fig. 3 as a function of λ for different methods to regulate the two-loop corrections. Obviously, the approach of neglecting electroweak D-terms fails for values of λ at which the masses entering the loop calculations become very light. However, for very large values of λ which are away from the poles, the agreement with the other calculations is also rather poor. In contrast, over the entire range of λ we see good agreement between the methods using regulator masses, if R = 10 −2 or 10 −3 is chosen, and the method of treating the Goldstones on-shell. It is interesting that for these values of R the minimum mass is √ R × M SUSY 100 GeV, i.e. logarithmic contributions involving the light scalars are being excised.
We note that the corrections from the consistent tadpole solution are small until λ becomes large, at which point we see significant deviations. However, as λ approaches 0.9 we see from Fig. 2 that the tree-level lightest Higgs mass approaches zero, so we expect our perturbative calculation of the "consistent tadpole solution" to break down and become unreliable.

Split SUSY
In Split SUSY scenarios [63][64][65][66][67][68], the SUSY scalars are much heavier than the gauginos and Higgsinos. Consequently, these models should be studied in an effective approach where all SUSY scalars are integrated out at some matching scale. The Lagrangian below this scale is given by where L SM is the Standard Model Lagrangian with Higgs potential (3.1). Because of the matching between the effective, non-supersymmetric model and the MSSM, the quartic Higgs coupling λ as well as the new Yukawa-like interactions g (1,2)(u,d) are not free parameters but fixed by the matching conditions at the scale M M . At tree level, the following relations hold: Here, g 1 and g 2 are the running gauge couplings of U (1) Y × SU (2) L and β is defined as the mixing angle of the two Higgs doublets in the MSSM (in contrast to the definition in the MSSM as a ratio of expectation values). There are important higher-order corrections to the matching conditions which are necessary to have a precise prediction for the Higgs mass at the low scale. In particular λ has been calculated including the two-loop SUSY corrections [69][70][71]. The numerical value of these corrections depends on the many SUSY parameters at the matching scale; however, a commonly taken useful approximation is to give the scalars a common mass M M , in which case the corrections can be given in terms of just this scale and the squark mixing. Moreover, in strict split SUSY where the fermion masses are protected by an R-symmetry (or another symmetry in Fake Split SUSY [72,73]) near the electroweak or TeV scale and well below M M , the squark mixing must by very small. In which case, the leading corrections to the Higgs quartic coupling are purely electroweak at one loop, and at two loops contain no logarithmic terms -meaning that they are very small (in particular since the strong gauge and top Yukawa couplings run to small values at higher scales), so using the tree-level relationship above can be good enough.
Below the scale M M , we must run to the scale of the fermion masses, before also integrating them out, and then running to the electroweak scale in the Standard Model. In some previous approaches, the running was performed all the way down to the electroweak scale, before calculating the Higgs mass in the full Lagrangian (5.1); however, it was found in [72] that in this approach it is necessary to include the three-loop leading logarithm involving the gluino mass to obtain good agreement between the two results -this is automatically resummed by the renormalisation group running in the former approach. In either case, the full contribution of the gauginos and Higgsinos to the matching conditions is only known in the literature to one-loop order [69].
Hence in this section we are interested in the effect of the two-loop corrections to the Higgs mass stemming from theg (1,2)(u,d) couplings which have not been studied in the literature before. They are expected to be small since they originate from electroweak interactions at the matching scale (and so, admittedly, one could argue that we should neglect them in the gaugeless limit). We shall not discuss the absolute value of the Higgs mass, for which we would need to include all higher-order corrections to the matching that have been calculated elsewhere, but only on the impact of the new two-loop corrections. The overall size of these corrections is rather insensitive to the exact matching conditions and we are using the above tree-level relations; but as we noted earlier, these should be a particularly good approximation for larger matching scales.
We make in addition the simplifying assumption that at M M the SUSY fermions are degenerate, i.e.   Fig. 4 for tan β = 1 and 10 for a calculation with and without the consistent tadpole solutions explained in Sect. 2.2. We show here results for M F up to 5 TeV. In order not to increase the theoretical uncertainty in the presence of new fermions in the multi-TeV range, we made use of the functionality in SARAH to perform the Higgs mass calculation in the effective SM [61]. For this purpose, a second matching is performed to extract λ at the renormalisation scale Q. The imposed matching condition is i.e. we perform a matching of the Higgs pole masses as suggested in [54], from which an effective λ is derived. λ is then evolved to m t using three-loop RGEs of the SM. At m t the Higgs mass is calculated within the SM at the two-loop level. The additional loop corrections discussed here enter the calculation of m Split h (M F ), i.e. the calculation of λ SM (M F ). We see that the additional corrections for SUSY fermions are always well below 1 GeV once the consistent solution to the tadpole equations are included. However, if those are not used, the misleading impression of sizeable corrections of a few GeV is given; it would be interesting to investigate this phenomenon further.

Two-Higgs-Doublet Model
The scalar potential of the CP-conserving 2HDM is defined in terms of scalar SU (2) L doublets in a basis {Φ 1 , Φ 2 }sometimes called the Z 2 basis -as One or both doublet(s) Φ 1 and Φ 2 may acquire VEVs if m 2 i j has one or two negative eigenvalues, and we write the doublets and their VEVs as We then define the angle β through the usual relation where v is defined by v 2 = v 2 1 + v 2 2 . The 2HDM hence has seven free parameters, which are λ i (for i ∈ {1, 2, 3, 4, 5}); m 2 12 ; tan β. (6.4) It is often more convenient to work in another basis -the so-called Higgs basis {H 1 , H 2 } -where the neutral component of the doublet H 1 is aligned in field space with the total VEV v, with a rotation of angle β We choose to write these two new doublets as In this new basis, following the notation of [74], the potential can be written as The CP-even physical states are eigenstates of the mass matrix which is diagonalised with an angle α, and they are given by The alignment limit is defined as the limit in which the neutral components of the Higgs-basis doublets are also mass eigenstates, or in other words, the limit in which one of the CP-even neutral scalar mass eigenstates is aligned with the VEV v. From Eq. (6.10) we see that this can be realised in two ways: 1. s β−α = 0 in which H carries the VEV and is identified with the SM-like Higgs. 2. c β−α = 0 which means that h is the SM-like Higgs.
We do not make any assumption on the size of the masses of the different scalars i.e. we do not suppose that we are in the decoupling limit as well. Consequently, at tree level we only require Z 6 v 2 → 0, and hence with the expression of Z 6 derived in [74], we have where λ 345 ≡ λ 3 + λ 4 + λ 5 . The simplest, and tan βindependent, way to fulfil this condition is to have which we will use in the following to constrain tree-level alignment. Also, we will require that the SM-like Higgs be the lightest mass eigenstate h (case (2) above), by ensuring that This implies that c β−α = 0, and thus, with the conventional choice that β ∈ [0, π 2 ] and |α| ≤ π 2 , we have 14) The constraint for tree-level alignment given in Eq. (6.12) reduces the number of free parameters of the model from seven to five, as two of the quartic couplings (e.g. λ 2 and λ 3 ) can be found as a function of the three other ones. For most scans and figures presented below, we worked in the type-I 2HDM if not indicated otherwise. However as the difference with type-II comes from the couplings of the scalars to the down-type quarks and to the leptons which are light and give much smaller contributions to the lightest Higgs mass than the top quark, we do not expect large effects on our results (even for large tan β, since the contributions typically involve the quark masses rather than just the couplings).

Renormalisation scale dependence of the Higgs mass computed with SPheno
The masses computed by SPheno are pole masses, which should in principle not depend on the renormalisation scale at which they are computed. Evaluating the variation of the masses with the scale Q hence provides a consistency check of our results and an estimate of the theoretical uncertainty as the variation of the two-loop masses with Q corresponds roughly to the three-loop corrections. For this purpose, we have tuned the λ i couplings to ensure a two-loop Higgs mass of 125.09 GeV, at scale Q = 160 GeV, together with tree-level alignment and find the following values (using HiggsBounds we have verified that this point in parameter space is not excluded by the current experimental constraints): At first we consider that these inputs are then given to SPheno as the value of the couplings at a scale Q that we vary in the range [100 GeV, 10,000 GeV], and we only consider the running of SM parameters; we find the results shown in Fig. 5 for the tree-level, one-loop and two-loop Higgs mass m h . Since phenomenological analyses typically supply the quartic couplings without reference to a higher-energy theory or the scale where they are determined, this plot shows the importance of the choice of that scale. We have verified that the renormalisation scale dependence of m h | tree is entirely due to the scale dependence of the Higgs VEV v, as the running of the quartic couplings is for the moment not applied. The renormalisation scale Q is seen to have only a limited effect on the two-loop value of m h which varies of about 2 GeV on the range of scales considered here, while the one-loop result varies by about 15 GeV. Since the twoloop curve is so flat, this shows that most of the variation in the calculation of the Higgs mass for the chosen quartics must come from variation of the Standard Model parameters, and that a two-loop calculation (rather than one-loop) is necessary not just for precision but also to ensure scale stability.
Using two-loop RGEs implemented in SARAH/SPheno, we can also include the evolution of the 2HDM parameters to obtain a more complete scale dependence of the masses, as shown in Fig. 6. Once more, the two-loop value of Higgs mass depends less on the renormalisation scale than the tree-level or one-loop values. This weaker dependence of the two-loop Higgs mass on Q, compared with the one-loop mass, even for choices of parameters that give large loop corrections is a first verification of the validity of our new two-loop routines. In the following we will therefore work at a fixed scale Q = m t , confident that the results will be for the most part independent of this choice. The relations defining the alignment limit in the beginning of this section are only valid at tree level and we expect them to receive corrections at one-and two-loop order, and in this section we will discuss the importance of these effects on the mixing angle of the neutral CP-even scalars α.
Scanning over the different free couplings of the model m 2 12 and λ i (i ∈ {3, 4, 5}) -we compare the values of the CP-even Higgs mixing angle α at tree level, one-loop and two-loop order, as shown in Fig. 7, and as expected, loop corrections cause deviations from the tree-level relation t α = −1/t β ⇔ c β−α = 0. The observations we can make from these plots are the following: 1. in the ranges of the parameters that we considered, the effect of loop corrections on the value of α is small, at most of the order of 1%; 2. the one-loop corrections to α show very little dependence on the quartic couplings λ i=3,4,5 ; 3. it appears that, for most parameter points, the two-loop corrections to α are of similar magnitude than the oneloop ones -although somewhat smaller when |λ i | 1; 4. for some parameter points, however, the two-loop corrections to α become significantly larger than the oneloop corrections; see the lower right plot in Fig. 7. We have verified that this happens when one of the quartic couplings λ i becomes large (typically |λ i | 1) -in the plot mentioned above of −1/t α as a function of λ 5 it is λ 4 that becomes smaller than −1. We may suspect the large two-loop effects are due to a loss of perturbativity: this will be discussed in more detail in the next section.

Perturbativity constraints
It is common in practice to use the physical scalar masses, the Z 2 breaking parameter m 12 as well as the angles α, β as input for the 2HDM in numerical studies. However, this input often hides that it corresponds to huge quartic couplings which spoil unitarity and the perturbative behaviour of the theory. Therefore, the constraints that all quartics must be smaller than 4π as well as the tree-level unitarity constraints [75][76][77] are applied to sort such points out. However, it was already shown in the SM that the limit of λ < 4π might be too weak [78]. We now have all the machinery at hand to impose another constraint on the 2HDM model namely that the radiative corrections to the Higgs mass converge. We show here in one example that this can be a much stronger constraint than tree-level unitarity, while a more detailed analysis of this constraint on the parameter space of 2HDM models is left for future work.
We Since the masses are treated as pole masses and only treelevel relations are used in the above work, no scale for the MS parameters is given. On the other side, it is usually checked that the translation of these masses into quartic couplings results in parameters which are allowed by tree-level unitarity. However, this treatment implicitly assumes that one can define at each loop level suitable counter-terms to renormalise the Higgs sector in a way that the masses can be kept constant, and that this renormalisation converges. This is, however, not the case for the parameter point defined by Eq. (6.16) as one can see as follows. The given input translates into the following set of quartic interactions using the tree-level relations 8 : These fulfil the tree-level unitarity constraints [75][76][77].
To check the perturbative behaviour, we show the scale dependence in Fig. 8. Here, we used the quartic couplings of  Eq. (6.17) as input and checked the scale dependence of the Higgs mass at different loop levels. For the evaluation of couplings to the considered scale, we used the two-loop RGEs calculated by SARAH. One sees that the scale dependence increases with increasing loop level. Of course, one might wonder if this is just an effect from our choice to define the quartic couplings at Q = m t as input. Therefore, we show in Fig. 9 the size of the loop corrections for different choices of our input scale Q. We see that the size of the loop corrections rapidly increases for Q > m t and the spread between one-and two-loop becomes even larger. Also choosing Q 160 GeV where the mass at one and two-loop level Fig. 9 The size of the one-and two-loop corrections of the lightest scalar mass as a function of the scale Q in at which the input masses of Eq. (6.16) are translated into quartic couplings seem to be roughly identical does not solve the problem: this is just a numerical coincidence and the scale dependence at two loops is even larger than at one loop.

Georgi-Machacek Model
The Georgi-Machacek Model [82] extends the SM by one real scalar SU (2) L -triplet η with Y = 0 and one complex scalar SU (2) L -triplet χ with Y = 1, which can be written as A very compact form to write the Lagrangian in a SU (2) L × SU (2) R invariant form is to express the doublet and triplet scalars as Here, φ are the components of the SM doublet. Using this notation, 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 [83]. The triplets obtain VEVs as where the custodial symmetry enforces v η = v χ ≡ v T , and there are no tree-level contributions to the ρ parameter. They further fulfil v 2 φ + 8v 2 T = v 2 , which allows us to define The free parameters of the model are then since μ 2 2 , μ 2 3 can be eliminated by the tadpole equations. The physical eigenstates can be organised into representations of the custodial symmetry as a five-plet (consisting of a doubly charged, singly charged a neutral CP-even scalar), a triplet (consisting of a singly charged and a CP-odd neutral scalar) and two CP-even singlets (where the Standard Model Higgslike boson is the lighter of the two). Expressions for the triplet mass m 3 , five-plet mass m 5 and singlet masses are given in, for example, [84]. m h , s H , m 5 seem to be a suitable choice for the input parameters and can be traded for λ 1 , λ 5 and v T . In the following we shall do this using tree-level relations derived from those in [84]. However, the choice to use masses instead of couplings as input can have the danger that one enters a nonperturbative regime without recognising it as we already have pointed out for the 2HDM. We will discuss the importance of higher-order corrections in general in this model in the following: in contrast for instance to the 2HDM, it is not possible to renormalise all mixing angles and masses onshell in this model. One reason for this is that the masses of the five-plet are only exactly degenerate at tree level but the custodial symmetry is not protected against loop effects [85]. Therefore, the number of mass parameters but also of rotation angles is extended at the loop level: one needs three instead of two angles to diagonalise the loop-corrected CPeven mass matrix, and also the CP-odd and charged Higgs mass matrix no longer share the same angle. Therefore, an MS renormalisation of the scalar sector is the natural option to check the impact of higher-order corrections to the masses Table 6 The scalar masses at tree and loop level for the parameter point λ 2 = λ 3 = λ 4 = 0, m 5 = 1 TeV and s H = 0. 75 and angles. We give in Table 6  We see in these numbers that not only a mass splitting between the components of the five-and seven-plets is induced at the one-loop level, but also that the loop corrections to the SM-like Higgs scalar can be huge. One can understand these large loop corrections for the chosen parameter point to some extent analytically: the one-loop corrections to the (1, 1)-element of the CP-even mass matrix are given in the effective potential in the limit m 5 v by Thus, for large values of m 5 and/or s H one can expect huge corrections to the mass. Note, there are additional loop corrections to the off-diagonal elements of the scalar mass matrix which can have a significant impact on the masses. Therefore, one needs a full numerical calculation already at the one-loop level to obtain an accurate number for the SM-like Higgs mass. Before we further investigate the loop corrections, we want to comment briefly on the choice for the renormalisation scale Q. In the SM, but also in other models like 2HDMs, it is suitable to set Q = m t to give a good convergence and ensure that there are no large logarithmic contributions from top loops. However, in the GM model the dominant loop corrections involve often scalar fields with masses ∝ m 5 . Therefore, the overall size of the loop corrections is usually smaller for Q = m 5 as one can see in Fig. 10.
We check now the Higgs mass in the (m 5 , s H ) plane proposed in [84] always using Q = m 5 . The other parameters are fixed in this plane to  The light Higgs mass at the one-and two-loop level is shown in Fig. 11. As expected, we see that the two-loop corrections are large for large s H and m 5 . In order to further demonstrate this, we show in Fig. 11 also explicitly the size of the oneand two-loop corrections for all three CP-even scalars. We see that in the upper right corner in the (s H , m 5 ) plane the two-loop corrections are much larger than the one-loop ones and the Higgs can even become tachyonic. For m 5 = 1 TeV, this already happens at s H > 0.5, while for m 5 = 1.5 TeV the upper limit of s H is as low as 0.25. For large m 5 this limit is much stronger than the one from perturbative unitarity of V V → V V scattering amplitudes which gives s H < 667 GeV m 5 [86]. Thus, even if it might still be possible to obtain the correct Higgs mass at two-loop level by adjusting the other input parameters or by absorbing finite corrections into counter-terms, the results in this parameter region should be taken with a lot of care. Most likely, they are meaningless. However, also for the other parameter regions with a reasonable hierarchy of the one-and two-loop corrections, one would need large adjustments in the input parameters to compensate for these loop corrections. These changes would then reflect in the couplings and some decay widths of the 125-GeV scalar will deviate for large s H and/or m 5 clearly from the tree-level expectation. Finally, one can also see in Fig. 12 that the loop corrections for the other scalars are sizeable and can shift the masses easily by tens to hundreds of GeV.

Conclusions
In this paper we have presented several varied results relating to the calculation of two-loop corrections to the Higgs mass in general models. Chief among these are: Higgs mass from consistently solving the tadpole equations to include more general minimisation conditions, in particular allowing fermion masses to be directly dependent on the parameters (such as μ in the MSSM) with the expressions given in Appendix B. 3. We compared our results with those available for the Standard Model. In particular, this allowed a comparison within the same code of calculations in two different gauges, and we also found that the electroweak corrections are negligible, while those from momentum dependence are very small. 4. We showed that our new computation does indeed remove the instabilities (sharp peaks in the Higgs mass for certain parameter choices) in the previous approach for supersymmetric models; however, the reader should be aware that there are still some limitations when scalar masses in the loops become small compared to the renormalisation scale. 5. We explored the corrections to the mixing angle in the alignment limit in the Two Higgs Doublet Model using Fig. 12 The size of the one-(left) and two-loop (right) corrections in the (s H , m 5 ) plane for the second (first row) and third (second row) CP-even scalar the MS couplings as inputs, and found that provided the quartic couplings are chosen to be small, the loop corrections are safely under control. 6. We explored the 2HDM and Georgi-Machacek models with masses as physical inputs and using tree-level relations to obtain MS couplings, as commonly done in the literature. We find that in most regions of the parameter space these lead to large quartic couplings, which rapidly lead to loss of control of the loop corrections. Perhaps surprisingly, this often occurs well before the couplings reach naive perturbativity bounds.
All of the shown results are available to the community with SARAH version 4.12.0, and we hope that this contributes to an efficient and more precise study of many extensions of the SM; this should open the avenue to much future work. It would be particularly interesting to explore more carefully the relationship between on-shell and MS calculations in non-supersymmetric models, to better understand how the divergent behaviour of the masses that we observe for the MS scheme translates into differences in physical couplings -or even possibly ruling out certain parameter regions of models as unphysical.
For the technical program of generic Higgs mass computations, it would be very interesting to compute the corrections to the electroweak VEV and top Yukawa coupling to the same precision that we can achieve for the Higgs mass from MS/DR inputs. It would also be interesting to complete the set of contributions with those stemming from electroweak couplings, even if we showed that these must be very small in the case of the Standard Model. from the ANR grant "HiggsAutomator" (ANR-15-CE31-0002). JB was supported by a scholarship from the Fondation CFM. FS is supported by ERC Recognition Award ERC-RA-0008 of the Helmholtz Association.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Loop functions
A.1Ṽ (x, 0, z, u) One of the key functions of the basis set is V (x, y, z, u). This is defined as It is singular as y → 0, so we define the regularised version (defined with one fewer explicit index to [38]): On the other hand, we require a slightly different regularised function: For the case x = 0, we can simply extract the result at vanishing external momentum: Then constructingṼ gives On the other hand, for x → 0 we cannot take the above limit. In principle we could start again from the expressions for V (x, y, z) given in the appendix of [38] and take the smooth limit x → 0. Instead, here we provide two direct derivations of a compact expression forṼ (0, z, u), with both general external momentum and then for our generalised effective potential limit. The starting point for the first derivation is the set of differential equations given in [39], in this case where the coefficients of the loop functions are themselves functions of s, x, y, z, u. However, here we encounter the problem that several of these coefficients are actually singular as y → 0 -so we cannot simply evaluate the right-hand side of the equation to determine V (x, 0, z, u)! However, we can obtain such a closed-form expression by using the ansatz y + Δ l log y + Δ 0 + · · · , (A.7) and substituting it into the above differential equation, to find f and f 1 : f log y + f 1 + · · · = (− 1 y + k 0 UU ) f 0 (s; x, z, u) + f (s; x, z, u)A(y) The The coefficients defined in the above are If we then make our generalised effective potential expansion, we find f 1 (s; 0, z, u) = − P SS (z, u) log(−s) s − log(−s) 2(u − z) 3 u 2 − z 2 − 2uz log u z + 1 4(u − z) 4 5(u + z) 3 + 8uz I (u, z, 0) + 2z log z(2u 2 − 11uz + z 2 ) + 2u log u(u 2 − 11uz + 2z 2 ) + 4uz(u + z) log u log z . (A.13) We do not need the limit when u = z = 0 because in that case we have couplings λ GGG . However, for z = 0 or u = 0 we do see that there is a smooth limit of the above.
Let us define If we now substitute in the standard expressions for I (z, u, 0) then we can simplify the above to We can also write it in a shorter but less symmetric form, . (A.23) We can then integrate this expression. For the case x → 0 we can simplify a little:

B Consistent solution of the tadpole equations with shifts to fermion masses
Here we give the two-loop shifts to the tadpoles and selfenergies due to shifts in fermion masses when we solve the tadpole equations consistently. We denote the undiagonalised fermion mass matrix as m I J . The mass-squared matrix is defined [58] as To illustrate this, consider the MSSM, where the tadpole equations read Solving for |μ| 2 we have This in turn will lead to a shift in the neutralino and chargino masses, which lead to a shift to the two-loop tadpoles.