One-loop Corrections to the Higgs Boson Invisible Decay in the Dark Doublet Phase of the N2HDM

The Higgs invisible decay width may soon become a powerful tool to probe extensions of the Standard Model with dark matter candidates at the Large Hadron Collider. In this work, we calculate the next-to-leading order (NLO) electroweak corrections to the 125 GeV Higgs decay width into two dark matter particles. The model is the next-to-minimal 2-Higgs-doublet model (N2HDM) in the dark doublet phase, that is, only one doublet and the singlet acquire vacuum expectation values. We show that the present measurement of the Higgs invisible branching ratio, BR$(H \to$ invisible $<0.11$), does not lead to constraints on the parameter space of the model at leading order. This is due to the very precise measurements of the Higgs couplings but could change in the near future. Furthermore, if NLO corrections are required not to be unphysically large, no limits on the parameter space can be extracted from the NLO results.


Introduction
Ever since the Higgs boson was discovered at the Large Hadron Collider (LHC) by the ATLAS [1] and CMS [2] collaborations, the measurement of the Higgs couplings to the remaining Standard Model (SM) particles became a powerful tool in constraining the parameter space of extensions of the SM. Another important ingredient when building extensions of the SM with dark matter (DM) candidates is the measurement of the invisible Higgs branching ratio. Very recently, a new result by the ATLAS collaboration combining 139 fb −1 of data at √ s = 13 TeV with the results obtained at √ s = 7 and 8 TeV was published. The observed upper limit on the SM-like Higgs (H SM ) to invisibles branching ratio (BR) is 0.11 [3], which is an improvement from the previous result with an invisible BR above 0.2. The results of the Higgs coupling measurements together with those of the invisible Higgs decay are our best tools at colliders to constrain extensions of the scalar sector of the SM with DM candidates.
We will focus on a specific phase of the next-to-minimal 2-Higgs-doublet model (N2HDM) with a scalar sector consisting of two complex doublets and one real singlet. Only one of the doublets and the singlet acquire vacuum expectation values (VEVs) and we end up with two possible DM candidates. This particular phase of the N2HDM is known as the dark doublet phase (DDP). The different phases of the N2HDM are described in detail in [4].
Our analysis will be performed by first imposing the most relevant theoretical and experimental constraints on the model. We then calculate the next-to-leading (NLO) electroweak corrections to the invisible decay of the SM-like Higgs boson, that is, the Higgs decaying into two DM candidates. The results will be presented for all allowed parameter space points, which will enable us to understand if NLO corrections can help to constrain the parameter space of the N2HDM. The NLO BR of Higgs to invisibles could be larger than the experimentally measured value for some regions of the parameter space.
As shown in a recent work [4], the constraints coming from the Higgs couplings to fermions and gauge bosons are enough to indirectly constrain the BR of the Higgs decay into invisibles to be below 0.1 in the N2HDM (DDP phase). So until recently, the Higgs BR to invisible was not a meaningful experimental result to constrain the parameter space. However, the new measurement by ATLAS, reaching now 0.11, is exactly at the frontier between the indirect bound coming from Higgs couplings and the direct one coming from the invisibles Higgs BR. So it is extremely timely to calculate the electroweak (EW) NLO corrections to the Higgs invisible BR.
The outline of the paper is as follows. In section 2, we will introduce the DDP phase of the N2HDM together with our notation. Section 3 is dedicated to the description of the different renormalization schemes used in this work. In section 4 we discuss the expressions of the Higgs invisible decay at leading order (LO) and at NLO. In section 5, the results are presented and discussed. Our conclusions are collected in section 6. There are two appendices where we discuss details of the renormalization procedure. as Z (1) We require the Z 2 symmetries to be exact, meaning that no soft breaking terms are introduced, and therefore the Higgs potential of the N2HDM is given by [4][5][6][7][8] where all parameters can be set real by rephasing Φ 1 or Φ 2 . In the N2HDM, there are four different minima, which break the SU (2) × U (1) Y symmetry into U EM (1), depending on the vacuum expectation values for the doublet fields and the singlet field, i.e. Φ 1 , Φ 2 and Φ S , respectively. The possible patterns are broken phase (BP): dark doublet phase (DDP): dark singlet phase (DSP): full dark phase (FDP): Φ 1 = 0 , Φ 2 = 0 , Φ S = 0 .
In the DDP, the components of the Higgs fields can be parameterized as where v = 246 GeV is the electroweak VEV and v S is the VEV of the singlet field. The doublet field Φ 1 corresponds the SM Higgs doublet, which contains the Nambu-Goldstone bosons G + and G 0 . Due to the unbroken Z 2 symmetry, the four dark scalars, H D , A D and H ± D do not mix, i.e., they are physical states. The lightest neutral dark scalar, which can be either H D or A D , is the DM candidate. On the other hand, the two CP-even Higgs fields ρ 1 and ρ S mix with each other. Together with the CP-even dark scalar H D , the mass eigenstates for the CP-even Higgs bosons can be expressed through a rotation matrix with the mixing angle α as, where by convention, we take m H 1 < m H 2 , and where we have introduced the short-hand notations c α ≡ cos α and s α ≡ sin α. Either H 1 or H 2 can be identified as the SM-like Higgs boson (H SM ) with a mass of 125 GeV. For later convenience, we define the rotation matrix as R, so that Eq. (9) can be rewritten by H i = R ij ρ j (i, j = 1, 2, 3), defining ρ 2 = H D and ρ 3 = ρ S .
The masses of the physical states can be written as Using these mass formulae together with the mixing angle α and the stationary conditions for Φ 1 and Φ S , i.e., the original parameters of the potential can be replaced by a new set to be used as input.
Together with the electroweak and singlet VEVs, we choose as our input the following 13 parameters, We assign Z 2 -even and Z 2 -even parity to the remaining SM fields and consequently only the Higgs doublet Φ 1 has Yukawa interactions. This in turn means that the Yukawa couplings are just the SM ones and that the dark scalars do not couple to the SM fermions. On the other hand, due to the kinetic term for Φ 2 , dark scalar-dark scalar-gauge boson type of vertices are allowed while the dark scalar-gauge boson-gauge boson type of vertices are forbidden by the Z (1) 2 symmetry. The Feynman rules for the vertices including two dark scalars and a gauge boson are written in Ref [4]. In particular, the trilinear scalar couplings that are relevant for the calculation of the invisible decay of the Higgs boson are given by (i = 1, 2) We note that in the case of cos α = 1 and λ 8 = 0 the expression of λ H i H D H D is exactly the IDM one. The other trilinear scalar couplings are also given in Ref. [4].

Renormalization
In this section, we discuss the renormalization scheme used in the calculation of the one-loop corrections to the Higgs boson H i (i = 1, 2) decay into a pair of DM particles, which we assume to be H D , unless otherwise stated, H i → H D H D . Thus, this section focuses on the renormalization of the scalar and gauge sectors. The renormalization of the fermion sector as well as any treatment of infrared divergence is not necessary for this particular process.
We perform the renormalization of the Higgs sector in the DDP of the N2HDM according to the procedure presented in Ref. [10] for the 2HDM and in Ref. [11] for the broken phase of the N2HDM. Although most of the parameters in the Higgs sector of the DDP are common with those of the broken phase, we describe the renormalization of all parameters in order to make the paper self-contained.

Gauge sector
The renormalization of all parameters and fields in the gauge sector is done using the on-shell (OS) scheme following Ref. [12]. As the three independent parameters in this sector, we choose the masses of weak gauge bosons and the electric charge, i.e., m W , m Z , and e, respectively. These parameters are shifted as Moreover, the bare fields for the gauge bosons in the mass basis are replaced by the renormalized ones as The OS conditions for these gauge fields are defined as where Σ tad,T W W and Σ tad,T ZZ denote the transverse part of the self-energies of the W and Z bosons. These contain the tadpole contributions due to our renormalization scheme choice. No "tad" superscript means that there is no contribution from the tadpole diagrams. The different tadpole schemes will be described below. The counterterm for the electric charge is determined from γeē in the Thomson limit and can be expressed as a function of the self-energies as This counterterm contains large logarithmic corrections arising from the small fermion masses, log m 2 f (f = t). We use the "G µ scheme" [13] in order to improve the perturbative behaviour. In this scheme, a large universal part of the O(α) corrections is absorbed in the leading order decay width by deriving the electromagnetic coupling constant α = e 2 /(4π) from the Fermi constant, G µ , as This allows us to take into account the running of the electromagnetic coupling constant α(Q 2 ), from Q 2 = 0 to the electroweak scale. In order to avoid double counting, the corrections that are absorbed in the LO decay width by using α Gµ have to be subtracted from the explicit O(α) corrections. This is achieved by subtracting the weak corrections to the muon decay, ∆r [12,14], from the corrections in the α(0) scheme. Hence, we redefine the charge renormalization constant as where (∆r) 1-loop is the one-loop expression for ∆r given by [12] (∆r Note that through the redefinition Eq. (26), the first term of δZ α(0) e in Eq. (24), which contains the large logarithmic corrections from the light fermion loops, cancels against the corresponding term in (∆r) 1-loop . The counterterms for the other EW parameters can be expressed in terms of those presented above. For example, the SU (2) L gauge coupling, g, is replaced by the tree Thus, the counterterm is given by

Higgs sector
In the Higgs sector, we have a total of 13 free parameters, given in Eq. (16), considering the two tadpoles T Φ and T S . We have to renormalize the scalar fields in the mass basis, H 1 , H 2 , H D , A D and H ± D . The counterterms are introduced via the shift of the input parameters, i.e, the masses of the scalar bosons, the mixing angle α of the CP-even Higgs bosons and the remaining original potential parameters that appear in the vertices of the processes under study, λ 8 and m 2 22 , where Φ denotes H 1 , H 2 , H D , A D , and H ± D . There is no need to renormalize λ 2 for this particular process. Apart from the tadpoles, the remaining two parameters are the VEVs. The electroweak VEV v is fixed by the W mass and the renormalization of v S will be discussed later.
The tadpole renormalization can be performed in different ways and we will discuss two approaches. These are designated by Standard Tadpole Scheme (STS) and Alternative Tadpole Scheme (ATS). The latter was originally proposed by Fleischer and Jegerlehner, in Ref. [15], for the SM. The ATS was also discussed in detail for the CP-conserving 2HDM in Ref. [10] and for the broken phase of the N2HDM in Ref. [11]. We will just briefly review the two schemes for completeness.
In the STS, the tree level tadpoles are replaced by and are chosen as the renormalization parameters. On the other hand, in the ATS, the VEVs are the renormalization parameters and are shifted as We use the ATS, which will now be explained in more detail. The reason to use this scheme is that, as shown by Fleischer and Jegerlehner, all renormalized parameters are gauge independent except for the wave function renormalization constants (or any parameter that depends on the wave function renormalizaton constants as is the case of the angle α in particular schemes).
Before moving to the discussion of the tadpole renormalization, we define the wave function renormalization constants of the scalar fields. The bare fields are replaced by the renormalized ones through Note that, for Eq. (32), the exact Z 2 symmetry ensures that the (3-k) and (k-3) components (k = 1, 2) are zero.

Tadpoles
In the ATS, the renormalized VEVs, which correspond to minima of the Higgs potential at loop-level, are regarded as the tree-level VEVs, namely, one imposes Using the tadpole conditions, one can derive expressions for δv and δv S , where the first term, T tree i , is zero from the stationary condition at the tree-level and the second term, f (δv, δv S ), denotes contributions from δv and δv S , which can be extracted by inserting Eq. (31) into the tree-level tadpole conditions. From Eq. (35) one obtains the following expressions for the VEV counterterms where R(α) denotes the 2 × 2 non-diagonal part of R (see Eq. (9)), and the one-loop tadpoles in the mass basis are given by T loop The left-hand side corresponds to the VEV counterterms δv H 1 and δv H 2 in the mass basis and the right-hand side coincides with the tadpole diagrams multiplied by the propagator for the Higgs bosons at zero momentum transfer, i.e. T loop . Therefore, Eq. (36) shows that δv H i can be regarded as the connected tadpole diagrams for H i . Once the counterterms for the VEVs are fixed, the shift is performed in all VEV terms in the Lagrangian. Hence in the ATS, one needs to insert tadpole diagrams in all amplitudes for which the original vertices contain one of the VEVs in addition to the usual one-particle irreducible diagrams. This general consequence is shown by focusing on specific amplitudes in Ref. [10] for the 2HDM and in Ref. [11] for the N2HDM.
Another important feature of the ATS is that, because the renormalized VEV is identified with the tree level VEV, the VEVs still have to be renormalized. For the EW VEV, the renormalized parameter is given by and the tree-level parameters g and m W are shifted as We have defined the shift of the tree level parameters related to the EW VEV as ∆v, which has no relation with δv. The same discussion holds for the singlet VEV v S . Once v S is related with some measurable quantity, a similar relation with Eq.(38) must exist, even if a physical process has to be used, and then ∆v S has to be introduced.

Mass and Wave Function Renormalization
The counterterms for the masses and the wave function renormalization constants (WFRCs) are determined by imposing the on-shell conditions for each scalar field. This yields the mass counterterms and the WFRCs where, to reiterate, Σ tad stands for the self-energy containing one particle irreducible (1PI) diagrams and tadpole contributions.

Mixing angle α
The renormalization of the mixing angle, α, requires special treatment since the gauge dependence of δα could result in a gauge-dependent physical process [16][17][18][19][20][21][22][23][24][25][26]. A gauge-independent amplitude can be obtained by starting with a gauge-independent definition of δα. One possible solution to avoid this gauge dependence is to apply the ATS for the renormalization of the tadpoles and to make use of the pinch technique [27,28], while keeping the on-shell renormalization for the mixing angle. This is the procedure that we adopt throughout this paper.
The expression for δα in the OS scheme can be derived by relating quantities in the gauge basis to the corresponding ones in the physical (mass) basis. This procedure is described in detail in [10,29,30]. Following [10] leads to the following expression for δα after adding the pinched terms, where Σ H 1 H 2 stands for pinched contributions to the H 1 − H 2 mixing self-energy, which remove the gauge-dependent part coming from the first two terms. They can be extracted from the expressions obtained for the broken phase of the N2HDM (see [11]) as The expression was obtained with the replacement (β − α 1 , α 2 , α 3 ) → (α, 0, 0) and the fact that the (A D , Z) and (H ± D , W ) loop contributions do not appear in our calculation due to the existence of an exact Z (1) 2 symmetry.

22
The counterterms for the quartic coupling δλ 8 and the invariant mass for the dark scalars δm 2 22 cannot be renormalized using OS conditions for the Higgs states. Hence we will renormalize these parameters using three different schemes: the MS scheme, a process-dependent scheme and one derivation of the latter that consists of taking the external momenta to be zero instead of taking them on-mass-shell.
In the MS scheme, the analytic expressions for the counterterms can be extracted from the beta functions at one-loop, yielding where ∆ div denotes the UV-divergent part, i.e., ∆ div = 1/ − γ E + log(4π). γ E is the Euler-Mascheroni constant and 1/ is the UV pole in dimensional regularisation. The beta functions are given in terms of the original potential parameters by β (1) where the Clebsch-Gordan coefficient C is given by C 2 = 5/3 and g 1 and g 2 denote the U (1) Y and SU (2) L gauge couplings, respectively. These expressions were derived using SARAH-4.14.2 [31][32][33][34][35].
For the process-dependent scheme, one can fix the counterterms δλ 8 and δm 2 22 , required in the one-loop decays H i → H D H D (i = 1, 2), by making use of the Higgs boson decays into a pair of dark CP-odd scalars. We choose as renormalization condition that the decay width for The counterterms δλ 8 and δm 2 22 defined by these conditions contain not only UV-divergent parts but also finite terms. The detailed explanation on how the counterterms are computed is given in Appendix A.
The process-dependent scheme takes all particles to be on-shell because it uses a physical process. This means, however, that the renormalization conditions Eq. (48) can only be used if the decay processes H i → A D A D are kinematically allowed. There is a way to circumvent this problem by not taking the particles on-shell.
The renormalization conditions Eq. (48) can be written as, because the tree-level amplitude is just a real constant. If instead we choose to use the same condition but with all external momenta equal to zero, we will not be restricting the parameter space of the model that can be probed. The third renormalization scheme is therefore defined by while using exactly the same two processes that were used for the process-dependent scheme.
Note that the problem in the on-shell case is related to the calculation of C 0 loop functions in forbidden kinematical regions [36]. We will refer to the two schemes as OS process-dependent and zero external momenta (ZEM) process-dependent in the following.

Determination of ∆v S
The quantity ∆v S , which is introduced by a similar relation to the one for the SM in Eq.
as given in Appendix B. We checked that by using Eq. (51) the one-loop amplitude for H 2 → H D H D is also UV-finite.

The Invisible Higgs Boson Decays at NLO EW
In this section, we calculate the one-loop corrections to the partial decay widths of the Higgs bosons decaying into a pair of DM particles. Hereafter, we regard the CP-even dark scalar H D as the DM candidate unless otherwise specified. We will therefore present the analytic expressions for the decay widths of H i → H D H D (i = 1, 2) at NLO.
The decay rate for H i → H D H D at LO is given by where the scalar coupling λ H i H D H D is given in Eq. (17). The 1PI diagrams contributing to the one-loop amplitude for the process H i → H D H D contain UV divergences that are absorbed by introducing the corresponding counterterms in the amplitude. Shifting all parameters in Eq. (17), we obtain the counterterms for the λ H i H D H D couplings, The counterterms δR i1 and δR i3 are those of the 3 × 3 mixing matrix for the neutral Higgs bosons, R (Eq. (9)). For instance, when i = 1, we obtain δR 13 = δs α = c α δα .
As previously discussed we have three options for the counterterms δλ 8 and δm 2 22 , namely, the MS scheme, the OS process-dependent scheme and the ZEM process-dependent scheme. The corresponding conditions and counterterms are given in Eq. (45), Eq. (49) and Eq.(50), respectively, together with Appendix A . In addition, performing the shift of the fields present in the tree-level Lagrangian for the H i H D H D vertices, we obtain Therefore, the counterterms for the one-loop amplitudes for H i → H D H D are given by With this counterterm, the renormalized one-loop amplitude for H i → H D H D is expressed as We can finally write the decay width at NLO as where the one-loop corrections are written as

Numerical Results
In this section, we analyze the impact of the one-loop corrections to the invisible decay of the SMlike Higgs boson. In Sec. 5.1, we start by discussing the behavior of the corrections to the partial decay width of H 1 → H D H D with the most relevant parameters of the model, namely the trilinear tree-level coupling of the 125 GeV Higgs with the two DM candidates and the mass difference between the two neutral scalars from the dark sector. We then perform a scan in the allowed parameter space, in Sec. 5.2, and present the results for the branching ratios for the invisible decays of the Higgs bosons H i (i = 1, 2). The calculations of the NLO corrections were performed using FeynRules 2.3.35 [37][38][39], FeynArts 3.10 [40,41] and FeynCalc 9.3.1 [42,43]. The same calculations were independently done using SARAH 4.14.2 [31][32][33][34][35], FeynArts 3.10 and FormCalc 9.8 [44]. Loop integrals were computed using LoopTools [44,45]. We have checked numerically that the results obtained with the two different procedures were in agreement.

Impact of the One-Loop Corrections on the Decay Rates
We start by analysing our model in the Inert Doublet Model (IDM) [9], which can be obtained as a limit of the DDP of the N2HDM by setting (in this order) λ 8 = 0, α = 0, v S → ∞. The parameters chosen take into account the bounds for the IDM presented in [46]. We will present numerical results for the one-loop corrected partial decay widths of the CP-even Higgs bosons to dark matter particles. For a particular choice of parameters, we will compare the three renormalization schemes for δλ 8 and δm 2 22 showing the one-loop corrections in the MS scheme, in the on-shell process-dependent scheme and in the ZEM process-dependent scheme. For this comparison we are not taking into account any theoretical constraints yet. The goal is to understand the theoretical behavior of the one-loop corrections. Numerical results considering the theoretical constraints as well as experimental constraints will be presented in the next section.
Among the 13 free parameters given in Eq. (16), the EW VEV is fixed by the input parameters from the gauge sector, which we choose to be m Z , m W and α Gµ . Using α Gµ allows us to resum large logarithms from the light fermion contributions. In this sense, our result for the decay width at LO does not correspond to the pure tree-level result as a large universal part of the O(α) corrections is already included at LO. The remaining 10 parameters, besides the two tadpoles T Φ = T S = 0, are set as follows: H 1 is the SM-like Higgs boson with m H 1 = 125.09 GeV, and the mass of the heavier Higgs boson H 2 is fixed as while the remaining mass parameters m A D and m 2 22 can be either scanned over or fixed in the following plots. We assume m A D > m H D , meaning that the dark scalar H D is the DM candidate. As previously stated we choose for λ 8 , α and v S , in that order, which is equivalent to take m 2 S , λ 6 , λ 7 and λ 8 equal to zero in the scalar potential in Eq. (3). This is in turn equivalent to the IDM potential. Hence, Eq. (63) gives the IDM limit in the DDP phase of the N2HDM and we should recover the IDM results. When the MS scheme is used to calculate δλ 8 and δm 2 22 , the one-loop amplitude for H 1 → H D H D depends on the renormalization scale µ and we set it as µ 2 = m 2 H 1 . In Fig. 1, we show the correlation between the tree-level coupling H 1 H D H D and the decay width for the corresponding process H 1 → H D H D at LO and NLO and for two different charged Higgs masses, m H ± D = 100 GeV and 500 GeV. In this plot, we set the mass of the CP-odd dark scalar to m A D = 62 GeV and vary m 2 22 in a range that forces the tree-level coupling to be |λ H 1 H D H D /(2v)| < 0.05 [46]. The upper bound for |λ H 1 H D H D /(2v)| corresponds to the current bounds for direct detection of DM from XENON1T [47]. From the left panel, in which the MS scheme results are shown, one can see a parabolic behaviour for the decay width at both LO and NLO, with the width vanishing at λ H 1 H D H D /(2v) = 0. The most important feature is that the NLO corrections strongly depend on the value of m H ± D and can be very large even for relatively small λ H 1 H D H D if the mass of the dark charged scalars is large (m H ± D = 500 GeV in the plot). In the right panel of Fig. 1, results for the two process-dependent schemes are shown. The behaviour of the results at NLO for m H ± D = 100 GeV is similar for all three renormalization schemes. In particular, we have confirmed that the result for the ZEM process-dependent scheme almost coincides with that for the OS process-dependent scheme. However, the NLO corrections for the decay width at m H ± D = 500 GeV are quite moderate in both process-dependent schemes, in contrast with the MS scheme.
In Fig. 2, we show the relative size of the NLO corrections  GeV is used because we want to compare the renormalization schemes in a region where they all can be applied. The SM-like Higgs decay into a pair of CP-odd scalars, H 1 → A D A D , has to be kinematically allowed so that Eq. (48) is applicable. The NLO corrections for the MS scheme, with m H ± D fixed to 100 GeV, are almost constant, i.e., they do not depend on the mass difference between the two dark neutral scalars. Nonetheless, as we have seen before, they do depend quite strongly on the charged Higgs mass. In both process-dependent schemes, the NLO corrections strongly depend on the mass difference, ∆m, but also on the value of the charged Higgs mass. For a low value of the charged Higgs mass, m H ± D = 100 GeV, the maximum value of the relative correction for the process-dependent schemes is ∆ NLO =4% at ∆m = 12 GeV, while the minimum is ∆ NLO ∼ 0%. These corrections increase for larger charged Higgs mass.
Considering m H ± D = 500 GeV, the value of Γ NLO /Γ LO − 1 has a minimum of about 4% (24%) for the OS (ZEM) case for ∆m = 0 and a maximum of about 40% (57%) for the OS (ZEM) case for ∆m = 12 GeV. This behaviour can be understood from the fact that there is a significant number of terms in M 1-loop H 1 →H D H D that are proportional to ∆m and, consequently, they have a large impact on the one-loop result. The latter is also proportional to the charged Higgs mass and, therefore, sizable corrections are found for m H ± D = 500 GeV. In the MS scheme, the NLO corrections for m H ± = 500 GeV are well above 100% in the entire mass range, ∆m ∈ [0, 12] GeV, and not shown in the plot.

Scan Analysis for the Branching Ratios
In this section, we will perform a scan over the allowed parameter space of the model. This will enable us to understand the overall behavior of the NLO corrections to the SM-like Higgs decays into a pair of DM particles. The evaluation of the branching ratio is performed using N2HDECAY [48] which is an extension of the original code HDECAY [49,50] to the N2HDM. The program computes the branching ratios and the total decay widths of the neutral Higgs bosons H 1 and H 2 , including the state-of-the art QCD corrections. Using the value of the partial widths evaluated by N2HDECAY, Γ N2HDECAY , we evaluate the branching ratios for H i → H D H D with the NLO EW corrections as where the correction factor δ EW H i →XX is defined by In Eq. (65), the total decay width is separated into the decays into the SM particles, Γ N2HDECAY H i →SM , and the decay into a pair of the scalar bosons Γ H i →ΦΦ , defined as where we include our computed EW corrections to the decays into neutral dark bosons, Here we highlight that, in the process-dependent scheme, δ EW disappears because of the renormalization condition Eq. (48).
We consider two different scenarios in our scan. In scenario 1, the lighter Higgs boson H 1 is identified as the SM-like Higgs boson and the other CP-even Higgs boson H 2 is heavier than the SM-like Higgs boson. In scenario 2, H 2 is the SM-like Higgs boson and the other Higgs boson H 1 is lighter than the SM-like Higgs boson. In both scenarios, the dark scalar H D is the DM candidate. The scan is performed for the two scenarios to examine the impact of the NLO corrections in the allowed parameter space. We use for both scenarios the following ranges for the parameters, 65 GeV < m H ± D < 1500 GeV , 10 −3 GeV 2 < m 2 22 < 5 · 10 5 GeV 2 , 1 GeV < v S < 5000 GeV , −π/2 < α < π/2 , We chose m H ± D to be above 65 GeV to prevent the SM-like Higgs boson decay into a pair of charged Higgs particles. Additionally, λ 2 is set positive due to the boundedness from below (BFB) conditions, see [8] for details on BFB of the N2HDM.
Since we focus on the case where H D is the DM particle, we assume m H D < m A D . Using ScannerS [51,52], we generate input parameter points that pass the most relevant theoretical and experimental constraints. For the theoretical constraints [4,8], ScannerS evaluates perturbative unitarity, boundedness from below and vacuum stability. The following experimental constraints are taken into account: electroweak precision data, Higgs measurements, Higgs exclusion limits, and DM constraints. These constraints are included in ScannerS via the interface with other high energy physics codes: HiggsBounds-5 [53] for the Higgs searches and HiggsSignals-2 [54] for the constraints of the SM-like Higgs boson measurements. For the DM constraints, the relic abundance and the nucleon-DM cross section for direct detection are calculated by MicroOMEGAs-5.2.4 [55][56][57]. The DM relic abundance has to be below the value measured by the Plank experiment [58] and the DM-nucleon cross section has to be within the bounds imposed by the XENON1T [47] results. All points presented in the plots have passed all the above constraints. In Fig. 3 we show two projections of the allowed parameter space in the planes (λ 8 , m 2 22 ) (left) and (sin α, m other ) (right), where m other is the mass of the non-SM-like Higgs boson. The red points are for scenario 1 and the blue points are for scenario 2. There are no particularly important features in the parameters λ 8 and m 2 22 that probe the dark sector as expected, except for theoretical constraints that limit the quartic couplings. As for sin α, due to the very SM-like behaviour of the discovered Higgs boson, sin α is either close to zero or close to ±1, depending on the considered scenario.
In Fig. 4, we show the correlation between the BR(H i → H D H D ) calculated at LO and at NLO in scenario 1 (left panel) and in scenario 2 (right panel). The red and blue points correspond to the calculations in the MS scheme and in the OS process-dependent scheme, respectively. This sample has points with m A D < 125/2 GeV. The first important thing to note is that in both scenarios the LO BR is always below 10%. The main reason for this to happen is the very precise measurements of the Higgs couplings to SM particles which indirectly limit the Higgs coupling to new particles.
The NLO corrections have a very different behaviour in the two renormalization schemes presented. For the MS scheme, the NLO corrections are not reliable with NLO BRs reaching 100% in both scenarios. Conversely, the OS process-dependent scheme is better behaved. This behaviour, in the OS scheme, can be traced back to the suppression of the NLO corrections by the mass difference between H D and A D , as explained in Sec. 5.1. In our analysis, the mass difference is in the range 0 GeV ∆m 6 GeV which leads to small corrections in the OS process-dependent scheme.
In Fig. 5, we show results for the ZEM process-dependent scheme. Again we display the correlation between the branching ratios at NLO and at LO in scenario 1 (left panel) and scenario 2 (right panel). We show results for two different samples of points, all calculated in the ZEM scheme. The grey points correspond to the previous sample where m A D < 125/2 GeV, while the blue points correspond to a range that is only allowed in the ZEM scheme, 125/2 GeV< m A D < 1500 GeV. The points for which m A D < 125/2 have an overall similar behaviour as the ones for the OS scheme, in the sense that the NLO BRs are all below 0.1. However, one can see that the corrections are much larger, even for this sample. When we look at the blue points the picture changes radically. This clearly shows that when the mass difference between H D and A D is large the corrections become unstable.
In order to understand to what extend these corrections depend on the renormalization schemes, we show, in  izontal lines corresponds to BR NLO /BR LO − 1 = ±50%. The plots clearly show that the OS process-dependent scheme is more stable with most corrections between -50% and 50%. In any case, the corrections in this scheme can still go up to 480%. As the corrections above 100% only occur for small values of the LO BRs, the NLO values of the BRs are still well below the experimental bound. The other two schemes are less stable and this is particularly true for the MS scheme.
A clearer picture of the results for the NLO corrections that can be trusted in terms of perturbation theory can be achieved by considering only the points for which the corrections are below 100%. In Fig. 7, we present the correlation between the branching ratios at NLO and at LO in scenario 1 (left panel) and scenario 2 (right panel), respectively. All points presented have NLO corrections below 100% and all points with NLO corrections above 100% were discarded. We conclude that the surviving points are all still below the current experimental limit for the Higgs invisible BR apart from a few grey points. One should keep in mind the theoretical uncertainties due to missing higher-order corrections. Additionally, the other decay channels do not include electroweak corrections, which has an impact on the branching ratios. Given these caveats, the points can still be considered compatible with the experimental results. 1 Therefore, no constraints on the parameter space come from including the NLO corrections. However, as the limit on the Higgs invisible BR improves, there are now ranges of allowed values for the NLO corrections that will certainly lead to constraints on the parameter space.
We end this section with a comment about the case where A D is the DM candidate. We have also calculated the NLO corrections to Higgs boson invisible decays in the case that A D is

Conclusions
In this work we have calculated the EW NLO corrections to the branching ratio of the SM-like Higgs boson invisible decay in the DDP of the N2HDM. We have analysed two different scenarios, one where the SM-like Higgs is the lighter of the visible two CP-even scalars and one where it is the heavier. There are, however, no significant differences between the two scenarios. The model has 13 input parameters from the scalar sector. Masses and wave function renormalization constants are renormalized on-shell. The rotation angle is renormalized by relating the fields in the gauge basis and in the mass basis and it is ultimately defined with the off-diagonal terms of the wave function renormalization constants. We apply a gauge-independent renormalization scheme based on the alternative tadpole scheme together with the pinch technique. Besides the EW VEV, that is renormalized exactly like in the SM, we still have four parameters left: λ 2 (which does not enter in any of the processes under study), λ 8 , m 2 22 and v S . Regarding λ 8 and m 2 22 the analysis was performed for three different renormalization schemes. The three schemes used for these parameters were the MS, the OS process-dependent and the ZEM process-dependent scheme. Only in the first one does v S need to be renormalized. The most stable scheme is the OS process-dependent one, but it can still lead to corrections above 100 % in some regions of the parameter space. It should be noted that one of the reasons for the OS scheme stability is that the mass difference between the two neutral dark scalars, m A D − m H D , is bounded to be below about 10 GeV. Note that the OS process-dependent scheme needs an allowed on-shell decay of the SM-like Higgs boson to both pairs of dark scalars.
With the LHC run 3 starting soon, the Higgs coupling and the invisible Higgs decay width measurements will become increasingly precise. It is clearly the time to understand what the NLO corrections can tell us about the models with more precise measurements. In fact, these experimental results can be the best, if not the only, tools available to probe the dark sectors postulated as extensions of the SM. One should stress that the parameters λ 8 and m 2 22 are only directly accessible through processes that involve the DM particles. The experimental sensitivity on the invisible decay width is now starting to become comparable to the limits imposed on the parameter space of the model from the coupling measurements.
We have found that the NLO corrections can be extremely large in some regions of the parameter space. Also, as we move to smaller values of the BR(H i → H D H D ), the corrections become larger and larger. This means that the more constrained the BR is the more unstable are the NLO corrections. As a perturbativity criteria, we rejected all points for which the NLO corrections relative to the LO results are above 100%. With this condition, the behaviour of NLO versus LO results is very much along the line BR NLO = BR LO . Still, if the experimental bound of BR(h → invisible) improves, for instance to 10 −2 , the NLO result would vary between ∼ 10 −4 to ∼ 2 × 10 −2 .

Appendix
A Determination of δλ 8 and δm 2 22 in the Process-Dependent Scheme In this appendix, we discuss how the counterterms δλ 8 and δm 2 22 are determined using the process-dependent scheme. As previously discussed, although our starting point is to force the amplitudes at LO and at NLO to be equal, we choose two different approaches. In one approach all particles are on-shell, which is equivalent to say that the condition is set on the actual physical process. In the other approach, the condition is set at the amplitude level taking all external momenta to be zero. The advantage of this second approach is not to curtail the allowed parameter space.
We will now describe in detail the renormalization procedure for the on-shell case and discuss the differences when the external momenta are set to zero at the end of this appendix. The onshell process-dependent renormalization condition is to impose the decay widths for H i → A D A D (i=1,2) calculated at NLO to be equal to the LO result, as expressed in Eq. (48) and repeated here for convenience, This results in two equations where δλ 8 and δm 2 22 are the only unknowns. The remaining renormalization constants are all fixed. Solving this set of equations, we can get expressions for these two counterterms. The renormalization conditions, given in Eq. (72), can be written as where dΦ 2 denotes the two-particle differential phase space volume and M tree/1-loop i ≡ M tree/1-loop H i →A D A D (i = 1, 2). The tree-level amplitude is given by where the scalar trilinear coupling λ H i A D A D is Taking into account that M tree i is a real constant, the renormalization condition simplifies to with the one-loop amplitude expressed as Hence the counterterm amplitudes can be written as Finally we obtain the following set of equations, which give us the expressions for δλ 8 and δm 2 22 . Note that the left-handed side of Eq. (81) corresponds to the linear combinations of δλ 8 and δm 2 22 that also appear in the counterterms for H i → H D H D .
The second process-dependent scheme, where all external momenta are set to zero, also starts from the same set of Eq. (81). The only difference is in the calculation of M 1PI i in which the external momenta are set to zero instead of on-shell. The two schemes are compared in Fig. 8, where we plot the ratio of the NLO corrections of the two process-dependent schemes, the zero external momenta over the on-shell scheme, in per-cent, as a function of the LO branching ratio. The left plot corresponds to the decay H 1 → H D H D while the right one corresponds to H 2 → H D H D . We conclude that the differences can be quite large. In fact, although we have cut the y-axis at 500 % for clarity, there are points where the corrections can go above 10 3 %, which, however, is not the case for the larger values of the LO branching ratios. The important point is that very large corrections only occur for the lower values of the BRs so that the NLO results for the larger values of the BRs are quite similar.

B Derivation of ∆v S
In this appendix, we derive the analytic expressions for ∆v S for the case where λ 8 and m 2 22 are renormalized in the MS scheme. As mentioned before, if these parameters are renormalized via a physical process there in no need to renormalize v S . We stated in Sec. 3.2.5, that δv S is determined such that the remaining UV divergence in the renormalized one-loop amplitude for  Therefore, they do not contribute to M CT H 1 →H D H D , which allows us to write