NLO electroweak corrections in general scalar singlet models

If no new physics signals are found, in the coming years, at the Large Hadron Collider Run-2, an increase in precision of the Higgs couplings measurements will shift the dicussion to the effects of higher order corrections. In Beyond the Standard Model (BSM) theories this may become the only tool to probe new physics. Extensions of the Standard Model (SM) with several scalar singlets may address several of its problems, namely to explain dark matter, the matter-antimatter asymmetry, or to improve the stability of the SM up to the Planck scale. In this work we propose a general framework to calculate one loop-corrections in BSM models with an arbitrary number of scalar singlets. We then apply our method to a real and to a complex scalar singlet models. We assess the importance of the one-loop radiative corrections first by computing them for a tree level mixing sum constraint, and then for the main Higgs production process $gg \to H$. We conclude that, for the currently allowed parameter space of these models, the corrections can be at most a few percent. Notably, a non-zero correction can survive when dark matter is present, in the SM-like limit of the Higgs couplings to other SM particles.


Introduction
With the start of Run 2, CERN's Large Hadron Collider (LHC) has entered the stage of precision measurements of the Higgs couplings to the Standard Model (SM) particles. Even though the particle physics community is focused on the search for direct signals of beyond the SM (BSM) physics, it may happen that no such signal is detected during Run 2. If this is the case, we need to take advantage of the precise determination of the relevant Higgs couplings to understand if any new physics contributions can be hidden behind those measurements. Scalar extensions of the SM have, in most cases, a decoupling limit where, if the new scalar states are heavy enough, the model can only be probed via radiative corrections. In fact, if we are faced with a situation where no direct hint of new physics is found, manifestations of BSM physics can only appear through deviations in the measured Higgs couplings.
In this work we will focus on extensions of the SM where an arbitrary number of singlets is added to the SM field content. These are the simplest extensions of the scalar sector that introduce a dark matter candidate [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. These models also allow for a strong first-order phase transition during the era of Electroweak Symmetry Breaking (EWSB) if the extension comprises at least two singlets [16][17][18][19][20]. Hence, at least two of the outstanding problems of the SM can be solved within the framework of these models, namely a candidate for dark matter and a solution to the matter antimatter asymmetry via electroweak baryogenesis. It should be noted that extensions with only one dark scalar singlet are basically excluded by the latest LUX results [21] when combined with the requirement that the dark matter relic density of the model matches the one obtained from the Cosmic Microwave background data. We have verified that this is not the case when at least two new singlets are added to the SM, one of them being dark and the other mixing with the SM-like Higgs.
If the LHC indeed does not find strong signs of new physics, such as new particle states, the scale for such new physics may be as large as the GUT or the Planck scales. This energy is unattainable by any current or planned collider experiments so we may have to work in a framework that is a good description of the fundamental interactions up to some high energy scale. Thus, any effective description that improves theoretical problems of the SM is an interesting candidate. In a previous work we have shown that the complex singlet extension of the SM also improves the stability of the SM. In fact, the presence of a heavier scalar state, which has to be heavier than about 140 GeV, can stabilise the SM up to the Planck scale [22].
In this article we focus on the issue of determining electroweak (EW) radiative corrections in general scalar SM extensions, with emphasis on the scalar singlet models framework. Our main goal is to find a general set of expressions that allows us to obtain next to leading order (NLO) electroweak corrections to the parameters of a model with any number of scalar singlet fields. We will go beyond the effective potential approach recently studied in [23], which is valid only when the new degrees of freedom are heavy. Thus, though we formulate our results to connect to that limit, they are valid for any external momentum scale (contrarily to the effective potential approximation which is valid for small external momenta). In our framework we obtain a set of conditions consistently truncated in an expansion in powers of which, once a number of consistent independent input parameters are chosen, deliver the NLO EW corrections for the remaining parameters. A special attention is given to the treatment of tadpoles and propagators and we provide a generic strategy to easily transform between different schemes. In connection with the effective potential approximation we also discuss, on general grounds, the issue of infrared divergences. We then apply our method to the real scalar singlet extension (RxSM) and to the complex scalar singlet extension (CxSM) of the SM. However, we note that the method is ready to be applied to SM extensions with an arbitrary number of singlet fields and that many of our formulas are also useful for other scalar extensions of the SM. In particular, our approach is especially suited for the automation of the computation of higher order corrections in general purpose numerical tools to scan the parameter space of scalar extensions of the SM [24,25].
Higher order corrections to singlet extensions of the SM have been performed in [26,27]. The corrections to the SM-like Higgs coupling to fermions and gauge bosons was shown to be of the order of 1 % [26]. Furthermore the corrections were maximal in the decoupling limit where the model becomes indistinguishable from the SM. Electroweak corrections to the decay H → hh were performed in [27]. With the main theoretical and experimental constraints taken into account, it was shown that corrections to the triple scalar vertex (Hhh) are of the order of a few percent. Calculations of higher order corrections in the complex singlet extension of the SM are still not available. With this work we will not only present a set of equations to renormalise the parameters of the theory at one loop but we will then also use them to calculate the electroweak correction to Higgs production via gluon fusion. This last calculation is performed near the decoupling limit with the main purpose to understand the contributions of the triple scalar couplings of the various scalars running in the loops at NLO. Clearly, with all the SM-like Higgs coupling close to the SM ones, the only large effects in the radiative corrections would have to come from such scalar-scalar interactions. The numerical analysis in our examples will be performed for three particular cases: the broken RxSM, with a new Higgs boson mixing with the SM-like one, and the broken and symmetric CxSM with, respectively, three mixing Higgs bosons, and two mixing Higgs bosons and a dark matter scalar. We will find that, consistently with earlier calculations for the NLO corrections to the decays, the corrections are very small, of the order of a few percent, also for production. Nevertheless, we will find that the presence of a dark matter particle can enhance the corrections, even very close to the SM-like limit, compared with the other models (though still in the few percent order).
The smallness of the electroweak corrections in the real singlet models calls for prudence in the claims of measurable differences relative to SM Higgs couplings. The interference effects for this kind of BSM scenarios was first addressed in [28] for the real singlet model, showing that interference effects to gg → h * , H ( * ) → ZZ → 4l can be important away from the non-SM scalar (H) peak region. Although the interference effects can be of up to order O(1) for the integrated cross sections for the 8 TeV LHC [29], judicious kinematical cuts can reduce the interference effects to O(10%). Interference effects at NLO QCD were discussed in [30] for the process gg → h * , H ( * ) → hh. It was shown that the double Higgs invariant mass can increase by up to 20% or decrease by up to 30% depending of the heavier Higgs mass. More importantly, interference effects can significantly distort the kinematic distribution around the resonant peak of the heavy Higgs. Recently the effects of higher order operators in the real singlet model [31] again showed that large cancellations can occur due to interference effects between the two sectors. In conclusion, if a significant deviation is found in Higgs couplings, the radiative corrections have to be combined with the interference for a proper interpretation of the results.
The structure of the paper is as follows. In the first two sections we start by defining our strategy. We present the Lagrangians and fix the notation in Sect. 2 and then, in Sect. 3, we obtain our master linear system that, given a choice of input parameters, provides as output the remaining renormalised parameters at NLO EW. The issue of infrared divergences in connection with the effective potential approximation is discussed in Sects. 3.1 and 3.2. In Sect. 4 we apply the procedure first to a general class of scalar singlet extensions of the SM and then specialise to the RxSM and to the CxSM, Sect. 4.4, for which we provide a numerical analysis in Sect 4.5. Our conclusions are summed up in Sect. 5 and several useful formulae/derivations are provided in the appendices.

Definitions and notation
To define a general four dimensional gauged Quantum Field Theory (QFT) Lagrangian we use the notation of [32] with a few adaptations [22,23]. We assume a decomposition of a general renormalisable Lagrangian where the gauge basis fields are such that: i) all scalar field multiplets are decomposed as N 0 canonically normalised real scalar fields, Φ i (i = 1, . . . , N 0 ), ii) all fermion multiplets are decomposed as a set of N 1/2 two-component Weyl fermions, Ψ I (I = 1, . . . , N 1/2 ) and iii) there are N 1 gauge bosons in the adjoint representation of the gauge group, i.e. A µ a (a = 1, . . . , N 1 ). We adopt the Einstein convention where repeated indices which are one up (superscript) and one down (subscript) are summed over. If the repeated indices are both down or both up they are not summed over. All (non-spacetime) latin indices are assumed to be in Euclidean space -they are lowered and raised with the identity matrix. The gauge basis interaction Lagrangian (i.e. suppressing kinetic terms) is then composed of the following terms: where the ghost fields are represented by ω a and c.c. denotes complex conjugation. We call this the L-basis following the nomenclature in [22,23], where the pure scalar, fermionic and gauge interaction coupling tensors are denoted, respectively by {L ... , Y ... , G ... } with . . . replaced by suitable sets of indices. Note that for a simple gauge group G abc is given by G abc = g f abc with g the gauge coupling constant and f abc the structure constants of the gauge group. For a direct product group we can still encode all information in G abc by requiring a block structure. This can be represented using sub-ranges for the indices a 1 = 1, . . . , n 1 , a 2 = n 1 +1, . . . , n 3 , etc... if components that have indices not all in the same sub-range are zero. More concretely we would have , etc... and, for example, G a 1 b 2 c 2 = 0.
A second form of the interaction terms is obtained after shifting the scalar fields by a general classical field configuration such that Φ i (x) = v i + φ i (x): Here we have introduced the notation Λ ... (T ) for the interaction coupling tensors containing a field of type T and scalar fields without derivatives where the type T runs over the three possible types of fields {S, F, G} (scalar, fermionic and gauge respectively) and . . . represents a set of indices. These couplings appear in the calculation of the effective potential which, at one loop in the Landau gauge, can be organised as a sum over field types -to be discussed in Sect. 3.1. For fermions the natural objects appearing in the Eq. (2.2) after the shift are the mass matrix M IJ and the Yukawa couplings Y IJk . However, we can also define a mass squared matrix and (effective) cubic and quartic couplings The latter can also be found in appendix A together with all other v i dependent parameters, Λ ...
(T ) -see also [23]. Finally, we can rotate all fields to their mass eigen-basis (named the λ-basis) through orthogonal or unitary transformations, for bosons and fermions respectively. Using the transformations (A.6) then we have Note that all couplings in the λ-basis, Eq. (2.3), whose indices are rotated according to the transformations induced by (A. 6), are now in lower case. For completeness, we provide in appendix A the relations between the various bases including the rotation matrices to obtain the mass eigenstates. We note that the latter can be represented collectively, for the field type T = {S, G, F }, by U (T ) (unitary or orthogonal) with the defining relation that all mass matrices are diagonalised: and reserve indices from the middle of the alphabet (i, j, k, . . .) for scalar field indices 1 . Note that whenever we use a matrix notation without explicit indices, Λ (T ) is assumed to represent Λ ab (T ) , i.e. the mass squared matrix. In this article we will need to compute the one-loop radiative corrections to the scalar mass eigenstates. These are determined by the poles of the radiatively corrected propagator, G ij , between scalar states i and j. It is well known that the Dyson re-summed inverse propagator G −1 ij can be expressed as [33] i where p µ is the external 4-momentum 2 . Furthermore where T j and Π ij , are the one-particle irreducible tadpole (1-point) and self-energy (2-point) functions. The tadpole term can in practice be set to zero by assuming an expansion of the theory around a minimum of the effective potential order by order in perturbation theory [34] if one works in Landau gauge (which we assume in this work). Finally we have defined ∆Σ ij p 2 ≡ Σ ij p 2 − Σ ij (0). The physical pole state x (labelling the N 0 physical scalar states) is defined to have an eigenvalue p 2 = M 2 x and an eigenvector E i x such that the pole conditions are obeyed: Here, we have suppressed the scalar indices inside the determinant so a matrix notation is used -see Eq. (2.5). In general, for unstable particles, the eigenvalue can have an imaginary part which relates to the width of the particle so one defines where m x is the physical mass and Γ x is the width of the particle. Finally, we will also need to extract the wave function renormalisation factor, Z x , from the inverse propagator. Projecting the inverse propagator along the normalised eigen-vector this is obtained from Thus, noting that the projection vectors do not depend on p 2 we obtain Z x from (2.9) 1 The latin indices from the beginning of the Latin alphabet will also be used later for the gauge field indices. However there is no danger of confusion because whenever we expand the expressions in the field type T all types of indices appear explicitly (scalar, fermionic and gauge indices). 2 We are using the metric signature convention (+ − −−).

General loop expansion and perturbative strategy
In this section we aim to organise the perturbative strategy to obtain the one loop corrections to observables. To do so we construct the loop expansions by consistently truncating in powers of . It will be convenient to use as expansion parameter the quantity ε ≡ /(4π) 2 . We start with some general considerations, assuming that a renormalisation scheme has been fixed and that, in principle, we have calculated a loop expansion for any observable as a function of a given set of input parameters for the theory. The number of such input parameters is equal to the number of running couplings appearing in the tree level Lagrangian. However, one may want to use a different set of input parameters, for example physical observables that are measured experimentally, to replace such Lagrangian parameters. With that choice then the Lagrangian parameters will be functions of such (observable) inputs. A more general and convenient approach may be to choose a mixture of observables and Lagrangian parameters, or even other theoretical parameters such as VEVs, as input.
To avoid choosing, a priori, a particular set of inputs, we adopt a perturbative strategy to obtain radiatively corrected relations among the various parameters. We formally loop expand all parameters to allow for a free choice of inputs. To specify which parameters are input we can then simply set their correction terms to zero. An advantage of this procedure is that we obtain a linear system of constraints that we can analyse to decide what are the available choices of sets of input parameters. 3 This flexibility is particularly well suited for the automation of the calculation of higher order corrections in general purpose tools to scan the parameter space of general scalar extensions of the SM [24,25].
To organise the system, let us denote collectively the set of Lagrangian parameters, VEVs and observables that we want to expand by Q A . For example A could run over the elements of the list C λ , v i , m 2 x , . . . with C λ representing the set of all Lagrangian parameters with 4 λ = 1, . . . , N λ ; v i the VEVs and m 2 x the physical masses. We assume that all quantities (for concreteness, renormalised in the MS-scheme) are loop expanded as A is the tree level quantity). Furthermore we assume that there is a set of constraints. Those can be, for example, the definition of observables (such as pole masses) or theoretical conditions such as the vacuum (or minimum) conditions. In general, the system of constraints is represented by where C Γ denotes a constraint, with the index Γ running over all the available constraints, and we suppress the index of Q A in arguments of functions for a lighter notation. Typically the constraints are also loop expanded as: Expanding up to linear order we obtain Equating order by order we obtain (for the first two orders) so once we solve the tree level constraints for Q (0) A we obtain a non-homogeneous linear system of constraints for the one-loop corrections Q A . For concreteness we now turn to the problem we want to solve. First, we wish to compute the one-loop corrections to the constraints that define the vacuum state and the physical masses in the scalar sector. With this we will obtain relations among the parameters of the scalar sector of the theory. We assume that the gauge and fermion sector couplings are input parameters that are known at the renormalisation scale µ. As already stated we formally allow all scalar couplings, all VEVs and all parameters defining the mass eigenstates to be loop expanded. Later we will decide, on a model by model basis, which parameters of the scalar sector are input. The set of constraints to impose are which are, respectively, the minimum conditions (tadpole equations) and the pole equations, which define the physical mass eigenstates -see also Eq. (2.7). The parameters of the theory that can, in principle, be expanded are where we have noted that the width is always zero at zeroth order so the Leading Order (LO) results appears only at first order. We also note that where V (n) is the n-loop effective potential and the self energy series only starts at first order. Using the general expansion, Eq. (3.5) and inserting the expansions, Eqs. (3.7) and (3.8), and assuming, without loss of generality, that we are in a field basis such that the tree level scalar eigen-states are aligned along each field direction, i.e. E (0)j x = δ j x , we obtain the tree level conditions and a linear system for the one-loop corrections (3.10) Here we use the notation [. . .] tree to denote a quantity evaluated with tree level arguments and ∆Σ (1) ij is evaluated at s = m (0)2 x . Taking the real and imaginary parts of the second constraint, and noting that, in the tree level basis, the tree level mass squared matrix is diagonal, we obtain the final result To this, we can add the conditions that the eigenstates are normalised which, up to one-loop order, translates to Now that we wrote the general system, Eq. (3.11), we discuss some possible choices of inputs. If, for example, one chooses to take as input parameters the set of running Lagrangian parameters C λ , we set C (i>0) λ ≡ 0. Then the one-loop shifts of all other quantities are computed using the system in Eq. (3.11). Within that choice, using the one-loop tadpole equation in the first line of Eq. (3.11) we obtain directly that, for consistency, for states that are massless at tree level, the first derivative of the one-loop effective potential evaluated at the tree level couplings must be zero and the corresponding VEV shift remains undetermined. Denoting the field space directions associated with the massless states by the sub-indices 1 , 2 , . . . and the ones associated with massive states with barred indices i,j . . ., then the remaining VEVs are obtained from (3.13) All that is left in this case is to solve the pole equations, Eq. (3.11) with C Other possible choices consist of perturbative inversions of the one-loop relations in Eq. (3.11). In our examples in Sect. 4 we will choose some input parameters to be physical quantities (such as masses) and others such as the Higgs VEV, mixing matrix elements and a few Lagrangian parameters. This is convenient, for example, to fix the one-loop Higgs mass to the experimental value of 125 GeV and the Higgs VEV to 246 GeV. Then the one-loop shifts of the remaining parameters are computed, ensuring that the relations among all parameters are correct to one-loop order.
A particularly choice that is useful, is one that decouples directly the corrections to the mass eigenstates states. Taking the anti-symmetric part of the second condition in Eq. (3.11) for i = x we obtain (3.14) So, assuming that the system allows the choice E we get the solution for all the corrections to the mass eigen-state expansions.
Finally, another set of quantities that we will need are the wave function renormalisation factors. Expanding Eq. (2.9) perturbatively we find

Coleman-Weinberg potential and self energies
The quantities we will need to evaluate in Eq. (3.11) are: the Coleman-Weinberg potential, its first and second derivatives and the variation in the self-energy functions. In the pole conditions one could, equivalently, simply compute the full self-energy functions. However, in Eq. (2.7) we have written the result in terms of the effective potential to separate out the p 2 -independent part and the p 2 -dependence. This is useful to connect to the p 2 → 0 approximation, which can be used if the dominant contributions to the radiative corrections are from heavy particles in the loops [23]. In that limit the effective potential encodes all the necessary information.
The general one-loop effective potential in the MS scheme is given by the Coleman-Weinberg potential. Recently [23] we have analysed the Coleman-Weinberg potential and obtained a closed form master formula for its one-loop derivatives with any number of external scalar field legs for a general theory as described in Sect. 2. Here we only review the two expressions that we need, i.e. the first and the second derivatives of the effective potential, respectively, and Here log (x) ≡ log(x/µ 2 ), with µ being the renormalisation scale; s T = 0, 1 2 or 1 is the spin of the field of type T and the λ ...
(T ) couplings have been defined in Eq. (2.3); S {ij} denotes symmetrisation of the indices; and Observe that the latin indices a, b, . . . in the λ ... (T ) tensors are to be replaced by scalar, fermionic or vector indices respectively according to T . The constant k T depends on the renormalisation scheme (for MS it is 3/2 for scalars and fermions and 5/6 for vector bosons).
Regarding the self-energies, they have been computed in [34]. Here we present the variation that we need, which is (s ≡ p 2 ) The various loop function variations can be obtained directly from the results in [34,35] and are provided in appendix B.

Infrared behaviour
It is well known [36,37] that the derivatives of second and higher orders of the (oneloop) Coleman-Weinberg potential can contain infrared divergences originating from the massless states running in the loops. However this is not a problem if all the p 2 dependent contributions are included because all such infrared divergences must cancel out.
In this section we verify this general cancellation explicitly and write the final result in a manifestly regular form suitable for numerical evaluation. First let us introduce an infrared regulator mass squared scale, . The second derivatives of the effective potential can be split as where the second term contains the contributions with internal sums over the indices a, b corresponding to two internal massless states -see Eq. (3.17). We recall that, in our notation, indices corresponding to massless state components are denoted by 1 , 2 , . . . when the type of the field, T , is not specified. Whenever T is specified, we use the indices With this notation, the IR-divergent piece is where, on the second line, we have series expanded in the cutoff, , and kept only the divergent and constant terms, respectively denoted by IR,div , we have explicitly expanded over spins. Moving on to the self-energies, we define a similar split Finally, using the fact that, in the mass-squared eigenbasis, m IJ only has non-zero elements between states I, J with the same mass, using Eq.

(3.24)
Therefore the divergent pieces cancel precisely. The final, explicitly finite, result for the full one loop contributions appearing in the pole equations is then IR,finite + ∆Σ (1) ij,finite (s) + ∆Σ (1) ij,IR,finite (s) . 4 Application to general scalar singlet extensions of the SM In this section we will apply the results of the previous section to the most General scalar singlet extension of the SM (GxSM) and then we specialise to a real (RxSM) and a complex (CxSM) singlet extensions. For simplicity we work in the MS scheme, so we replace the k T by their numerical values.

Definition of the GxSM
The most general scalar singlet extension of the SM is obtained by adding to the Lagrangian a set S k (k = 1, . . . , N S ) of real scalar hypercharge zero singlet fields with a general renormalisable scalar potential. The Lagrangian density of the model for the interaction terms is then where ∆(S) and V (S) are polynomials that are, respectively, up to quadratic and quartic in the fields S k (without constant terms) and H is the SM Higgs doublet. In this framework the full scalar potential is then One or more of the singlet fields can mix with the Higgs boson provided that ∂ k ∆(v i ) = 0 for at least one value of k at the electroweak symmetry breaking vacuum with a choice of VEVs v, v k such that Here h is the SM Higgs field fluctuation, G 0 , G + are the Goldstones and s k are the singlet field fluctuations around the vacuum. The new scalar singlet fields, S k , do not couple directly to other SM fields. As a consequence, the tree level coupling of the scalar mass eigenstates that are a mixture of singlet field fluctuations, s k , with the Higgs boson fluctuation, h, to the other SM particles is simply scaled by a mixing factor (compared with the Higgs couplings in the SM). We note that, at tree level, the Higgs field fluctuation is decomposed in terms of scalar mass eigenstates as (see Eq. (A.6)) where we have ordered the set of scalar fields after the VEVs shift as Here κ j is the scaling factor to apply to the SM coupling of an SM-like Higgs of the same mass as the state R j , to obtain the coupling of that state in the GxSM. Due to the orthogonality of the mixing matrix we have that, at tree level, which means that the SM-like coupling is shared among the Higgs like states [38]. As a consequence the one-loop radiative corrections to the scalar mass eigenstates will contain some SM-like contributions suitably suppressed by the dilution factors κ j and also contributions exclusively due to the new scalar sector.

NLO parameter shifts
In this section we compute the NLO shifts of the parameters in the GxSM.

One-loop tadpoles
The contributions to the tadpole conditions, Eqs. (3.16) all contain a coupling factor λ (T )aai . For fermions and vector bosons we know that, for the GxSM, the couplings to the massive scalars are simply scaled by a κ j factor and the couplings to Goldstone bosons are precisely the same as in the SM. Thus, noting that from now on we no longer deal with space-time indices, if we use Greek indices α = 1, . . . , 4 to denote the four scalar degrees of freedom in the SM and define a dilution tensor D α i (n s is the number of non-Goldstone real scalars) where we have denoted the SM couplings on the right hand side with the superscript SM. Then one can check that (see also appendix D) where, in the last line, we have evaluated the result keeping only the dominant top quark contribution in the fermion sector and the electroweak vector boson contributions -see appendix D.

One-loop poles and wave function renormalisations
From the pole equations, Eq. (3.11) and using Eq. (4.7) one can show that where, again, the barred indices run only over eigenstates with a non-zero mass. In practice we will be interested in the components i, x such that s = m 2 x so we have The other quantity that we will use are the one-loop wave function renormalisation factors for the massive states given by which, we can check, also involves ∂ s B s . The term that will be most relevant is x . With all these ingredients, we will specialise these formulas in Sect. 4.5 to obtain the parameter shifts in particular scalar singlet models.

Corrections to mixing sums
In the GxSM, generically, there is a mixing of the SM Higgs field fluctuation with singlet fields. This typically results, at tree level, in a block n by n mixing matrix with n the number of non-dark scalar mass eigenstates with the tree level sum rule, Eq. (4.5), for the suppression factors of each scalar x to other SM particles. If we denote the tree level suppression factor by κ (0) x , at one loop using Eq. (3.7), the one loop mass eigenstates in the gauge basis are (with j, x running over the mass eigenstates) where now we denote by κ x the one loop corrected mixing factor. From this, and using Eq. (3.14), we obtain the order ε correction to this sum, which is given by This result is in fact independent of our choice of one-loop input parameters within the choice we made to normalise the mass eigenstates such that E (1) xx = 0. Thus, though we will evaluate Eq. (4.15) in a specific scheme, the result is fixed within our class of schemes.
In our numerical results we will be interested in assessing the importance of the one loop corrections to the parameters of the theory. The quantity in Eq. (4.15) is a good one to test the importance of these corrections as it is a shift of a tree level value that is 1, and that does not depend on other choices within our class of schemes. Finally, this quantity is also of interest because it will contribute to the NLO sum rule that is expected to exist among the effective couplings of the mixing Higgs bosons to SM particles such as to preserve unitarity. A complete computation of such effective couplings is, however, beyond the scope of this study.

NLO gluon fusion cross section
In Sect. 4.5, we will evaluate the NLO electroweak corrections to the SM-like Higgs production cross-section in the gluon fusion channel that are due to the new scalar sector couplings. The current collider data already sets the suppression factor for the SM-like Higgs to be very close to unity and the suppression factors of the other new Higgs bosons to be very close to zero, i.e. κ 2 h ∼ 1 and κ 2 i =h 1. Therefore we will focus on this limit by systematically dropping terms that are suppressed by κ i =h or higher powers. Furthermore, using the standard assumption of factorisation of the QCD higher order corrections (see for example [39]), we focus only on the NLO electroweak corrections. The NLO amplitude for gluon fusion in the GxSM is of the form: (4.17) Squaring the amplitude and replacing κ 2 h = 1 (we are keeping here terms linear in κ h because the sign can be important depending on the convention) we obtain where the σ where (4.22) Here we denote SM limit quantities with the superscript SM so that we can separate out δ SM , which are the NLO correction as in the SM, and the new contributions due to the extended scalar sector δ GxSM . The factor C hhf was computed 5 in [40,41], so we use the value C hhf = 0.0066, which is independent of the centre of mass energy. Note that the κ h dependence is important because both signs are allowed in the conventions used in our scans.
We are interested in observing if there are scenarios where the new scalar contributions can correct considerably the cross-section. It is well known [39] that the SM corrections are small, δ SM 5%. Our aim is to compute the new factor δ GxSM . Finally, note that the SM-like parts in δZ h − δZ SM h , for T = S cancel out exactly in the difference, in the limit we are using, so we only have to evaluate the difference over the scalar contributions without Goldstones δZ h − δZ SM h scalars using Eq. (4.9). One can check that, at one loop, δZ h = (1 − Z −1 h ), which we can obtain from Eq. (4.12). In particular we get that in the SM limit the corresponding contribution is in agreement with [40]. Another important issue relates to the crossing of thresholds. In particular when a threshold for the Higgs boson to decay to a pair of a lighter bosons is crossed the wave function renormalisation contains a singularity at the threshold if evaluated with a real pole mass [42,43]. This problem can be cured by working with the full complex pole mass where the width of the Higgs boson in included in the imaginary part of p 2 [42]. In our calculations -see e.g. Eq.(3.15) -we expand around the real tree level masses. Thus, to avoid these unphysically large deviations, we will ignore scenarios close to these thresholds, within a 5 GeV mass window. Away from these thresholds this approximation is known to work well -see for example Figs. 6 and 7 of [43] for the EW corrections to gluon fusion around the W W threshold in the SM. The near threshold cases have the potential to produce extra enhancements. This limit is beyond the scope of our study and will be left for future work.
In our numerical analysis we use the following values for the relevant SM parameters (consistently with the code sHDECAY used to generate the tree level samples [44]): Note also that, in our normalisation, the SM-like Higgs triple coupling is given by

Particular models
In Sect. 4.5 we analyse samples for two benchmark models to illustrate the size of the EW NLO corrections. Here we provide a brief summary of the models.

The real singlet model (RxSM)
The potential for the model is Here the (real) couplings of the theory are m, λ, λ HS , m S and λ S and S is a real singlet field with a Z 2 symmetry (S → −S). In this model S is a dark matter candidate if v S ≡ S = 0, or it is a new scalar mixing with the Higgs if v S = 0. We will focus on the latter because the former seems to be very close to being ruled out except in the region around m h 125 /2 and for very large dark matter masses (see for instance [45][46][47]). The model has five independent input parameters which we choose to be {α, m 1 , m 2 , v, v S } both at tree level and at one loop. Here m 1 < m 2 are, respectively, the masses of the scalar eigenstates h 1 and h 2 decomposed as In this way, we guarantee that there is an SM-like Higgs boson in our one-loop corrected samples that is compatible with the observed Higgs boson. In sum, the one loop shifted parameters are then all the Lagrangian parameters, {m 2 , λ, λ HS , m 2 S , λ S }, and the eigenstate h i = h 125 . The VEVs were also chosen to be input, so their shift is set to zero.

The complex singlet model (CxSM)
The potential for the model is Here S = (S + iA)/ √ 2 with S, A real fields. This model has a dark phase if v A ≡ A = 0 and v S = 0, with A the dark matter candidate and two Higgs bosons mixing. In this phase the A → −A symmetry is preserved. In the broken phase v A = 0 and we have three Higgs bosons which may be visible at colliders. We will investigate both phases of this model.

Dark phase:
In the dark phase, v A = 0, the set of inputs is similar to the RxSM. We choose {α, v, v S , a 1 , m 1 , m 2 , m D } where m 1 < m 2 are, respectively, the masses of the scalar eigenstates h 1 and h 2 , and m D is the mass of the dark matter particle. The mass eigenstates are now decomposed as where we note that, again, we can choose to shift only the state that is not the SM like Higgs boson. After analysing the one loop system, Eq. (3.11), one concludes that the set of shifted parameters can be given again by all the Lagrangian parameters except a 1 , i.e. {m 2 , λ, δ 2 , b 2 , d 2 , b 1 }, and the eigenstate h i = h 125 . The VEVs were also chosen to be input, so their shift is zero.
Broken phase: Regarding the broken phase of the model, when v A = 0, the input parameters are now chosen to be {α 1 , α 2 , α 3 , v, v S , m 1 , m 3 }. The three mass eigenstates h 1 , h 2 , h 3 and their masses m 1 < m 2 < m 3 are decomposed as where α i ∈ [−π/2, π/2] and now we do not display explicitly the choice of one loop corrections to the mass eigenstates for brevity. In this case the one loop system, Eq. (3.11), forces us to introduce shifts for more than one mass eigenstate. Nevertheless, it still allows us to fix one of the mass eigenstate so, again, we choose to keep the SM-like Higgs mass and eigenstate as input. In addition we choose three corrections to the other mass eigenstates to be non-zero, such that we can solve Eq.
Both the RxSM and the CxSM were recently analysed in light of the LHC run-1 data, and compared with the NMSSM to determine if they could be distinguished from the latter in the Higgs sector [22,44]. We will use the samples that were generated in [44], which are compatible with all the latest theoretical and phenomenological constraints. For the case with dark matter we have also applied the latest dark matter direct detection bounds from the LUX experiment [21]. For further details on the constraints applied in the samples we refer the reader to [44].

Numerical results
In this section we use the specific models presented in the previous section to illustrate the importance of the NLO EW corrections for the the parameters of the theory and for the gluon fusion production cross-section of the observed 125 GeV SM-like Higgs boson. In Sect. 1 we have noted that, according to previous studies [26,27], the NLO corrections to decay widths in the RxSM model are small, typically in the order of a few percent. Here we go beyond the simplest scalar singlet extension to compare the various models and, at a phenomenological level, focus on the production rather than the decay.
In Fig. 2 we first present the corrections to the tree level mixing sum relation. We have observed, Eq. (4.15), that the correction to this tree level relation is, in fact, a good proxy to evaluate the importance of the NLO corrections to the input parameters because it does not depend on further choices specific to the scheme (other than the normalisation condition of the mass eigenstates). We start by discussing the RxSM model (top left panel). The simplicity of this model allows us to interpret the distribution of points in the scan more straightforwardly. It is also useful to interpret the distribution of points for the other models because it is a limiting case of those models. We separate the two scenarios where the SM-like Higgs is the lighter (grey) and the heavier (purple) of the two mixing Higgs bosons. We observe that the magnitude of the correction to the tree level value of 1 for the mixing sum, ranges approximately between +4% and −0.5%. The various features of the distribution of points are due to the combination of theoretical and phenomenological constraints. One can see that for small masses the mixing sum correction is small and rises sharply from ∼ m h 125 /2 GeV to ∼ 100 GeV as the decay channel h 125 → h 1 + h 1 closes. This is due to stronger constraints from negative searches at colliders which force the model to be more SM-like 6 so that the new scalar is more decoupled (hence it also contibutes less in the loops). For masses larger than ∼ 125 GeV, the collider constraints become progressively more restrictive up to the opening of the decay of the heavy Higgs to a pair of SM-like Higgs bosons (h 125 ) where there is again a sharp rise. The top boundary for large masses is due to the electroweak precision data constraints. For the lower boundary, we observe that negative corrections become possible at about ∼ 350 GeV, i.e. the threshold for decay to a pair of top quarks.
For the CxSM dark phase (right panel), we observe a distribution of points that is similar to the RxSM with the difference that larger magnitude negative values are allowed. The colour scale, which represents the dark matter mass, shows that this is possible in scenarios where the dark matter particle running in the loops is lighter than about ∼ 500 GeV. Observing the yellow points, corresponding to large dark matter masses, we recover the lower boundary observed for the RxSM (see yellow region). This is expected because the contributions from the dark matter loop propagators are suppressed for large masses.
Regarding the CxSM broken phase (bottom panels) the distribution of points is more complicated because we have three mixing Higgs bosons. Nevertheless, from the vertical axes, we see that the magnitude of the upper and lower ranges of the mixing sum correction is similar. In the two panels we indicate the scenario where the SM-like Higgs is the lightest, next to lightest and heaviest of the three mixing Higgs bosons, respectively with grey, purple and yellow points. In the left panel we have the smallest mass in the horizontal axis and in the right panel the largest mass. In the two scenarios of the left panel (yellow and purple points) we observe that the points distribute similarly to the RxSM. This is consistent with the observation, in the RxSM, that a light scalar in the loop corrections produces positive corrections. The only difference is that the yellow points in the peak region for masses larger than ∼ 62 GeV are more suppressed. This is consistent with the fact that the SM coupling is diluted over two small mixing factors for each of the two light scalars, which further suppresses each contribution. The yellow points stop at ∼ 118 GeV because we have applied a cut to avoid degenerate scenarios with a distance of 3.5 GeV between any two scalar states. This is why the purple layer (h 125 ≡ h 2 ) stops at ∼ 121 GeV in the left panel. In the right panel we can observe the scenario h 125 ≡ h 1 in the grey points. There we see that for scenarios where h 3 is heavier than about 350 GeV we can have a negative correction. This is consistent with the RxSM observation that negative corrections are possible in this scenario when the heavy scalar is above this mass. The main difference with the RxSM is that the lower boundary is not so sharply defined. This is simply due to loss in density in the scan for the CxSM, which has more parameters making it harder to collect large amounts of points. Finally, note that for the purple points, where the SMlike Higgs is the next to lightest, it is not possible to obtain a negative correction, which indicates that the positive contributions from one lighter Higgs are enough to push the correction to positive values. Also note that, contrarily to the RxSM plot, we cannot see any sharp rise at m h 3 ∼ 2m h 125 because it is always possible, for fixed m h 3 , to have decays involving a lighter non-SM like Higgs with a different mass. This extra free parameter in the scan dilutes such a boundary. Now we turn to the discussion of the NLO EW scalar contributions to the corrections to gluon fusion in the limit, already discussed in Sect. 4.3, where the mixing factor of the SM-like Higgs boson is close to the SM limit κ 2 h 125 → 1. In our numerical analysis we have evaluated, for each model, the quantity δ GxSM , in Eq. (4.22) specialised to the RxSM and CxSM models using the samples already discussed. In addition, we have selected points within 10%, 5% and 2% of the limiting case κ 2 h 125 → 1. Close to this limit our approximation is then reliable. Furthermore, applying this increasingly tighter constraint simulates the experimental situation where the SM-like Higgs boson couplings are measured to be SM-like with an increasingly higher accuracy. Understanding how large the new scalar corrections are allowed to be, when such tight experimental constraints become available, may then provide improved bounds on each model. Finally, in all plots, we have applied a 5 GeV mass window exclusion around thresholds for the opening of SM-like Higgs to scalars decays to avoid the singular behaviour of the wave function renormalisation discussed at the end of Sect. 4.3.
In Fig. 3 we can first observe, by inspecting the vertical axes for each model, that the NLO EW scalar contributions to the corrections for the two models without dark matter (top left and bottom panels) are always negative, whereas for the CxSM dark (top right) positive values are possible as we move away from the limit κ 2 h 125 → 1 (blue and red points). Focusing first on the RxSM (top left) we observe that the distribution of points away from the SM-limit and the various peaks and thresholds follow closely the same patterns already observed in Fig. 2. For masses of the new Higgs in the range ∼ 70 GeV to ∼ 100 GeV we observe that the corrections can deviate away from the SM by about −1.6%, −1.1% and −0.8% for the red, blue and green points respectively. For larger masses, above the SM-like Higgs mass, we observe well defined boundaries for the three layers and the correction can be at most about −1%, −0.5% and −0.2%, respectively for the red blue and green points. For the broken phase of the CxSM (bottom panels), the ranges of values for the corrections are similar. In the bottom left panel, with the lightest Higgs boson mass on the horizontal axis, we can see more clearly the exclusion mass window we have applied around m h 125 /2 to avoid dealing with the singular limit. In the bottom right panel we represent the same points as a function of the largest mass. By observing the two bottom panels we do not recover directly the light and the heavy region of the RxSM plot since the points spread down not only for small values of the masses. This is because we have more than one Higgs boson coexisting with the SM-like Higgs. For the bottom right panel, we can have one scenario where h 1 ≡ h 125 and another one where h 2 ≡ h 125 . In the latter scenario we have checked that we can have h 1 in the mass region around 100 GeV to contribute to the correction with a larger negative value in a way similar to the RxSM, whereas the other scenario is very similar to the RxSM. This mixture of scenarios and the fact that the CxSM has more free parameters explains the distribution of points and the absence of sharply defined boundaries in the scatter plots.
Finally, we discuss the dark phase of the CxSM which is the one that allows for larger CxSM -broken phase deviations even for the scenarios closest to the SM-like limit. In the top right panel of Fig. 3 we see that the green points can spread down to negative values as negative as the other two layers. Furthermore there are also points where the correction can be positive. We have checked that these deviation, even in the limit very close to the SM, are possible due to contributions from light to intermediate mass dark matter in the region m D 500 GeV similarly to Fig. 2 (top panel). For heavier dark matter scenarios we recover the RxSM distribution of points.

Conclusions
In this work we have developed a general framework for the calculation of NLO EW corrections in scalar extensions of the SM with focus on models with an arbitrary number of scalar singlet fields added to the SM field content. The derived set of equations can be applied for a wide class of choices of renormalisation scheme with the only restriction that the chosen input parameters consistently allow for a solution of the truncated linear system -see also Eq. (3.11) and related discussion. Our results go beyond the effective potential approach since they are valid for any external momentum in the inverse propagators. We also point out that our general procedure is well suited for the automation of such corrections within a general purpose numerical tool such as ScannerS [24,25]. We then applied our method to specific models: the real scalar singlet extension and the complex scalar singlet extension of the SM. First, in order to assess the importance of the NLO EW corrections to the parameters of the theories, we calculated the NLO correction to the mixing sum x κ 2 x , which is one at tree-level. Other than the normalisation condition for the one-loop mass eigenstates this quantity is scheme independent. Thus it provides a good measure of the general trend for the magnitude of the one-loop corrections (in this case to a tree level mixing sum relation). The correction we found for this sum was at most 4%, already hinting that, for physical processes, the NLO EW corrections are small.
Thus, we then moved on to the evaluation of the NLO EW corrections for a physical process: gluon fusion production of the SM-like Higgs. We worked in the limit where the SM-like Higgs has couplings to SM particles close to the SM values. This is the preferred region of the allowed parameter space given the latest Higgs signal measurements from the LHC. We separated out the new contributions due to the extended scalar sector from the fixed SM-like NLO EW contributions (∼ 5%). We found that, similarly to earlier results for the decays in the real singlet model, the NLO EW corrections to gluon fusion production are of the order of a few percent in all models. In fact, even though we have examined three different scenarios, one in the RxSM and two in the CxSM, the general conclusion is similar: the new scalar sector corrections range between about ∼ −2% and ∼ 0.1% (or between about ∼ 3% and ∼ 5% if we sum the SM-like contribution). In all scenarios without dark matter, we also observe that the more we approach the SM-like limit, the smaller the new scalar sector corrections, thus consistently recovering the SM NLO EW correction. The exception is the dark matter phase of the CxSM where the dark matter particle can produce a non-zero correction even in that limit. Nevertheless, in this model the overall effect is to shift the SM correction from 5% to about 4%. The main conclusion of this study is then that future improvements in the measurements of gluon fusion production of the SM-like Higgs in these singlet models will, likely, not be able to probe radiative effects due to the new scalars. Combined with previous results in real singlet studies for the decay, where small corrections were also found, and with the fact that interference effects can be large, such minimal scalar singlet extensions will not be easily probed in future Higgs boson precision measurements. Thus, these minimal scalar singlet extended models have to be probed through direct searches for the new particles in their spectrum.
Despite the smallness of the corrections that we have found there are a few open questions. With a new dark scalar in the spectrum we have found that the corrections can deviate from the SM. In scenarios with multi-singlet dark matter, could the corrections be large enough to shift the Higgs boson properties visibly within the precision of future measurements? On the other hand we have not studied the mass region near the threshold where the SM-like Higgs can decay to a pair of scalars. Can we have a considerable enhancement of the corrections close to this threshold? This could be particularly interesting because the corresponding Higgs decay channel to a pair of scalars is kinematically suppressed near the threshold making it unlikely to be observed directly. But in the dark case this is precisely the region where the Higgs-dark-dark coupling is allowed to be larger by the dark matter relic density constraints 7 . These, and other questions, will be left for future investigations.

A Relations between bases
The L-basis relates to the Λ-basis as follows. The field dependent scalar couplings are and the field dependent gauge couplings are We have defined the mass-squared matrices, Λ ij (S) , and Λ ab We also define fermionic cubic and quartic effective vertices, involving two scalars and two fermions, which are hermitian with respect to fermionic indices: The matrices used to rotate from the Λ-basis to the λ-basis are defined through

B Loop functions
Here we present a summary of the loop functions that we have used, which can be obtained from [35]. The basic scalar loop function that we use is 8 where ∆ ≡ s 2 + x 2 + y 2 − 2(sx + sy + xy) s 2 (B.2) The various limits that are necessary are: B s (0, 0) = 2 − log s + iπ (B.8) 8 Here the factor i is an infinitesimal quantity to define the integration contour on the complex plane. Finally, the derivatives of the loop functions that are necessary to obtain the wave function renormalisation factors can all be expressed in terms of the following derivative

C Some useful identities
In this section we prove a few useful identities. We first want to relate a contraction between the effective cubic fermion-fermion-scalar couplings, with a contraction of the Yukawa couplings with the mass matrices for massless fermion states as follows: where we have used Here t L , b L and t R are the left handed Weyl fermions which give, respectively, the left handed part of the top and bottom quarks, and the right handed part of the top quark.
Since each fermion is a triplet of colour there is an extra colour contraction between each fermion/anti-fermion pair. If we organise the three sets of left handed Weyl fermions in a vector ψ I → (ψ 1 , . . . , ψ 9 ) = (t 1 L , t 2 L , t 3 and the real scalars in a another vector 4) and note that in the SM this decomposition is already in the mass eigenbasis, we can read m IJ and y IJk directly from Eq. (D.1) with the definitions in Eqs. (2.3). We can also obtain λ (F )IJk and λ (F )IJkm using, Eqs. (A.5).
The one-loop corrections from the y t couplings to the pole equations in the SM that are used in the text are written in terms of the following vector, which depends on the momentum-squared scale s: