Electroweak corrections to Higgs boson decays in a Complex Singlet extension of the SM and their phenomenological impact

The complex singlet extension CxSM of the Standard Model (SM) is a simple extension of the SM with two visible Higgs bosons in the spectrum and a Dark Matter (DM) candidate. In this paper we complete the computation of the next-to-leading (NLO) electroweak (EW) corrections to on-shell and non-loop-induced Higgs decays. Our calculations are implemented in the code EWsHDECAY which also includes the relevant QCD corrections. Performing an extensive parameter scan in the model and including all relevant theoretical and experimental single- and di-Higgs as well as DM constraints, we obtain a viable parameter sample. We find that current DM constraints are able to test the model in DM mass regions where collider searches are not sensitive. The relative EW corrections turn out to be large for scenarios with relatively large couplings, threshold effects or small leading-order (LO) widths. Otherwise, they are of typical EW size and can amount up to about 20–25%. The theory uncertainty derived from the change of the renormalization scheme dependence then is of a few per cent. While the NLO corrections applied in the constraints due to single- and di-Higgs searches impact the validity of specific parameter points, the overall shape of the allowed parameter region is not yet sensitive to the EW corrections. This picture will change with further increased experimental precision in the future and necessitates precise predictions on the theory side as presented in this paper.


Introduction
After the discovery of the Higgs boson by the Large Hadron Collider (LHC) experiments AT-LAS [1] and CMS [2] open puzzles of the Standard Model (SM) of particle physics are still awaiting their solution.One of the most prominent ones is the question for the nature of Dark Matter (DM).While there has been no direct discovery of new physics so far the precise investigation of the properties of the Higgs boson can advance our knowledge on beyond-the-SM (BSM) physics.In [3] we investigated the Complex Singlet extension of the SM (CxSM) where a complex scalar singlet field is added to the model.Imposing a Z 2 symmetry on one of the additional scalar degrees of freedom leads to a DM candidate which only couples to the Higgs boson and which can be tested at the LHC in Higgs-to-invisible decays.We computed the nextto-leading order (NLO) electroweak (EW) corrections to the Higgs decay into two DM particles and investigated the impact on the interplay between the allowed parameter space of the model and the constraints from the LHC experiments on the branching ratio of Higgs to invisible, which has been bounded to below 11% by ATLAS [4].
In this paper we further advance our precision predictions for the Higgs-to-invisible branching ratio and other observables by completing the NLO EW corrections to all Higgs boson decays which are on-shell and not loop induced.The calculated NLO EW corrections are implemented in the extension sHDECAY [5] of the program HDECAY [6,7], which computes the partial decay widths and branching ratios of the CxSM Higgs boson sector.This allows us to take into account the state-of-the-art higher-order QCD corrections and combine them consistently with our newly computed EW corrections so that we get the most precise predictions for the CxSM Higgs boson decay widths and branching ratios available at present.The program code has been made publicly available and can be downloaded at the url: https://github.com/fegle/ewshdecay/For the computation of the EW corrections we provide the possibility to choose between different renormalization schemes to cancel the UV divergences of the higher-order corrections.The application and comparison of different renormalization schemes allows us to estimate the uncertainty in the decay width predictions due to missing higher-order corrections.In order to investigate the impact of our improved predictions we first perform a scan in the parameter space of the model and keep only those data points that are in accordance with all relevant theoretical and experimental constraints.For this parameter sample we investigate the sizes of the obtained EW corrections in the various decay widths for different model set-ups, depending on which of the visible Higgs bosons behaves SM-like.We analyze the impact of the renormalization scheme choice.We investigate in detail the origin of possibly large corrections and if and how they can be tamed.This allows us to get a better understanding of the impact of EW corrections and their related uncertainties.
We then move on to the investigation of the phenomenological impact of our increased precision.For the obtained allowed parameter sample we will study the change of the Higgs-toinvisible branching ratio with respect to the given limits by the LHC searches.We will further investigate how this relates to other DM observables like direct detection and the relic density.We will also use our predictions to estimate the impact of higher-order EW corrections on di-Higgs production.The process is important for the determination of the trilinear Higgs selfinteraction and the experimental verification of the Higgs mechanism.Our improved predictions will help us to get a better understanding of the BSM physics landscape and possible DM candidates.
Our paper is organized as follows.In section 2, we give a short introduction of the CxSM and set our notation.In Section 3, we present the renormalization of the CxSM and describe the computation of the NLO decay widths.We furthermore present and discuss the calculation of the next-to-next-to-leading order Higgs-to-Higgs decay width that becomes relevant for parameter configurations with vanishing LO widths and hence also vanishing NLO widths.Section 4 describes the implementation of our corrections in sHDECAY.In Sec. 5 we detail our numerical scan together with the applied theoretical and experimental constraints before we move on to Sec. 6 which is dedicated to the numerical analysis.We first present the allowed parameter regions and investigate the Dark Matter observables in our model.We then analyze the obtained sizes of the EW corrections and investigate remaining theoretical uncertainties estimated by applying different renormalization schemes.Finally, we discuss the phenomenological impact of the newly computed higher-order corrections.Our conclusions are given in section 7.

The Model
The CxSM is obtained by adding a complex singlet field to the SM Higgs sector.We parametrize the scalar doublet field Φ and the new singlet field S as where H, S and A denote real scalar fields and G 0 and G + the neutral and charged Goldstone bosons for the Z and W + bosons, respectively.The vacuum expectation values v, v A and v S of the corresponding fields can in general all be non-zero so that all three scalar fields mix with each other.In our version of the model we impose invariance of the potential under two separate Z 2 symmetries acting on S and A, under which S → −S and A → −A, and additionally require v A to be zero, so that A is stable and becomes the DM candidate.The renormalizable potential is given by with all parameters chosen to be real.We choose v S = 0 so that the other Z 2 symmetry is broken and S and H mix with each other. 1 The mass eigenstates of the CP-even field h i (i = 1, 2) are obtained through the rotation with The mass matrix M in the gauge basis (H, S) reads 1 Compared to the CxSM with a DM candidate that was discussed in [5] and implemented in the code sHDECAY [8], we here choose a different model by imposing an extra Z2 symmetry on S (which is subsequently broken).The additional term a1S in the potential of [5] does not appear then in our potential Eq. (2), i.e. a1 = 0 here.
with the tadpole parameters T 1 and T 2 defined via the minimization conditions, At tree level we have T i = 0 (i = 1, 2).The mass of the DM candidate A is given by and the remaining mass values are the eigenvalues of the mass matrix M, The scalar spectrum of the CxSM consists of two visible Higgs bosons h 1 and h 2 , where by definition m h 1 < m h 2 , and a DM scalar A. One of the h i (i = 1, 2) is the SM-like 125 GeV Higgs boson.The mixing of the two scalars leads to a modification of their couplings to the SM particles given by the factor k i , where g H SM SM SM denotes the SM coupling between the SM Higgs and the SM particle SM .As input parameters of our model we choose in terms of which the parameters of the potential are given by Two of the input parameters are fixed; the mass of the SM-like Higgs boson has to be equal to 125.09 GeV, and the doublet VEV v is given by v = 1/ √ 2G F ≈ 246.22 GeV, where G F denotes the Fermi constant, which we choose in the following as input parameter instead of v.

Electroweak Corrections
In addition to the already calculated EW corrections to the Higgs decays into two DM particles we present in this paper our new calculation of the EW corrections to the remaining non-loopinduced on-shell decays of the visible h i .These are, if kinematically allowed, the Higgs-to-Higgs decays h 2 → h 1 h 1 , the decays into massive gauge bosons, h i → V V (V = W ± , Z) and the decays into fermions h i → f f .Note that we do not include EW corrections to the loop-induced decays into photons or gluons, as they would be of two-loop order.We furthermore do not include EW corrections to off-shell decays, so that for the SM Higgs we do not consider corrections to off-shell decays into massive gauge bosons.And also for heavier Higgs bosons below the top-pair threshold we do not include EW corrections into off-shell tops.
The higher-order corrections involve UV divergences that have to be cancelled through the process of renormalization.We replace the bare fields and parameters of our Lagrangian by the renormalized ones and their corresponding counterterms.For the isolation of the divergences we work in D = 4− dimensions so that the divergences appear as poles in .The finite parts of the counterterms are determined by the chosen renormalization scheme.We offer several schemes which fulfill the following requirements2 : -We require on-shell (OS) renormalization conditions wherever possible.
-The chosen renormalization schemes preserve gauge-parameter-independent relations between the input parameters and the computed observables.
-If possible, renormalization schemes that lead to unnaturally large corrections are avoided.We call good renormalization schemes in this context "numerically stable".
-If possible, process-dependent renormalization schemes, i.e. renormalization schemes that depend on a physical process, are avoided.
The reason for the latter condition is the exclusion of parameter scenarios where the chosen process is kinematically not allowed.Numerically stable conditions reflect a good convergence of the higher-order corrections.Gauge-parameter-independent relations allow to relate different observables to each other.On-shell conditions can make use of measured experimental parameters like the masses of the particles.

Renormalization of the CxSM
In the following we present the renormalization of the CxSM.Since the renormalization conditions have already been presented in the literature or can be taken over from other models, we restrict us here to the minimum and refer for further details to the literature.Our main goal here is setting our notation for the various renormalization schemes.

The Scalar Sector
In [3] we computed the NLO EW corrections to the Higgs boson decays into DM pairs, h i → AA and introduced the renormalization of the scalar sector.Here we only list the main ingredients and refer to [3] for further details.Electroweak corrections to the scalar sector requires the field and mass renormalization of h 1 , h 2 and A, the renormalization of the tadpoles T 1 and T 2 , of the mixing angle α and of the singlet VEV v S .We apply the following renormalization conditions: • Mass and field renormalization of h 1 , h 2 , A: For this, we choose OS conditions.
• Tadpole renormalization: The tadpole renormalization is related to the way we choose the VEVs at 1-loop order so that the minimum conditions hold.We follow the scheme proposed by Fleischer and Jegerlehner [21] for the SM.In this way, all counterterms related to physical quantities will become gauge independent.The obtained VEV is the true VEV of the theory.Note that in this scheme the self-energies contain additional tadpole contributions, and in the virtual vertex corrections additional tadpole contributions have to be taken into account if the resulting coupling that comes along with it exists in the CxSM (cf. the appendix of [9] for a discussion in the 2HDM).
• Renormalization of the mixing angle α: We use the OS-pinched and the p * -pinched scheme that was introduced in [9] for the 2HDM and in [13] for the N2HDM.It is based on the scheme proposed in [22,23] which relates the mixing angle counterterm to the field renormalization matrix constant.Combined with the Fleischer-Jegerlehner scheme for the treatment of the tadpoles and the electroweak VEV and the pinch technique [24,25] to extract the gauge-independent part unambiguously, the two schemes introduce a gaugeparameter-independent counterterm for α.The two schemes differ only in the choice of the value for the squared external momenta in the pinched self-energies, which is either the mass squared of the incoming scalar or the squared mean of the masses of the incoming and outgoing scalar, respectively.
• Renormalization of the singlet VEV v S : For the renormalization of the singlet VEV v S we offer the choice between the process-dependent scheme and the zero external momentum (ZEM) scheme [26].In the process-dependent scheme, we choose either h 1 → AA or h 2 → AA for the renormalization of v S as they contain in their h 1/2 AA coupling the parameter v S .In order to be used for the renormalization, they have to be kinematically allowed and must be different from the process that we want to renormalize.The counterterm δv S of v s is then extracted from the process by demanding that the NLO amplitude is equal to the leading-order (LO) amplitude.We call the renormalization scheme based on the h 1 → AA decay 'OSproc1' and the one based on h 2 → AA 'OSproc2'.Since it is not guaranteed for each valid parameter point of the model that these processes are kinematically allowed, we also offer the ZEM scheme introduced in [26] for both process choices, called 'ZEMproc1' and 'ZEMproc2' in the following.It avoids the kinematic restrictions of the OS decays into OS final states by setting the squared external momenta to zero in the decay process that is used for the renormalization.We ensure this counterterm to be gauge independent by using the pinched versions of the self-energies in the wave function renormalization constants that occur in the ZEM counterterm of v S .

The Gauge and the Fermion Sector
The mass and field renormalization counterterm constants, δm V and δZ V of the massive gauge bosons V = W ± , Z (we do not need to renormalize the photon) are obtained in the OS scheme.
The counterterm δZ e for the electric charge is determined from the photon-electron-positron (γeē) vertex in the Thomson limit.In the computation of the decay widths we use the G µ scheme [27] in order to improve their perturbative behavior by absorbing a large universal part of the EW corrections in the LO decay width.To avoid double counting we have to subtract the corresponding NLO part from the explicit EW corrections of our NLO calculation.We achieve this by redefining the charge renormalization constant accordingly.For details, we refer to [26].The mass and field renormalization constants δm f and δZ f are chosen in the OS scheme.For details, see [9].
In Table 1 we summarize for convenience the various renormalization schemes used in the computation of the EW-corrected decay widths.

The EW-Corrected Decay Widths at NLO
The NLO decay width Γ NLO for the decay of the scalar Higgs h i into two final state particles XX can be written as the sum of the LO width Γ LO and the one-loop corrected decay width Γ (1) , The one-loop corrected Γ (1) is obtained from the interference of the LO and the NLO amplitudes M LO and M NLO , respectively, We get for the individual LO decay widths Γ LO h i →XX of the h i Higgs decays into a lighter Higgs pair, two massive gauge bosons and a fermion pair, respectively, with h i = h 1,2 and φ a = h 1 , A, where V = Z, W ± and δ s = 1 (2) for V = W ± (Z), and with the color factor N c = 3(1) for quarks (leptons).The couplings are normalized as Since we consider only on-shell decays, the corresponding mass values must be such that m h i ≥ 2m X (X = φ a , V, f ).The h i coupling factors g h i SM SM to the SM particles SM ≡ f, V were given in Eq. (9).The trilinear couplings between the scalar particles that we need for the computation of the NLO decay widths, are given by, The one-loop correction Γ (1) consists of the virtual corrections, the counterterm contributions and -if applicable -the real corrections.The counterterms cancel the UV divergences and the real corrections the infrared (IR) divergences, if they are encountered in the virtual corrections.This happens if a massless particle is running in the loop.For example in Higgs decays into charged W ± bosons a photon can be exchanged in the loop diagrams.The IR divergences are cancelled by the real corrections, that contain bremsstrahlung contributions where a photon is radiated from the charged initial and final state particles, and of diagrams that involve a four particle vertex with a photon.For details, see e.g.Ref. [9] which describes the procedure for the 2HDM that can be easily translated to our model.The virtual corrections are built up by the pure vertex corrections and the external leg corrections.The latter vanish due to the chosen OS renormalization of the external particles.The vertex corrections comprise all possible 1particle irreducible diagrams.The counterterm contribution contains the involved wave function renormalization constants and parameter (couplings, masses, mixing angles) counterterms.

EW-Corrected Decay Width at NNLO
It turns out that for certain parameter configurations the LO decay width vanishes.This can happen for the decay h 2 → h 1 h 1 .The LO decay amplitude is given by As expected, this amplitude vanishes in the SM-like limit α = π/2(0) where the portal coupling δ 2 vanishes, the h 1 (h 2 ) decouples and h 2 (h 1 ) is the SM Higgs boson.The amplitude vanishes, however, also for parameter configurations with Since the NLO decay width is proportional to the LO decay width it also vanishes.For In the SM-like limit also the NLO amplitude M NLO h 2 →h 1 h 1 is zero.However, for tan α = v S v we can have non-vanishing contributions to M NLO h 2 →h 1 h 1 = 0 (cf.Fig. 1) so that we obtain a non-vanishing decay width which is of next-to-next-to-leading order (NNLO) and given by In contrast to the computation of the NLO decay width we now also have to take into account the imaginary parts of the NLO amplitude M NLO h 2 →h 1 h 1 .While the amplitude is still UV-finite, we have to ensure that the imaginary parts do not destroy the gauge-parameter independence of the NNLO amplitude.This is achieved by taking into account also the imaginary part of the wave function renormalization constant.It cancels the gauge-parameter dependence of the imaginary part of the leg contribution that is left over after applying our renormalization conditions.We note, that the consideration of the imaginary part of the wave function renormalization constant to get the NNLO amplitude gauge-parameter independent does not have any effect on our previous calculation of the NLO widths, as here we always take the real part of the NLO amplitude.
In order to calculate the full NNLO decay width when we move away from tan α = v S /v we would also need to take into account the term 2Re(M LO h 2 →h 1 h 1 M NNLO † h 2 →h 1 h 1 ) in Eq. (27).For this, we would need to calculate the NNLO decay amplitude M NNLO h 2 →h 1 h 1 , which is beyond the scope of this work.Instead, we only include the approximate NNLO width, given by Γ NNLO, approx Sufficiently away from tan α = v S /v, the NLO width should dominate about the NNLO one, and the incomplete NNLO calculation should not add much to the theoretical uncertainty.Sufficiently close to tan should not contribute much to the NNLO width as the LO amplitude is close to zero so that the approximation Eq. ( 28) should be good enough.In the intermediate region, however, only the complete NNLO computation would give information on the relative importance of the missing NNLO contribution, 2Re(M LO h 2 →h 1 h 1 M NNLO † h 2 →h 1 h 1 ), compared to the one taken into account by us, |M NLO h 2 →h 1 h 1 | 2 , and would hence give information on the associated theoretical uncertainty due to the incomplete NNLO calculation.
To illustrate this, we show in Fig. 2 for the benchmark point and varying mixing angle α the decay width Γ(h 2 → h 1 h 1 ) at LO, NLO and approximate NNLO given by Eq. ( 28).For α = arctan(v S /v) = 1.33 the LO and NLO widths vanish, and the result is given by the NNLO decay width, which is exact here and amounts to the value of 1.8 × 10 −5 GeV.For α > 1.33 the LO and NLO widths become non-zero with the NLO width gaining in importance compared to the approximate NNLO result3 .Only at α = 1.388, however, the NNLO width amounts to less than 10% of the NLO result.In the relatively large transition region α ∈ [1.33, 1.388] it is totally unclear, however, if the large difference between NNLO and NLO is due to incomplete cancellations in the approximate NNLO width and/or due to a small NLO width close to its vanishing point at α = 1.33.Note, finally, that as expected all amplitudes vanish in the SM-like limit which is obtained for α = π/2.

Implementation in sHDECAY: EWsHDECAY
We implemented the EW one-loop corrections to the Higgs decays derived in this work in the singlet extension sHDECAY of the code HDECAY [6,7] where we updated the underlying HDECAY version to version 6.61.The resulting code is called EWsHDECAY.Note, that we impose an extra Z 2 symmetry on S in contrast to the CxSM version implemented in sHDECAY [8].In order to use our version of the CxSM therefore the input parameter a 1 has to be set to zero in the input file, cf. also footnote 1.The Fortran code HDECAY and thereby its extension sHDECAY include the state-of-the-art QCD corrections to the partial decay widths.For the consistent combination of our EW corrections with HDECAY which uses the Fermi constant G F as input parameter, we use the G µ , respectively the G F scheme, in the definition of our counterterm for the electric charge, cf.Subsection 3.1.2.
Note, that sHDECAY also calculates the off-shell decays into final states with an off-shell topquark t * , φ → t * t (φ = h 1 , h 2 ) and into final states with off-shell gauge bosons V * (V = W, Z), φ → V * V * .The EW corrections are only computed for on-shell decays, however.This means that if only off-shell decays are kinematically allowed for the mentioned final states, then the EW corrections are not included and only the LO decay width is calculated.Be reminded also, that we did not include EW corrections to the loop-induced decays into photonic and gluonic final states.We furthermore assume that the QCD and EW corrections factorize.The relative QCD corrections δ QCD are defined relative to the LO width Γ CxSM,LO calculated by sHDECAY, which contains (where applicable) running quark masses in order to improve the perturbative behavior.The relative EW corrections δ EW on the other hand, are obtained by normalizing to the LO width with OS particle masses.The QCD and EW corrected decay width into a specific OS and non-loop-induced final state, Γ QCD&EW is hence obtained as And the branching ratio of the higher-order (HO) corrected decay width h i (i = 1, 2) into a specific final state is calculated by with the total width given by, As mentioned above, the QCD corrections to the decays into colored final states are those as implemented in the original HDECAY version, for details see [7].As mentioned above, the EW corrections are only included for non-loop-induced final states and for on-shell decays.Otherwise the off-shell tree-level decays are used if implemented in sHDECAY, as is the case for the tt, ZZ and W W final states, hence and Note, that for the fermionic final states we included the EW corrections only of the third generation.As the decay widths in fermion pairs of the first two generations are much smaller, the effect of not including the EW corrections here is negligible.The partial widths Γ HO (h i → XX) in Eq. ( 31) are given as defined in Eqs.(30) to (34).
The The model for which we compute the EW corrections is the CxSM in the DM phase, so that the user has to choose 'icxSM=4' after setting 'isinglet=1'.The input values are then given in the block named 'complex singlet dark matter phase'.We consider a specific version of the model where the parameter a 1 does not appear in the model, so that it has to be set to 0 for the computation of the EW corrections.Therefore, if the parameter is chosen to be non-zero, no EW corrections will be calculated.Note, that m 1,2 denote the lighter and heavier of the two visible Higgs boson masses, and m 3 is the mass value of the DM particle.For the inclusion of the EW corrections the input file hdecay.in of sHDECAY has been extended by the following lines: The warning reminds the user that our computation of the EW corrections only applies to the CxSM in the DM phase with the a 1 = 0 in the potential.By setting the input parameter 'ielwcxsm' equal to 1 (0), the EW corrections will be computed (omitted).The next three parameters choose the applied renormalization scheme.Setting 'vsscheme' to 1, the singlet VEV counterterm will be computed in the process-dependent scheme using the on-shell decay h 1 → AA (h 2 → AA) by choosing 'pdprocess =1 (2)' (called 'OSproc1' and 'OSproc2', respectively, in Tab. 1).When 'vsscheme' is set to 2, then the process-dependent renormalization scheme for v s is evaluated at zero external momenta for the chosen decay process (called 'ZEMproc1' and 'ZEMproc2', respectively in Tab. 1).With 'ralph mix=1 (2)' the OS-pinched (p * -pinched) scheme is chosen for the renormalization of the mixing angle α.
The next parameter setting chooses with 'DeltaE' the detector resolution needed in the computation of the IR corrections, which we set by default equal to 10 GeV, cf.[9].The following two parameters concern the computation of the NNLO corrections to the decay width h 2 → h 1 h 1 .If 'NNLOapp' is set equal to 1 then the NNLO width is computed for parameter configurations with tan α values in the vicinity of v S /v, namely tan α ∈ {v s /v−deltaNNLO, v s /v+deltaNNLO}.Be aware, however, that for tan α = v s /v, the NNLO computation is incomplete, as discussed above.
The last four parameters refer to the conversion of the input parameters.The change of the renormalization scheme allows for an estimate of the uncertainty in the NLO EW corrections due to the missing higher-order corrections.For such an estimate to be meaningful also the input parameters have to be changed consistently.For the conversion of the input parameters, when going from one to the other renormalization scheme, our code uses an approximate formula based on the bare parameter p 0 which is independent of the scheme so that, where p x and δp x denote the renormalized parameter p and its counterterm δp in the standard scheme and in the user specified scheme, respectively.Setting the flag ' Paramcon = 1', the parameter conversion is applied.In this case, the user can then also choose which renormalization scheme is the 'standard' scheme, by setting the three parameters '(stdvs,stdproc,stdalpha)'.
If the user specified renormalization scheme '(vsscheme,pdprocess,ralph)' and the standard scheme '(stdvs,stdproc,stdalpha)' are identical, the input parameters v S and α remain unchanged.Otherwise they are converted from their values understood to be given in the standard scheme '(stdvs,stdproc,stdalpha)' to the values they take in the user specified scheme '(vsscheme,pdprocess,ralph)'.The converted parameter values are given out in the file 'Paramconversion.txt'.If, however, the user chooses 'Paramcon = 0' the parameters are not converted when the specified renormalization scheme differs from the standard scheme.In this case a comparison of the results for different renormalization schemes to estimate the theoretical uncertainty would not, however, be meaningful.
Note that we do not give out the NLO value of a specific decay width if it becomes negative 4 .In this case, and also if a chosen renormalization scheme cannot be applied because it is kinematically not allowed, the decay width is computed and given out at LO in the EW corrections together with a warning.Remark, finally, that when we take the SM limit of our model and compare the results for the EW corrections to those obtained from HDECAY for the SM, we differ in the decay into gauge boson final states as HDECAY includes the EW corrections also to the off-shell SM Higgs decays into W W/ZZ and to the loop-induced decays into gluon and photon pairs.The code can be downloaded at the url: https://github.com/fegle/ewshdecayApart from short explanations and user instructions, we also provide sample input and output files on this webpage.

Parameter Scan
In order to obtain viable parameter points for our numerical analysis we performed a scan in the parameter space of the model and kept only those points which fulfill the relevant theoretical and experimental constraints described below.They are checked for by ScannerS [42,43] which we used for the parameter scan, with the scan ranges summarized in Tab.The theoretical constraints that are taken into account are boundedness from below, perturbative unitarity 5 and stability of the vacuum.For detailed information, we refer to [3].
As for the experimental constraints, we first note that in our model the ρ parameter is equal to 1 at tree level and there are no tree-level flavour-changing neutral currents, as the gauge singlet does not couple to fermions and gauge bosons in the gauge basis.We check for compatibility with EW precision data by requiring the S, T, U parameters [47] to be consistent with the measured quantities at 95% confidence level.Compatibility with the LHC Higgs data and exclusion bounds is checked through the ScannerS link with HiggsSignals [48,49] and HiggsBounds [50,51].We require the signal rates of our SM-like Higgs boson to be in agreement with the experimental data at the 2σ level.Through the link with MicrOmegas [52] we ensure that the DM relic density of our model does not exceed the measured value.Concerning DM direct detection, the tree-level cross section is negligible in this model [53,54], but the loopcorrected DM-nucleon spin-independent cross section [55,56] has to be checked to be below the experimental bounds [57][58][59] with the most stringent bounds given by the LUX-ZEPLIN experiment [59].Further detailed information on the experimental constraints that we applied is given in [3].
In our sample we take off parameter scenarios where the deviation of any other neutral scalar mass from the h 125 mass is below 2.5 GeV.This suppresses interfering di-Higgs signals which would require a further special treatment to correctly consider theory and experiment beyond the scope and focus of the present work.
We also used the program BSMPT [60,61] to check if the thus obtained parameter sample leads to vacuum states compatible with an EW VEV v = 246 GeV also at NLO, i.e. after including higher-order corrections to the effective potential.Provided this to be the case, we checked if we have points providing a strong first-order EW phase transition (SFOEWPT), one of the Sakharov conditions [62] to be fulfilled for EW baryogenesis.We found that none of the parameter points with a viable NLO EW vacuum provides an SFOEWPT.
Finally, we applied constraints from the resonant di-Higgs searches performed by the LHC experiments.We proceed here as described in [63] to which we also refer for a detailed description, giving here only the most important points.In order to apply the constraints, we calculated for parameter points of scenario I, where m h 2 > 2m h 1 ≡ 2m h 125 the production cross section σ(h 2 ) at next-to-next-to-leading-order QCD using the code SusHi v1.6.1 [64][65][66] and multiplied it with the branching ration BR(h 2 → h 1 h 1 ).We then checked if the thus obtained rate is below the experimental values given in Refs.[67][68][69] for the 4b, [70][71][72] for the (2b)(2τ ), [73,74] for the (2b)(2γ), [75] for the (2b)(2W ), [76] for the (2b)(ZZ), [77] for the (2W )(2γ) and [78] for the 4W final states.There are also experimental limits from non-resonant searches.They do not constrain our model so far, however, as will be discussed below.
The parameter points which are found to be compatible with all applied constraints can be divided up into two samples, depending on which of the visible Higgs particles is the SM-like one, called h 125 in the following.For the two samples, we will adopt the following notation, 6 Numerical Analysis We start our numerical analysis by investigating what is the impact of our extended scan compared to [3] on the allowed parameter regions of the model.We move on to the discussion of the overall size of the computed EW corrections before investigating the remaining theoretical uncertainty due to missing higher-order corrections.We then discuss the phenomenological impact of these corrections w.r.t the Higgs decays into DM particles and the impact on the allowed parameter regions.Finally, we investigate Higgs pair production in the context of our model and how it is affected by our corrections.

Allowed Parameter Regions and DM Observables
In our new scan we go up to m A values of 1 TeV, as compared to [3], we now also allow for scenarios where the decay of the SM-like Higgs boson into DM particles is kinematically closed.This leads to an enlarged parameter space, a larger variation in the particle couplings and masses and hence effects in the NLO EW corrections that would not appear in the more restricted sample.We therefore show the corresponding plot to Fig. 1 in [3]. Figure 3 (upper) displays the allowed combinations of α and v S for scenario I, i.e. h 1 = h 125 , (left) and for scenario II, i.e. h 2 = h 125 , (right). 6The lower plot shows the corresponding values of non-125-GeV Higgs boson mass m S vs. v S .Compared to [3], we see that α and v S are not linked any more as for the large DM matter masses, that are now included in our scan, the kinematic constraints inducing this relation are not applicable any more.This leads also to a larger allowed maximum range for |α| which has increased in scenario I from about 0.27 to about 0.34 for almost all parameter points of scenario I.This is the maximal allowed value compatible with the Higgs data except for parameter points where the visible Higgs bosons are close in mass and the Higgs signal is built up by the two Higgs bosons.These are the outliers with larger mixing angles up to about 0.37.If we dropped the exclusion of parameter points with non-SM-Higgs masses closer than 2.5 GeV to the SM-like Higgs mass then the mixing could be even larger. 7Also in scenario II larger deviations from α = ±π/2 are possible for the increased parameter scan. 8The lower plots, apart from becoming more dense, did not change the shape.The perturbative unitarity constraints are reflected in the linear relation between m S and v S , cf. the discussion in [3].
Figure 4 displays the allowed values of the DM mass m A versus the non-SM-like Higgs mass mass m S for scenario I (red points) and scenario II (blue points).For m A values below 62.5 GeV, only mass values in the vicinity of m h 125 /2 or m S /2 are allowed, as then efficient DM annihilation via h 125 or the non-SM-like Higgs particle can take place such that the DM relic density constraints can be met (cf.our discussion in [3]).For large m A values where annihilation via other SM particles can take place all mass values up to the upper limit of our scan are allowed.The limits from XENON1T [57] and LUX-ZEPLIN [59] are given by the dashed and full line, respectively.The grey shaded region corresponds to the neutrino floor.The color code represents the non-SM-like Higgs mass value in GeV.
In Fig. 5 we show the effective spin-independent nucleon DM detection cross section as a function of the DM mass m A for scenario II (left) and scenario I (right).The multiplication of the spin-independent cross section SI σ with the factor accounts for the fact that the relic density (Ωh 2 ) A of our DM candidate A might not account for the entire observed DM relic density (Ωh 2 ) obs DM = 0.120 ± 0.001 [79]. 9In our model the direct detection cross section at LO is negligible due to a cancellation [53,54] so that we present the one-loop result calculated in [55,56].In [3], we made a mistake that we correct here so that now both the exclusion limits from Xenon1T [57] (dashed) and LUX-ZEPLIN [59] (full) as well as the neutrino floor (grey region) move down by one order of magnitude.The effect is that for both scenarios I and II we have now parameter points above the neutrino floor that can hence be tested by direct detection experiments.Moreover, we find that in the parameter region m A ≥ m h 125 /2, where we cannot probe DM at the LHC, the LUX-ZEPLIN experiment [59] is sensitive to the model in the region 66 GeV ∼ < m A ∼ < 78 GeV.Future increased precision in the direct detection experiments will allow to test the model in (large) parts of the still allowed m A range in scenario I (scenario II).
Figure 6 displays the relic density for all points passing our constraints for scenario II (left) and scenario I (right).As can be inferred from the plot for the whole allowed region up to m A = 1 TeV, where our scan ends, there are parameter scenarios that saturate the relic density.The color code represents the non-SM-like Higgs mass value in GeV.

Size of NLO EW Corrections
In Tab. 4 we show the relative NLO corrections to the partial widths Γ of the decay of h i (i = 1, 2) into the various possible final states XX (X = τ, b, t, W, Z, A, h 1 ), defined as for scenario I (Tab.4) and scenario II (Tab.5) and different renormalization schemes.These are for the renormalization of the singlet VEV v S the process-dependent and the ZEM scheme, using either the h 1 → AA or h 2 → AA decay and for the renormalization of α the OS-or the p -pinched scheme.The notation for the schemes is summarized in Tab. 1.We also note that we discard all parameter points where the large relative corrections are beyond -100%, in order to not encounter unphysical results with negative decay widths at NLO.We furthermore emphasize that the results given in a specific renormalization scheme assume the input parameters to be given in this specific renormalization scheme.We hence do not start from a certain scheme a and then move on to a different scheme b by consistently converting the input parameters of scheme a to scheme b.The results given in the various schemes hence have to be looked at and discussed each one by one.They are meant to give an overview of what sizes of corrections can be expected in the various schemes.
As can be inferred from Tab. 4, for h 1 = h 125 and v S renormalization through the on-shell decay processes h 1,2 → AA the relative NLO EW corrections are of typical size with at most 25%, apart from the corrections to the decay h 2 → h 1 h 1 in case the decay h 2 → AA is chosen for the renormalization of v S .The reason for this is two-fold.Large relative corrections appear due to a very small LO width resulting from parameter points where tan α ≈ v S /v.Since the case h 1 = h 125 requires small α values (below about 0.37) in order to comply with the experimental constraints, this requires v s ∼ < 90GeV.4: Scenario I: Relative NLO corrections δ EW (h i → XX) for all possible decay processes in scenario I, i.e. m h 1 = 125.09GeV, and different renormalization schemes ren v S -ren α .For δv S we have OS renormalization and ZEM renormalization in the two possible processes, ren v S = OSproc1, OSproc2, ZEMproc1, ZEMproc2.For δα we have the two possible schemes ren α = OS, p * .The renormalization schemes and their notation are summarized in Tab. 1. trilinear couplings g h 1 h 2 h 2 , g h 2 h 2 h 2 and g h 2 AA (remind that cos α ≈ 1 for h 1 = h 125 ). 10 In case we use h 1 → AA for the renormalization we are kinematically more constrained and also have to comply with the DM constraints so that points with vanishing LO coupling can no longer occur in the sample and we do not have large corrections here.
When we look at the corrections for scenario I using ZEM renormalization of v S , given in Tab. 4, we see that the relative corrections to the h 2 decays apart from h 2 → AA in the ZEMproc2 scheme become large.The reasons can be as above small LO widths combined with large couplings.Additionally, threshold effects in the derivative of the B 0 function [80] can appear at 2m A − m h 2 ∼ > 0. Such kinematic scenarios can only appear in the ZEM scheme where for the renormalization, we are not restricted any more to scenarios with on-shell h 1,2 → AA decays.These threshold effects increase the diagonal wave function renormalization constants leading to large corrections.In the decays h 2 → AA we cannot have these threshold effects, as h 2 must decay on-shell into AA.'proc2' renormalization scheme, we see that in the latter they are of moderate size, whereas for 'proc1' renormalization they become large.The reason is that the counterterm δv S in the ZEM scheme can become large when there is a large mass difference between the initial and final state particles in the process used for renormalization in the ZEM scheme.
Note that in the ZEM scheme for all but the Higgs-to-Higgs decays the sizes of the relative corrections are the same independently of the used process for the renormalization of v S since these decays do not necessitate the renormalization of v S .In the OS case for the renormalization of v S , cf.Tab.4,they differ, however, because due to the kinematic restriction we have different parameter samples depending on which process we use for the renormalization of v S .
In scenario II, where h 2 = h 125 , we infer from Tab. 5 (upper) that now also the decay h 2 → h 1 h 1 can get large relative corrections in case h 1 → AA is used for renormalization.The reason for this enhancement is a vanishing LO decay width for tan α = v S /v.Since in scenario II α ≈ π/2, this does not require v S to be small and hence such parameter scenarios are not in conflict with the DM constraints.Note finally, that the h 2 decays into W W, ZZ and tt are kinematically closed in scenario II, and we do not consider off-shell decays when we compute the EW corrections.
If we use h 2 → AA for the renormalization of v S then the h 1 → bb and τ τ decays can get large relative corrections.This can be due to threshold effects for 2m A − m h 1 ∼ > 0 or due to large couplings between scalars.For the h 1 → AA decay, the kinematics required for this process combined with the experimental constraints leads to no strong enhancement of the relative corrections.The reasons for large relative corrections to the h 2 → h 1 h 1 decays are a vanishing LO width or large couplings (which in contrast to scenario I are not correlated with a vanishing LO width any more) or threshold effects for 2m A − m h 1 ∼ > 0.
If we choose ZEM renormalization, cf.Tab. 5, then h 1 → bb and h 1 → τ τ get large corrections mainly due to threshold effects in the derivative of the B0 function for 2m A − m h 1 ∼ > 0. This is independent of the chosen process for the v S renormalization as it does not enter in these processes.We also have a parameter point here where the large corrections are due to large couplings between scalar particles.The h 2 decays only show large corrections in the decays into scalar particles.For the h 2 → h 1 h 1 decays the reasons can be a small LO width, large couplings or threshold effects.In the case of the h 2 → AA decays we cannot have threshold effects from 2m A − m h 2 ∼ > 0 as they are kinematically not possible.If we use h 1 → AA for the renormalization of v S threshold effects for 2m A − m h 1 ∼ > 0 can appear and lead to large corrections.Large couplings between scalar particles can additionally enhance the corrections independently of the used renormalization process for v S .

Size of NLO EW Corrections after Cuts
To confirm the above observations we excluded in a next step from our parameter sample all parameter points where 0 GeV ≤ 2m A − m h i ≤ 9 GeV (i = 1, 2) to avoid large threshold corrections, and points with λ ijkl > 4π, where λ ijkl stands for all possible quartic couplings between the scalars of our model.This eliminates possibly large couplings (which were not eliminated yet by the constraints from perturbative unitarity).The resulting relative NLO corrections are summarized in Tab.6 for scenario I and in Tab.7 for scenario II.7: Scenario II after cuts: Relative NLO corrections δ EW (h i → XX) for all possible decay processes in scenario II, i.e. m h 2 = 125.09GeV, and different renormalization schemes ren v S -ren α after applying the cuts described in the text.small, at most 29%.Only the relative corrections to the h 2 → h 1 h 1 width in the OSproc2 and the ZEM schemes and the h 2 → AA decays in the ZEMproc1 scheme can become large.The large relative corrections for h 2 → h 1 h 1 are due to a small LO width that also entails large couplings between the scalar particles (see discussion above).In the OSproc1 scheme the relative corrections are small also for this decay channel as here the additional kinematic constraint 2m A < m h 1 = 125.09GeV allows for DM decays so that additional experimental constraints have to be considered such that the conditions for the large corrections are not met.Additionally, in both decay channels we can have large corrections in the ZEMproc1 scheme because of the counterterm δv S which can become large when there is a large mass difference between the initial and final state particles in the process h 1 → AA used for its renormalization.
In scenario II, cf.Tab. 7, all relative corrections remain below 19% apart from those to the h 2 → h 1 h 1 decays.They are due to small LO widths.Additionally, in the ZEMproc1 and ZEMproc2 scheme the NLO corrections can become large again due to the counterterm δv S that can become large when there is a large mass difference between the initial and final state particles in the process used for its renormalization (here h 1 → AA or h 2 → AA).
We can summarize: For h 1 = h 125 , the relative corrections to both h 1 and h 2 decays are in general decent being at most 20 to 30% as long as OS conditions are applied in the processdependent renormalization of v S .The exception are the relative corrections to the decay h 2 → h 1 h 1 which can become large due to small LO width entailing also large couplings among the scalars.For h 2 = h 125 also the corrections to h 1 decays into bb and τ τ can get large due to threshold effects or large couplings between the scalars.If we use the ZEM scheme for the process-dependent renormalization of v S , in scenario I all h 2 decays get large corrections, in scenario II the h 1 decays into bb and τ τ and the h 2 decays into scalar pairs get large corrections.For a better perturbative convergence, it is hence advisable to use the OS scheme in the process-dependent renormalization of v S .However, this also restricts the possible parameter scenarios that can be used, as the kinematic constraints for the OS decays used for renormalization have to be met.We therefore use the ZEM scheme, more specifically the ZEMproc2 scheme, in the following as our standard scheme for the NLO corrections.

Theoretical Uncertainty
We can get an estimate of the theoretical uncertainty due to missing higher-order corrections by investigating the NLO results in different renormalization schemes.For this comparison to be meaningful, we have to consistently convert the input parameters of scheme a to scheme b when moving from scheme a to scheme b.For definiteness, in this investigation our starting scheme a is ZEMproc2-OS with the input parameters assumed to be given in this scheme.We then convert the input parameters to the other schemes under investigation and compute the higher-order corrections in these schemes.For the conversion of the input parameters we use an approximate formula based on the bare parameter p 0 which was given in Eq. (35).
In the following, we give the results of the EW corrected decay widths for two sample benchmark points.We define the relative uncertainty on the investigated Higgs decay width Γ due to missing higher-order corrections, estimated based on the renormalization scheme choice as (b=OSproc1-OS,OSproc1-p * , OSproc2-OS, OSproc2-p * , ZEMproc1-OS, ZEMproc1-p * , ZEMproc2-p * ) The first benchmark point is chosen from the scenario I sample, and is defined by the following input parameters (m h 1 = 125.09GeV) BP1: m h 2 = 590.48GeV, m A = 61.93GeV , v s = 446.14GeV, α = 0.1654 .
Table 8 displays for BP1 the relative EW corrections δ EW Eq. ( 38) and the relative uncertainty δ ren Eq. (39) due to missing higher-order corrections based on a renormalization scheme change.
The table also contains the input values which change with the renormalization scheme, namely v S and α.We show the results for all applied renormalization schemes.First of all note, that δ ren = 0 for ZEMproc2-OS for all decays as this is our original renormalization scheme.For all other schemes, we see that α barely changes and v S changes by at most 2% when the renormalization scheme is altered.In line with this observation, we see that the change of the EW corrected decay widths with the renormalization scheme is at most 3.1%.This is as expected for relative EW corrections δ EW that are found to be of small to moderate size for this parameter point, not exceeding 15%.
We investigate a second benchmark point characterized by larger EW corrections and theoretical uncertainties.It is chosen from the scenario II sample, and is defined by the following input parameters (m The relative EW corrections and uncertainties in the corrected decay widths for BP2 are shown in Table 9.Since for these parameter values an OS decay h 1 → AA is kinematically not possible we cannot use this process for renormalization so that the EW corrections cannot be calculated in the OSproc1-OS and OSproc1-p * schemes.In the OSproc2-OS and OSproc2-p * schemes the decay width h 2 → AA is used for the renormalization so that the LO decay width Relative corrections δ EW and relative uncertainties δ ren due to the renormalization scheme change at NLO EW for the various decay widths of benchmark point BP1 and for all applied renormalization schemes.The input scheme is ZEMproc2-OS.Also given are the input parameters α and v S , which change with the renormalization scheme. Γ(h 2 → AA) is compared to the NLO decay width, computed in the ZEMproc2-OS scheme, for the derivation of the theoretical uncertainty, leading to a relatively large value of δ ren which is, however, artificial, due to the LO-NLO comparison.We see, however, that compared to BP1 we have larger theoretical uncertainties in the process h 2 → AA in the ZEM schemes, ranging up to ∼ 19%.This can be understood by looking at the input parameters.The large difference between the non-SM-like Higgs mass m h 1 and v S leads to relatively large couplings involved in the wave function renormalization constants, increasing the counterterms, the corrections and the relative uncertainty.The uncertainties of the remaining NLO decays are of small to moderate size.
Overall, we observe that the uncertainties in the EW corrections are small or of moderate size, apart from scenarios where large scalar couplings are involved.Here the corrections and the related uncertainties can become significant to large.

Higgs-to-Invisible Decays
We first checked if the approximate NLO branching ratio of h 125 into a DM pair, BR NLO,approx (h 125 → AA), defined in [3], where we included the NLO EW corrections only in the decay width Γ(h 125 → AA) was sufficiently good, by comparing it with the BR NLO (h 125 → AA) ≡ BR(h 125 → AA) defined in Eq. ( 31), where we included the EW corrections to all on-shell, non-loop-induced decays in the width.Figure shows the relative δ approx between the two approaches, defined as as a function of the non-SM-like scalar mass m S for the four different renormalization schemes OSproc1/2-OS, OSproc1/2-p * , ZEMproc1/2-OS, and ZEMproc1/2-p * .The relative difference can be written as where δ = 1 (0) for m h 125 ∼ > 2m S ( ∼ < 2m s ) and k i has been defined in Eq. ( 9).Note that in the analysis of [3], we used a sample of parameter points, where h 125 decays into DM particles are kinematically possible, hence m A ∼ < m h 125 /2.This means that for scenario I (red points in the figure) where h 1 ≡ h 125 , the non-SM-like decay h 2 → AA is always possible and can be used for the renormalization of v S .In the other case, h 2 ≡ h 125 , not all parameter points necessarily fulfill the condition where h 1 → AA can be used for renormalization.Then no EW corrections are calculated and both our new and the previous results should agree so that δ approx = 0.These are the light blue points in the upper left and upper right plots of Fig. 7 which show the results for the case that the OS condition is used for the renormalization of v S .Using ZEM, we do have this kinematic constraint so that δ approx = 0 for all parameter points.We remark, that for the light blue points δ approx is not exactly zero.This is due to a small off-set in the total widths (below 3%) used in [3] and in this paper which stems from using an older HDECAY version in [3].For the blue points (scenario II) with m S ≤ m h 125 /2, the h 125 → h 1 h 1 decay is open and its NLO corrections that were not included in [3] play a role in δ approx .We see that their effect is larger if the ZEM scheme is applied as we already learned from our investigations in Subsec.6.2.Nevertheless, overall the deviations remain below ±4%.In all remaining points, the deviations in δ approx are due to the inclusion of the EW corrections to the fermionic decays.They remain below 1.5%.The plot demonstrates that for the allowed parameter points in both scenarios I and II the approximation used in [3] is quite good and the deviation is well below the experimental precision.The application of the approximation in [3] for this decay process was hence justified.

Impact on Allowed Parameter Regions
In order to investigate the phenomenological impact of the EW corrections we generated two parameter samples: Sample 1 is the parameter sample which we also used in the previous sections.It is based on the check of the single-and di-Higgs constraints using the branching ratios at LO in the EW corrections.Sample 2 is the parameter sample, where we include the EW corrections in the branching ratios as defined in Eq. (31).As renormalization scheme, we use the ZEM scheme, more specifically ZEMproc2-OS, as it allows to use the largest amount of parameter scenarios11 .We do not apply any mass cuts or cuts on the couplings to suppress large NLO corrections.We only make sure to exclude parameter scenarios where the EW corrections lead to unphysical negative decay widths.
We then investigate if the allowed parameter regions of sample 2 change with respect to sample 1.We find that more points are rejected when we include the NLO corrections compared to the case where we only take LO decay widths (see also discussion below).Unfortunately, however, the shape of the allowed parameter regions overall did not change.The NLO corrections hence do not have a direct impact on the exclusion of certain parameter regions of our model yet.With increasing experimental precision in the future, it is evident, however, that higher-order EW corrections have to be considered.

Higgs Pair Production
Higgs pair production is one of the most prominent processes investigated at the LHC.Its measurement allows the extraction of the trilinear Higgs self-interaction [81] and thereby the ultimate experimental test of the Higgs mechanism [82,83].At the LHC, the dominant Higgs pair production process is gluon fusion into Higgs pairs [81,84,85] which at LO is mediated by heavy quark triangle and box diagrams [86][87][88].The experiments search for Higgs pairs both in resonant and non-resonant Higgs pair production.The limits on the Higgs self-coupling between three SM-like Higgs bosons in terms of the SM trilinear Higgs coupling are at 95% CL κ λ ∈ [−0.4,6.3] given by ATLAS [89] and κ λ ∈ [−1.24, 6.49] given by CMS [90].After applying all constraints described in Sec.
So there is still some room left to deviate from the SM value, values equal to zero are excluded, however.When compared to other models as e.g.those discussed in [63], we see that the trilinear self-coupling of the SM-like Higgs in the CxSM, when h 1 ≡ h 125 , is in general more constrained than in the CP-conserving (2HDM) and CP-violating 2-Higgs-Doublet Model (C2HDM), the next-to-2HDM (N2HDM) or the next-to-minimal supersymmetric extension (NMSSM) where for some of the models the coupling can still be zero and also take negative values.In case of scenario II, the lower limit is below the one in the (C)2HDM and the N2HDM.It is interesting to note, that in the case that we do not cut out degenerate Higgs scenarios with large mixing angles where the discovered Higgs signal is built up by two scalar resonances close to each other, larger deviations in the trilinear couplings would be possible, with negative or zero trilinear coupling values.After application of the NLO EW corrections and hence using the sample 2, we see that the Higgs self-couplings are not further constrained compared to sample 1.
While the range of allowed trilinear Higgs self-couplings is not sensitive yet to the EW corrections on the Higgs decays, the di-Higgs production cross sections are sensitive.More specifically, the maximum allowed di-Higgs cross section depends on whether or not EW corrections are included.When we use sample 1 where we calculate the resonant di-Higgs rate by multiplying the NNLO QCD h 2 production cross section obtained from SuShi with the LO branching ratio taken with caution, however, due to remaining theory uncertainties.And we remind that we do not apply the resonant di-Higgs constraints when the narrow-width approximation is not justified any more.Since the total width can be smaller or larger depending on the parameter point and the inclusion or not of the EW corrections in the decay widths this of course also has an impact on the excluded region.

Points with Vanishing Coupling λ h 1 h 1 h 2
The discussion so far has shown that the NLO corrections do not impact the shape of the parameter regions that are still allowed but they affect individual parameter points.We want to discuss now the interesting situation of parameter configurations where the trilinear coupling λ h 1 h 1 h 2 is close to zero.Here, the impact of the NLO corrections could become important.We found that this is indeed the case.For some of these parameter points the NLO corrections to h 2 → h 1 h 1 shifted the point above the allowed Higgs data constraints.As an example we look at the benchmark point BP5, defined as ( We are close to the case tan α = v S /v where the coupling λ h 1 h 1 h 2 vanishes so that the partial decay width Γ(h 2 → h 1 h 1 ) is small.The inclusion of the NLO corrections can considerably change this.And indeed we find, that the relative correction is (in the ZEMproc2-OS scheme) This shift is so large that the parameter point is not allowed any more when NLO corrections are included.

Conclusions
In this paper we computed the EW corrections to the decay widths of the visible Higgs bosons of the CxSM, the SM extended by a complex singlet field.Applying two separate Z 2 symmetries and requiring one VEV to be zero, the model contains one stable DM candidate.The EW corrections complete a previous computation of our group where the EW corrections to the Higgs decays into a DM pair were computed, by calculating the EW corrections to the decays into on-shell fermionic, massive gauge boson and Higgs pair final states.We do not compute corrections to off-shell decays nor to the loop-induced decays into gluon or photon pairs.We apply two different renormalization schemes for the renormalization of the mixing angle and four different schemes for the renormalization of the singlet VEV.The latter is renormalized process-dependent both on-shell and in the so-called ZEM scheme that we introduced in a previous calculation.The renormalization of the mixing angle is gauge-parameter independent by applying the alternative tadpole scheme and extracting the gauge-independent self-energies through the pinch technique.Our computation has been implemented in the code EWsHDECAY which has been publicly made available.The user has the options to choose between different renormalization schemes, including also the possibility to convert the input parameters, which has to be applied for a meaningful interpretation of the change of the renormalization scheme.
For the numerical analysis we performed an extensive scan in the CxSM parameter space including theoretical and experimental constraints as well as DM constraints.We also considered the constraints from di-Higgs searches at the LHC.We found that our parameter points are able to saturate the DM relic density, and that direct search limits are already sensitive to parts of the model for parts of the DM mass region where the collider searches in invisible SM-like Higgs decays are not sensitive due to kinematic constraints.
The analysis of the impact of the relative EW corrections showed that they can be large in case of relatively large involved couplings, due to threshold effects or because the LO decay width is small.Here, the NLO width can become important and exclude parameter points that would be allowed at LO.For parameter configurations with exactly vanishing LO, and hence also, NLO width, we included the option to calculate the next-to-next-to-leading-order (NNLO) width.The applied formula is, however, only approximate when moving away from the exact zero LO value.If we exclude parameter configurations with large couplings or large threshold effects, the relative corrections are of typical EW size of up to 25%.The estimated theory uncertainty due to missing higher-order corrections based on the change of the renormalization scheme with converted input parameters is of a few per cent.The phenomenological impact of the computed EW corrections is that specific parameter points that are allowed at LO would be excluded at NLO.The overall distribution of the allowed parameter points, however, is not yet affected by the EW corrections.With increasing future precision in the experiments, this will change, however, and the here computed corrections provide an important contribution to reliably constrain the allowed parameter space of the model.
We also investigated the phenomenology of di-Higgs production.Our parameter scan revealed that deviations of the SM-like trilinear Higgs self-coupling from the SM value are still allowed.However, vanishing or negative trilinear coupling values are not allowed any more.The allowed range of the trilinear coupling barely changes when we include the NLO instead of the LO branching ratio in the Higgs-to-Higgs decay in resonant di-Higgs production and compare with the experimental results on resonant di-Higgs searches.For the investigated parameter points with maximum di-Higgs cross sections at LO, respectively at NLO, we found that the relative corrections on the di-Higgs cross sections are small and below the theory uncertainty.This picture may change, however, when the full EW corrections to the complete di-Higgs production cross section are considered.

Figure 1 :
Figure 1: Sample diagrams contributing to the NLO decay width h 2 → h 1 h 1 for tan α = v S /v.

Table 2 :
The scan ranges used for the generation of parameter points with ScannerS.The mass m s denotes the scalar mass of the non-125 GeV Higgs boson.the SM input parameters needed for our calculation are given in Tab. 3.These are the values suggested by the LHC Higgs Working Group [44].

3 :
The SM parameter values used in the numerical evaluation taken from[45] and following the recommendations of[44].Here, m b denotes the MS bottom quark mass which is taken at the scale of the bottom quark mass, m b , and h 125 denotes the SM-like Higgs boson with a mass of 125.09GeV.

Figure 4 :
Figure 4: Allowed DM mass values m A versus non-SM-like Higgs mass values m S in GeV for scenario II (h 2 = h 125 , blue points) and scenario I (h 1 = h 125 , red points).

Figure 5 :
Figure5: Effective spin-independent nucleon DM cross section as a function of the DM mass for scenario II (left) and scenario I (right).The limits from XENON1T[57] and LUX-ZEPLIN[59] are given by the dashed and full line, respectively.The grey shaded region corresponds to the neutrino floor.The color code represents the non-SM-like Higgs mass value in GeV.

Figure 6 :
Figure 6: Relic density as a function of the DM mass m A for both II (left) and scenario I (right).The color code represents the non-SM-like Higgs mass value in GeV.

Figure 7 :
Figure 7: Relative uncertainty δ approx in the EW NLO BR(h 125 → AA) due to using the EW correction in the partial decay width for h 125 → AA only, as a function of the mass m S of the non-SM-like Higgs boson for scenario I (red) and II (blue points) in the four possible different renormalization schemes.