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 ℤ2 odd particles — scalar S, pseudo-scalar A, and charged Higgs H±. In such a case, the co-annihilation processes play a significant role to reduce DM relic density. As long as a co-annihilation governs the total interaction rate in the early universe, a small annihilation rate is expected to reach a correct DM relic density and its coupling λS between DM pair and Higgs boson shall be tiny. Consequently, a negligible DM-nucleon elastic scattering cross section is predicted 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 quartic self-coupling λ2 between ℤ2 odd particles indeed contributes the one-loop quantum correction and behaves non-trivially for the co-annihilation scenario. Interestingly, the parameter space, which is 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 antiproton excess. The parameter space can be further probed at the future high luminosity LHC.


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 except for the gravitational force. To understand the particle nature of DM, several detection methods have been developed in the past decades such as DM direct detection (DD), indirect detection (ID) and colliders. Although some of DD and ID analyses have reported anomalies [1][2][3][4][5], DM signal is still absent in the LHC searches [6] so that the DM JHEP06(2020)033 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 due to the couplings between DM and SM particles are too weak to be detected under the current sensitivities. Hence, an upgrade design of instrument and 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 search strategies for the weakly interacting massive particles (WIMPs) are still hard to work [9]. The second possibility comes from the compressed mass spectrum models in which the next lightest dark particle and DM have a small mass splitting [10]. The signals from such compressed mass spectra are usually predicted with soft objects and then vetoed or polluted from SM backgrounds at colliders [11]. Due to the small mass splitting, the next lightest dark particle can be long lived. Interestingly, the interaction couplings between dark sector and SM fields may be not suppressed for this possibility and they can be testable in the DD or ID searches [12].
The inert two Higgs doublet model (i2HDM) [13][14][15][16] is a 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 difficult 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 within 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 DM-nucleon 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 quartic self-coupling λ 2 between Z 2 odd bosons and the mass splittings, JHEP06(2020)033 ∆ 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) contributions 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 model parameter space, particularly for the co-annihilation scenario, the DM spinindependent 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 section 2, 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 section 3, we consider both the theoretical and experimental constraints used in our likelihood functions. In section 4, 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 section 5.

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 section 2.3.

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 vacuum expectation value (VEV). These two doublets can be given as .

JHEP06(2020)033
Here, G ± and G 0 are charged and neutral Goldstone bosons respectively. The symmetry breaking pattern for the doublets are H T 1 = 0, v/ √ 2 and H T 2 = (0, 0), where v ≈ 246 GeV. In the end, we have five physical mass eigenstates: two CP-even neutral scalar h and S, one CP-odd neutral scalar A, and a pair of charged scalars H ± .
Before going to the detailed calculation, let us briefly recap the main features of the i2HDM. First, Z 2 -odd particles S, A and H ± are not directly coupled to SM fermions while Z 2 -even Higgs h is identified as the SM Higgs with mass ∼ 125 GeV. Second, owing to the exact Z 2 symmetry, there is no tree-level flavor changing neutral current. Finally, the DM candidate can be either S or A depending on their masses, but it is hard to phenomenologically distinguish one from the other [18]. Here, we restrict ourselves to focus on the CP-even scalar S as the DM candidate rather than the CP-odd pseudo scalar A.
Unlike the general two Higgs doublet model, since the mixing term −µ 2 12 (H † 1 H 2 + h.c.) is forbidden by the exact Z 2 symmetry, 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. Because two parameters the VEV v and Higgs mass can be fixed by the experimental observations and one parameter can be eliminated by 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 fourpoints 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 Assuming m S < m A , we can see that λ S is always smaller than λ A at 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 {m S , ∆ 0 , ∆ ± , λ 2 , λ S }. (2.6)

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. (2.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π 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. (2.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. (2.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 This extra mass splitting can be at most about O(100) MeV [17,23] and it is usually smaller than the one from RGEs. 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. On the other hand, 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 have a correction of several hundred MeV from the oneloop contributions. However, the effects from the loop corrections can be safely ignored in this analysis, as we require ∆ ± ≥ 1 GeV. Indeed, if one takes the one-loop corrections on ∆ 0 and ∆ ± into account, the parameter space can only be slightly shifted but our result remains unchanged. Thus, we do not include these corrections in this analysis.

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 given by 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 s W stands for sin θ W with θ W being the weak mixing angle. Similarly, the couplings c jk V and c jk A for lepton sectors can be represented as and for quark sectors where j (k) runs over up-type (down-type) fermions and V CKM is the Cabbibo-Kobayashi-Maskawa matrix. We can apply the similar expression for the decay mode H ± → Aff . The eq. (2.11) and eq. (2.12) show that the lifetimes of A and H ± are sensitive to ∆ 0 and ∆ ± , respectively. For example, if ∆ 0 < 2.6 GeV, the decay width Γ(A → Sff ) < 2 × 10 −10 GeV implies that the lifetime of A is longer than the long-lived particle criterion at the LHC, ∼ 1 µm/c. 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.

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 the perturbativity, stability, and unitarity will be discussed. For the current experimental constraints, we will consider the collider, relic density, DM direct detection and DM indirect detection constraints.

Theoretical constraints
Once the extra Higgs doublet has been introduced, the theoretical constraints of Higgs potential in i2HDM, such as the 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, in this analysis, we use the physical mass basis as our inputs except for λ 2 and λ S = (λ 3 + λ 4 + λ 5 )/2. 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 the 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], and these three parameters are correlated to each other. Following by refs. [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].

Scalar bosons production at the LEP
Generally speaking, new scalar bosons (S, A and H ± ) can be produced either singly or doubly at the colliders. However, due to the protection of the extra Z 2 symmetry in the i2HDM, all of these new scalar bosons can only be produced doubly. Therefore, those searches of single new scalar boson production in LEP, Tevatron and LHC cannot be applied in the i2HDM case. We review the searches for new scalar boson pair productions at the LEP in the following. First, if new scalar bosons are lighter than W or Z boson, the decay channels such as W ± → {SH ± , AH ± } and/or Z → {SA, H + H − } can be detected by LEP. Utilizing the null signal detections reported by LEP [28], one can obtain With the above criteria, the W and Z bosons cannot directly decay into these new scalar bosons.

JHEP06(2020)033
Second, taking S as the DM candidate, the CP-odd A can decay into SZ ( * ) , while the charged Higgs boson H ± can decay into W ±( * ) S. If H ± is heavier than A, the decay channel H ± → W ±( * ) A → W ±( * ) SZ ( * ) 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 bosons. To certain extents, the signatures for charged Higgs searches in i2HDM can be similar to the supersymmetry searches for charginos at the e + e − and hadron colliders [29][30][31]. Nevertheless, the cross sections for fermion and scalar boson pair productions 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 the 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 ± [32]. In order to properly take the differences into account, we veto the parameter space based on the 95% OPAL exclusion [33]. The exclusion has been recast and projected on (m H ± , m A ) plane presented in figure 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 [34]. 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. [35]. In our approach, we use the exact exclusion region on (m S , m A ) plane as given in figure 7 of ref. [35] to veto the parameter space.
Since reconstructing these three LEP constraints with a precise likelihood will 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 a half of the SM-like Higgs boson h, the Higgs exotic decays h → SS/AA/H + H − can be opened. These exotic Higgs decays can modify the total decay width of the SM-like Higgs boson as well as the SM decay branching ratios which can be constrained by the current Higgs boson measurements and further tested by the future Higgs boson precision experiments. For the compressed mass spectrum scenario with the DM candidate S, the final states for h → SS are invisible while h → AA/H + H − are missing energy plus very soft jets or leptons. In the case that these jets or leptons are too soft to be detected at the LHC, the signatures of h → AA/H + H − are identical to h → SS. Recently, both ATLAS and CMS have reported their updated limits on the branching ratio of Higgs invisible decays [36,37]. 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 [36]. Similarly, the CMS collaboration has also reported an upper limit on the invisible branching ratio BR(h → inv.) < 0.19 at 95% confidence level [37] by the combining searches for Higgs-strahlung, VBF, and also gluon fusion (ggH) processe.s 1 On   1 The CMS collaboration has reported their first search for the Higgs invisible decays via tth production channel at s = √ 13 TeV [38], but the constraint BR(h → inv.) < 0.46 at 95% confidence level is much JHEP06(2020)033 the other hand, a recent global-fit analysis on the SM-like Higgs boson measurements using 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 [39].
In the near future, BR(h → inv.) is expected to reach the limit less than about 5% at the HL-LHC [40]. For the sake of conservation, we only use the result from CMS [37] in this analysis.

Diphoton signal strength R γγ in the i2HDM
Beside exotic Higgs decays, the rate of the SM-like Higgs boson decaying into diphoton can also be modified. In particular, the new contribution adding to the SM one is the charged Higgs triangle loop. Since, at the leading order, the couplings between the SM-like Higgs boson and SM particles are unchanged, the production cross section of the Higgs boson will be the same as the SM one. Hence, we can obtain the diphoton signal strength in the i2HDM by normalized to the SM value: The exact formula for the partial decay width of h → γγ in the i2HDM can be found in refs. [41,42] and BR(h → γγ) SM = 2.27 × 10 −3 are taken from PDG data [28]. We apply the public code micrOMEGAs [43] by using the effective operators as implemented in ref. [44] to calculate BR(h → γγ) i2HDM in this study. Recently, both ATLAS and CMS collaborations have reported their searches for the Higgs diphoton signal strength [45,46]. In particular, a combined measurements of the Higgs boson production from ATLAS [45] gives R γγ = 1.08 +0.13 −0.12 . On the other hand, the measurements of Higgs boson production via ggH and VBF from CMS [46] give R γγ = 1.15 +0. 15 −0.15 and 0.8 +0.4 −0.3 , respectively. One can see that all of these measurements are in agreement with the SM prediction. In this analysis, we only use the latest ATLAS result [45] to constrain the model parameter space.

3.2.5
Mono-X and compressed mass spectra searches at the LHC 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 [47] and lepton [48]. 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 has a small mass splitting with the DM and decays into very soft and undetectable particles, it can also contribute to 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, for m S < m h /2 case, the missing transverse energy is usually low that the mono-jet constraint is less efficient.
We recast the current ATLAS mono-jet search [47] by using Madgraph 5 [49] and Madanalysis 5 [50]. It turns out that the current search excludes λ S > 3 × 10 −2 for the case of m S < m h /2, and λ S > 5.0 for higher DM masses. We will see later that these limits are much weaker than other DM constraints.
weaker than the combined one in ref. [37].

JHEP06(2020)033
Mono-lepton: the mono-lepton signal in this model is raised from the process pp → SH ± with H ± → Sl + ν. However, the current mono-lepton search from ATLAS [48] is not really sensitive to the small mass splitting ∆ ± . Indeed, the signal efficiency is too low to be detected because the transverse mass distribution of the lepton and missing transverse momenta in the final state is not large enough. 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 from CMS [51] and ATLAS [52,53] Collaborations. These typical signatures are sensitive to any model with compressed mass spectra if the production cross section is large enough. 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. The stransverse mass m T 2 is sensitive to the mass-splitting ∆ ± .
We recast the ATLAS SUSY compressed mass spectra search [53]. The matrix element generator Madgraph 5 [49] is used to generate the signal events at leading order which are then interfaced with Pythia 8 [54] for showering and hadronization, and Delphes 3 [55] for the detector simulations. The Madanalysis 5 package [50] is used to recast the experiment results. We apply the same preselection requirements and signal regions selection cuts as in ref. [53]. Two opposite-charged muons in the final states are chosen for recasting in this study. In our parameter space of interest, a suitable signal region is the one labeled as SR-E-low in ref. [53] with the muon pair 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.

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 ∼ e −m/T . Under the thermal DM scenario, if the mass splittings between S, A, and H ± are as small as the case considering in this analysis, the number densities of three particles at T are comparable to each other and their co-annihilating processes play an important role of reducing the relic density. In this subsection, we summarize the dominant channels for annihilation and co-annihilation in the i2HDM.
Depending on the specific DM mass range, the DM annihilation is dominated by different channels. For the Higgs resonance regions, i.e. m S m h /2, the annihilations of SS to SM fermions via the Higgs boson exchange are dominated, especially for the process JHEP06(2020)033 It is easy to see that the annihilation cross section from eq. (3.4) is dramatically decreased when m S > m h /2. For the region of m S > m h /2, the annihilation of SS to W + W −( * ) via the four-vertex interaction, the s-channel h exchange, and the t and u-channel with H ± exchange, becomes the dominant channel. The contribution from four-vertex diagram is usually dominant while the one from charged Higgs is typically small. The one from h exchange becomes important only at around the Higgs resonance region. However, the total SS → W + W −( * ) annihilation can be significantly reduced by the cancellations among these contributions. First, the cancellation between the four-vertex diagram and h exchange contribution can take place at λ S = −(s − m 2 h )/(2v 2 ) [56]. Second, if the masses of S and H ± are nearly degenerate, the cancellation between the four-vertex diagram and charged Higgs contribution can also occur. In figure 1, we show the cross section times relative velocity (σv rel ) for the process SS → W W * → W lν as a function of ∆ ± with various values of JHEP06(2020)033 λ S . The cross section is computed by using MadGraph5. The DM mass is fixed to be m S = 64 GeV and its relative velocity in the early universe is taken as v rel /c = 0.02 in order to realize a cancellation between the four-vertex diagram and h exchange contribution for λ S = −(s − m 2 h )/(2v 2 ) −4.2 × 10 −3 . We also present the cyan line (λ = 0 with v rel = c/3) as a reference. Three benchmark values of λ are selected based on the allowed region given in a previous study [18]. For a positive value of λ S = 10 −3 (black solid line), the contribution from the four-vertex diagram and s-channel Higgs exchange are dominant. When the s-channel Higgs exchange is removed by setting λ S = 0 (blue dashed line), σv rel is reduced but still sizable. For λ S = −4.2×10 −3 (red dashed-dotted line), one can see that σv rel is the smallest and even can be neglected if the charged Higgs mass becomes heavy. This is due to the maximal cancellation between the four-vertex diagram and h exchange contribution. The remaining σv rel is mainly from the charged Higgs contribution. One can also see that from the black solid and blue dashed lines, σv rel is slightly dropped if the mass splitting ∆ ± decreases from ∼ 250 GeV down to ∼ 1 GeV. This is the result of the cancellation between the four-vertex diagram and t/u channel charged Higgs contributions.
Unlike the most parameter space of the i2HDM, the difficulty of a compressed mass spectrum scenario is that additional co-annihilation channels in this scenario make underabundant relic density. Ideally, the annihilation has to be properly switched off, but it is particularly hard for SS → W W ( * ) because the coupling of four-vertex diagram SSW W is the electroweak coupling. Therefore, the only way to reduce the annihilation cross section of SS → W W * is via cancellation as shown in figure 1. Since the mass splitting ∆ ± in this work is always greater than 1 GeV, it is interesting to see the role of the t/uchannel with H ± exchange playing in the cancellation. By engaging with micrOMEGAs package, 2 we have checked that if we change ∆ ± from 20 GeV to 500 GeV, the relic density is slightly changed by a few percent. Hence, the t/u-channel with H ± exchange are not the leading contribution to suppress the contribution from SS → W W * annihilation in the co-annihilation scenario.
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 splitting scenario, the most dominant 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 [43]. The numerical result of the relic density which has been taken into account the annihilation JHEP06(2020)033 and co-annihilation contributions, is required to be in agreement with the recent PLANCK measurement [20]: Ω CDM h 2 = 0.1199 ± 0.0027. (3.5) 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 nextto-minimal i2HDM is indeed interesting but beyond the scope of our current study. Here, we only consider the one component DM scenario.

DM direct detection
As indicated by several Higgs portal DM models [57][58][59], 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 [60].
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 predicts a huge cross section but it has already excluded by the latest XENON1T result [61,62]. The inelastic interactions are inefficient when the ∆ 0 is larger than the momentum exchange O(200 keV) [61]. Hence, it is safe to ignore the inelastic interactions when requiring ∆ 0 ≥ 10 −3 GeV. We not 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 which depends 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 . Its value can be either positive or negative. Therefore, we can introduce a factor R to illustrate the loop-induced effects, The correction parameter δλ is computed by using LoopTools code [63]. We consider two scenarios in this work: the XENON1T likelihood (Poisson distribution) is obtained with and without R by using DDCalc code [64]. 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.

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. Thanks to a better measured dSphs kinematics which gives smaller systematic uncertainties than other DM indirect detection, the most reliable limit on the DM annihilation cross section at the DM mass ∼ O(100 GeV) comes from Fermi dSphs gamma ray measurements so far [65].
On the other hand, several groups have found out some anomalies such as GCE [3,[66][67][68] 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-factors as implemented in LikeDM [69]. 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. In this paper, we only focus on the region of m S < 100 GeV. For the DM indirect detection at the region of m S > 100 GeV, the future Cherenkov Telescope Array may give a severe limit [70].

Numerical method
In a similar procedure developed in our previous works [18,[71][72][73][74][75], we use the likelihood distribution given in the we perform 35 Markov chains in the five dimensional parameter space, 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 the EWPT data (for a detailed analysis, see ref. [77]). 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 cannot be 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 [78]. 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.

Co-annihilation and mixed scenarios
As mentioned in section 3.3, the compressed mass spectrum scenario guarantees the number density of S, A and H ± at the same temperature before freeze-out are comparable. Nevertheless, the thermal averaged cross section of co-annihilation could be very different with annihilation one. Therefore, the interaction rate can describe the effects from both the thermal averaged cross section and their mass splitting. If the interaction rate is below the Hubble expansion rate, the freeze-out mechanism occurs. Additionally, it is hard to cease 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 JHEP06(2020)033 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 ) planes are shown in the two upper panels of figure 2. We found that the co-annihilation rate (Γ tot. − Γ ann ) is always larger than 45% in our parameter space of interest. 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 efficient to reduce JHEP06(2020)033 the relic density. On the other hand, ∆ 0 < 8.7 GeV is caused by the LEP-II constraints. We note that there is no upper limit of ∆ ± from neither OPAL exclusion nor PLANCK measurement. The small m S and ∆ ± regions are excluded by the OPAL as shown in the right bottom panel of figure 2.
Once the phase space of SS → W + W −( * ) process is suppressed at the early universe, the co-annihilation rate dominates Γ tot before freeze-out. As a result, the co-annihilation scenario f ann < 0.15 only located at the region m S < 70 GeV. The annihilation process SS → W + W −( * ) is still at least 5% contribution to Γ tot .
There are three edges of the allowed region in (m S , ∆ 0 ) plane. The left-bottom and right-bottom corner are excluded by too little and too much relic density, respectively. The upper limit on the mass splitting, ∆ 0 < ∼ 8.8 GeV, is due to the LEP limit and this also leads to that f ann 55%. Therefore, the expected features of annihilation are hidden in the region of the Mixed scenario. For example, near the Higgs resonance region (m S 62 GeV), one can see a small kink at the edge ∆ 0 ≈ 7.6 − 7.7 GeV of the Mixed scenario region.

Results
In this section, we present the 2σ allowed region based on the total likelihoods, as shown in table 1. Those constraints have been discussed in section 3 and they are referred to the phrase "all constraints", unless indicated otherwise. Figure 3 shows the 2σ allowed region by taking into account 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 figure 2, 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 JHEP06(2020)033 for understanding the different features of mixed scenario (blue crosses) and co-annihilation scenario (red squares), we label some specific regions shown in the left panel of figure 3 and discuss some features 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. One has to keep in mind that the coupling λ S can be either very small or a negative value in this region. It depends on whether SS → W W * is close (a tiny λ S ) or opened (a negative λ S ).

Region (b).
Comparing with the lower panel in figure 2, 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 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 relic abundance. However, this enhances σ SI p to be 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 figure 3.

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 to suppress co-annihilation and its lower bound is varied with respect to m S . The mixed scenario can reach a larger DM mass region as compared with the co-annihilation scenario. On the other hand, the lower limit of ∆ ± for the mixed scenario from the OPAL exclusion is more released. As shown in the right panel of figure 3, the mass splitting ∆ ± > ∼ 19 GeV for the mixed scenario.
Region (e). The LEP-II exclusion is presented. Together with the current XENON1T constraint, they yield 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 at the LHC Run 3, especially for the compressed mass spectra searches, can probe both co-annihilation and mixed scenarios. On the left panel of figure 3, we 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 with the muons invariant mass window: 3.2 GeV ≤ m µ + µ − ≤ 5 GeV [53]. The significance is given by z = S/ B + δ 2 B , where S is the number of signal event, B is the number of background event and δ B is the background uncertainty which is assumed to be 10% of the number of background event at the future JHEP06(2020)033  LHC. 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. [53], will certainly improve the significance, however, it is beyond the scope of this work. On the other hand, the mass-splitting ∆ ± is sensitive to the stransverse mass m T 2 , similar recasting method can be done according to ref. [53]. We will return to these two parts in a future work.
As shown in eq. (3.6), the NLO effects on the DM-proton scattering cross section can be understood by comparing the λ S and λ eff S couplings. In figure 4, we show the 2σ allowed region taking into account all constraints. The elastic scattering cross section σ SI p in the two left panels are based on the tree level coupling λ S , and the two right panels are based on . The scatter plots in 2σ allowed region on 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. process in the early universe, the λ S of the first region is required to be small to fulfill the relic density constraint, however 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 the one in the annihilation (Higgs resonance) region. However, the NLO corrections can enhance σ SI p to be more detectable in the near future direct detection experiments, and thus the effective coupling |λ eff S | can reach about 0.003 as seen in the upper right panel. Note that λ eff S can be generated by gauge boson loops as the g 2 g 2 term shown in eq. (2.7) from one-loop RGEs. This reveals that some fine-tuning of the parameters is needed in order to reach a tiny |λ eff S | which is much less than O(10 −3 ) in figure 4.
Regarding the second region m S m h /2, the four-vertex diagram SSW W process becomes too sufficient to reduce the relic density. Hence, λ S needs to be negative so that a cancellation between this diagram and SS → h * → W W * can occur. We can see that the exact cancellation happens at the strip region, λ S ≈ −(s − m 2 h )/(2v 2 ), 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 .
The two lower panels of figure 4 show the 2σ distribution of λ S (left panel) and λ eff Co-annihilation Mixed Figure 6. The scatter plots on the (m S , σ SI p ) plane. Allowed points are in 2σ agreement with all constrains except the XENON1T measurements. The DM-proton scattering cross section is computed at LO in the left panel while the one taking into account the one loop correction is presented in the right panel. The color scheme for the scatter points is the same as figure 2. The solid black line represents the limit from XENON1T [60]. The projected sensitivities of LZ [79] and DARWIN [80] are dashed green and dashed magenta lines, respectively. The orange region is the neutrino background [81].
kinematically opens, the asymmetry between positive and negative λ S can be found in the co-annihilation scenario. 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 corrections can generally increase the size of |λ eff S |, especially for the positive λ S region which has been excluded if only tree level is considered. Note that the value of |λ eff S | can be smaller than |λ S | if the loop correction parameter δλ and λ S are opposite signs as shown in appendix B.
In figure 5, the NLO effects are more distinguishable on the plane of (λ S × 10 3 , λ 2 ). The coupling λ 2 is not involved in LO processes, but it causes some interesting behaviors once the NLO corrections are included. There are two main features. The first feature is that some cancellations between the tree level diagrams and loop level diagrams are taking place. Particularly, these tree-loop cancellations can be found in the region of λ 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 can be saved when including NLO corrections. The second feature is that the NLO contribution generally increases σ SI p . Apart from the tree-loop cancellation region, more parameter space have been ruled out comparing with only LO computation. Figure 6 shows the DM-proton scattering cross section as a function of DM mass, m S , at tree-level (left panel) and one-loop level (right panel). Note that the DM direct detection constraints are not included in the scatter point regions to demonstrate its exclusion power. Again, the inelastic scattering can be neglected at the region ∆ 0 O(200 keV) [61]. Hence, JHEP06(2020)033  [82]. The green error bar is the 1σ signal region for the antiproton excess [83].
the only significant contribution to the detection rate is the elastic scattering via t-channel Higgs exchange.
Because both scenarios are generally required a small λ S to fulfill the relic density constraints, making a drawback to detect the co-annihilation and Higgs resonance region. In the left panel, we can see that the elastic scattering cross section in the resonance region, m S < ∼ m h /2, is overall lower than the future projected sensitivities from LZ [79] and DARWIN [80]. 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 in the Higgs resonance and the co-annihilation regions are both significantly enhanced. Thanks to the loop contributions, especially from the large |λ 2 |, co-annihilation region can be testable in the future direct detection searches.
Lastly, we would like to discuss the detection of the DM annihilation at the present. In figure 7, we show the 2σ distribution in (m S , σv ) plane, allowed by all constraints in table 1. With a conservative treatment, we have taken only the Fermi 15 dSphs gamma ray data into our total likelihood. However, we also show the current most stringent limit (solid black line) which is obtained by combing the latest data from Fermi-LAT, HAWC, HESS, MAGIC and VERITAS [82]. To illustrate the antiproton anomaly, we also present the signal region for bb final state [83] as a comparison.
For the region of m S < m h /2, the dominant annihilation channel is SS → h → bb where λ S is suppressed due to the relic density constraint. Thus, the velocity averaged cross section at this region is a relatively small value, σv ∼ 10 −27 cm 3 s −1 . The Higgs JHEP06(2020)033 resonance presents at the region of m S m h /2, but the large cross section part can be excluded if one takes into account the combined limit (black line). The DM annihilation at the present universe with mass m S > m h /2 are almost entirely to W W * final state and its spectrum dN/dE is similar to the bb final state. Except a small part of blue crosses (mixed scenario), most of parameter space in this region are still survived. Interestingly, the viable regions of the model including the mixed and co-annihilation scenarios are in agreement with the antiproton anomaly within 3σ [83]. Hence, 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.

Summary and future prospect
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 the 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 the 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. 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 among 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, 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.
Motivated by the non-trivial correlations between compressed mass spectra, coannihilations, and NLO corrections of σ SI p , we conduct a global scan to comprehensively explore the parameter space of compressed mass spectra in the i2HDM, including 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 JHEP06(2020)033 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 spaces for the co-annihilation scenario are located at the 60 GeV m S 66 GeV, 7.4 GeV ∆ 0 8.7 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 −10 −3 δλ 0.0 when the contribution from ∆ ± dominates the δλ.
Next, the correlation between λ S and λ 2 is non-trivial at the NLO level. For λ 2 > 1, 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. Therefore, the cross section σ SI p for co-annihilation scenario can be significantly enhanced but ruled out by XENON1T. It becomes testable in the future direct detection experiments.
Finally, the 95% allowed region predicted in the model 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 figure 3.
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,  Table 2. The relevant equation numbers to calculate Feynman diagrams with effective interaction 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 the W boson contribution inside the one-loop box type diagram, for example.
The higher twist gluon operators have been neglected in the above effective Lagrangian. Based on eq. (A.1), the scattering amplitude and spin-independent cross section of dark matter and nucleon can be written as, (2)) , (A.3) 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. (A.3) and can be given as, At the next leading order, we closely follow the calculations in ref. [21] with the following classifications, • 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 2.
Here (B.5) means equation (B.5) in ref. [21] with the 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 2 with the following form, 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 [43]. Therefore, the loop-induced effects can be simply written as the form in eq. (3.6) with four input parameters: m S , ∆ 0 , ∆ ± and λ 2 inside δλ.

B Supplemental figures
In this appendix, we show some plots which are useful to understand the details of NLO effects in the parameter space.
The 2σ distribution for (∆ 0 , ∆ ± ) planes are shown in figure 8. The black dashed line represents 2σ contour for the LO calculation while the scatter points correspond to the one taking into account the NLO effects. From the left panel of figure 8, one can see that the NLO contribution takes effect only to the small ∆ 0 or ∆ ± region, but it does not change the LO result significantly.
In figure 9, 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 figure 9, we discuss the cases of λ S > 0 (upper-left panel) and λ S < 0 (upper-right panel) on the (λ 2 , δλ) plane, separately. 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 an opposite sign of δλ, while a small value of |λ S | can either accompany with a same sign 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 figure 9. The color codes indicate the value of ∆ ± (lower-left panel) and the JHEP06(2020)033 . The 2σ allowed region on the plane (λ 2 , δλ). Here, σ SI p is calculated at next leading order.
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 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 co-annihilation 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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.