Confronting dark matter co-annihilation of Inert two Higgs Doublet Model with a compressed mass spectrum

We perform a comprehensive analysis for the light scalar dark matter (DM) in the Inert two Higgs doublet model (i2HDM) with compressed mass spectra, small mass splittings among three $\mathbb{Z}_2$ odd particles---scalar $S$, pseudo-scalar $A$, and charged Higgs $H^\pm$. In such a case, the co-annihilation processes play a significant role to reduce the DM relic density. On the other hand, a co-annihilation rate dominated in the early universe implies that a small annihilation rate is needed to reach a correct DM relic density. This results a tiny coupling $\lambda_S$ between DM and standard model Higgs boson and a negligible DM-nucleon elastic scattering cross section at the tree-level. In this work, we include the one-loop quantum corrections of the DM-nucleon elastic scattering cross section. We found that the DM self-coupling $\lambda_2$ indeed contributes a sizeable one-loop quantum correction and behaves non-trivially for the co-annihilation scenario. Interestingly, the parameter space, which are allowed by the current constraints considered in this study, can predict the DM mass and annihilation cross section at the present compatible with the AMS-02 anti-proton excess. The parameter space can be further probed at the future high luminosity LHC.


I. INTRODUCTION
Roughly a quarter of the Universe is made of Dark Matter (DM), but many experimental results also reveal that DM is weakly or even not interacting with the Standard Model (SM) sector. To understand the particle nature of DM, people have developed several detection methods in the past decades such as DM direct detection (DD), indirect detection (ID), and colliders. Although some of DD and ID analysis have reported anomalies [1][2][3][4][5], DM signal is still absent in the LHC searches [6] so that the DM properties (e.g., spin, mass, and couplings) are still not able to be determined. By considering interactions between DM and SM particles within the detector energy threshold, two possibilities arise from null signal detection. The first possibility is that couplings between DM and SM particles are too weak to be detected under the current sensitivities. Hence, an upgrade design of instrument or longer time of exposure may be needed in order to catch the DM signal [7,8].
Nevertheless, if a very tiny coupling between DM and SM particles would be favored by future measurements, the current strategies for the weakly interacting massive particles (WIMPs) are hard to work [9]. The second possibility of DM escaping from detectors of collider experiments may owe to a compressed mass spectrum between the next lightest dark particle and DM [10]. The next lightest dark particle can be long lived because of a small mass splitting. On the other hand, the signal of compressed mass spectrum is usually vetoed or polluted from soft objects of SM backgrounds at colliders [11]. Interestingly, the interaction couplings between dark sector and SM fields may be not suppressed for the second possibility and they may be testable in the DD or ID searches [12].
The inert two Higgs doublet model (i2HDM) [13][14][15][16] is the simplest spin-zero DM model within the framework of two Higgs doublets, and it can naturally realize the compressed mass spectrum [17]. There are three Z 2 -odd scalar bosons: scalar S, pseudo-scalar A, and charged Higgs H ± . Either S or A can play the role of a DM candidate, but it is hard to distinguish one from the other in the phenomenological point of view [18]. In this paper, we only discuss the scalar S playing the role of the DM candidate. The discrete Z 2 symmetry in i2HDM can be taken as an accident symmetry after the symmetry breaking of a larger continuous symmetry. These residual approximated symmetries of the scalar potential can force the mass spectrum of these exotic scalars to be compressed. Given the fact that the SM Higgs doublet H 1 is Z 2 -even and the second doublet H 2 is Z 2 -odd, one can have m S = m A at the leading order as long as the term λ 5 (H † 1 H 2 ) 2 + h.c. in the scalar potential is vanished or omitted [15]. Similarly, if further removing the term λ 4 (H † 2 H 1 )(H † 2 H 1 ) in the scalar potential, all three of Z 2 -odd scalars are degenerated at the leading order [19]. Hence, the compressed mass spectrum can be naturally realized in the i2HDM.
As the thermal DM scenario, a compressed mass spectrum usually results a sufficient co-annihilation in the early universe because the lightest and the next-lightest Z 2 -odd particles are with considerable number density before freeze-out [10]. Requiring the correct relic density to be in agreement with the PLANCK collaboration [20], the sum of WIMPs annihilation and co-annihilation rates shall be a certain range. If the co-annihilation rate overwhelms the total interaction and the annihilation is inactive, the coupling λ S between DM pair and SM Higgs boson will be very small which may lead to undetectable signals at the current DD and collider searches. In particular, due to the smallness of λ S , the DMnucleon elastic cross section is suppressed at tree level so that some regions of parameter space cannot be probed by the current DDs. As reported in Ref. [21], the elastic scattering cross section can be enhanced at the loop level which non-trivially depends on the size of λ 2 and the mass splittings, ∆ 0 = m A − m S and ∆ ± = m H ± − m S . Comparing with Ref. [21] in which only the Higgs resonance regions are discussed, we further investigate the impact of the loop corrections on the co-annihilation scenario in this work.
In this paper, a global analysis of the compressed mass spectrum scenario within the framework of the i2HDM is performed under the combined constraints from theoretical conditions, collider searches, relic density, XENON1T, and Fermi dSphs gamma ray data.
In order to highlight the role of compressed mass spectrum in the early universe, we further divide the allowed region into two groups: the co-annihilation and mixed scenarios. This classification is based on the DM relic density reduction before freeze-out which is mainly governed by the annihilation or co-annihilation processes. The leading order (LO) and next-leading order (NLO) contribution to the DM-nucleon elastic scattering cross section are also computed and compared in the context of these two scenarios. In some regions of the parameter space in the model particularly for the co-annihilation scenario, the DM spin-independent cross section at the NLO can be significantly enhanced and hence probed by the present XENON1T result. Including this NLO enhancement, one can pin down the parameter space of m S , λ 2 , ∆ 0 and ∆ ± . The surviving parameter space is also compatible with AMS-02 antiproton anomaly and can be further probed at the future DD and the high luminosity LHC (HL-LHC) searches for the compressed mass spectrum.
The remaining sections of the paper are organized as follows. In Sec. II, we briefly revisit the i2HDM model, the corrections of mass spectrum beyond the tree-level, and the decay width of heavier Z 2 -odd particles. In Sec. III, we discuss all the constraints from both theory and experiments used in our likelihood functions. In Sec. IV, we present our numerical analysis and the 2σ allowed regions with and without the loop corrections to the DM-nucleon scattering calculation. Finally, we conclude in Sec. V.

II. INERT TWO HIGGS DOUBLET MODEL
In this section, we first review the structure of i2HDM and its model parameters. We then discuss the possible one-loop contributions of ∆ 0 and ∆ ± , including renormalization group equations (RGEs) and electroweak symmetry breaking (EWSB). Finally, the decay widths of A and H ± in the case of compressed mass spectrum are given in Sec. II C.
A. Parameterization of the i2HDM scalar potential The i2HDM [13] is the simplest version of DM model within the two Higgs doublets framework. Compared with the single scalar doublet in the SM, the i2HDM has two scalar doublets H 1 and H 2 under a discrete Z 2 symmetry, H 1 → H 1 and H 2 → −H 2 which is introduced to maintain the stability of DM. The Z 2 symmetry cannot be spontaneously broken, so that H 2 never develops a VEV. These two doublets are (1) Here, G ± and G 0 are charged and neutral Goldstone bosons respectively.  [18]. Hence, we restrict ourselves to focus on the CP-even scalar S as DM candidate rather than CP-odd pseudo scalar A.
Unlike the general two Higgs doublet, the mixing term −µ 2 12 (H † 1 H 2 + h.c.) in the scalar potential is forbidden by the exact Z 2 symmetry and the scalar potential of i2HDM has a simpler form, After the electroweak symmetry breaking, there remains eight real parameters: five λs, µ 1 , µ 2 and the VEV v for the scalar potential. The v is fixed to be 246 GeV from observed weak gauge boson masses and the mass of h is fixed to 125 GeV. Together with the Higgs potential minimum condition, only five real parameters (µ 2 2 , λ 2 , λ 3 , λ 4 and λ 5 ) are inputs of the model. Note that the quartic coupling λ 2 is only involved in the four-points interaction of Z 2 -odd scalar bosons (|H 2 | 4 ), which is a phenomenologically invisible interaction at the tree-level. Nevertheless, the role of λ 2 is important to the calculations of the DM-nuclei elastic scattering cross section at the one-loop level [21].
Conventionally, it is more intuitive to adopt the physical mass basis as inputs where we denote Under the condition of m S < m A , we can see that λ S is always smaller than λ A in the tree-level parameterization. Reversely, the quartic couplings λ i=1,3,4,5 in terms of these 4 physical scalar masses and µ 2 2 are given by For the scenario with the compressed mass spectra, the mass splitting parameters ∆ 0 = m A − m S and ∆ ± = m H ± − m S instead of m A and m H ± are more useful. Hence, our input parameters are B. Scalar mass splittings beyond the tree level Considering a compressed mass spectrum in the i2HDM, namely small ∆ 0 and ∆ ± , the couplings λ 4 and λ 5 are naturally small as shown in Eq. (5). However, possible modifications to λ 4 and λ 5 from renormalization group equations (RGEs) beyond the tree-level may not be ignored. The one-loop RGEs of λ 4 and λ 5 are represented as [17,22] (4π) 2 dλ 4 dt = −3λ 4 (3g 2 + g 2 ) + 4λ 4 (λ 1 + λ 2 + 2λ 3 + λ 4 ) +2λ 4 (3y 2 t + 3y 2 b + y 2 τ ) + 3g 2 g 2 + 8λ 2 5 , (4π) 2 dλ 5 dt = −3λ 5 (3g 2 + g 2 ) + 4λ 5 (λ 1 + λ 2 + 2λ 3 + 3λ 4 ) where t = ln(µ/µ 0 ) with renormalization scale µ divided by the electroweak scale µ 0 = 100 GeV. Since all terms in the right hand side of Eq. (8) are proportional to λ 5 , once we set λ 5 = 0 at any reference scale, its value does not change in the one-loop RGEs. Still, there is one term proportional to g 2 g 2 in the right hand side of Eq. (7); Even if we set λ 4 = λ 5 = 0 at a specific reference scale, the value of λ 4 can be modified by the one-loop RGEs. Therefore, unlike the neutral mass splitting ∆ 0 , the charged mass splitting ∆ ± cannot be extremely small. More details can be found in Ref. [17].
Additionally, there is a finite contribution to ∆ ± from EWSB at one-loop level [23], and it is given by where f (x) is defined as As reported in Ref. [17,23], this extra mass splitting can be at most about O(100) MeV.
Therefore, the correction on ∆ ± from RGEs is usually larger than the EWSB correction at one-loop level.
As aforementioned, ∆ 0 can be very small if the discrete Z 2 symmetry on H 2 coming from global U (1) symmetry [15]. The possible one-loop contributions for ∆ 0 can be neglected once λ 5 ∼ 0 at any reference scale. Even if we set λ 4 = λ 5 = 0 at a specific reference scale based on global SU (2) symmetry or custodial symmetry on H 2 in the i2HDM [19], ∆ ± may still receive some non-negligible one-loop contributions of several hundred MeV. Note that we do not take these one-loop corrections on ∆ 0 and ∆ ± into account of our analysis. We require ∆ ± ≥ 1 GeV such that the loop corrections to ∆ ± is not important. However, once these one-loop corrections are considered, the ∆ 0 and ∆ ± in our parameter space can only be slightly shifted but our result remains unchanged.
C. Decay widths of A and H ± in compressed mass spectra In the compressed mass spectra of i2HDM, the dominant decay modes for A and H ± are A → SZ * and H ± → SW ± * , with off-shell W/Z bosons. After integrating out W and Z bosons, the decay widths for A → Sff and H ± → Sff channels can be approximately written as (12) where N i(j) c is the color factor of the i(j)-th species and f , f are SM fermions. The step function Θ comes from the four-momentum conservation. The couplings a i V and a i A can be expressed as where i runs over all SM fermion species, Q i (T 3 i ) is the charge (third component of isospin) for the i-th species, and the sine of weak mixing angle is s W = sin θ W . Similarly, the couplings c jk V and c jk A for lepton sectors can be represented as and for quark sectors where j runs over up-type fermions and k over down-type, and V CKM is the Cabbibo-Kobayashi-Maskawa matrix. We can apply the similar expression for the decay mode The Eqs. (11) and (12) show that the lifetimes of A and H ± are sensitive to ∆ 0 and ∆ ± , respectively. For example, if ∆ 0 < 2.6 GeV, the Γ(A → Sff ) < 2 × 10 −10 GeV implies that the lifetime of A, cτ A > 1 µm which is a criterion for a long-lived particle at the colliders.
In this analysis, nevertheless, both ∆ 0 and ∆ ± are required to be larger than 5 GeV due to the current constraints. Therefore, A and H ± cannot be long-lived particles.

III. CONSTRAINTS
In this section, we summarize the theoretical and experimental constraints used in our analysis. First, the theoretical constraints for the i2HDM Higgs potential such as perturbativity, stability, and unitarity will be discussed. For the current experimental constraints, we will consider collider, relic density, DM direct detection and DM indirect detection constraints.

A. Theoretical constraints
Once the extra Higgs doublet has been introduced, the theoretical constraints of Higgs potential in i2HDM, such as perturbativity, stability, and tree-level unitarity, have to be properly taken into account. As studied in the literature [24], these theoretical constraints are generically implemented in the Higgs basis λ i parameters. However, we use the physical mass basis as our inputs except for λ 2 and λ S = (λ 3 + λ 4 + λ 5 )/2. In this work, we employ the mass spectrum calculator 2HDMC [24] to make the conversion between these two bases and take care of the Higgs potential theoretical constraints. We collect those parameter points which have passed perturbativity, stability, and tree-level unitarity constraints.

Electroweak precision tests
In the i2HDM, electroweak precision test (EWPT) is sensitive to the mass splitting among these Z 2 odd scalar bosons [15]. In addition, those data can be parametrized through the electroweak oblique parameters S, T , and U [25], but these three parameters are correlated to each other. Following by Ref. [26,27], we can write down the form of χ 2 EWPT as where the covariance matrix and bases ∆ S , ∆ T , and ∆ U are given by PDG data [28]. detections reported by LEP [28], one can obtain With the above criteria, the W and Z cannot directly decay to these new scalar bosons.
Second, taking S as DM candidate, the CP-odd A can decay to SZ ( * ) , while the charged Higgs boson H ± can decay to W ±( * ) S. If H ± is heavier than A, the decay channel can also be opened. Therefore, the final states of the two production processes e + e − → H + H − and e + e − → SA can be the signatures of missing energy together with multi-leptons or multi-jets, depending on the decay products of W ± and Z. To certain extents, the signatures for charged Higgs searches in i2HDM can be similar to the supersymmetry searches for charginos at the e + e − or hadron colliders [29,30]. Nevertheless, the cross sections for fermion and scalar boson pair production are scaled by β 1/2 and β 3/2 respectively, where β is the velocity of the final state particle in the center-of-mass frame. Hence, one can expect production of the scalar pair is suppressed by an extra factor of β compared with the fermionic case. The limits for a fermion pair (chargino-neutralino) production cannot directly applied on the scalar boson pair production such as H ± H ∓ and SH ± [31]. In order to properly take the differences into account, we veto the parameter space based on the 95% OPAL exclusion [32]. The exclusion has been recast and projected on (m H ± , m A ) plane presented in Fig. 5 of Ref. [17].
Finally, for SA production mode, we can mimic the neutralino searches at LEP-II via e + e − → χ 0 1 χ 0 2 followed by χ 0 2 → χ 0 1 ff [33]. The process e + e − → SA followed by the cascade A → SZ ( * ) → Sff can give similar signature and the detail analysis had been carefully done in Ref. [34]. In our approach, we use the exact exclusion region on (m S , m A ) plane as given in Fig. 7 of Ref. [34] to veto the parameter space.
Since reconstructing these three LEP constraints with a precise likelihood can cost a lot of CPU-consuming computation, we only use hard-cuts to implement them into our analysis.

Exotic Higgs decays
Once these Z 2 odd scalar bosons are lighter than half of SM-like Higgs boson h, the Higgs  [35,36]. Including the Higgs-strahlung pp → ZH/W H and the vector boson fusion (VBF) processes, the ATLAS collaboration has reported an upper limit on the invisible branching ratio BR(h → inv.) < 0.26 at 95% confidence level [35]. Similarly, the CMS collaboration has also reported an upper limit on the invisible branching ratio BR(h → inv.) < 0.19 at 95% confidence level [36] by the combining searches for Higgs-strahlung, VBF, and also gluon fusion (ggH) processes 1 . On the other hand, the recent works about global fits of the SM-like Higgs boson searches from ATLAS and CMS data suggested a more aggressive constraint on the branching ratio for nonstandard decays of the Higgs boson to be less than 8.4% at the 95% confidence level [38].
For the sake of conservation, we only apply the constraints from the direct searches of Higgs invisible decays in our analysis based on Ref. [36]. In the near future, BR(h → inv.) will reach to the limit less than about 5% at the HL-LHC [39].

Diphoton signal strength R γγ in the i2HDM
Beside exotic Higgs decays, the SM-like Higgs boson decays into diphoton rate can also be modified in the i2HDM. The new contribution adding to SM one is the charged Higgs triangle loop. Since the couplings between the SM-like Higgs boson and the SM particles are the same as the SM at the leading order, the production cross sections of the Higgs boson will be the same as the SM ones. Hence, we can obtain the diphoton signal strength in the i2HDM by normalized to the SM value which is given by The exact formula for the partial decay width of h → γγ in the i2HDM can be found in Ref. [40,41] and BR(h → γγ) SM = 2.27 × 10 −3 are taken from PDG data [28]. We apply the public code micrOMEGAs [42] by using the effective operators as implemented in Ref. [43] to calculate BR(h → γγ) i2HDM in this study. Recently, both ATLAS and CMS collaborations have reported their searches for the Higgs diphoton signal strength [44,45]. In particular, 1 The CMS collaboration has reported their first search for the Higgs invisible decays via tth production channel at s = √ 13 TeV [37], but the constraint BR(h → inv.) < 0.46 at 95% confidence level is much weaker than the combined one in Ref. [36].
after the combined measurements of Higgs boson production in the diphoton decay channel, the measured R γγ is 1.08 +0. 13 −0.12 from the ATLAS [44]. On the other hand, the measurements of Higgs boson production via ggH and VBF from the CMS [45] give R γγ to be 1.15 +0. 15 −0.15 and 0.8 +0.4 −0.3 , respectively. One can find all of these measurements are in agreement with the SM prediction. In this work, we only use the latest ATLAS result R γγ = 1.08 +0. 13 −0.12 to constrain the i2HDM parameter space.

Mono-jet:
One possible way to search for DM at the LHC is looking at final states with a large missing transverse energy associated with a visible particle such as jet [46] and lepton [47]. In the i2HDM, the mono-jet signal is a pair of DM produced by the Higgs boson and accompanied with at least one energetic jet. If the pseudo-scalar A degenerates enough with the DM and decays into very soft and undetectable particles, it can contribute the mono-jet signature.
For the case of m S > m H /2, since the DM pairs are produced through an off-shell Higgs, the cross section of the mono-jet process is suppressed. On the other hand, a light DM produces smaller missing transverse energy which makes the mono-jet constraint be less efficient.
We recast the current mono-jet search at ATLAS [46] by using Madgraph 5 [48] and Madanalysis 5 [49]. It turns out that the current search excludes the value of λ S > 3×10 −2 for the case of m S < m H /2, and λ S > 5.0 for the higher mass of the DM. We will see later that these limits are much weaker than other DM constraints.

Mono-lepton:
The mono-lepton signal in this model arises from pp → SH ± (H ± → Sl + ν) process. However, the current mono-lepton search in Ref. [47] is not sensitive to small mass splitting ∆ ± . This is because the transverse mass of the lepton in the final state is not large enough which results the signal efficiency too low to be detected. Therefore, this constraint cannot be applied in this work.

Compressed mass spectra search:
Searches for events with missing transverse energy and two same-flavor, opposite-charge, low transverse momentum leptons have been carried out at CMS [50] and ATLAS [51,52].
These typical signatures are sensitive to any model with compressed mass spectra. In the i2HDM, the pairs of AS, AH ± and H ∓ H ± can be produced at the LHC via the qq fusion and VBF processes. The heavier scalar A then can decay into a dilepton pair via an off-shell Z boson, such that the dilepton invariant mass (m ll ) is sensitive to the mass-splitting ∆ 0 .
On the other hand, the charged Higgs H ± can decay into a lepton and a neutrino via an off-shell W boson, such that the transverse mass m T 2 is sensitive to the mass-splitting ∆ ± .
Here, we recast the SUSY compressed mass spectra search at the ATLAS [52]. We use Madgraph 5 [48] to generate the signal events at leading order which were then interfaced with Pythia 8 [53] for showering and hadronization, and Delphes 3 [54] for the detector simulations. The Madanalysis 5 package [49] is used to analyze the experiment results. We apply the same preselection requirements and signal regions selection cuts as in Ref. [52].
Two opposite-charged muons in the final states is a good signature. In our parameter space of interest, a suitable signal region is the one labeled as SR-E-low in Ref. [52] with the muons invariant mass window: 3.2 GeV ≤ m µ + µ − ≤ 5 GeV. Due to the small production cross section of AS, AH ± and H ∓ H ± , the current data at the LHC cannot probe the parameter space of interest, but the sensitivity at future HL-LHC may be expected to reach it.

C. Relic density
Assuming a standard thermal history of our Universe, the number density of a particle with mass m at the temperature T can be simply presented by a Boltzmann distribution It is easy to see that the annihilation cross section from Eq. (19) is dramatically decreased For the region of m S > m h /2, the annihilation of SS to W + W −( * ) via the four-points interaction, the s-channel h exchange, and the t and u-channel with H ± exchange, become the dominant channels. The contribution from t and u channels is small and negligible. The squared amplitude of the four-points interaction and the h exchange channels can be given by [55] where g is the electroweak coupling constant. Eq. (20) shows that if λ S is negative, the cancellation between the four-points interaction and the h exchange channels can be occurred. Because the annihilation of the four-points interaction at the early Universe cannot be neglected at the region of m S < 500 GeV [18], the correct relic density within the coannihilation scenario requires a negative λ S to cancel the annihilation contribution as shown in Eq. (20).
The co-annihilation contribution to the relic density is more complicated than the annihilation process. In particular, it depends on the size of mass splitting ∆ 0 and ∆ ± . Since we are focusing on the compressed mass spectrum scenario, the co-annihilation happens naturally and cannot be expelled from the full mass regions. In such a small mass splittings scenario, the most sdominant co-annihilation channels are: • SA → ff for a small ∆ 0 .
Unlike other subdominant co-annihilation channels, these two co-annihilation cross sections are only involved with the SM couplings. Naively speaking, the relic density in the coannihilation dominant region is essentially controlled by the two mass splittings ∆ 0 and ∆ ± .
We evolve the Boltzmann equation by using the public code MicrOMEGAs [42]. The numerical result of the relic density which has been taken into account the annihilation and the co-annihilation contributions, are required to be in agreement with the recent PLANCK measurement [20]: We would like to comment on the multi-component DM within the framework of the i2HDM. If there exists more than one DM particle in the Universe, the DM S can be only a fraction of the relic density and the DM local density. The DM constraints from relic density, direct detection, and indirect detection can be somewhat released. However, an important question followed by adding more new particles to the Lagrangian is whether the Higgs potential is altered and theoretical constraints are still validated. Such a next-tominimal i2HDM is indeed interesting but beyond the scope of our current study. Here, we only consider the one component DM scenario.

D. DM direct detection
As indicated by several Higgs portal DM models [56][57][58], the most stringent constraint on DM-SM interaction currently comes from the DM direct detection. This is also true in the case of the i2HDM, if we only consider the DM-quark/gluon elastic scattering via t-channel Higgs exchange at the leading order. Simply speaking, one would expect that the size of S − S − h effective coupling (λ S at tree-level) can be directly constrained by the latest XENON1T experiment [59].
However, for a highly mass degenerated scenario ∆ 0 0, the DM-quark inelastic scattering Sq → Aq can be described by an unsuppressed coupling Z − S − A whose size is fixed by the electroweak gauge coupling. Indeed, such the DM inelastic scattering scenario results in a huge cross section and hence already excluded by the latest XENON1T result [60,61].
However, the inelastic interactions are inefficient when the ∆ 0 is larger than the momentum exchange MeV [60]. Hence, it is safe to ignore the inelastic interactions when requiring ∆ 0 ≥ 10 −3 GeV. Note that ∆ 0 has to be greater than 5 GeV when the DM relic density constraint is considered.
Going beyond the leading order calculation, we calculate the corrections at the next leading order in the i2HDM. We fold all the next leading order corrections into the effective coupling λ eff S and their values are dependent on the relative energy scale we have set. Once the input scales of our scan parameters are fixed at the EW scale, the effective coupling λ eff S will be modified at the low energy scale where the recoil energy of DM-quark/gluon scattering is located. As shown in the Appendix A, the one-loop correction δλ is a function of m S , ∆ 0 , ∆ ± , and λ 2 . Their values can be either positive or negative. Therefore, we can introduce a factor R to illustrate the loop-induced effects, The effective coupling λ eff S is introduced here and we compute the correction δλ by using LoopTools code [62]. We consider two scenarios in this work: the XENON1T likelihood (Poisson distribution) is obtained with and without R by using DDCalc code [63]. Note that both the tree and loop level are isospin conserving. Hence, the value of R for DM-proton and DM-neutron scattering are approximately the same.

E. DM indirect detection
In addition to the singlet Higgs DM whose annihilation is only via SM Higgs exchange, the i2HDM at the present universe can be also dominated by the four-points interaction SSW + W − . If the DM annihilation is considerable, e.g. at dwarf spheroidal galaxies (dSphs) or galactic center where DM density is expected to be rich, some additional photons or antimatter produced by W bosons or SM fermion pair would be detected by the DM indirect detection. Unfortunately, none of the indirect detection experiments reports a positive DM annihilation signal but gives a sever limit on the DM annihilation cross section. Because dSphs kinematics is more precisely measured and it reduces the systematic uncertainties, the most reliable limits at the DM mass ∼ O(100 GeV) for DM annihilation cross section are, so far, from Fermi dSphs gamma ray measurement [64].
On the other hand, several theoretical groups based on the current astrophysical knowledge found out some anomalies such as GCE [3,[65][66][67] and AMS02 antiproton excess [4,5] which might be able to be explained by the DM annihilation. Interestingly, these two anomalies are located at the DM mass ∼ 50 − 100 GeV region where coincides with the mass region discussed in this paper. Hence, we adopt a strategy in this work that only Fermi dSphs gamma-ray constraints are included in the scan level, but our allowed parameter  space is compared with antiproton anomaly.
The differential gamma-ray flux due to the DM annihilation at the dSphs halo is given by The J-factor is J = dldΩρ(l) 2 , where the integral is taken along the line of sight from the detector with the open angle Ω and the DM density distribution ρ. We adopt 15 dSphs and their J-factor as implemented in LikeDM [68]. We sum over all the DM annihilation channels ch. The annihilation branching ratio BR(ch) and energy spectra dN ch γ /dE γ are computed by using micrOMEGAs in which the three-body final states (e.g. SS → W W * → W + l − ν) are properly taken into account.

A. Numerical method
In a similar procedure developed in our previous works [18,[69][70][71][72][73], we use the likelihood distribution given in the Table I Here, we only choose ∆ ± up to 30 GeV but one can freely extend it to a larger value until m H ± ∼ 250 GeV disfavored by EWPT data. On the other hand, without the helps from the charged Higgs co-annihilations, we have checked that the total co-annihilation contribution to the relic density is no longer larger than 85%.
In order to scan the parameter space more efficiently, we set the range of λ 2 up to 4.2 allowed by the unitarity constraint [18]. We compute the mass spectrum, theoretical conditions and oblique parameters by using 2HDMC. All survived parameter space points are then passed to micrOMEGAs to compute the relic density and the tree-level DM-nucleon elastic scattering cross section. The XENON1T statistics test is computed by using DDCalc.

The loop corrections of DM-nucleon elastic scattering cross section is calculated by using
LoopTools. In the end, we have collected more than 2.5 million data points and our achieved coverage of the parameter space is good enough to pin down the contours by using "Profile Likelihood" method [75]. Under the assumption that all uncertainties follow the approximate Gaussian distributions, confidence intervals are calculated from the tabulated values of ∆χ 2 ≡ −2 ln(L/L max ). Thus, for a two dimension plot, the 95% confidence (2σ) region is defined by ∆χ 2 ≤ 5.99.

B. Co-annihilation and mixed scenarios
As mentioned in Sec. III C, the compressed mass spectrum scenario guarantees the number density of S, A and H ± at the same temperature before freeze-out are comparable, but the thermal averaged cross section of co-annihilation is not necessary to be a similar size of annihilation. Therefore, we have to consider simultaneously the impacts from both the thermal averaged cross section and their mass splitting in order to figure out the contribution from co-annihilation processes, namely interaction rate. If the interaction rate is below the The red square and blue cross are the co-annihilation (f ann < 15%) and mixed regions (f ann >

15%), respectively
Hubble expansion rate, the freeze-out mechanism occurs. Additionally, it is hard to cease the contributions from the annihilation processes, particularly from four-points interactions whose the contributions are always significant at the mass region m h /2 < m S 500 GeV.
Here, we introduce a new parameter f ann. to account for the fraction of interaction rate attributable to the annihilation. Following the convention from micrOMEGAs, the fraction f ann. can be given by where Γ ann. and Γ tot. are the annihilation and total interaction rate before freeze-out, respectively. In this analysis, we assume that the co-annihilation domination before freeze-out acquires f ann. < 0.15.
The allowed regions of theoretical conditions, LEP, and PLANCK 2σ projected on (∆ 0 , f ann. ) and (∆ ± , f ann. ) plane are shown in the two upper panels of Fig. 1. We first found that the co-annihilation rate (Γ tot. − Γ ann. ) is always larger than 0.45 in our parameter space.
Therefore, there is no pure annihilation contribution in this compressed mass spectrum scenario. We also see that the relic density constraints can give a lower limit for both ∆ 0 and ∆ ± because the co-annihilation processes are too sufficient to reduce the relic density.
On the other hand, ∆ 0 < 8.7 GeV is caused by the LEP-II constraints. We note that there are no upper limit of ∆ ± from neither OPAL exclusion nor PLANCK measurement.
Once the phase space of SSW W four points interaction process is so suppressed at the early universe, the co-annihilation rate usually dominates Γ tot before freeze-out. This leads that co-annihilation scenario f ann. < 0.15 only located at the region m S < 70 GeV. Even if the phase space is suppressed, the annihilation channel SS → W + W − is still at least 5% contribution to Γ tot .
There are three edges of the allowed region in (m S , ∆ 0 ) plane. As mentioned before, the upper limits of ∆ 0 is due to the LEP-II limit; the left-bottom and right-bottom corner are excluded by too little and too much relic density, respectively. This also results in an upper limit on m S . The small m S and ∆ ± regions are excluded by the OPAL as shown in the right bottom panel of Fig. 1.

C. Results
In this section, we present the 2σ allowed region based on the total likelihoods, as shown in Table I. Those constraints have been discussed in Sec. III and they are referred to the phrase "all constraints", unless indicated otherwise. Fig. 2 shows the 2σ allowed region by applying all constraints. The left panel is on the (m S , ∆ 0 ) plane and the right panel is on (m S , ∆ ± ) plane. Here, the σ SI p is computed at LO. By comparing with the lower panels in Fig. 1, we can easily see that the allowed range of the DM mass is shrunk to 60 < m S / GeV < 72 after applying all constraints. Since the correlations between m S , ∆ 0 , and ∆ ± are non-trivial but important for understanding the different features of mixed scenario (blue crosses) and co-annihilation scenario (red squares), we label some specific regions as shown in the left panel of Fig. 2  for these regions as follows: Region (a), the small ∆ 0 can provide an efficient co-annihilation to reach the correct relic density even without the annihilation. Generally speaking, once the annihilation process SS → W ± W ∓ is open, it can be too efficient to reduce the relic density. The only way to get rid of annihilation is via a cancellation between four points interaction and s-channel Higgs exchange. Therefore, the coupling λ S shall be negative. On the other hand, the λ S can be very small in this region as long as SS → W ± W ∓ is close.
Region (b), comparing with the lower panel in Fig. 1, involving the Higgs invisible decay constraint lifts the lower limit on DM mass to 60 GeV.
Region (c), the co-annihilation scenario cannot reach the large m S regions, particularly m S < 64 GeV for ∆ 0 = 8.7 GeV and m S < 66 GeV for ∆ 0 = 7.4 GeV. This is due to the current DM direct detection constraint. In particular, a more negative value of the coupling between the DM and Higgs boson is needed in the larger DM mass region so that the cancellation between SS → h * → W W * and the four points interaction SSW W can occur to satisfy the correction relic density. However, this gives rise to the DM-proton scattering cross section and hence excluded by the XENON1T measurements. We also note that, due to the OPAL exclusion, a smaller DM mass region results in a larger ∆ ± value. For the co-annihilation scenario, one can see that ∆ ± > 23 GeV as shown in the right panel of Fig. 2.
Region (d), it is totally opposite to the region (a). The relic density at this region mainly comes from SS → W ± W ∓( * ) annihilation. Therefore, ∆ 0 shall be large enough in order to suppress co-annihilation contributions. The lower bound of ∆ 0 is varied with respect to m S . The mixed scenario can reach a larger DM mass region as compared with the co-annihilation scenario, hence, due to the OPAL exclusion, the lower limit on ∆ ± can be weaker. As shown in the right panel of Fig. 2, the mass splitting ∆ ± > 19 GeV for the mixed scenario.
Region (e), the LEP-II exclusion is presented. Together with the current XENON1T constraint, we found two important upper limits: ∆ 0 < 8.8 GeV and m S < 72 GeV.
As mentioned in the previous section, the current searches at the LHC is not yet to bite the parameter space. However, we find that the extended searches for LHC Run 3, especially for the compressed mass spectra searches, can make both co-annihilation and mixed scenarios to be testable. On the left panel of Fig. 2, we also show the future prospect contours of 2σ significance from the LHC compressed mass spectra searches with the integrated luminosity of 250 fb −1 (black line), 500 fb −1 (green line), and 750 fb −1 (orange line). Here, we fix λ S = 0, λ 2 = 1 and ∆ ± = 28 GeV as a benchmark point and only recast the dimuons final state of SR-E-low signal region [52]. One can see that both scenarios can be partly probed at the LHC Run 3 with L = 250 fb −1 . The parameter space in the co-annihilation scenario can be mostly covered if the integrated luminosity reach 500 fb −1 , while one needs the luminosity about 750 fb −1 to probe the whole parameter space in the mixed scenario.
A combined analysis of various signal regions, as the strategy presented in Ref. [52], will certainly improve the significance, however, it is beyond the scope of this work. We will return to this in the future.
In Fig. 3, we discuss the loop correction parameter δλ as a function of ∆ ± , λ S , and λ 2 .
Although the ∆ 0 can alter the size of δλ, it is not significant due to its smallness required by the relic density constraint. In the two upper frames of Fig. 3, we separately discuss the cases of λ S > 0 (upper-left panel) and λ S < 0 (upper-right panel) on the (λ 2 , δλ) plane. The cross section σ SI p is proportional to the effective coupling squared (λ eff S ) 2 which is sensitive to the relative signs of |λ S | and |δλ| if they are in the same order. In particular, the cross section at NLO can even be lower than LO if the cancellation between λ S and δλ occurs.
A large value of |λ S | tends to accompany with a opposite sign of δλ, while a small value of |λ S | can either accompany with a same or opposite sign of δλ. Such a configuration results a small σ SI p at NLO as well as escapes from the XENON1T constraint. The scatter plots of all allowed points on the (λ 2 , δλ) plane are shown in the lower frames of Fig. 3. The color codes indicate the value of ∆ ± (lower-left panel) and the governed relic density channels (lower-right panel). For the region of λ 2 > 0.2, the loop correction parameter δλ highly depends on λ 2 . Especially, due to the cancellation between loop diagrams, δλ drops down at λ 2 ∼ 1.4. However, it significantly increases for a larger value of λ 2 . On the other hand, ∆ ± plays a significant role in δλ for the region of λ 2 < ∼ 0.2. In particular, the lower limit on ∆ ± which is due to the OPAL exclusion, can give a lower limit on the loop correction |δλ|. Furthermore, in this region, |δλ| is restricted to be small (|δλ| < 1.1 × 10 −3 ) due to our choice of ∆ ± < 30 GeV. As aforementioned, the OPAL exclusion is more stringent constraint on the co-annihilation than the mixed scenario, hence a stronger lower limit on |δλ| is set for the co-annihilation scenario. This results in the blank strip for the coannihilation scenario (red crosses) appears at around the region −10 −3 < δλ < 0. We note that the loop correction |δλ| is slightly larger for the mixed scenario in the region λ 2 > ∼ 1.4. As shown in Eq. (22), we can understand the NLO correction by comparing λ S and λ eff S in function of m S and ∆ 0 . The former one presents the hSS coupling at LO while the later one includes the NLO corrections. In Fig. 4, we show the 2σ allowed region including all constraints. The elastic scattering cross section σ SI p in the two left panels are based on the tree level coupling λ S , but the two right panels are based on λ eff S . Let us start with the two upper panels of Fig. 4. We can see two interesting regions: i) m S m h /2, and ii) m S m h /2.
Because of Higgs resonance or co-annihilation process in the early universe, the λ S of the first region is required to be small to fulfill the relic density constraint, but such a small coupling makes an undetectable σ SI p in the present XENON1T experiment. Particularly, the λ S in the co-annihilation region can be even smaller than annihilation (Higgs resonance) region.
However, the NLO correction can enhance σ SI p to be more detectable in the near future direct detection experiments, and thus the effective coupling |λ eff S | can reach to about 0.003 as seen in the upper right panel. Regarding the second region m S m h /2, the process with four points interaction SSW W becomes too sufficient to reduce the relic density. Hence, the negative λ S is required to make a cancellation between the process with four points interaction and SS → h * → W W * . We can see that the exact cancellation happens at the strip region where the mixed scenario is absent because the universe is over abundant. Once the co-annihilation mechanism is triggered, the correct relic density can be still obtained in this region and nearby. For co-annihilation scenario, the above cancellation also makes an upper bound of DM mass with respect to the size of λ S . However, the effective coupling |λ eff S | in this region can be larger than λ S due to the loop contributions from Eq. (22) and the Appendix A, especially for the small λ S region.
The two lower panels of Fig. 4 show the 2σ distribution of λ S (left panel) and |λ eff S | (right panel) in function of ∆ 0 . The maximum value of ∆ 0 for co-annihilation scenario is due to the combined constraints from LEP-II and relic density. Interestingly, the asymmetry between positive and negative λ S is mainly in the co-annihilation scenario. As mentioned in the previous paragraph, the positive λ S in the co-annihilation scenario is disfavored because the four points interaction SSW W annihilation process opens. In addition, the co-annihilation gradually losses the power if splitting ∆ 0 increases and therefore the positive λ S has to be decreased once co-annihilation dominates at the early universe. However, the negative λ S is rather favorable because of the cancellation between diagrams of the four points interaction and SS → h * → W W * . Again, we can see from the lower-right panel that the NLO correction can generally increase the size of |λ eff S |, especially for the positive λ S region which has been excluded if only tree level is considered. the plane of (λ S ×10 3 , λ 2 ). The scattered points represent σ SI p at NLO while the region inside black dashed contours indicates that σ SI p is computed at LO. The color scheme is the same as Fig. 1.
σ SI p used in likelihood while the scatters are based on the NLO computation. First, the NLO contribution is only sensitive to a small ∆ 0 or ∆ ± region, but the larger lower limit ∆ 0 ∼ 7.5 GeV is required by relic density constraints. From the distribution in the (∆ 0 , ∆ ± ) plane, we can see that it is hard to make a significant difference between LO and NLO computations in the case of compressed mass spectra.
Those blue cross points present the mixed scenario, but the relic densities of red square points are governed by co-annihilation. Clearly, we find that lower limits of ∆ ± for the mixed scenario are lower than co-annihilation scenario. For the mixed scenario, the lower limit of ∆ ± can be easily understood from the right panel of Fig. 2. The small value of ∆ ± for the mixed scenario is still allowed by the OPAL constraint. Namely, the ∆ ± 24 GeV for co-annihilation scenario located at the m S 66 GeV has been ruled out. Hence, those points in the mixed scenario with ∆ ± < 24 GeV are located at the m S > 66 GeV.
In the right panel of Fig. 5, the NLO effects are more sizable. The coupling λ 2 is not involved in LO processes, but it has some interesting behaviors once the NLO corrections are included. There are two main features. The first feature is some cancellations between the tree level diagrams with a negative λ S and loop level diagrams with large λ 2 . This cancellations can be found in the region λ S < −4.8 × 10 −4 and λ 2 > 1 as well as the region λ S > 4 × 10 −4 where the parameter space were excluded from LO computation, but they are saved when including NLO corrections. The second feature is that the NLO contribution generally increases σ SI p . Apart from the cancellation region, more parameter space have been ruled out comparing with only LO computation. In Fig. 6, we show the difference between σ SI p with tree-level computation (left panel) and next-leading order included (right panel). In order to present the exclusion power of XENON1T, we do not include the XENON1T constraint in Fig. 6. As the previous discussion, the inelastic scattering is not sufficient at the ∆ 0 > O( MeV) region and the dominant contribution of elastic scattering is still t-channel Higgs exchange. This makes a drawback to detect the co-annihilation and Higgs resonance region because both scenarios are generally required a small λ S to fulfill the relic density constraints. In the left panel, we can see that the resonance region m S < m h /2 is overall lower than the future DARWIN projection. Furthermore, the co-annihilation region at m S < m h /2 is even below the neutrino floor and hard to be detected under current strategies of DM direct detection.
Strikingly, once the next-leading order correction for σ SI p is considered, the elastic scattering cross section σ SI p of the Higgs resonance region and the co-annihilation region are both enhanced. Thanks to the loop contributions, especially from the large |λ 2 |, co-annihilation region becomes more visible in the future direct detection. The boost of σ SI p can be even with the order 10 3 for the region near the resonance.  [79]. The green error bar is the 1σ signal region for anti-proton excess [80].
Lastly, we would like to discuss the detection of the DM annihilation at the present. In Fig. 7, we show the 2σ distribution in (m S , σv ) plane, allowed by all constraints in Table I if the antiproton anomaly would be confirmed in the future, one would expect to see a signal at this DM mass region from the future HL-LHC compressed mass spectrum searches.

V. SUMMARY AND FUTURE PROSEPECT
The compressed mass spectrum in the i2HDM is an interesting topic. It can escape from most of current collider constraints and its relic density in this region can be mainly governed by S and A co-annihilations. To obtain the correct relic density and minimize the contribution from annihilation channels, there are two main features of the co-annihilation scenario. The first one is the requirement for a negligible |λ S | near the Higgs resonance region. As a type of Higgs portal DM model, one can intuitively reduce λ S to minimize annihilation contributions. However, this method only works at the region where SS → W + W − * is not open. Once this channel is kinematically-allowed at the early universe, the SSW + W − four points interaction is very sufficient to reduce the relic density and λ S is not relevant. To lower down the contribution from the SSW + W − four points interaction, one has to choose a negative sign of λ S to make a cancellation between four points interaction and s-channel Higgs exchange. Therefore, the second feature is that a negative λ S is needed if SS → W + W − * channel is kinematically-allowed at the early universe.
Because of a small or negative λ S , the NLO corrections for DM-nucleon elastic scattering become more important. We first fix the input scales of our scan parameters at the EW scale, then the effective coupling λ eff S will be modified from quantum corrections at the low energy DM direct detection scale. At the tree-level, a small value of |λ S | predicts a cross section σ SI p smaller than the neutrino floor. However, it can be testable if σ SI p includes the NLO corrections. If the value of |λ S | is large and λ S is negative, a cancellation between LO and NLO contribution may be needed in order to escape from the present XENON1T constraint. Such a NLO contribution is sensitive to not only the mass splitting ∆ ± but also the coupling λ 2 . Interestingly, the coupling λ 2 plays no role at the tree-level phenomenology, neither Higgs nor DM. Hence, motivated by the non-trivial correlations between compressed mass spectra, co-annihilations, and NLO correction of σ SI p , we conduct a global scan to comprehensively explore the parameter space of compressed mass spectra in the i2HDM.
In this paper, we have performed a global statistical analysis of the i2HDM with five parameters (m S , ∆ 0 , ∆ ± , λ 2 , λ S ) at the EW scale. Particularly, the parameters ∆ 0 and ∆ ± are adopted for searching the compressed mass spectrum and co-annihilation scenario. By using the profile likelihood method, the survived parameter space were subjected to constraints from the theoretical conditions (perturbativity, stability, and tree-level unitarity), the collider limits (electroweak precision tests, LEP, and LHC), the relic density as measured by PLANCK, the σ SI p limit from XENON1T, and the σv limit from Fermi dSphs gamma ray data. For the computation of σ SI p , the LO and NLO calculations are considered separately in the likelihood. We have shown the 2σ allowed points grouped by co-annihilation scenario (f ann < 15%) and mixed scenario (f ann > 15%) in two dimensional projections of the parameter space.
We found that the viable parameter space of co-annihilation scenario is located at the 60 GeV m S 66 GeV, 7.4 GeV ∆ 0 8.8 GeV, and ∆ ± 23.0 GeV. The NLO correction δλ is sensitive to λ 2 when it is greater than 0.2 while the contribution from ∆ ± dominates the δλ if λ 2 < 0.2. Due to the interplay between OPAL exclusion and PLANCK relic density constraint, the co-annihilation scenario cannot be realized with the condition 0.0 δλ 0.1 when the contribution from ∆ ± dominates the δλ.
Next, the correlation between λ S and λ 2 is non-trivial at the NLO level. For λ 2 > 1, the lower bound of λ S is shifted from −4.5×10 −3 to −6×10 −3 due to the cancellation between λ S and δλ. However, the same sign between λ S and δλ can enhance σ SI p so that the parameter space with a large value of λ S and λ 2 is excluded. For λ 2 < 1, the ∆ ± plays the role to NLO corrections. The allowed λ S for a small λ 2 is shifted from −4.5 × 10 −3 λ S 3.8 × 10 −3 (σ SI p at LO) to −3.8 × 10 −3 λ S 5.0 × 10 −3 (σ SI p at NLO). Therefore, the cross section σ SI p for co-annihilation scenario can be significantly boosted and ruled out by XENON1T. It is testable in the future direct detection experiments.
Finally, the DM with compressed mass spectrum mainly annihilates to W W * . The i2HDM predicts the region which coincides with the AMS-02 antiproton anomaly within 3σ range. This region can be also tested by the compressed mass spectra searches at the LHC. Although we found that the region is not sensitive to the current LHC searches, it can be partially probed with future luminosity 250 fb −1 and mostly probed with luminosity 750 fb −1 as shown in the left panel of Fig. 2.

Acknowledgments
We thank Tomohiro Abe and Ryosuke Sato for some useful discussions and suggestions in the part of spin-independent cross section at the next leading order in i2HDM. Y. Appendix A: Spin-independent cross section at the next leading order In this Appendix, we outline the calculations of spin-independent cross section at the next leading order from Ref. [21]. We have checked the consistency of our numerical results.
First, the effective interaction of the dark matter and quark/gluon can be represented as where O q µν is the quark twist-2 operator with the following form, The higher twist gluon operators have been neglected in the above effective Lagrangian.
Based on Eq. (A1), the scattering amplitude and spin-independent cross section of dark matter and nucleon can be written as, where µ is the reduced mass with the form µ ≡ m S m N /(m S + m N ). m S and m N are dark matter and nucleon masses, respectively.
At the leading order, only Γ q and Γ G are non-zero in Eq. (A3) and can be given as, At the next leading order, we closely follow the calculations in Ref. [21] with the following classifications,  tion of the dark matter and quark/gluon at the next leading order from Ref. [21]. Here (B.5) means equation (B.5) in Ref. [21] with W boson contribution inside the one-loop box type diagram, for example.
• One-loop box type diagrams • One-loop Higgs vertex correction diagrams • Two-loop gluon contribution diagrams We list relevant equation numbers to calculate above diagrams from Ref. [21] in the Table II.
Here (B.5) means equation (B.5) in Ref. [21] with W boson contribution inside the one-loop box type diagram, for example.
We can further define Γ q and Γ G at the next leading order as, Finally, according to Eq.(3.21) − (3.24) in Ref. [21] and arguments therein, the one-loop correction δλ are combination of functions in the Table II with the following form, (2)). (A7) where f q , f g , q(2) andq(2) are matrix elements for nucleons and f N ≡ 2 9 + 7 9 q f q . Their exact values for proton and neutron are taken from the default values of MicrOMEGAs [42].
Therefore, the loop-induced effects can be simply written as the form in Eq. (22) with four input parameters : m S , ∆ 0 , ∆ ± and λ 2 inside δλ.