Thermal Relic of Self-Interacting Dark Matter with Retarded Decay of Mediator

The existence of a light mediator is beneficial to some phenomena in astroparticle physics, such as the core-cusp problem and diversity problem. It can decouple from Standard Model to avoid direct detection constraints, generally realized by retard decay of the mediator. Their out-of-equilibrium decay process changes the dark matter (DM) freeze-out via temperature discrepancy. This type of hidden sector (HS) typically requires a precision calculation of the freeze-out process considering HS temperature evolution and the thermal average of the cross-section. If the mediator is light sufficiently, we can not ignore the s-wave radiative bound state formation process from the perspective of CMB ionization and Sommerfeld enhancement. We put large mass splitting between DM and mediator, different temperature evolution on the same theoretical footing, discussing the implication for DM relic density in this HS. We study this model and illustrate its property by considering the general Higgs-portal dark matter scenario, which includes all the relevant constraints and signals. It shows that the combination of BBN and CMB constraint favors the not-too-hot HS, $r_{\mathrm{inf}}<10^2$, for the positive cubic interaction of mediator scenario. On the other hand, the negative cubic interaction is ruled out except for our proposed blind spot scenario.


Introduction
Most of the dark matter (DM) community have considered the situation in which DM mass is comparable with mediator mass [1], e.g. 100 GeV neutralino with Higgs or W/Z bosons being the mediators. As a result, these processes belong to contact short-range interaction [1,2]. However, with the development of indirect detection and DM N-body simulation, long-range interaction attracts lots of attention such as Positron Excess [3]. Besides, although cold dark matter (CDM) in ΛCDM model is very successful at explaining the current observations on large scales [4], there are various discrepancies between N-body simulations of collisionless CDM and astrophysical observations on galactic scales. Self-interacting dark matter (SIDM) was thus proposed to solve the core-cusp [5,6] and missing satellites [7][8][9][10][11] problems on galactic and smaller scales. The energy transfer from self-interaction in the central regions of DM halos provides a method to heat DM particles so that an isothermal core can be created easily. In SIDM [12][13][14][15][16][17][18][19], the cross section within long-range interaction is velocity-dependent in order to escape bullet-cluster constraint at relatively high velocity and obtains large cross section at low velocity. Furthermore, Sommerfeld enhancement (SE) on the annihilation cross sections can significantly deplete DM relic density. This is mainly because DM pair annihilation occurs at small non-relativistic velocity at freeze-out temperature.
The non-observational evidence of WIMP in both direct detection and indirect detection challenges the traditional form of SIDM, which leads to the study of viable generalizations of SIDM to overcome stringent limitations such as pseudo-Dirac SIDM [20]. The basic strategy for model building is to find a mechanism that realizes correct DM abundance in the early universe without causing a large elastic scattering rate with nucleon. Hidden sector (HS) DM [21][22][23][24][25] is one of the most attractive options following the strategy, where DM freezes out of thermal equilibrium entirely within its HS.
In this paper, for simplicity, the decoupled HS only contains two particles: DM and mediator. The large mass hierarchy between DM and mediator offers long-range interaction. As a result, SIDM can be naturally realized in this class of models. In addition, radiative bound state formation (BSF) and SE are inevitable, which lead to large annihilation cross section. To avoid large elastic scattering between DM and nucleon, we hypothesize the coupling between the mediator and visible sector (VS) particles is ultra-tiny so that the HS is completely decoupled from the VS thermal bath. As a consequence, HS have its own temperature, which may also increase the annihilation cross section. The mediator provides a portal between the HS and VS, and can decay into VS particles and connect the two sectors. Unlike other models, we can achieve correct velocity dependence at late times without assuming particle-antiparticle asymmetry in the HS and can get the full relic abundance. The purpose of this paper is thus to perform a precise calculation on the HS relic density with large mass splitting by illustrating General Higgs Portal Dark Matter (GHPDM) as our benchmark model of the HS.
The structure of this paper is organized as follows: In Sec. 2 We derive the temperature ratio between the HS and the VS in terms of entropy conservation, and implement it to the Boltzmann equation. And then briefly review the general contents of our benchmark model, GHPDM, in which SE and formation of bound state can be realized. The constraints from BBN an CMB are also calculated. We give the semi-analytical formula to calculate the self-interaction cross-section in Sec. 3. A simple chi-squared analysis shows that mass ratio m φ /m mχ smaller than 10 −4 is favored. In Sec. 4, we show the analytical results of SE and radiative BSF cross-section. The implication of relevant DM properties such as relic abundance involving SE effect and BSF are shown explicitly. We finally conclude in Sec. 5.

Hidden sector temperature and Boltzmann equation
When DM and mediator as a system effectively decoupled from SM thermal bath, the two sectors undergo different temperature evolution, which is convenient to define the ratio of the HS and VS temperatures, r = T h /T . Here quantities with subscript or superscript h stand for HS. We assume both sectors are populated after inflation (reheating). Therefore the ratio of their initial temperature r inf = T h inf /T inf can be regarded as a free parameter and initial condition. With the expansion of the universe, the evolution of temperature is determined by entropy densities namely s = 2π 2 /45 g * s (T )T 3 and s h = 2π 2 /45 g h * s (T h ) T 3 h , where g * s (T ) and g h * s (T h ) are the effective relativistic entropy degrees-of-freedoms (d.o.f.) of VS and HS. Even though the temperatures of the two sectors are independent of each other, the conservation of comoving entropy densities, d(sa 3 )/dt = 0, determines the evolution of ratio as follows [24,26] r = g * s (T ) g * s,inf where g * s,inf and g h * s,inf are relativistic entropy d.o.f. at inflation for VS and HS. To solve for the temperature ratio r in Eq. (2.1), we need to present the evolution of the d.o.f. g * and set up the initial condition r inf . Different initial conditions lead to different DM phenomenology which will be discussed in Sec. 4.
Assume that the temperature at inflation is much higher than the mass of any particle in two sectors, then d.o.f. at inflation equal (g s + 7 8 g f ), where g s and g f are the intrinsic d.o.f. of scalars and fermions. For example, if there are one complex scalar φ and one Majorana fermion χ in HS, then g h * s,inf = 3.75. If VS only consists of SM particles, then g * s,inf 106.75. Evolution of g * s (T ) for SM particles is given in [27]. Assuming the chemical potential are negligible, effective entropy d.o.f. g * s as a function of T for a particle with a mass m and intrinsic d.o.f. g can be obtained by where z = m/T and positive sign in the formula stands for fermion and negative for scalar particle. The Eq. (2.1) can be solved after determining the variation of g * s (T ) and g h * s (T h ). In Fig. 1, we show the evolution of the r with respective to the x = m χ /T for the different choice of parameters. Here, as an example, we adopt minimal supersymmetric standard model (MSSM) with all supersymmetric particle masses of 50 TeV in the VS, and its d.o.f at inflation is g * s,inf 228.75. The first drop of r values around x = 10 −3 are the decoupling of supersymmetric particles. The ratio r in blue and orange line start to increase after x = 1, since DM becomes non-relativistic. In green line, m φ as heavy as DM mass m χ = 100 GeV, they became non-relativistic simultaneously. After the decoupling of lightest particles in HS, which is the mediator φ in our case, the comoving number density is conserved, equally phase space density f ∼ e −E/T h is constant. From E ∼ a −2 and T ∼ a −1 , we know that and it can be seen from the tail of the curves in Fig. 1.
Within the temperature ratio evolution r known, the Boltzmann equation describing HS χ, φ before entropy injection from mediator is straightforward to obtaiṅ If χ is not self-conjugated, we will replace σv by σv /2. We make some general comments on the Boltzmann equation: • When the DM mass is almost degenerate with mediator mass, Boltzmann equations in (2.4) reduce to Co-decaying DM [28]. However, Self-interacting DM requires the mediator φ much lighter than DM χ. The large hierarchy is the main difference between our model and Co-decaying DM. Due to the existence of long-range interaction, the cross-section σv χχ→φφ includes both SE and BSF effects.
• When the decay of φ into SM particles is prompt, the HS remains thermal equilibrium with the VS. In that sense, only one temperature in the Boltzmann equation T h = T appears. However, we assume the mediator decay is too tiny to have direct detection constraints. The HS and VS undergo different temperature evolution described by the ratio r. For the varying values of r, the corresponding solution becomes different.
Since the portal coupling between mediator φ and SM sector is ultra-tiny, mediator decay is retarded decay process. The out-of-equilibrium decay injects more entropy to the thermal bath and violates entropy conservation. Therefore our equation is only valid up to the time that the mediator decay process dominates. Furthermore, additional entropy injection with the early matter-dominated universe will dilute the pre-existing relic density of DM. In our setup, the dilution factor ∆ is not sizeable because of the small mediator mass.
• Since DM and mediator annihilate with each other efficiently to maintain the kinetic equilibrium, the number density of mediator n φ tracks with thermal equilibrium n φ,eq . We will calculate the elastic scattering process χφ → χφ in Sec. 4 to prove it. Then, the Boltzmann equation (2.4) reduces to conventional WIMP type equatioṅ There are some modifications that we need to consider in Eq. (2.6). One is that we should perform the thermal average of cross section σv in T h rather than T. Furthermore, thermal equilibrium is also modified to be evaluated at hidden temperature T h , where K 2 is the Bessel function of the second kind. Despite conventional WIMP scenario, the Hubble constant includes the contribution of the HS, GeV is the plank mass, g * and g h * are the energy d.o.f of the VS and HS respectively.

Benchmark model
From here, for the specific calculation of some observables, we implement GHPDM as the benchmark model. When the HS contains only a Majorana fermion DM χ and a complex scalar mediator s, the relevant Lagrangian is written as follows where h u and h d stand for the up and down-type Higgs doublets respectively. The Higgs portal of the Standard Model provides the opportunity for coupling to a very light scalar field s via the renormalizable operator |s| 2 (h 2 u + h 2 d ) and super-renormalizable operator sh u h d . It allows for the existence of large direct detection cross section and Higgs invisible decay width. To avoid these two constraints, we assume the portal coupling λ and A λ to be ultra-tiny. The κ quantifies the interaction between DM and mediator, while A κ describes the self-coupling between singlets. The GHPDM can be viewed as the low energy effective theory of general NMSSM. The UV complete interaction is determined by the superpotential and soft susy breaking term: (2.11) Since we are only interested in the DM phenomenology, there is no need to consider the UV completion-supersymmetry in detail. The complex singlet mediator can be further divided into CP-even and odd part, Therefore our HS is described in the following lagrangian.
The input parameters of this model are thus {m χ , m φ , α = κ 2 /4π, , A κ }. We should mention that this HS is not limited to supersymmetry but UV completion of the Higgsportal model. The CP-odd singlet a could provide a long-range interaction. One main reason for giving up this option is that the long-range potential mediated by the CP-odd singlet is the function of the spin structure which becomes negligible after averaging spin configuration [3]. As a result, we only focus on Yukawa potential mediated by CP-even singlet φ, (2.14)

Constraints from BBN and ∆N eff
To avoid overclosing the universe, the mediator φ and a ultimately decay into Standard Model particles. The most important constraint for its lifetime comes from Big Bang Nucleosynthesis (BBN) Γ φ > 10 −25 GeV. The decay width of the mediator comes from mass-insertion calculation [28]. The allowed decay channels are almost the same as Higgs decay if the kinematics channel is open such as φ → ff , gg, γγ. Here we only list φ → gg for brevity.
Constraints from BBN are shown in Fig. 2. The grayed shaded region corresponds to the decays occurring after BBN. An additional requirement of HS places upper bound on the coupling . The orange shaded regions are excluded assuming that the HS remains thermal equilibrium with SM thermal bath via the condition Γ ∼ 0.1H eff (T h = m φ ).
The energy density of the mediator affects the measurement of the effective number of relativistic neutrino species N eff . Current CMB experiment thus can probe our HS via 10 -4 10 -3 10 -2 10 -1 10 0 10 -11 It is easy to find that the maximal ∆N eff for our model is 0.1 at the limit of a massless mediator. So we can safely ignore ∆N eff constraint. Note that if mediator dominantly decays into a neutrino, the constraint re-appears. However, the branching ratio into neutrino is tiny in the Higgs portal.

Self-Interacting Dark Matter Realization
There are three main approaches to realize velocity-dependent self-interacting cross-sections. One is the existence of a light mediator to produce velocity-dependent cross-section by nonperturbative resummation. The second approach borrows the similarity of neutron-proton scattering where confinement in the HS adapts finite-size effect [29] to decrease the crosssection at low velocity. Finally, resonant SIDM [30] can naturally host a velocity-dependent behavior. In our paper, we focus on the first approach.
When the momentum of DM k is much less than the mass of mediator m φ , the crosssection is mainly from the S-wave. However, for k m φ , higher angular momentum > 0 states become more important when it goes to classical limit k/m φ → ∞. Commonly used momentum transfer cross section [31][32][33] is more suitable than total cross-section σ = dΩ(dσ/dΩ) due to the scattering velocity and angle dependence of dσ/dΩ. However, problems will emerge when the DM is the identical particle. Therefore, viscosity cross section [14] σ is adopted from heat conductivity [34]. Viscosity cross-section may better characterize DM halo dynamics since DM scattering is described by heat conductivity in fluid formulations of self-interacting DM [35,36].
In the partial wave analysis, the differential cross section is achieved by determining the phase shift δ via solving the Schrodinger equation Thus, the approximated viscosity cross section is obtained For identical particles, spatial part of the wave equation requires to be even (odd) for symmetric (anti-symmetric) spin parts. Then the differential cross section is given by where the positive (negative) sign correspond to an even (odd) spatial wave function, respectively. From Eq. (3.6) and P (−x) = (−1) P (x), we know that, for even (odd) case contribution from odd (even) phase shifts vanishes. The viscosity cross sections for even and odd partial wave cases are For realistic unpolarized scattering particles, contribution to the total cross section comes both from even and odd viscosity cross section. For example, cross section for selfinteracting Majorana DM is The last step to get useful cross section is to find analytical solution of the phase shift δ represented in terms of which correspond to the dimensionless momentum and strength of the potential relative to the kinetic energy, respectively. Viscosity cross section for attractive potential is [37] σ Semi−Classical where K 0 , K 1 and K 2 are the Bessel functions of the second kind. These formula is valid in semi-classical limit i.e. η > 1. When η is smaller than 0.4, it belongs to quantum regime where Hulthen potential can be used to solve Schrodinger equation where the phase δ is determined as follows δ = arg i Γ(l + + l − − 2) Γ(l + )Γ(l − ) , l + = 1 + η 1.6 (i + i 3.2β − 1), (3.14) Within the regime 0.4 < η < 1, we adapt the interpolation method σ Interpolation in the overlap region. In terms of these three cross-sections, we can cover most of interesting parameter space of SIDM.
Calculation of velocity averaged viscosity cross-section do not simply follow the procedure of multiplying the σ V by velocity and taking average, instead energy transfer rate related to the viscosity cross-section is defined by [37] where v rel and v 0 are DM relative velocity and velocity dispersion of v which obeys Maxwell-Boltzmann distribution. Averaged cross section per unit massσ v rel /m χ of fermionic DM as a function of v rel is presented in Fig. 3 for the benchmark point m χ = 200 GeV, m φ = 5 MeV, α = 0.3. As is shown in the Fig. 3, our model is suitable as a solution to the small scale problems, since the resulting curve is compatible to the data points from five clusters, seven low surface brightness spiral galaxies and six dwarf galaxies. There are two relevant parameters in fitting particle input into data (η, β) →σ v /m χ . However η and β contain the velocity which is not directly relevant to projection on HS. Thus we choose the mass ratio m φ /m χ and α as input parameters and perform best fit on this two-parameter model with m χ = 1 TeV. Figure 4 displays the chi-squared distribution. The value of α has tiny effect on chi-squared which means correct self-interaction cross section is almost independent of α once we fix the value of m χ . The results of the correlation analysis shows that the ratio m φ /m χ should be smaller than 10 −4 . It is compatible with the semi-classical regime. Throughout the paper we require m φ /m χ < 10 −4 .

Cross-Section and Implication on Relic Density
In this section we consider χχ → φφ, aa, φa processes to compute the relic density. It encodes the underlying story on the SE, radiative BSF, and HS freeze-out. Typically, HS freeze-out process includes only born level calculation σ eff v rel = σ tree v rel . However, non-perturbative corrections from SE and BSF are not negligible when requiring SIDM in Eq. (2.14). The effective cross-section in our model is thus written as follows Evaluating the tree-level Feynman diagrams yield the born level cross sections [39] σ ann (χχ → φφ) v rel 17 256π (4. 2) The first two processes are p-wave contribution σ p v rel and the third one is s-wave crosssection σ s v rel . The expression is valid only in large m χ /m φ,a limit. The long-range effect is re-summed to solve the Schrodinger equation with Yukawa potential V (r). The short distance annihilation cross-section is encoded in the absorptive part of forward scattering amplitude called Wilson coefficient f [ 2s+1 L j ]. The schematic form of s-wave Sommerfeld enhanced cross-section is written as σ eff v rel = σ tree v rel S 0 . In the limit of Hulthen potential, we can present SE in s-wave and p-wave by and respectively. Here ε v ≡ v rel / (2α) and ε φ ≡ m φ / (αm χ ). As a result, the thermal averaged cross-section with the SE is obtained by In most parameter spaces, the contributions of the s-wave predominate over the p-wave, although it suffers a dangerous constraint from CMB ionization. However, we propose a blind spot scenario A κ = −3κm χ , where the s-wave cross-section vanishes automatically. The CMB constraint is thus alleviated. However, the BSF cross-section remains s-wave to re-introduce CMB constraint, even though it is small.
The cross section of radiative BSF in ground state, i.e. n = 1, l = m = 0 is [40] exp(−4ζ arccot(ζ)), (4.6) where the characteristic parameter ζ = α/v rel represents the ratio between Bohr momentum and momentum exchange of unbound particles, while ξ = µα/0.84m φ is the ratio of interaction range and Bohr radius. Both are larger than 1 to allow for the existence of the bound state, existing as phase space suppression in the following formula. The effectively BSF cross-section is also the function of the ionization and decay process where ionization rate of bound states is The binding energy E B determines the relative strength of the ionization rate. The annihilation decay is responsible for the decay of the bound state where we replace the Sommerfeld factor by the squared of bound state wave-function The Eq. (4.9) is only valid at the s-wave dominated case, since the conservation of angular momentum forbids the p-wave decay process. Instead, the only allowed decay channel is 4φ process B → φφφφ, (4.10) The Boltzmann equation in Eq. (2.4) can be numerically solved for the cross section in Eq. (4.1) to obtain DM relic density. In Fig. 5, we give the solid lines which satisfy correct relic density, constraints from various experiments are also presented. We can obtain the correct relic density in blue lines where the initial condition on reheating is r inf = 1. To show the deviation of parameter space, we also consider different values of r inf = 10 −1 , 10 1 with purple and orange lines. In the bottom panel, we find larger r inf requires larger interaction strength α and is ruled out by CMB constraint. Large r inf induces a large effective Hubble parameter to obtain correct relic density, a larger annihilation cross-section is required that favors larger α. The panel on the top left shows that positive A κ does not affect relic density to a large extent, since cross-section is almost the function of α and m χ . For negative A κ on the top right panel, the cancellation between A κ and m χ will make the cross-section smaller than the required one. As a result, light DM compensates for the cancellation effect. When negative A κ becomes smaller than m χ , it will be the same as positive A κ . In addition, the CMB constraint rules out the r inf larger than 10 for the positive A κ . While negative A κ scenario almost does not have a viable parameter space after implementing the CMB and the BBN constraints. The observation data from the  • Green-shaded region for dwarf. To solve small-scale issues, a large self-interaction cross-section in dwarf galaxies is required. However, it can not be arbitrarily large. Instead astrophysical data puts an upper bound on self-interaction cross-section on dwarf scale i.e. σ/m χ < 50 cm 2 g −1 .
The bottom panel displays the weak dependence of self-interaction cross-section on α and m χ . While A κ plays no role in SIDM, the vertical line in the top panel indicates the sensitivity of SIDM on m χ .
• Magenta-shaded region for cluster. On the cluster scale, according to the astrophysical observations, the self-interacting cross-section must reduce to be smaller than 1 cm 2 g −1 .
• Gray shaded region for BBN. The SIDM works well in the HS itself, which is not sensitive to different temperature evolution. While BBN is the function of decay width of mediator, it is sensitive to the r inf . The decoupling from kinetic equilibrium can be realized by selecting = 10 −8 for r inf = 10 −1 , 10 0 and = 6 × 10 −9 for r inf = 10. Since we choose the ratio m φ /m χ = 10 −5 in both figures, the constraint on a mediator from BBN can be interpreted as the constraint on DM mass. Therefore we have two vertical lines for the BBN boundary. It is easy to find that the BBN constraint is the most strict one for light DM.
• Blue-shaded region for CMB The DM annihilation products inject energy into the CMB and lead to anisotropies. It provides a sensitive probe of DM annihilation during the dark ages. The s-wave cross-section is bounded from above [41], in order not to distort the CMB spectrum, which is roughly [42] lim v→0 (σv) < 3 × 10 −24 cm 3 sec −1 × m χ TeV . (4.11) Even though the negative A κ leads to null parameter space for our model, there is still a blind spot, in which the s-wave annihilation automatically vanishes and the p-wave annihilation process dominates. Therefore, the CMB constraint becomes relevant only if the BSF cross-section becomes large. The task, for now, is to compute the limit of v → 0 for BSF cross-section. The SE in Eq. (4.6) is where a ± = 1 + iw(1 ± 1 − x/ω), ω = m χ v rel /2m φ and x = 2α/v rel . When the velocity is smaller than m φ /m χ , it has a non-vanishing limit (4.13) By substituting into the s-wave BSF cross-section, we obtain the analytical cross-section (σ BSF ) v rel →0 = 4608απ 3 1 − (16m 2 φ /α 4 m 2 χ ) e 4 m φ m χ sin 2 (π αm χ /m φ ) . (4.14) Generally, it is larger than the CMB bound in Eq. (4.11). The loophole is that the BSF cross-section vanishes when the phase space is not open, i.e. α < 6.32 × 10 −3 , which we use it to constrain the coupling in Fig. 6. Finally, we show the benchmark points that satisfy both relic density and constraints in Tab. 1. Even though, large cancellation between negative A κ and m χ can increase the fraction of BSF cross-section, BBN prevents it to be larger than 5%.
lived mediator. Unfortunately, these observational consequences would not, as many other DM models, differentiate our model from others which predict same signal detection.

Conclusion
The purpose of the paper is to explore the long-range potential effect in the hidden sector. It provides the correct self-interaction cross-section and a new s-wave annihilation contribution from the bound state formation. Our study quantifies the influence of different temperature ratios r inf on relic density in General Higgs portal dark matter model. Even though the bound state formation cross-section is not sufficient to affect relic density, the existence of the new s-wave contribution plays a crucial role in CMB constraint. For the positive A κ , the combination of CMB and BBN constraint favors the heavy dark matter and not-too-hot hidden sector. For the negative A κ , no parameter space for the model except for the blind spot scenario that we propose.