Electroweak Corrections to Dark Matter Direct Detection in the Dark Singlet Phase of the N2HDM

Direct detection experiments are the only way to obtain indisputable evidence of the existence of dark matter (DM) in the form of a particle. These experiments have been used to probe many extensions of the Standard Model (SM) that provide DM candidates. Experimental results like the latest ones from XENON1T lead to severe constraints in the parameter space of many of the proposed models. In a simple extension of the SM, the addition of a complex singlet to the SM content, one-loop corrections need to be taken into account because the tree-level cross section is proportional to the DM velocity, and therefore negligible. In this work we study the case of a DM particle with origin in a singlet but in a larger framework of an extension by an extra doublet together with the extra singlet providing the DM candidate. We show that in the region of interest of the present and future direct detection experiments, electroweak corrections are quite stable with a $K$-factor very close to one.


Introduction
The only way to unmistakably identify a dark matter (DM) particle is in direct detection experiments. In the mass region of the so-called Weekly Interacting Massive Particles (WIMPs) the latest and most restrictive constraints were obtained by the XENON1T collaboration [1,2]. In this type of experiments, when DM interacts with XENON it creates light and electric charge and theses signals provide information about the energy and location of the initial collision. Since direct detection experiments play the major role in probing the WIMP region it is important to understand in great detail the DM-nucleon cross sections in the different models. There is a particularly interesting case, the one of the extension of the SM by a complex singlet that leads to a tree-level DM-nucleon cross section proportional to the DM velocity and therefore to a negligible rate [3]. The calculation of the electroweak corrections to DM-nucleon scattering in this model was performed in [4][5][6] and shown to be several orders of magnitude above the tree-level result.
In previous works we have also calculated the electroweak corrections [7,8] in a vector DM model [9]. In this case the tree-level cross section is not negligible and electroweak corrections, in the region not excluded by XENON1T, are quite stable with a K-factor close to 1 (K = σ NLO /σ LO -the ratio of the next-to-leading order to the leading order cross section). In this letter we discuss a scenario where the DM candidate originates from a singlet but now within the larger framework of the Dark Singlet Phase (DSP) [10,11] of the next to 2-Higgs Doublet Model (N2HDM). The first point to note is that in this case there is no tree-level cancelation. Hence, the leading order cross section is not negligible. The main question we would like to answer is if the corrections are still stable and not too large when the parameter space of the visible sector is enlarged which is the case of the DSP os the N2HDM. The DM candidate is singlet-like but the visible sector is now a Z 2 symmetric 2HDM, with a new set of parameters and extra contributions to the electroweak corrections. As we have discussed in great detail all the steps of the calculations in our previous works [6][7][8] and also because there are no major changes in the methodology we will whenever possible refer the reader to those works and will just focus on the differences for the model under study.
The outline of the letter is as follows. In section 2, we will introduce the DSP of the N2HDM together with our notation. Section 3 contains a brief description of the renormalization procedure used in this work. In section 4 we calculate the electroweak corrections to the spinindependent direct detection cross section. In section 5, the results are presented and discussed. Finally, we present our conclusions in section 6.

The Dark Singlet Phase of the N2HDM
The model considered in this work is the DSP of the N2HDM [12][13][14]. The Higgs sector of the SM is extended by one complex SU (2) L doublet with hypercharge +1, and one real SU (2) L singlet with hypercharge 0. We focus on a particular phase of the four possible dark phases of the N2HDM, the DSP, where the singlet field has a vanishing vacuum expectation value (VEV) and does not couple to the SM fields, making it a DM candidate. A detailed discussion of the different dark matter phases of the N2HDM can be found in [10,15]. The Yukawa version of the model is type I meaning that that all quarks and leptons couple to only one of the doublets. The Higgs potential is simplified by requiring invariance under the two Z 2 symmetries, which allows us to write the most general CP-conserving and renormalizable scalar potential invariant under these Z 2 symmetries as with three real mass dimension parameters m 2 11 , m 2 22 , m 2 s and eight real dimensionless parameters λ 1 · · · λ 8 . The symmetry Z (1) 2 is spontaneously broken by the doublet VEV. In the DSP, both doublets acquire VEVs but the singlet VEV vanishes which leaves the symmetry Z (2) 2 unbroken. After electroweak symmetry breaking (EWSB) the doublet and singlet fields can be parametrized in terms of the VEVs v 1 and v 2 , and component fields as where φ + 1 , φ + 2 are complex charged fields, ρ 1 , ρ 2 , ρ s η 1 and η 2 are neutral fields. After EWSB, the CP-even fields ρ 1 and ρ 2 mix and give rise to the CP-even mass eigenstates h 1 and h 2 defined such that m h 1 ≤ m h 2 , and either h 1 or h 2 can be identified with the 125 GeV SM Higgs boson. Similarly η 1 and η 2 mix to give a pseudoscalar mass eigenstate A and the neutral Goldstone boson G 0 . Finally, φ + 1 , φ + 2 mix to give a charged Higgs H + and the charged Golstone boson G + . The singlet field ρ s does not mix with any of the doublet fields, nor does it couple to any SM particles. Moreover, the unbroken Z (2) 2 symmetry gives rise to a dark parity, such that ρ s ≡ χ with mass m χ emerges as a DM candidate in the model.
The mass eigenstates can be expressed in terms of the gauge eigenstates via rotation matrices as follows, where the rotation matrices are parametrized as 1 The VEVs of the two doublets are related to the SM VEV (v ≈ 246.22 GeV) as with v = 2m W /g, where m W is the mass of the W ± boson. Further details on the minimization conditions can be found in [10,11]. The final set of independent input parameters chosen is v , α , tan β , m χ , m h 1 , m h 2 , m A , m H ± , λ 6 , λ 7 , λ 8 .

Renormalization of the Model
In the following, we present the renormalization of the N2HDM DSP in order to calculate the electroweak (EW) corrections to the scattering process of the scalar DM particle χ with a nucleon. We follow the prescription presented in [16] by adapting the renormalization of the unbroken phase of the N2HDM to our scenario. This is done by taking the limit v s = α 2 = α 3 = 0 and by additionally renormalizing the parameters λ 7 and λ 8 , which are the ones that enter our calculation.
The bare input parameters p 0 defined in Eq. (8) are expressed in terms of the renormalized parameters p and their respective counterterms δp as whereas bare fields φ 0 are expressed in terms of the renormalized fields φ and the wave-function renormalization factors (WFRs) Z φ as where Z φ is a matrix in case the fields mix at one-loop. The renormalization conditions define the finite parts of the counterterms. In this work, we will use the on-shell (OS) renormalization scheme to fix the renormalization constants for the masses and fields. The tadpoles are treated in the Fleischer and Jegerlehner (FJ) [17] scheme. A detailed description of the scheme and its consequences for gauge independence can be found in [16,[18][19][20]. In the following sections we just present a brief description of the renormalization of the model sector by sector, giving the expressions required for the renormalization of the direct detection process χp → χp.

Scalar Sector
In the DSP of the N2HDM, after EWSB there are four neutral scalars (two CP-even, one CPodd and the DM candidate χ) and one charged Higgs pair. The OS conditions for the physical Higgs states result in the following mass counterterms, where Σ tad (p 2 ) are the self-energies containing all tadpole topologies. Since the tadpoles are absorbed into the self-energies, explicit tadpole counterterms do not appear in the mass counterterms [16]. The fields are renormalized in terms of the WFR constants where , and the other two cases are obtained by replacing {h 1 , h 2 } by {G, A} and {G ± , H ± }. Finally the field strength renormalization for the DM particle χ is expressed in terms of its self-energy Σ χχ (p 2 ) as

Gauge Sector
The gauge sector of the model is renormalized through OS conditions. The masses, couplings and fields are expressed in terms of their counterterms as e → e + δZ e , g → g + δg , where m W and m Z are the W and Z boson masses, respectively, e is the electric charge and g is the weak SU (2) L coupling. The OS conditions for the masses give rise to the counterterms where the superscript T indicates the transverse part of the self-energy, which also includes tadpole contributions. The counterterm for the electric charge is fixed in the Thomson limit as in the SM and is expressed in terms of the Weinberg angle θ W as Using the above expression we can then fix the counterterm δg as Finally, the WFR constants for the gauge fields are given by

Quark Sector
In the quark sector the OS scheme is applied for each quark. The renormalized quark fields are expressed in terms of their left-and right-handed components, with counterterms for each component as follows with q ∈ {u, d, s, c, b, t}. In order to fix the counterterms we need to define the structure of the quark self-energies, where the self-energy superscripts L, R and S respectively correspond to the left-handed, righthanded and scalar parts of the quark self-energies, and P L,R are the left-and right-handed projectors. Using the above expression the quark WFR constants and mass counterterms in terms of the self-energies containing the tadpole topologies, are defined as

Renormalization of the Mixing Angles
Following the renormalization prescription for mixing angles in the 2HDM, the angles α and β are renormalized as proposed in [16,21,22]. The scheme connects δα and δβ to the off-diagonal WFR constants of the scalar sector. Following again [16] the angle counterterms are The counterterm for δβ can be derived either from the charged sector or the CP-odd sector using the same steps, and therefore we have two possible expressions for δβ given by and 3.5 Renormalization of λ 7 and λ 8 We are left with the parameters λ 7 and λ 8 to complete the renormalization of the model. We will use a process-dependent scheme with the on-shell decays h i → χχ (i = 1, 2). The NLO amplitude A N LO consists of the LO decay amplitude A LO , the vertex corrections A V C and the counterterm amplitude A CT , where the index h i denotes the decaying particle. The renormalization condition is such that we force the LO decay width to be equal to NLO decay width. With the Higgs coupling C χχh i between h i and two DM particles given by this condition gives rise to a system of equations for δλ 7 and δλ 8 such that and this concludes our renormalization programme. We can now proceed to the calculation of the EW corrections.

Electroweak Corrections to the SI Cross Section
The spin-independent (SI) DM-nucleon cross section can be written in terms of an effective coupling, α n , such that α n χ χ n n = iA n = iα n u n u n = i · 2m n α n , where u n u n = 2m n (m n is the nucleon mass) because we assume that the velocity of the DM particle is negligibly small. With this definition the DM-nucleon cross section takes the form where m χ is the DM mass. As the nucleon is a bound state, the DM-nucleon coupling receives contributions both from valence quarks (q = u, d, s) and from the gluons. The SI DMnucleon cross section is calculated using a parton basis, with the operators considered in the non-relativistic limit. Its most general form is given by [23] with the operators which are built with the DM field χ, the quark spinor q and the gluon field strength tensor G a µν ; α s is the strong coupling constant. The quark-DM interaction is encoded in the operator O q S while the gluon-DM interaction is encoded in O g S . Finally, the twist-2 operator O q µν also contributes to the SI cross section. The expectation values of the operators in Eq. (36) are written as [24][25][26] n| m qq q |n ≡ m n f n q , where the nucleon matrix elements f n q and f n g are determined by lattice calculations. Their numerical values are given in Appendix A. The QCD trace anomaly relates the heavy quark Q = b, c, t operators with the gluon field strength tensor [26] corresponding to the Feynman diagram in Fig. 1 and can therefore be determined by first calculating the (heavy) external quark process and then using Eq. (38) to determine the effective gluon interaction. These amplitudes are then used to the extract the Wilson coefficients C g S . The DM nucleon cross section σ n is given at NLO by where, according to the previous discussion, the LO and NLO form factors are given as We will neglect the term proportional to C g S because the matching in Eq. (38) cannot be used at EW NLO. As discussed in [6], taking into account the EW NLO corrections of the gluon contributions would require a proper mixed QCD-EW matching of the QCD trace anomaly. Therefore the NLO contributions considered are (43) Here f i q are the Wilson coefficients from the upper vertex corrections, lower vertex corrections, mediator corrections and the box contributions. The g box q are the box contributions proportional to the second momenta of the quarks q n (2). The values for q n (2) are also given in Appendix A.

Upper vertex corrections
For the extraction of the Wilson coefficients of the upper vertex the one-loop corrections to the coupling χχh i (i = 1, 2) need to be calculated. For this purpose χ is taken on-shell and it is assumed that the momentum transfer goes to zero. Thus the incoming momentum p in of the dark matter particle equals the outgoing momentum p out ≡ p. The NLO amplitude consists of the LO amplitude A LO , the virtual vertex corrections A V C and the counterterm A CT . In the limit of zero momentum transfer the LO amplitudes for the upper vertex topology read (i = 1, 2) with C χχh i given in Eq. (30), the Higgs h i coupling to a quark pair given by and u(ū) denoting the spinor of the (anti-)quark. The counterterm amplitudes read where the counterterm of the coupling is obtained by varying the trilinear coupling C χχh i , Therefore the NLO amplitude at zero momentum transfer is given by

Mediator corrections
To obtain the Wilson coefficient from the mediator corrections the one-loop corrections to the propagator and the corresponding counterterms need to be determined. They can be expressed in terms of the renormalized one-loop propagator (i, j = 1, 2) with the renormalised self-energy matrix with δT denoting the tadpole counterterm matrix and M ≡ m 2 h i δ ij the diagonal mass matrix. The details of the calculation can be found in [6]. Here we just note that the self-energies receive extra contributions relative to the complex singlet extension in [6] because there are new scalars in the loop from the second doublet.

Lower vertex corrections
The lower vertex corrections are also calculated exactly as in [6]. The difference is again the contribution of the new scalar particles in the loop. However from the point of view of the renormalization procedure nothing changes. A very detailed discussion on the different problems arising in this calculation is presented in our Ref. [8]. Of particular importance is the treatment of infrared divergences and the discussion on the heavy quark contributions to the process.

Box diagrams
Finally, the calculation of the box corrections is also discussed in detail in our previous works, Refs. [6][7][8] and again was shown to be one order of magnitude smaller than the main contribution.

Results and discussion
One of the Higgs bosons, either h 1 or h 2 , is the SM-like Higgs boson with a mass of 125.09 GeV [27]. The other CP-even Higgs can be lighter or heavier than the SM-like Higgs boson. The points presented in the scatter plots were generated using ScannerS [28,29] where the most relevant experimental and theoretical constraints were taken into account. ScannerS checks if the potential is bounded from below, that there is a global minimum and that perturbative unitarity holds. Agreement with the electroweak precision measurements at the 2σ level is enforced using the S, T, U [30,31] parameters. Collider bounds from Tevatron, LEP and LHC are encoded in HiggsBounds 5.6.0 [32] and HiggsSignals 2.3.1 [33]. We ask for a 95% confidence level agreement using the exclusion limits for all available searches for non-standard Higgs bosons, including Higgs invisible decays. Branching ratios are calculated using AnyHdecay 1.1.0 [29]. The code includes the Higgs decay widths for the N2HDM N2HDECAY [34], with stateof-the art higher-order QCD corrections. The code N2HDECAY is based on the implementation of the N2HDM in HDECAY [35,36]. All EW radiative corrections in HDECAY are turned off for consistency.
The DM relic abundance is calculated using MicrOMEGAs [37], and a bound on its value is set by the current experimental result (Ωh 2 ) obs DM = 0.1186±0.002 from the Planck Collaboration [38]. We require the calculated relic abundance to be equal or below its experimental central value plus 2σ, that is, we allow the DM not to saturate the relic density and therefore define a DM fraction where (Ωh 2 ) χ stands for the calculated DM relic abundance in our model. As for direct detection, the XENON1T [1,2] experiment provides the most stringent upper bound on the spinindependent DM nucleon scattering. The ranges of the input parameters for the scan are shown in Table 1. m Φ denotes the masses of Φ = h i , A, H ± , where h i is the non-125 GeV Higgs. Note that the constraints are applied and therefore the allowed parameter space will be a small fraction of the initial space. In Fig. 2 we present a scatter plot of the K-factor, defined as K = σ N LO/ σ LO , as a function of the σ LO DM-nucleon spin-independent cross section. In the left panel we show the behaviour with λ 8 in the colour bar while in the right panel the DM mass is shown in the colour bar. The XENON1T experiment sets an upper limit on the cross section of 10 −45 cm 2 valid for any value of the DM mass (the limit is stronger for smaller DM masses). Therefore, in the region of interest it is clear that the correction are very stable with the bulk of the points just below 0.94. The largest value of the correction yields a K = 1.12 and the lowest value of K is just below 0.75. Hence the correction are stable and there is in general a slight decrease of the cross section at NLO. We have also checked that in the region of interest none of the free parameters play a special role in the K-factor values. In the left panel we show the behaviour with λ8 in the colour bar while in the right panel the DM mass is now shown in the colour bar.
In Fig. 3 we now show the DM-nucleon cross section including the correction factor f χχ , at NLO (left) and LO (right) compared to the XENON1T limit (orange line), as a function of the DM mass. The points shown are such that they are all below the XENON1T line at NLO as can be seen in the left plot. In the right plot we show that some of the points would not comply with XENON1T limit if calculated at LO.

Conclusions
We have calculated the spin-independent DM-nucleon scattering cross section for the DSP of the N2HDM including higher-order corrections. One of the main goals of this work was to check if the parameter space of the model compatible with the most important theoretical and experimental constraints would give rise to large and/or unstable higher-order corrections. We found that in this model the corrections are stable with a K-factor close to one. The reason behind this result is that the main corrections come from the upper vertex, mediator and lower vertex, where the diagrams are similar to the ones in the CxSM. In the CxSM, although the LO cross section turns out to be negligible due to a peculiar Feynman diagram cancellation, the NLO corrections are quite stable as discussed in [4][5][6]. The same is true for the vector dark matter (VDM) model discussed in [7,8]. For the VDM the LO cross section is not negligible and the K-factor is quite stable and close to one, except for large values of the dark gauge coupling.
Still, in the model discussed in this work, there are new particles in the loops, like for instance in the case of the CP-even self-energies which contribute to the mediator correction, one could in principle expect sizeable corrections which is not the case. The masses and couplings in this model are already very constrained by the LHC results but a such a stable result could not be anticipated.
From the phenomenological point of view the overall conclusions are the following. The NLO corrections can increase the LO results to values where the XENON1T experiment becomes sensitive to the model, or to values where the model is even excluded due to cross sections values above the XENON1T limit. But the reverse is also true even if not as common. Parameter points that might be rejected at LO may render the model viable when NLO corrections are included. We conclude that as a first approximation the LO cross section is a very good approximation but if a DM candidate is detected NLO corrections should be taken into account in order to either validate or exclude the model.