Strong Dark Matter Self-Interaction from a Stable Scalar Mediator

In face of the small-scale structure problems of the collisionless cold dark matter (DM) paradigm, a popular remedy is to introduce a strong DM self-interaction which can be generated nonperturbatively by a MeV-scale light mediator. However, if such the mediator is unstable and decays into SM particles, the model is severely constrained by the DM direct and indirect detection experiments. In the present paper, we study a model of a self-interacting fermionic DM, endowed with a light stable scalar mediator. In this model, the DM relic abundance is dominated by the fermionic DM particle which is generated mainly via the freeze-out of its annihilations to the stable mediator. Since this channel is invisible, the DM indirect detection constraints should be greatly relaxed. Furthermore, the direct detection signals are suppressed to an unobservable level since fermionic DM scatterings with a nucleon appear at one-loop level. By further studying the bounds from the CMB and BBN on the visible channels involving the dark sector, we show that there is a large parameter space which can generate appropriate DM self-interactions at dwarf galaxy scales, while remaining compatible with other experimental constraints. ar X iv :1 91 0. 01 23 8v 3 [ he pph ] 8 O ct 2 01 9


Introduction
Despite great theoretical and experimental efforts in searching for dark matter (DM) during the last several decades [1][2][3][4], we still know little about the nature of DM. Until very recently, the most popular DM candidate has been the collisionless cold DM. However, by comparing the N -body simulations of such DM with astrophysical observations, a number of discrepancies have been found at the scale of dwarf galaxies, such as the cusp-vs-core problem [5][6][7][8] and the too-big-to-fail problem [9,10]. One intriguing possibility to solve these small-scale structure problems is to introduce large enough DM self-interaction [11][12][13][14][15][16][17][18], even though there are also some other more conventional solutions [19,20]. Concretely, the desired DM self-interaction should satisfy 0.1 cm 2 /g < σ T /m DM < 10 cm 2 /g, where σ T denotes the so-called momentum transfer cross section and m DM is the mass of the DM particle. On the other hand, DM self-interaction is severely constrained by observations at the galaxy cluster scale to be in the range σ T /m DM < 1 cm 2 /g [21][22][23][24][25][26]. Note that one important difference between DM particles in these two systems are their average velocities v, with v = 30 km/s in dwarfs and v = 1000 km/s in clusters. Therefore, the observations favor a velocity-dependent DM self-interaction cross section [17].
Note that in this scenario the DM self-interaction cross section increases as the DM velocity becomes small, which is helpful to evade the aforementioned cluster-scale constraint. Furthermore, the observed DM relic density can be naturally obtained by freeze-out of annihilations of DM particles into light mediators [37][38][39][40]. However, was the mediator unstable and decaying into SM particles, this scenario would face constraints of direct and indirect DM searches. Simplest realizations of the above scenario are models with weakscale fermionic DM and MeV-scale vector or scalar mediator. For the case with a vector mediator, it was shown in Ref. [41,42] that the model was ruled out by indirect detection constraints from Cosmic Microwave Background (CMB) [43][44][45][46][47], Fermi-LAT [48] and AMS-02 [49][50][51][52][53]. On the other hand, in the model with a scalar mediator, even though the DM indirect detection constraints can be avoided since the fermionic DM annihilation is p-wave dominated [54], the DM direct detection upper bounds [55] imply the longevity of the scalar mediator, which would modify the primordial abundances of light elements during Big-Bang nucleosynthesis (BBN) [56][57][58]. Though, there have already been many possible solutions proposed to avoid such tensions for models with scalar and vector mediators [54,[59][60][61][62][63][64][65][66][67][68][69].
In the present paper, inspired by Refs. [65,68,69], we construct and study a DM model in which the dominant DM component is a fermionic particle with a stable light scalar mediator. The stable mediator constitutes a subdominant DM by the freeze-out of its annihilation into additional light scalar particles. In this model, the observed DM relic density is mainly obtained by the freeze-out of the annihilation of a fermionic DM pair into the scalar mediators, which cannot be probed by DM indirect detection. Moreover, the fermionic DM is also free from any DM direct search constraints, since its nuclear recoils appear at one-loop level. Thus, the model is expected to be less constrained compared to the counterpart with an unstable scalar mediator, and has the potential to reconcile the aforementioned conflict among different experiments. However, as will be shown below, this model might be experimentally tested and constrained by observations of BBN [70,71] and CMB [44][45][46][47], since the processes involving dark-sector particles can still leave their tracks in modifications of the primordial abundances of light elements and it can change the CMB power spectrum. Therefore, the main question in the following is whether we can find the parameter space which can produce the desired strong DM self-interactions, while still be consistent with the current experimental bounds.
The paper is organized as follows. In Sec. 2, we briefly introduce our model and clarify notation and conventions. Sec. 3 is devoted to the calculation of the DM relic density. In Sec. 4, the constraint from the SM-like Higgs boson invisible decay is investigated. The constraints from CMB and BBN on the dark Higgs properties are discussed in Sec. 5. In Sec. 6, we show the analytic expressions of nuclear recoil cross sections for both DM components in the limit of zero momentum transfer, showing that they are effectively invisible under current experimental status. In Sec. 7, we discuss DM indirect detection constraints for visible channels involving dark sector particles. Discussion of the calculation of DM self-interactions is given in Sec. 8. Then we show our numerical studies in Sec. 9. Finally, we present our conclusions in Sec. 10.

The Model
Our model is a simple extension of the SM by including a Dirac fermion χ and two real scalars S and φ. We firstly impose the following stabilizing symmetry on the model [72]: with other fields being neutral. The symmetry implies the following new scalar potential terms and the following terms involving χ Since we assume S = 0, the stabilizing symmetry remains unbroken. Note also that the above Lagrangian is invariant under its subgroup of this symmetry: χ → −χ [72]. Therefore, when mass of the scalar S is below half of the Dirac fermion χ mass, both particles are stable and can contribute to the final DM relic density. Since our intention is to enhance the fermionic DM self-interactions by an exchange of the scalar S, we will consider the case with χ much heavier than S. Concretely, χ is assumed to be at the electroweak scale with its mass m χ ∼ O(1 GeV ÷ 100 GeV), while the mass of S will be of O(MeV). However, for such a small mass of S, the annihilation of S into SM particles are typically inefficient to reduce its relic abundance to be subdominant. Thus, we need to introduce an additional real scalar φ so that there is an extra annihilation channel for S to deplete its abundance. In order to avoid unnecessary parameters in the Lagrangian, we impose an extra Z 2 symmetry: φ → −φ, which implies the following potential terms Note that the latter Z 2 symmetry would be broken by the VEV of φ = v φ so that the perturbation ϕ = φ − v φ can mix with the neutral component h of the Higgs doublet which is defined as H ≡ (0, (v H + h)/ √ 2) T in the unitary gauge with v H = 246 GeV. By minimizing the total scalar potential, we can determine non-zero VEVs of H and φ solving the following two equations: We further expand the potential written in terms of perturbation fields S, h, and ϕ up to the second order, determining the mass squared of S as m 2 S = µ 2 S + (κ HS v 2 H + κ Sφ v 2 φ )/2 and the following mass squared matrix for h and ϕ: (2.6) We can diagonalize the matrix above defining the mass eigenstate h 1,2 in terms of the perturbations and mixing angle θ as follows where s θ ≡ sin θ and c θ ≡ cos θ. Denoting the masses of h 1,2 by m 1,2 , the following relations between parameters hold If we assume that h 1 is the SM-like Higgs state, then m 1 = 125 GeV and the model can be parametrized by the following 9 free parameters v φ , m 2 , m S , m χ , s θ , λ S , κ HS , κ Sφ , g Y . (2.9)

Dark Matter Relic Density
In the present model, there are two stable particles: the heavy DM candidate χ and the stable mediator S, both of which will contribute to the final DM relic abundance. As mentioned before, the stable mediator S should be hierarchically lighter than the fermionic DM χ, i.e., m S m χ , so that the self-interactions of χ can be non-perturbatively enhanced to the level large enough to solve the small-scale structure problems. In this section, we would like to discuss how the DM relic density can be obtained via the freeze-out mechanism in this scenario.
The number densities of the two DM components, n χ and n S , can be yielded by solving the following two coupled Boltzmann equations: where n χ is the total number density of both χ andχ, while n S that of the light scalar mediator S. On the right-hand side of both equations, we have only shown the most important processes for the determination of the DM relic densities, including the dominant annihilations for both components (χχ → SS for χ and SS → h 2 h 2 for S), the semiannihilation (SA) process χχ(χχ) → Sh 2 , and the conversion process χS →χh 2 (χS → χh 2 ) [72]. The corresponding Feynman diagrams are listed in Fig. 1. Other processes are suppressed either by the small mixing angle s θ like χS →χh 1 or by the loop factors such as the χχ annihilations into SM particles. Note that when the freeze-out of the heavy fermionic DM χ takes place, cross sections for χχ, χχ andχχ annihilations are dominated by the p-wave contribution [54]. By further taking into the multiple exchange effects of the light mediator S, we should also multiply them by the following p-wave Sommerfeld enhancement factor [33,[73][74][75][76] where the factor S s is the Sommerfeld factor for s-wave annihilations with its explicit form given by .
Here, we have defined a ≡ v/(2α X ) and c ≡ 6α X m χ /(π 2 m S ), with α X ≡ g 2 Y /(4π) and v is the relative velocity between χ andχ. However, we find that this effect is so small that it can be neglected in the final prediction of relic densities of both DM components.
Furthermore, as will be discussed in Sec. 5, the dark sector should decouple from the visible SM sector at least before the QCD phase transition at T QCD ≈ 100 ÷ 200 MeV, here after we will adopt T dec ≈ 500 MeV. After decoupling the two sectors evolve independently with their own temperatures [69]. Under the assumption of the entropy conservation in each sector, the temperature ratio between the dark and SM sectors can be estimated as a function of the SM plasma temperature T where g * S (T ) and g * S (T D ) are the relativistic degrees of freedom relevant for the entropy in the visible and dark sectors, respectively. In the present model with the desired DM masses, the stable scalar mediator S always freezes out after the decoupling, while the relic density of the fermionic DM χ can be also affected as long as its mass is lighter than 12 GeV [69]. In order to account for the impact of this thermal decoupling of the dark sector on the DM relic abundances quantitatively, we follow Ref. [37] simply multiplying the relic abundances of S and/or χ obtained assuming equal temperatures in both sectors by the correction factor ξ(T f ), where T f denotes the freeze-out temperatures of the corresponding stable particles. In our work, we numerically solve the coupled Boltzmann equations in Eq. (3.1) using the modified MicrOMEGAs v4.3.5 code [77,78] which takes into account the aforementioned Sommerfeld enhancement and early decoupling effects. As noted in Ref. [69], when calculating the final relic abundance of the subdominant component S, the χχ annihilation into a pair of S and the semi-annihilation process play crucial roles, even though they already froze out and cannot further modify the abundances of χ andχ.

The Constraint from Higgs Boson Invisible Decay
Since the light scalar mediator S is assumed to be below 1 GeV, its interactions are subject of constraints from the SM-like Higgs boson invisible decay. The SS final state 1 contributes to the total invisible width as follows where in the second equality we have assumed that m S m 1 and s θ 1. In order to satisfy the current LHC upper limit Br inv < 0.24 [1] one has to assume that κ HS < 1.65 × 10 −2 , which implies that the annihilation of S into the SM particles is insufficient to deplete the stable mediator S when it freezes out. Therefore, S freezes-out mainly via the annihilation SS → h 2 h 2 , and does not depend much on the value κ HS . In the following, in order to satisfy the constraint from h 1 invisible decays we adopt κ HS = 10 −2 .

Constraints on h 2 Decays from CMB and BBN
Since the dark sector is assumed to be in thermal equilibrium with the SM particles at early times, the dark Higgs boson h 2 would exist abundantly in the early Universe. In the following, we are going to consider light dark Higgs boson h 2 with O(MeV) mass, so the only decay products are e + e − and γγ. In particular, when m 2 > 2m e , h 2 mainly decays into e + e − , while for the h 2 mass below this threshold the γγ mode induced at one-loop order would become dominant. However, these two channels are strongly constrained by the CMB observations. If the respective decays are effective at redshifts z 2 × 10 6 , the produced electromagnetic energy cannot fully thermalize with the background plasma and thus induces the spectral distortion in the CMB [79][80][81]. It is shown in Ref. [47] that this constraint can be transformed into the limit on the h 2 lifetime τ h 2 < 10 5 s. This indicates that we can only have the e + e − as the dominant decay channel of h 2 , since one-loop induced γγ channel always gives the lifetime much larger than this bound. Therefore, we take the h 2 mass to be m 2 = 1.5 MeV as in Ref. [69]. The decay rate of h 2 → e + e − is given by (5.1) The limit on the lifetime implies s θ > 2.55 × 10 −7 . Hereafter, we fix s θ = 5 × 10 −7 for simplicity. The presence of h 2 also affects the yields of light elements predicted by the theory of BBN, with the following two possible effects [1,82]: (i) When the lifetime of h 2 exceeds 10 4 s, the electromagnetic products of h 2 decays can destroy the formation of light nuclei, such as deuterium and helium, via the photo-disintegration process.
(ii) The states h 2 and S at O(MeV) mass might be in thermal equilibrium with the SM sector particles until temperature T 10 MeV and, as an extra relativistic degrees of freedom (dofs), may boost the cosmological expansion rate before BBN. This contribution of dark radiation to the energy density in the Universe is conventionally parametrized by the effective number of neutrino species ∆N eff .
It was argued in Ref. [69] that the constraint (i) can be easily evaded for light h 2 , e.g. for m 2 = 1.5 MeV, since the electromagnetic decay products, mainly the e + e − pairs, cannot generate photons whose energy could reach the photo-disintegration threshold E dis = 2.2 MeV for deuterium [83,84].
As for the constraint (ii), the current upper bound on the dark radiation during the BBN is set to be ∆N eff < 0.2(0.36) at the 2σ (3σ) level by exploiting the updated baryonto-photon ratio from the CMB and the latest nuclear reaction rates [1,82]. If h 2 and/or S remain in equilibrium with the SM plasma during BBN, this constraint is obviously violated. The simplest way to avoid the dark radiation bound is to decouple the dark and visible sectors at earlier times well before the temperature T 10 MeV so that the dark sector is much cooler than the SM one. It could be shown that it is enough to decouple the dark sector before the QCD phase transition occurring around T QCD 100 ÷ 200 MeV, during which the SM relativistic dofs are greatly reduced due to the confinement of light quarks and gluons in QCD. Quantitatively, if both h 2 and S are lighter than 10 MeV and are decoupled from the SM sector at T dec 500 MeV, their contributions to the dark radiation in terms of the equivalent neutrino numbers can be calculated as follows where we have taken the relativistic dofs in the SM before the QCD phase transition to be g SM (T dec ) = g SM (500 MeV) = 61.75 and g SM (10 MeV) = 10.75, as usual. If h 2 is the only relativistic dof in the dark sector, ∆N eff would be reduced by half so that the constraint would be even weaker. Hence, no matter what particle content is in the dark sector, the dark radiation constraint can be easily satisfied as long as the SM and dark sectors cease to be in the thermal contact before QCD phase transition. Note that the thermal decoupling between the visible and dark sectors is mostly controlled by the total rate of the h 2 pair annihilation into SM particles Γ h 2 h 2 →SM SM , which determines the chemical equilibrium between the two sectors. The decoupling happens while Γ h 2 h 2 →SM SM (T dec ) = H(T dec ) [85], where H(T dec ) denotes the Hubble expansion rate at the decoupling temperature T dec . Below T dec , the h 2 annihilation cannot catch up with the cosmic expansion so that it is effectively turned off. It is shown in Fig. 1 of Ref. [69] that when the Higgs portal coupling κ Hφ obeys the condition κ Hφ ≤ 4 × 10 −4 , the decoupling of the dark sector particles from the SM plasma can be achieved prior to the QCD phase transition. In the following, we will take κ Hφ = 4 × 10 −4 as in Ref. [69] in order to maximize the allowed parameter space. Together with the previously fixed s θ = 5 × 10 −7 , we can obtain the VEV of φ as v φ = 79.4 MeV.

Dark Matter Direct Detection
As mentioned before, in the strongly self-interacting fermionic DM model with an unstable O(MeV) scalar mediator, the upper bound from the current DM direct detection experiments requires a too long mediator's lifetime that is not consistent with the BBN [56][57][58]. This problem is one of the main motivations to study the present model with a stable scalar mediator. Note that both χ and S contribute to the observed DM relic density. Thus, their scatterings with nucleons would cause recoils observable in DM underground detectors, which could severely constrain the model. Note that both the scalar-nucleon SN and fermion-nucleon χN scatterings are mediated by the two Higgs particles h 1,2 . As h 2 is assumed to be light, we need to take into account its light mediator effect in the DM-nucleon interactions [86,87]. The scattering of the fermionic DM χ against nucleon is induced at one-loop order as shown in Fig. 2(a). The nuclear recoil cross section is given by with the effective Higgs-nucleon coupling f N 0.3 [85,88,89]. Here the factor (4µ 2 χN v 2 + m 2 2 ) in the denominator represents the possible light h 2 mediator effects [86,87], where v ≈ 220 km/s denotes the typical DM velocity in the Milky Way relative to the Solar system, and the loop fucntion F (t) is defined as follows which is valid when 0 < t < 4. If we take the benchmark values for the model parameters as g Y ∼ 0.1, m χ ∼ 100 GeV, and m S ∼ 10 MeV, the above formula gives the χ-nucleon cross section to be of O(10 −48 cm 2 ), which is well below the best experimental bound from XENON1T experiment [55]. Hence, we can ignore the DM direct detection constraints on the heavy DM χ properties. It turns out that the contribution from the stable mediator S is also subdominant. Its tree-level scattering cross section on nucleon (see Fig. 2 with µ SN ≡ m S m N /(m S + m N ) being the reduced mass in the S-N system. In the above formula we have ignored the momentum transfer for h 1 and even for h 2 . Note that the mass of S is expected to be of O(1 ∼ 100 MeV). In this range, the nuclear recoil energy is well below the lowest thresholds obtained in the direct detection experiments. This problem can be overcome if one considers electron recoils [90,91], however in the discussed model the Higgs-mediated DM-electron interactions are strongly suppressed. Furthermore, it is also possible to utilize the highly energetic mediators S arising from χ annihilations [92,93] or collisions with cosmic rays [94], but these bounds are to weak to constrain the model. Concluding, the current DM direct searches cannot put any useful constraints from either the heavy DM χ or the light stable mediator S scattering. Therefore we are able to avoid the conflict faced by the self-interacting fermionic DM models with the unstable mediator investigated in Refs. [56][57][58].

Dark Matter Indirect Detections
It is well known that one of the most relevant difficulties to construct a viable model with sufficient DM self-interactions originates from DM indirect detection constraints [4]. This is because almost the same intermediate DM bound states that generate the DM selfinteractions would also induce the Sommerfeld enhancement of the DM annihilation [41,42], which is the main target for DM indirect searches. At present, the most important DM indirect detection constraints comes from the Planck measurement of CMB spectral distortion [43], a gamma-ray probe of the dwarf spheroidal galaxies by Fermi-LAT [48], and the AMS-02 search of high energy positrons in the Milky Way [49,50]. As shown in Refs. [68,69], if the light vector mediator which helps to form the DM bound states is stable, such constraints can be efficiently avoided, since most annihilations of dominant DM particles go into invisible light mediators. A similar mechanism takes place in the present model with a stable light scalar mediator, since the dominant fermionic DM density is provided by the freeze-out of its unobservable annihilation channel χχ → SS. Nevertheless, there are still several channels in the dark sector which could be probed by the indirect detections.
One important channel is the semi-annihilation process χχ(χχ) → Sh 1,2 followed by the h 1,2 decays. However, this process can be proven to be dominated by the p-wave cross section, which approaches zero in the non-relativistic limit with the relative velocity v rel → 0, in spite of a potential significant p-wave Sommerfeld enhancement. We have explicitly computed the cross sections with the DM velocities of dwarf galaxies and the Milky Way, and compared the results with the constraints given by AMS-02 and Fermi-LAT [95]. As a result, this process cannot impose any relevant constraints. This conclusion agrees with that given in Ref. [54].
Another possible signal stems from the stable scalar mediator S annihilation into the visible particles via the Higgs portal interactions. Since the coupling of S with the Higgs doublet is constrained to be very small as it has been discussed in Sec. 4, its dominant annihilation goes into the h 2 pairs which subsequently decay into the electrons and positrons. Since S is assumed to be lighter than 1 GeV, the relevant constraint can only be given by the CMB observation, while the energy of produced e + /e − is always too small to be constrained by Fermi-LAT or AMS-02. The explicit upper bound from CMB is given as follows where σv rel 4e ± CMB (m S ) represents the upper bound given in Ref. [45] on the S annihilation cross section with two pairs of electrons and positrons in the final state. The factor (Ω S h 2 /Ω DM h 2 ) 2 accounts for the suppression from the S energy density fraction, with Ω DM h 2 0.12 the total measured DM relic density [1]. Moreover, the model is further constrained by the DM indirect searches for the process χS →χh 2 andχS → χh 2 , in which the signal also originates from the h 2 decays. For simplicity, we will use χS →χh 2 to denote both processes in the following formulas and plots. In the non-relativistic limit relevant to the DM indirect searches, the energy of the produced h 2 can be estimated to be the mass of S. By considering the suppression from respective density fractions of the fermionic DM and the stable mediator, we can write down the following constraint for this process where the factor 1/2 on the left hand side accounts for the fact that only one pair of e ± is generated in this reaction. Note that in contrast to the annihilation of fermionic DM particles, the above two processes are not subject to the Sommerfeld enhancement, and the corresponding cross sections are s-wave dominated so that they approach non-zero values at the non-relativstic limit v rel → 0.

Dark Matter Self-Interactions
The main motivation for the introduction of the DM self-interactions is to solve the smallscale structure problems in our Universe, such as the cusp-vs-core problem and the toobig-to-fail problem at the scale of the dwarf galaxies. By fitting the data, the required DM self-interaction per unit DM mass should be 0.1 cm 2 /g < σ T /m χ < 10 cm 2 /g for dwarf galaxies with the typical DM velocity at v 30 km/s [11][12][13][14][15][16][17][18], where σ T denotes the momentum transfer cross sections in the DM scatterings. On the other hand, such a strong DM self-scattering can also leads to the observable effects at the galaxy cluster scale with v 1000 km/s. The absence of these effects in the the galaxy cluster data strongly constrains the DM self-interactions at this cosmological scale, with the conservative upper limit given by σ T /m χ < 1 cm 2 /g [21][22][23][24][25].
With the mass hierarchy between the scalar S and the Dirac fermion χ, the DM selfinteractions can be greatly enhanced by the formation of the χχ or χχ bound states via the mediation of S [28,29,32,32,33,33]. Since the incoming and outgoing states are different, it is difficult in analyzing the DM self-interactions in this form of the Yukawa couplings. In order to overcome this problem, we define the following two Majorana fermions With this transformation of fermion basis, the corresponding Lagrangian is modified to Accordingly, the original two Z 2 symmetries are transformed into the following three ones Thus, the scalar S and two Majorana fermions χ ± are all DM candidates. With the Lagrangian in Eq. (8.2), the Yukawa terms now are diagonalized and the incoming and outgoing fermions keep same, which simplifies our following discussion. In particular, it is easy to see that the potential in the χ + χ + and χ − χ − systems are attractive, given by while χ + χ − scatterings are repulsive with the following potential: Accordingly, the momentum transfer cross section σ T , which characterizes the DM self-interactions, is defined in the present model as follows [54] where σ ±± T and σ ±∓ T are the cross sections for the χ ± χ ± and χ ± χ ∓ scatterings, respectively, and, as argued before, correspond to the ones induced by the attractive and repulsive potentials, i.e., σ a T and σ r T . In our computation of the momentum transfer cross sections, we follow Ref. [33,54] to numerically solve the Schrödinger equations with the potentials given in Eqs. (8.4) and (8.5), which effectively sums over multiple exchanges of S in the interactions between two heavy fermionic DM particles as illustrated in Fig. 3. This is particularly useful when we work in the parameter space of the so-called quantum resonant regime in which the attractive potential leads to the formation of resonances of χ bound states. In our calculation, we also take into account effects arising from the quantum indistinguishability of identical particles in the initial and final states following the procedure given in Ref. [54].
Note that when g 2 m χ /(4πm S ) 1, the non-perturbative momentum transfer cross sections would reduce to the analytical formula given in Ref. [54] for the Born approximation. In another limit where m χ v/m S > 10, we are entering the classical regime in which it is very difficult to solve the Schrödinger equations and the approximate expressions for the momentum transfer cross sections given in Ref. [34] can be applied. In our following numerical calculations, we shall adopt both analytical expressions in these two parameter regimes.

Numerical Results
In this section, we present our numerical results in relevant regions of the parameter space. We begin by presenting the results in Fig. 4 spanned by the portal coupling κ Sφ and the light stable mediator mass m S for different values of the fermionic DM mass m χ = 2, 10, 100, 1000 GeV from top left to bottom right, respectively. Following the discussion in Sec. 5, we have fixed the mass of the dark Higgs boson h 2 to be m h 2 = 1.5 MeV to avoid the CMB and BBN constraints for the late-time h 2 decays as much as possible. Note that to make the annihilation channel SS → h 2 h 2 kinematically allowed in order to deplete the mass density of S, we choose the lower limit of S mass to be m S = 2 MeV in all panels of Fig. 4. Then we fix κ Hφ = 4×10 −4 , κ HS = 10 −2 , s θ = 5×10 −7 , and v φ = 79.4 MeV, so that we are left with only four free parameters. Therefore, the Yukawa coupling g Y in Eq. (2.3) can be determined for each set of given values of (m χ , m S , κ Sφ ) by the requirement that χ and S constitute all of the measured DM relic density with Ω DM h 2 = Ω χ h 2 + Ω S h 2 0.12. We also show, in Fig. 4, the contours of constant relic density fraction Ω S /Ω DM as black dashed lines with the three curves in each panel corresponding to Ω S /Ω DM = 0.1, 10 −4 , and 10 −7 , from left to right respectively. In Fig. 4, this fraction is found to rise with decreasing portal coupling κ Sφ and increasing m S , what implies that the cross section of SS → h 2 h 2 process is reduced. At some point, the cross section of SS → h 2 h 2 becomes too small to eliminate enough S so that, no matter what value of g we choose, the relic density of S alone would overclose the Universe. This part of parameter space should be excluded as the orange shaded regions indicate.
The main purpose of this paper is to provide a viable DM model in which the strong DM self-interactions can resolve the structure problems at the dwarf galaxy scale. In Fig. 4, we show the light and dark green regions of the parameter space where the momentum transfer cross section of χ is generated in the ranges of 0.1 cm 2 /g < σ T /m χ < 1 cm 2 /g and 1 cm 2 /g < σ T /m χ < 10 cm 2 /g, respectively. We also display the bound on the DM self-interactions σ T /m χ 1 cm 2 /g from the galaxy cluster scale as the yellow shaded area. Note that this constraint is relevant only when the dominant DM particle χ is as light as m χ = 2 GeV. For heavier χ, it is always satisfied for the whole parameter regions of interest. Moreover, when m χ = 100 GeV, the region which can explain the cosmological  In each panel the second (dark) Higgs boson mass is fixed to be m 2 = 1.5 MeV, while the dark Yukawa coupling g Y is obtained by requiring Ω χ h 2 + Ω S h 2 0.12 at each point. The orange shaded areas are excluded, because DM relic density exceeds the measured value for any chosen g Y . The light and dark green regions represent the parameter space which can generate a momentum transfer cross section σ T of the fermionic DM self-interaction at the scale of dwarf galaxies in the range of 0.1 cm 2 /g < σ T /m χ < 1 cm 2 /g and 1 cm 2 /g < σ T /m χ < 10 cm 2 /g, respectively, while the yellow region is excluded by the bound σ T /m χ 1 cm 2 /g on the scale of galaxy clusters. The red and blue shaded regions indicate parameters which are ruled out by the CMB constraints on the late energy injections from the processes SS → h 2 h 2 and χS →χh 2 (χS → χh 2 ).
small-scale problems appears in the form of many spikes that are a characteristic feature of the quantum resonant regime discussed in Sec. 8. On the other hand, the relevant parameter regions for light fermionic DM masses m χ = 2 and 10 GeV correspond to the Born regime, while when χ becomes as heavy as O(1000 GeV), we enter the non-perturbative classical regime to generate large enough DM self-interactions.
Furthermore, it is interesting to see that in all of the panels in Fig. 4, the signal regions with Ω S 0.1Ω DM are represented as straight bands, indicating that the DM self-interaction cross section is insensitive to κ Sφ . This fact can be understood as follows.
When Ω S 0.1Ω DM , the observed DM relic density is dominated by that of the fermionic DM, with its annihilation channel χχ → SS. Thus, with a fixed fermionic DM mass in each plot, the Yukawa coupling g Y remains invariant as a consequence of a nearly constant σv χχ→SS , so the fermionic DM self-interaction is only the function of m S . On the other hand, when S occupies a fraction 0.1 of the DM relic abundance, the fermionic DM abundance should be reduced and the value of g Y varies accordingly. This is reflected in the lower left panel of Fig. 4 by the bands which turn upward as they go into the parameter regions with Ω S 0.1Ω DM .
We also show the constraints from DM indirect detection in Fig. 4. As discussed in Sec. 7, the constraints originates mainly from the CMB upper bounds on energy injections during the recombination era due to the processes SS → h 2 h 2 and χS →χh 2 (χS → χh 2 ), which are represented as red and blue shaded areas in Fig. 4. From these four plots, it is easy to see that the annihilation of SS into a pair of h 2 strongly constrains the parameter space of a large S relic density fraction Ω S /Ω DM > 0.1. On the other hand, the energy injection induced by the process χS →χh 2 can exclude additional regions. Firstly, it is found that the allowed regions in each plot of Fig. 4 have upper boundaries, which reflect the fact that the cross section of this process is suppressed too much when m S goes up beyond the boundary. Secondly, below this S mass limit, this process excludes more parameter space than SS → h 2 h 2 , extending the exclusion region deeply into the part with even smaller S density fraction. This can be understood as follows: the scattering rate between S and χ(χ) is proportional only to the first power of the density fraction of S, while the rate for the S self-annihilation is further suppressed by the two powers. Lastly, when the mass of χ becomes to be of O(100 ∼ 1000 GeV), the region with lighter S is excluded by the χS →χh 2 signal, which is the result of large coupling κ Sφ that overcompensate the suppressed S density when computing the rate for this process.
It is clearly seen from Fig. 4 that, when the dominant fermionic DM mass is below several hundred GeV, there is always a large parameter space which can accommodate DM relic density and explain the cosmological small-scale problems while satisfying the constraints of DM indirect detections. For the remaining parameter regions, the DM density fraction from the stable light mediator S is always subdominant, which is constrained by the late-time electromagnetic energy injection during the recombination from the processes SS → h 2 h 2 and χS →χh 2 . However, in the bottom right panel of Fig. 4, the required DM self-interaction contours are excluded by the above two processes during the CMB formation, therefore this model disfavors the fermionic DM candidate as heavy as several TeV.
We also present the same data points in the m S -m χ planes in Fig. 5, with the four panels corresponding to the cases with κ Sφ = 10 −4 , 10 −3 , 10 −2 , and 10 −1 from upper left to lower right. The color coding in each panel is the same as that in Fig. 4. From the upper left panel with κ Sφ = 10 −4 , it follows that for a sufficiently small κ Sφ O(10 −4 ), most of the parameter space which can give rise to the interesting DM self-interaction cross section at the scale of dwarf galaxies is excluded by the combinations of constraints from DM selfinteractions at galaxy clusters and from CMB observations. Only the corner with a light fermionic DM of m χ < 2 GeV and a small mediator mass m S < 4 MeV is allowed by the data. On the other hand, with the increase of κ Sφ , the upper bounds from energy injections at recombination becomes more and more relaxed, and larger part of the parameter space is allowed to account for the small-scale structure problems. It is seen from Fig. 5 that the CMB constraint from SS → h 2 h 2 is sensitive to the S density fraction and excludes at least regions with Ω S /Ω DM > 0.1, while the process χS →χh 2 implies a strong limit for the region with large fermionic DM mass m χ . Especially, the latter excludes all DM self-interaction signal regions with m χ > 1 TeV, which is consistent with the results shown in Fig. 4.

Conclusions
Sufficiently strong DM self-interactions provide a possible solution to the small-scale structure problems arising in the collisionless cold DM paradigm. In order to generate such a large DM self-scattering, one popular strategy is to introduce a light mediator to enhance the cross section nonperturbatively. Furthermore, this leads to the velocity-dependent DM self-interactions which help to evade the strong constraint at the scale of galaxy clusters. The simplest realization of this scenario is the weak-scale fermionic DM model with a MeV-scale light scalar or vector mediator, in which the DM relic density is controlled by the freeze-out of the DM annihilation into mediators. However, in the case with a light scalar mediator which can decay into visible SM particles like photons or electrons, the constraints from the DM relic abundance, DM direct detections and DM indirect searches strongly limit this scenario. In order to reconcile this conflict, inspired by Ref. [65,68,69], we have studied a model with a dominant fermionic DM candidate χ and a stable light mediator S, where the density of the latter state is depleted due to its efficient annihilation into a new scalar particle h 2 . Consequently, the model has two advantages: first, the dominant DM χ relic density can be obtained by the conventional freeze-out mechanism where the dominant annihilation channel χχ → SS becomes invisible and can escape the DM indirect searches. Secondly, the model is even free from DM direct detection constraints since scattering of the main DM component χ with a nucleon appears at one-loop level and the light mediator is too light to have any observable signals.
Following Ref. [69], we have fixed κ HS = 10 −2 by imposing experimental constraints from the SM-like Higgs boson invisible decays, and set s θ = 5 × 10 −7 , m 2 = 1.5 MeV, κ Hφ = 4 × 10 −4 and v φ = 79.4 MeV to avoid the limits on the properties of the dark Higgs h 2 from the CMB spectral distortion and BBN observations. Note that even though the main annihilation channel for the freeze-out of the dominant DM particle χ is invisible, the model is still constrained by the CMB bounds on late-time electromagnetic energy injection during recombination from other visible processes such as SS → h 2 h 2 and χS →χh 2 (χS → χh 2 ) followed by h 2 decays. Imposing all experimental constraints and scanning relevant regions of the parameter space, we have found (see Figs. 4 and 5) that there remains a large parameter region which can accommodate DM self-interactions at the level sufficient to explain the structure problems at dwarf galaxies scale while agreeing with bounds from galaxy clusters. In particular, it is found that the fraction of S in the DM relic abundance is constrained to be smaller than 0.1 in this model. Also, in order to explain the cosmological small-scale structure problems, the fermionic DM mass should be lighter than O(TeV), while a relatively large coupling κ Sφ > 10 −4 is favored.