Thermodynamic Evolution of Secluded Vector Dark Matter: Conventional WIMPs and Nonconventional WIMPs

TL;DR

This paper investigates the thermodynamic evolution of secluded vector dark matter, demonstrating how a long-lived mediator can boost annihilation cross sections and fit galactic gamma-ray excess data, with implications for future detection efforts.

Contribution

It provides a first-principles thermodynamic analysis of secluded vector dark matter, including coupled Boltzmann equations, and explores parameter space consistent with gamma-ray excess and experimental bounds.

Findings

Long-lived mediators can enhance dark matter annihilation cross sections.

Parameter space fits Galactic center gamma-ray excess data.

Future dwarf spheroidal galaxy observations can probe relevant dark matter models.

Abstract

The secluded dark matter resides within a hidden sector and self-annihilates into lighter mediators which subsequently decay to the Standard Model (SM) particles. Depending on the coupling strength of the mediator to the SM, the hidden sector can be kinetically decoupled from the SM bath when the temperature drops below the mediator's mass, and the dark matter annihilation cross section at freeze-out is thus possible to be boosted above the conventional value of weak interacting massive particles. We present a comprehensive study on thermodynamic evolution of the hidden sector from the first principle, using the simplest secluded vector dark matter model. Motivated by the observation of Galactic center gamma-ray excess, we take two mass sets $\sim{\cal O}(80\, \text{GeV})$ for the dark matter and mediator as examples to illustrate the thermodynamics. The coupled Boltzmann moment equations for number densities and temperature evolutions of the hidden sector are numerically solved. The formalism can be easily extended to a general secluded dark matter model. We show that a long-lived mediator can result in a boosted dark matter annihilation cross section to account for the relic abundance. We further show the parameter space which provides a good fit to the Galactic center excess data and is compatible with the current bounds and LUX-ZEPLIN projected sensitivity. We find that the future observations of dwarf spheroidal galaxies offer promising reach to probe the most relic allowed parameter space relevant to the boosted dark matter annihilation cross section.

Equations

$${\mathcal{L}}_{\text{kinetic}}\supset -\frac{1}{4}{X}_{\mu \nu }{X}^{\mu \nu }+(D_{\mu }{\Phi }_S{)}^{†}(D^{\mu }{\Phi }_S),$$

$${\mathcal{L}}_{\text{scalar}}=-{\mu }_H^2|{\Phi }_H{|}^2-{\mu }_S^2|{\Phi }_S{|}^2-\frac{{\lambda }_H}{2}({\Phi }_H^{†}{\Phi }_H{)}^2-\frac{{\lambda }_S}{2}({\Phi }_S^{†}{\Phi }_S{)}^2-{\lambda }_{HS}({\Phi }_H^{†}{\Phi }_H)({\Phi }_S^{†}{\Phi }_S),$$

$${\Phi }_H=\frac{1}{\sqrt{2}}(v_H+{\varphi }_h+i{\sigma }_h),{\Phi }_S=\frac{1}{\sqrt{2}}(v_S+{\varphi }_s+i{\sigma }_s),$$

$$\displaystyle\left(\begin{array}[]{c}\phi_{h}\\ \phi_{s}\end{array}\right)$$

$$\left(\begin{array}[]{cc}\lambda_{H}v_{H}^{2}&\lambda_{HS}v_{S}v_{H}\\ \lambda_{HS}v_{S}v_{H}&\lambda_{S}v_{S}^{2}\end{array}\right)=\left(\begin{array}[]{cc}m_{h}^{2}c_{\alpha}^{2}+m_{S}^{2}s_{\alpha}^{2}&(m_{h}^{2}-m_{S}^{2})s_{\alpha}c_{\alpha}\\ (m_{h}^{2}-m_{S}^{2})s_{\alpha}c_{\alpha}&m_{S}^{2}c_{\alpha}^{2}+m_{h}^{2}s_{\alpha}^{2}\end{array}\right)\,.$$

$$H={\left(\frac{8\pi G}{3}{\rho }_t\right)}^{1/2},$$

$${\rho }_t=\frac{{\pi }^2}{30}\left(g_{\mathrm{eff},\mathrm{SM}}(T)T^4+g_{\mathrm{eff},\mathrm{h}}(T_h)T_h^4\right).$$

$${\sigma }_N=\frac{{\mu }_{XN}^2{m}_N^2{f}_N^2{g}_{\mathrm{dm}}^2}{4\pi }\frac{{\sin }^{2}2\alpha }{v_H^2}{\left(\frac{1}{m_S^2}-\frac{1}{m_h^2}\right)}^2,$$

$$\frac{\partial f_h}{\partial t}-Hp\frac{\partial f_h}{\partial p}=C[f_h],$$

$$f_h(E_h,T_h)=\frac{1}{e^{(E_h-{\mu }_h)/T_h}-1},$$

$$C[f_{h_1}]=$$

$$\times (2\pi {)}^4{\delta }^{(4)}(p_{h_1}+p_{h_2}+\cdots +p_{b_1}+p_{b_2}+\cdots -p_{h_1^{\prime }}-p_{h_2^{\prime }}-\cdots -p_{b_1^{\prime }}-p_{b_2^{\prime }}-\cdots )|M{|}^2$$

$$\displaystyle\times\frac{\Delta-\Delta^{\prime}}{S_{y}\,S_{y}^{\prime}}\Big{[}f_{h_{1}^{\prime}}f_{h_{2}^{\prime}}\cdots f_{b_{1}^{\prime}}f_{b_{2}^{\prime}}\cdots(1+f_{h_{1}})(1+f_{h_{2}})\cdots(1\pm f_{b_{1}})(1\pm f_{b_{2}})$$

$$\displaystyle\ -f_{h_{1}}f_{h_{2}}\cdots f_{b_{1}}f_{b_{2}}\cdots(1+f_{h_{1}^{\prime}})(1+f_{h_{2}^{\prime}})\cdots(1\pm f_{b_{1}^{\prime}})(1\pm f_{b_{2}^{\prime}})\Big{]}\,,$$

$$d{\Pi }_i\equiv \frac{d^3{p}_i}{(2\pi {)}^{3}2E_i},$$

$$f_i=e^{-(E_i-{\mu }_i)/T_i}(1\pm f_i)\simeq e^{-(E_i-{\mu }_i)/T_i}.$$

$$C[f_{h_1}]$$

$$\times$$

$$n_X(T_X)=g_X\int \frac{d^3{p}_X}{(2\pi {)}^3}{f}_X(T_X),n_S(T_S)=g_S\int \frac{d^3{p}_S}{(2\pi {)}^3}{f}_S(T_S).$$

$$\frac{d{n}_X}{dt}+3H{n}_X=$$

$$\displaystyle-\langle\sigma v^{2}\rangle_{XXS\to SS}\bigg{(}n_{X}^{2}n_{S}-(n_{X}^{\text{eq}})^{2}\frac{n_{S}^{2}}{(n_{S}^{\text{eq}})}\bigg{)}+\frac{1}{3}\langle\sigma v^{2}\rangle_{SSS\to XX}\bigg{(}n_{S}^{3}-(n_{S}^{\text{eq}})^{3}\frac{n_{X}^{2}}{(n_{X}^{\text{eq}})^{2}}\bigg{)}\,,$$

$$\frac{d{n}_S}{dt}+3H{n}_S=$$

$$\displaystyle-\bigg{(}\langle\sigma v\rangle_{SS\to\sum_{ij}{\rm SM}_{i}\,{\rm SM}_{j}}n_{S}^{2}-\langle\sigma v\rangle_{SS\to\sum_{ij}{\rm SM}_{i}\,{\rm SM}_{j}}(T)(n_{S}^{\text{eq}}(T))^{2}\bigg{)}$$

$$\displaystyle+\langle\sigma v\rangle_{XX\to SS}\bigg{(}n_{X}^{2}-(n_{X}^{\text{eq}})^{2}\frac{n_{S}^{2}}{(n_{S}^{\text{eq}})^{2}}\bigg{)}+\frac{1}{6}\langle\sigma v^{2}\rangle_{XXX\to XS}\bigg{(}n_{X}^{3}-n_{X}n_{S}\frac{(n_{X}^{\text{eq}})^{2}}{n_{S}^{\text{eq}}}\bigg{)}$$

$$\displaystyle+\frac{1}{2}\langle\sigma v^{2}\rangle_{XXS\to SS}\bigg{(}n_{X}^{2}n_{S}-(n_{X}^{\text{eq}})^{2}\frac{n_{S}^{2}}{(n_{S}^{\text{eq}})}\bigg{)}-\frac{1}{2}\langle\sigma v^{2}\rangle_{XSS\to XS}\bigg{(}n_{X}n_{S}^{2}-n_{X}n_{S}n_{S}^{\text{eq}}\bigg{)}$$

$$\displaystyle-\frac{1}{2}\langle\sigma v^{2}\rangle_{SSS\to XX}\bigg{(}n_{S}^{3}-(n_{S}^{\text{eq}})^{3}\frac{n_{X}^{2}}{(n_{X}^{\text{eq}})^{2}}\bigg{)}-\frac{1}{6}\langle\sigma v^{2}\rangle_{SSS\to SS}\bigg{(}n_{S}^{3}-n_{S}^{2}n_{S}^{\text{eq}}\bigg{)}\,,$$

$$y_X(x_X;x)$$

$$y_S(x_S;x)$$

$$\displaystyle\frac{d}{dt}=\frac{m_{X}^{2}}{x}\bigg{(}\frac{8\pi^{3}G}{90}\bigg{)}^{1/2}\frac{h_{\rm eff}}{g_{*}^{1/2}}\frac{d}{dx}\,,$$

$${\tilde{h}}_{\mathrm{eff}}\equiv h_{\mathrm{eff}}[1+(1/3)(d\ln h_{\mathrm{eff}}/d\ln T)],$$

$$s(T)=\frac{2{\pi }^2}{45}{h}_{\mathrm{eff}}(T)T^3,\rho (T)=\frac{{\pi }^2}{30}{g}_{\mathrm{eff}}(T)T^4.$$

$$\frac{d{y}_X}{dx}=$$

$$\displaystyle+\frac{\delta_{\text{dof}}}{x^{5}}\frac{\pi}{\sqrt{90}}\frac{h_{\text{eff}}(T)}{g_{*}^{1/2}(m_{X})}\frac{m_{X}^{2}}{M_{\rm pl}}\Bigg{[}-\frac{1}{3}\frac{\langle\sigma v^{2}\rangle_{XXX\to XS}}{(\langle\sigma v\rangle^{(0)}_{XX\to SS})^{2}}\bigg{(}y_{X}^{3}-y_{X}(y_{X}^{\text{eq}})^{2}\frac{y_{S}}{y_{S}^{\text{eq}}}\bigg{)}$$

$$\displaystyle+\frac{1}{3}\frac{\langle\sigma v^{2}\rangle_{SSS\to XX}}{(\langle\sigma v\rangle^{(0)}_{XX\to SS})^{2}}\bigg{(}y_{S}^{3}-\frac{y_{X}^{2}}{(y_{X}^{\text{eq}})^{2}}(y_{S}^{\text{eq}})^{3}\bigg{)}-\frac{\langle\sigma v^{2}\rangle_{XXS\to SS}}{(\langle\sigma v\rangle^{(0)}_{XX\to SS})^{2}}\bigg{(}y_{X}^{2}y_{S}-\frac{y_{S}^{2}}{y_{S}^{\text{eq}}}(y_{X}^{\text{eq}})^{2}\bigg{)}\Bigg{]}\,,$$

$$\frac{d{y}_S}{dx}=$$

$$\displaystyle-\frac{\delta_{\text{dof}}}{x^{2}\langle\sigma v\rangle^{(0)}_{XX\to SS}}\Big{[}\langle\sigma v\rangle_{SS\to\sum_{i}{\rm SM}_{i}{\rm SM}_{i}}y_{S}^{2}-\langle\sigma v\rangle_{SS\to\sum_{i}{\rm SM}_{i}{\rm SM}_{i}}(x)\,(y_{S}^{\text{eq}}(x))^{2}\Big{]}$$

$$\displaystyle-\frac{\delta_{\text{XS}}\delta_{\text{dof}}}{x^{2}}\bigg{[}\frac{(y_{X}^{\text{eq}})^{2}}{(y_{S}^{\text{eq}})^{2}}y_{S}^{2}-y_{X}^{2}\bigg{]}+\frac{\delta_{\text{dof}}}{x^{5}}\frac{\pi}{\sqrt{90}}\frac{h_{\text{eff}}(T)}{g_{*}^{1/2}(m_{X})}\frac{m_{X}^{2}}{M_{\rm pl}}\Bigg{[}\frac{1}{6}\frac{\langle\sigma v^{2}\rangle_{XXX\to XS}}{(\langle\sigma v\rangle^{(0)}_{XX\to SS})^{2}}\bigg{(}y_{X}^{3}-y_{X}(y_{X}^{\text{eq}})^{2}\frac{y_{S}}{y_{S}^{\text{eq}}}\bigg{)}$$

Full text

I Introduction

Motivated by particle physics, the theoretical studies and experimental searches have for many decades focused on the popular class of the dark matter (DM) candidates, called the weakly interacting massive particles (WIMPs). In the WIMP scenario, when the dark matter becomes nonrelativistic, its comoving number density is exponentially depleted through Boltzmann suppression and keeps the thermal equilibrium with the bath until freeze-out. The resulting DM with the weak scale interaction and mass can provide the correct relic abundance today.

Many DM experiments are thus motivated by the WIMP scenario. Nevertheless, no conclusive observations have been made by the direct detection searches, Large Hadron Collider (LHC), and other collider experiments. Several groups have reported the GeV gamma-ray excess around the Galactic center (GC) Goodenough:2009gk ; Hooper:2010mq ; Hooper:2011ti ; Abazajian:2012pn ; Gordon:2013vta ; Huang:2013pda ; Daylan:2014rsa ; Calore:2014xka ; Calore:2014nla ; Karwin:2016tsw ; TheFermi-LAT:2017vmf , for which, however, the allowed WIMP dark matter models have been also severely constrained by the current null results of the direct detection Aprile:2017iyp ; Cui:2017nnn ; Akerib:2016vxi ; Akerib:2018lyp and collider experiments. In light of these measurements, an interesting paradigm that goes beyond the “conventional” WIMP scenario and becomes more and more popular is known as “secluded (WIMP) dark matter”. In this paradigm, the dark matter candidate may reside within one of the hidden sectors and communicates with the visible sector through a lighter metastable mediator, which weakly couples the standard model (SM) to the WIMP. As such, the DM signals, suppressed at the direct detection and colliders, could be observable in indirect measurements Pospelov:2007mp ; Ko:2014gha ; Berlin:2014pya ; Escudero:2017yia ; Ko:2014loa ; Abdullah:2014lla ; Martin:2014sxa ; Kim:2016csm ; Yang:2017zor ; Profumo:2017obk .

The mechanism for the secluded WIMP dark matter was discussed by Pospelov, Ritz, and Voloshin Pospelov:2007mp . In this mechanism, the WIMP can still be a thermal relic, and the dominant DM annihilation channel is into a pair of unstable mediators which ultimately decay into SM particles. Basically, for this model, as long as the mediator decays before the beginning of the big bang nucleosynthesis (BBN), the effective number of neutrino species and abundance of helium and deuterium will not be modified, as compared with the standard BBN, so that the result can be easily compatible with the current Planck measurement Ade:2015xua .

As for building secluded DM models, many people restricted their works to the parameter space relevant to the WIMP scenario where the hidden sector is in chemical and thermal equilibrium with the bath prior to freeze-out  Ko:2014gha ; Berlin:2014pya ; Escudero:2017yia ; Ko:2014loa ; Abdullah:2014lla ; Martin:2014sxa ; Kim:2016csm ; Yang:2017zor . However, for the case that the dark sector has kinetically decoupled from the bath, due to its weak couplings to the SM particles, before it becomes nonrelativistic, if the secluded DM annihilates into nearly degenerate mediators which later decay out-of-equilibrium with the bath, the DM density will be exponentially depleted through the decay process of the mediator, instead of following Boltzmann suppression Dror:2016rxc . Moreover, during the period of time in which the dark sector is out of thermal equilibrium with the bath, if the $3\to 2$ number changing interactions are allowed and efficiently active, the hidden sector can first undergo an epoch called “cannibalism”. See the related discussions in Refs. Farina:2016llk ; Pappadopulo:2016pkp ; Yang:2018fje ; Berlin:2016gtr . Alternatively, $3\to 2$ DM annihilation mechanism is also relevant to the strongly-interacting massive particles (or called SIMP) Hochberg:2014dra and elastically decoupling relic (or called ELDER) Kuflik:2015isi ; Kuflik:2017iqs scenarios.

In this paper, to have a thorough understanding about the thermodynamics of the secluded dark matter from the first principle, we will study the simplest secluded vector dark matter model, taken as an example in which the vector dark matter and the mediator within the hidden sector are in thermal equilibrium with each other before freeze-out, but may be kinetically decoupled from the SM bath at temperature $T\sim m_{X,S}$, depending on the couplings to the SM, where $m_{X}$ and $m_{S}$ are the masses of the DM and hidden scalar, respectively.

We separately obtain the evolution equations of number densities and temperatures for the hidden species, by taking suitable moments of the Boltzmann equation. We will give a detailed result of describing chemical and kinetic decouplings of the hidden sector from the thermal bath. We will show that, depending on the coupling strength of the mediator to the SM, the relic annihilation cross section is likely to be boosted above the conventional WIMP value. The present study can be easily generalized to a generic case.

Using two mass sets: (i) $m_{X}=80$ GeV, $m_{S}=0.8m_{X}$, and (ii) $m_{X}=80$ GeV, $m_{S}=0.99m_{X}$, we numerically solve the thermodynamic evolution of the hidden sector, which can be either in thermal equilibrium or out of equilibrium with the bath before the DM freezes out, and moreover, is secluded from the visible sector with small interaction rates compatible with colliders and direct detection bounds. Use of the present mass sets of the hidden sector is motivated by the observed GC gamma-ray excess which can be accounted for by this model via one-step cascade annihilation Ko:2014gha ; Escudero:2017yia . More detailed discussions about the GC allowed region, which are constrained by the astrophysical and cosmological measurements as well as the LUX-ZEPLIN projected sensitivity Akerib:2018lyp , will be presented in Sec. VI.

The rest of this paper is organized as follows. In Sec. II, we start with an introduction of the vector DM model which is UV-complete. In this model, the hidden sector contains an abelian vector dark matter and a complex scalar. The former is a gauge boson associated with a dark (hidden) gauge symmetry $U_{X}(1)$, while the latter is charged under $U_{X}(1)$. In Sec. III, the model parameters constrained by direct detection and collider experiments will be described first. In Sec. IV, we present a general description of Boltzmann equation in the framework of an expanding Universe which is homogeneous and isotropic. We further consider the moments of Boltzmann equations that are relevant to the evolutions of the number densities and temperatures for the hidden species. In Sec. V, two sets of mass parameters which can account for the GC gamma-ray excess are used in the numerical analyses. The results are given and discussed. The parameter space relevant to the GC gamma-ray excess and concerning the current limits and prospects are further discussed in Sec. VI. In Sec. VII, we draw the conclusions. All technical derivations are collected in Appendices.

II The Model

The simplest secluded vector dark matter can be made of the abelian gauge bosons, $X_{\mu}$’s, which get the mass from the vacuum expectation value (VEV) of the hidden complex scalar field $\Phi_{S}$ due to the spontaneously dark gauge symmetry $U_{X}(1)$ breaking, where the $Z_{2}$ symmetry, $X_{\mu}\to-X_{\mu}$ and $\Phi_{S}\to\Phi_{S}^{*}$, is imposed to stabilize the dark matter Ko:2014gha . The relevant kinetic Lagrangian (${\cal L}_{\text{kinetic}}$) and the scalar potential (${\cal L}_{\text{scalar}})$ are given by

[TABLE]

where $\Phi_{H}=(H^{+},H^{0})^{\rm T}$ is the SM Higgs doublet, $X_{\mu\nu}=\partial_{\mu}X_{\nu}-\partial_{\nu}X_{\mu}$, and the covariant derivative is defined as $D_{\mu}\Phi_{S}=(\partial_{\mu}+ig_{\rm dm}Q_{\Phi_{S}}X_{\mu})\Phi_{S}$. Here $g_{\rm dm}$ is the gauge coupling and $Q_{\Phi_{S}}$ is the $U_{\rm dm}(1)$ charge of $\Phi_{S}$. After spontaneous symmetry breaking, the Higgs fields develop non-zero VEV’s,

[TABLE]

where the CP-odd states, $\sigma_{h}$ and $\sigma_{s}$, respectively becomes the longitudinal components of the $Z$ boson and $X_{\mu}$; the dark matter thus obtain a mass, $m_{X}=g_{\rm dm}Q_{\Phi_{S}}v_{S}$. In the present paper, we will simply use $Q_{\Phi_{S}}=1$.

The scalar fields $(\phi_{h},\phi_{s})$ can be expressed in terms of mass eigenstates of physical Higgses $(h,S)$ as

[TABLE]

and the mass squared matrix in the former basis can be parametrized in terms of masses of the latter and the mixing angle $\alpha$,

[TABLE]

Here and throughout the paper, we adopt the abbreviations: $s_{\alpha}\equiv\sin\alpha$ and $c_{\alpha}\equiv\cos\alpha$. Using $v_{H}\simeq 246$ GeV and $m_{h}=125.18$ GeV pdg2018 , we will take $m_{X},m_{S}$, $g_{\rm dm}$ and $\alpha$ as the independent parameters in the following analysis.

The branching ratios of the hidden scalar, $S$, with a mass of $m_{S}\lesssim 2m_{h}$, are depicted in the left panel of Fig. 1, where, in the range giving a good fit to the GC gamma-excess data, the scalar mass satisfies $m_{S}\lesssim m_{X}\lesssim 130$ GeV. The related partial widths of the hidden scalar $S$ are summarized and discussed in Appendix A, where the results are relevant to the studies of the relic abundance and indirect detection searches. In the right panel of Fig. 1, we show several values of the $S$ decay width, $\Gamma_{S}$, compared with the evolution of $2H$ (the inverse time interval of the radiation dominated epoch), where $H$ is the Hubble rate which is given by the Friedmann equation,

[TABLE]

with the total energy density being

[TABLE]

Here, $g_{\rm eff,SM}$ and $g_{\rm eff,h}$ are the effective relativistic degrees of freedom of the SM and hidden sector at the temperatures $T$ and $T_{h}$, respectively. In Fig. 1, we have simply adopted $T_{h}=T$. As will be shown in Eqs. (37) and (64), and discussed in Sec. V(iv), if the nonrelativistic scalar $S$ is kept in kinetic equilibrium with the thermal bath via its inverse decay $\text{SM}~{}\text{SM}\to S$, then this kinetic energy injection rate to $S$ will be larger than the Hubble cooling rate, i.e., roughly $\Gamma_{S}\gtrsim 2H$. In the right panel of Fig. 1, the $S$ particles with the width corresponding to $\alpha=1\times 10^{-5}$ or $1\times 10^{-6}$ can be in thermal equilibrium with the bath when $m_{X}/T\gtrsim 0.4$ or 3 (with $m_{X}$=80 GeV) (see Sec. V(iv), where a more precise estimation is given). However, for the hidden scalar with a much smaller mixing angle $\alpha=5\times 10^{-7}$ or $1\times 10^{-7}$, because the ratio of the Hubble cooling rate to heating rate, \sim 2Hn_{S}(T_{S})T_{S}/\big{(}\Gamma_{S}n_{S}^{\rm eq}(T)\,T\big{)}, is much larger than 1 due to the fact that $n_{S}(T_{S})T_{S}\gg n_{S}^{\rm eq}(T)T$ for the nonrelativistic $S$ (see Figs. 4 and 5), where $n_{S}(T_{S})$ is the number density of $S$ at its temperature $T_{S}$, and $n_{S}^{\rm eq}(T)$ is the equilibrium number density of $S$ at the corresponding bath temperature $T$, the hidden scalar thus starts to undergo out-of-equilibrium decay at the cosmological time $\simeq$ the $S$ lifetime, $(2H)^{-1}\simeq\Gamma_{S}^{-1}$. As shown in the right panel of Fig. 1, the corresponding out-of-equilibrium temperature is about $m_{X}/T\sim$ 6 or 30 for $\alpha=5\times 10^{-7}$ or $1\times 10^{-7}$. The underlying physics and a more precise estimation will be given in Sec. V(vi).

III Direct detection and LHC constraints

In this paper, we will use two sets of the masses for the dark matter and mediator: (i) $m_{X}=80$ GeV, $m_{S}=0.8m_{X}=64$ GeV, and (ii) $m_{X}=80$ GeV, $m_{S}=0.99m_{X}=79.2$ GeV, to study the thermal evolution of the hidden sector. These two sets can provide a good fit to the GC gamma-ray excess data. For the first set, when the hidden sector with a sizable mass gap undergoes the cannibal process, the down-scattering rate, $XX\to SS$, can be significantly larger than the up-scattering one, $SS\to XX$. For the second set, the hidden sector is nearly degenerate, and can be further constrained by the gamma-line searches at the indirect detection. Moreover, because the low-velocity DM annihilation cross section is zero in $m_{S}\to m_{X}$ limit, a larger $X$-$S$ coupling is needed to account for the GC data and the DM relic abundance.

In the secluded DM model, the direct detection measurements and colliders weakly constrain the parameter region allowed by the GC excess result. The spin-independent cross section for a vector dark matter particle scattering off a single nucleon via a scalar mediator $S$ exchange is given by

[TABLE]

where $\mu_{XN}=m_{X}m_{N}/(m_{X}+m_{N})$ is the reduced mass of the dark matter ($X$) and nucleon ($N$), and $f_{N}=\sum_{q}\langle N|\bar{q}q|N\rangle m_{q}/m_{N}\simeq 0.3$ Cline:2013gha . The parameter space constrained by XENON1T Aprile:2017iyp is shown in Fig. 2, where in the right panel the bound by the LUX-ZEPLIN (LZ) projected sensitivity Akerib:2018lyp is given. In the left panel of Fig. 2, the allowed parameter region on the ($m_{X},m_{S}$) plane for a given value of $\alpha$ is above the corresponding dashed line, where we have limited $\alpha\leq\pi/4$ which is suitable for the case with a small $\alpha$. As for $\pi/4\leq\alpha\leq\pi/2$, the bound is same as that with a mixing angle $=\pi/2-\alpha$, because $\sin 2\alpha=\sin 2(\pi/2-\alpha)$.

The paramter constraint from the invisible Higgs decay, which is less than 25% at the 95% CL pdg2018 , is much weaker than that from direct detection. Meanwhile, for the present case, $h\to S\,S$ is kinematically forbidden.

IV Thermal evolution of the nonrelativistic hidden particles

The evolution of the phase space distribution $f_{h}$ (with $h\equiv X\text{ or }S$) of the hidden sector particles in the homogeneous isotropic Friedmann-Robertson-Walker Universe is described by the Boltzmann equation,

[TABLE]

where $H$ is the Hubble expansion parameter, $p=\sqrt{E_{h}^{2}-m_{h}^{2}}$ is the momentum of the hidden particle, and $C[f_{h}]$ is the collision term. During the process of the thermal evolution, the distribution of the hidden sector particles follows Bose-Einstein statistics,

[TABLE]

with $\mu_{h}$ the chemical potential of the particle species $h$. In the present case, we consider $m_{X}\gtrsim m_{S}\sim{\cal O}(10-100~{}\text{GeV})$ and the thermal evolution that the elastic scattering $XS\leftrightarrow XS$ can keep the $X$ and $S$ particles in thermal equilibrium ($T_{X}=T_{S}$) until kinetic decoupling temperature $T_{X}^{\rm kd}\equiv T_{X}(a_{X}^{\rm kd})$, below that we have $T_{X}(a)\simeq T_{X}^{\rm kd}\cdot(a_{X}^{\rm kd}/a)^{2}$, where $a$ is the cosmic scale factor and $a_{X}^{\rm kd}$ is its corresponding value at $T_{X}^{\rm kd}$.

For a hidden sector particle, $h_{1}$, the generic form of one of the collision terms described by the interaction “$h_{1}h_{2}\cdots b_{1}b_{2}\cdots\leftrightarrow h_{1}^{\prime}h_{2}^{\prime}\cdots b_{1}^{\prime}b_{2}^{\prime}\cdots$” can be written as

[TABLE]

where

[TABLE]

$h_{i}^{(\prime)}$ is the particle of the hidden sector with temperature $T_{X}$ for $X$ or $T_{S}$ for $S$, $b_{i}^{(\prime)}$ is the relativistic SM particle with temperature $T$, $|M|^{2}$ invariant under times reversal and reflection is the square of the amplitude summed over the internal degrees of freedom (dof), $g_{i}$, of all the initial and final particles, $S_{y}$ and $S_{y}^{\prime}$ are the symmetric factors in the initial and final states, respectively, $\Delta$ and $\Delta^{\prime}$ are the numbers of the species which are the same as $h_{1}$ and participate in the interaction in the initial and final states, respectively; note that if the particle composition in the initial state is exactly the same as that in the final state (e.g. elastic scattering $S+\text{SM}\leftrightarrow S+\text{SM}$, and elastic self-scattering $SS\leftrightarrow SS$), the moment result can be non-vanishing (see Eq. (39) for instance), and an additional factor “1/2” needs to be added in $C[f_{h_{1}}]$ to avoid double-counting. Taking $SSS\leftrightarrow XX$ as an example, we have $S_{y}\equiv 3!$, and $S_{y}^{\prime}\equiv 2!$. Moreover, we have $\Delta\equiv 3,\Delta^{\prime}\equiv 0$ for considering the Boltzmann equation of the $S$ particles, while $\Delta\equiv 0,\Delta^{\prime}\equiv 2$ for $X$. Here the $1\pm f_{i}$ terms with plus and minus signs encode the influence due to Bose enhancement and Pauli blocking, respectively.

We are interested in reactions dominated by the phase space region where the average number of particles in a single-particle state is much less than 1, i.e., $1\pm f_{i}\simeq 1$, and thus approximate the distributions as

[TABLE]

Basically, this is a good approximation even for high or low temperature. The collision term can be then given as the following form,

[TABLE]

where $f_{i}(T_{i})\equiv f_{i}(E_{i},T_{i})$. It should be noted that the relativistic SM, $X$ and $S$ are defined by the different temperatures, $T,T_{X}$ and $T_{S}$, respectively. In our case, the dark matter and mediator are in thermal equilibrium, i.e. $T_{X}=T_{S}$, until their kinetic decoupling; we will further discuss this point in the following sections.

IV.1 The Boltzmann moment equation for the number densities of hidden sector particles

To get the coupled Boltzmann equations for the number densities of $X$ and $S$, we form the moment by multiplying Eq. (15) with “1” and integrating over the momentum space,

[TABLE]

In our case, the sufficiently large interactions within the hidden sector can maintain thermal equilibrium among the hidden sector particles, i.e., $T_{X}=T_{S}$, before the DM freezes out. The resulting equations of number densities are given by

[TABLE]

[TABLE]

where $K_{i}$ is the modified Bessel function of the second kind with $x_{S}\equiv m_{X}/T_{S}$ and $x\equiv m_{X}/T$, $n_{i}^{\rm eq}$ is the equilibrium number density with vanishing chemical potential, $\Gamma_{S}$ is the total decay width of $S$ into SM final states, and $\langle\sigma v\rangle_{i}$ and $\langle\sigma v^{2}\rangle_{i}$ are respectively the thermally averaged cross sections for $2\to 2$ and $3\to 2$ annihilation processes denoted by the subscript “i”; the details for these results are given in Appendices B and D. Note that only the terms involving $n_{S}^{\text{eq}}$ and, meanwhile, relevant to $S\to{\rm SM}_{i}~{}{\rm SM}_{i}$ and $SS\to{\rm SM}_{i}~{}{\rm SM}_{i}$ on the right hand side (RHS) of Eq. (23) are functions of the bath temperature, “$T$”, which is explicitly indicated, whereas the remaining ones appearing in Eqs. (22) and (23) are functions of “$T_{S}$” or “$T_{X}$”. Here and in the following analysis, we will use $T_{X}=T_{S}$ due to the fact that the hidden sector particles keep thermal equilibrium before the DM freeze-out.

Because the comoving number density of dark matter is conserved after freeze-out, we introduce the normalized yields for the hidden sector,

[TABLE]

with $x_{i}\equiv m_{X}/T_{i}$ and $x\equiv m_{X}/T$ being the temperature variables of the hidden sector particles and thermal bath, respectively, the yields $Y_{i}\equiv n_{i}(T_{i})/s(T)$ being the number density normalized by the total entropy density, and $g_{*}$ being the effectively total number of relativistic dof; see below for definition. Here $\langle\sigma v\rangle^{(0)}_{XX\to SS}$ is the leading approximation of $\langle\sigma v\rangle_{XX\to SS}$ which is s-wave. In the following discussion, we will simply use $y_{i}(x)\equiv y_{i}(x_{i};x)$. We use $x\equiv m_{X}/T$ as the evolution variable and trade the time derivative in the Boltzmann equations to be

[TABLE]

where the relativistic degrees of freedom, $h_{\rm eff}$ and $g_{*}^{1/2}\equiv\tilde{h}_{\rm eff}/g_{\rm eff}^{1/2}$ with

[TABLE]

are defined via

[TABLE]

Thus, we can further rewrite the Boltzmann equations to be

[TABLE]

[TABLE]

where $M_{\rm pl}\equiv(8\pi G)^{-1/2}=2.44\times 10^{18}$ GeV is the reduced Planck mass,

[TABLE]

and the equilibrium value of $y_{i}$ is given by

[TABLE]

Again, it should be noted that in Eq. (30) the terms involving $y_{S}^{\text{eq}}$’s and relevant to $S\to{\rm SM}_{i}~{}{\rm SM}_{i}$ and $SS\to{\rm SM}_{i}~{}{\rm SM}_{i}$ are functions only of “$x$”, as shown explicitly, while other quantities appearing Eqs. (29) and (30) are instead defined as functions of “$x_{S}$”, which are not shown explicitly, before freeze-out. In the following section, we will exhibit the evolution of $T_{S}/T$ as a function of $x$, i.e., a function of the bath temperature $T$.

The relic abundance is found to be

[TABLE]

which can be determined by matching the present-day DM relic abundance $\Omega_{\rm{DM}}=(0.1198\pm 0.0026)/h^{2}$ pdg2018 ; Ade:2013zuv , where $Y_{X}^{\infty}$ is related to $y_{X}^{\infty}=y_{X}(x\to\infty)$ (see also Eq. (24)), $s_{0}=2891~{}\text{cm}^{-3}$ is the visible entropy density today, $\rho_{c}=3H_{0}^{2}/(8\pi G)$ is the critical energy density, and $h\simeq 0.678$ is the Hubble constant $H_{0}$ of the present day in units of $100~{}\text{km}\,s^{-1}\text{Mpc}^{-1}$. $y_{X}^{\infty}$ is related to $x=x_{f}(\equiv m_{X}/T_{f})$ with $T_{f}$ the freeze-out temperature, and can be understood as follows. Well after DM freeze-out, which occurs at $x=x_{f}$, the Boltzmann equation in Eq. (29) can be approximated as

[TABLE]

Solving the equation, we get

[TABLE]

where we use the fact that the value of $y_{X}$ at $x=x_{f}$ is significantly larger than $y_{X}^{\infty}$, and we can approximate $T_{X}\approx T_{S}$ in the calculation (see Figs. 3, 4, and 5 for the temperature dependence in the next section).

IV.2 The Boltzmann moment equation for $T_{S}/T$

We consider the case that the DM can be kept in thermal equilibrium with the hidden scalar before DM freeze out, but may be highly decoupled from the SM thermal bath. Here we focus on the study about the temperature evolution of the hidden scalar $S$, and then discuss the DM temperature evolution after freeze out. The temperature evolution of nonrelativistic $S$ is affected by the following elastic scattering and (species) number changing interactions — (i) annihilation: $SS\leftrightarrow\text{SM~{}SM}$, (ii) elastic scattering: $S+\text{SM}\leftrightarrow S+\text{SM}$, (iii) cannibalization including $SSS\leftrightarrow SS$, $XSS\leftrightarrow XS$, $XXS\leftrightarrow SS$, $SSS\leftrightarrow XX$, and $XXX\leftrightarrow XS$, and (iv) decay: $S\leftrightarrow\text{SM~{}SM}$.

We adopt the definition of the temperature,

[TABLE]

which is suitable not only for the nonrelativistic case at low temperatures, $T<m_{S}$, but also for relativistic case at high temperature, $T\lesssim m_{S}/0.01$. The Boltzmann moment equation of the hidden scalar’s temperature can be formed by multiplying Eq. (15) with ${\bf p}_{S}^{2}/(3E_{S})$ and then integrating over the momentum space. Thus, we arrive at the form of the temperature evolution equation,

[TABLE]

where

[TABLE]

which is approximately to be “0” for nonrelativistic particles or “1” for ultra-relativistic ones. Here we have denoted the collision term as $[f_{S}\cdot{\bf p}_{S}^{2}/(3E_{S})]$ which is related to $C[f_{S}]$ (see Eq. (17) or (20)) with the replacement

[TABLE]

corresponding to a process $i(\text{initial~{}state})\leftrightarrow f(\text{final~{}state})$ with a prime for the final state. In the following, the collision term due to various interactions will be discussed term by term in details.

After the hidden sector is chemically decoupled from the bath, i.e., its number density production rate is less the expanding rate of the Universe, we have $n_{h}a^{3}={\rm constant}$ from Eqs. (22) and (23) if cannibalization can be neglected. Moreover, well after the cannibal epoch, if the nonrelativistic hidden sector is out of thermal equilibrium with the bath, we have $T_{h}a^{2}={\rm constant}$ as read from Eq. (37).

IV.2.1 The collision term due to $SS\leftrightarrow\text{SM}_{1}~{}\text{SM}_{2}$

The collision term resulting from $S(p_{S,1})S(p_{S,2})\leftrightarrow\text{SM}_{1}(p_{1^{\prime}})~{}\text{SM}_{2}(p_{2^{\prime}})$ is given by

[TABLE]

where $m!=2$ for identical final state particles or 1 otherwise, and we have used the energy conservation $E_{1^{\prime}}+E_{2^{\prime}}=E_{S,1}+E_{S,2}$, and the relation,

[TABLE]

Here the thermal average, for which the detailed description is provided in Appendix C, is given by

[TABLE]

where $v_{\rm lab}$ is the relative velocity measured in the laboratory frame where one of the incoming particles is at rest. For the present $S$-wave annihilation, the ratio of $\langle\sigma v\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle/(\langle\sigma v\rangle T)$ is $0.75\sim 1.05$ for $x(\equiv m_{X}/T)>2$ with $m_{X}=80$ GeV, and becomes unity in the limit $x\to\infty$. For a typical case with $\sigma v$ being constant, the ratio is equal to one. See also discussions in Appendix C.

IV.2.2 The collision term due to $S\leftrightarrow\text{SM}_{1}~{}\text{SM}_{2}$

The collision term resulting from $S(p_{S})\leftrightarrow\text{SM}_{1}(p_{1^{\prime}})~{}\text{SM}_{2}(p_{2^{\prime}})$ is given by

[TABLE]

where $m!=2$ for identical final state particles or 1 otherwise, $\Gamma_{S\to\text{SM}_{1}\text{SM}_{2}}$ is $S\to\text{SM}_{1}\text{SM}_{2}$ decay width, and

[TABLE]

which approaches to $n_{S}^{\rm eq}(T)$ in a nonrelativistic limit. The result of Eq. (43) is also correct if replacing $\Gamma_{S\to\text{SM}_{1}\text{SM}_{2}}$ with the relevant three- (or more-) body decay mode.

IV.2.3 The collision term due to cannibal annihilations among the hidden sector particles

The collision term arising from cannibal annihilations among the hidden sector particles contains the following processes: $SSS\leftrightarrow SS$, $XSS\leftrightarrow XS$, $XXS\leftrightarrow SS$, $SSS\leftrightarrow XX$, and $XXX\leftrightarrow XS$. Here, the result for $XSS\leftrightarrow XS$ will be shown, while for the others can be derived in a similar way. When the temperature of the hidden sector drops below $m_{X,S}$, the role of the cannibalization becomes important. If the hidden sector kinetically decouples from the bath at $T\lesssim m_{X,S}$, its temperature will decrease logarithmically with the cosmic scale factor during cannibalization (see Eq. (71) for discussion). The collision term resulting from $X(p_{X})S(p_{S,1})S(p_{S,2})\leftrightarrow X(p_{X^{\prime}})S(p_{S^{\prime}})$ is given by

[TABLE]

where we have used $p_{S^{(\prime)},i}^{\mu}=(E_{S^{(\prime)},i},{\bf p}_{S^{(\prime)},i})$ for the $S$ and $p_{X^{(\prime)}}^{\mu}=(E_{X^{(\prime)}},{\bf p}_{X^{(\prime)}})$ for the $X$, and have approximated three initial hidden particles ($XSS$) that annihilate or are produced in the nonrelativistic limit, i.e., $E_{S,i}\approx m_{S},E_{X}\approx m_{X}$, such that

[TABLE]

IV.2.4 The collision term due to elastic scattering: $S+\text{SM}\leftrightarrow S+\text{SM}$

Here we consider the elastic scattering, $S(p_{S})+{\rm SM}(k)\leftrightarrow S(p_{S^{\prime}})+{\rm SM}(k^{\prime})$, where $p_{S^{(\prime)}}^{\mu}=(E_{S^{(\prime)}},{\bf p}_{S^{(\prime)}})$, $k^{(\prime)\mu}=(\omega^{(\prime)},{\bf k}^{(\prime)})$ and “SM” stands for one of the relativistic SM particles that can participate the interaction. The collision term for this elastic scattering takes the following form,

[TABLE]

where the hidden scalar scattering with all relativistic SM fermions is taken into account. Under the typical condition $m_{S}\gg T\sim\omega$, this term can further reduce to a semi-relativistic Fokker-Planck-type equation Bringmann:2006mu ; Bringmann:2009vf ; Gondolo:2012vh ; Visinelli:2015eka ; Binder:2016pnr ; Binder:2017rgn ,

[TABLE]

where the momentum relaxation rate is given by

[TABLE]

for which the sum runs over all relevant relativistic SM species, and the differential elastic scattering cross section is

[TABLE]

with $|M_{Sf\to Sf}|^{2}$ the square of the scattering amplitude summed over initial and final spin states. Note that in Eq. (49), we have followed the approach given in Ref. Gondolo:2012vh to adopt the $t$-average matrix $(8{\bf k}^{4})^{-1}\int_{-4{\bf k}^{2}}^{0}dt(-t)d\sigma/dt$ due to the fact that the scattering amplitude squared in our case vanishes around $t=0$ in the relativistic limit $m_{f}\to 0$; therefore, it is unsuitable to take the result at zero momentum transfer of the $t$-channel as done in Ref. Bringmann:2006mu .

Taking into account the elastic scattering $Sf\to Sf$ which is dominated by the amplitudes with the SM Higgs or hidden scalar mediated in the $t$-channel, we find the amplitude squared to be

[TABLE]

with $N_{c}^{f}\equiv 3\,(1)$ for quarks (leptons), and the couplings shown in Eqs. (90), (93), (94), and (95). Averaging over $t$ for the scattering amplitude squared, we get the momentum relaxation rate to be

[TABLE]

Here the transferred momentum $t$ in the denominator of the amplitude squared is neglected in the calculation consistent with the requirement $-t<4|{\bf k}|^{2}\sim T^{2}\ll m_{S}^{2}$. Therefore, for the case with the resulting elastic decoupling temperature $\sim m_{X}$ as shown in the left panel of Figs. 3, 4, and 5, the kinetic transition rate should be overestimated, i.e., the true value of $x_{\rm el}(\equiv m_{X}/T_{\rm el})$ should be less than what is shown (see (iii) in Sec. V for the definition of $x_{\rm el}$). However, such overestimation does not affect our conclusions.

IV.2.5 The temperature evolution equation for the hidden scalar

After including all interaction terms, we arrive at the Boltzmann moment equation for the temperature of the hidden scalar,

[TABLE]

where

[TABLE]

which approaches 1 in a nonrelativistic $S$ limit. If considering the temperature below which the DM and hidden scalar are kinetically decoupled, i.e., $T_{X}\not=T_{S}$, we need to further include the following two terms to the RHS of Eq. (53),

[TABLE]

where the first and second terms are the kinetic energy-transfer rates by elastic scattering ($XS\leftrightarrow XS$) and by annihilation ($XX\leftrightarrow SS$), respectively. This impact will be discussed in (iv) of Sec. V. Some related results are collected in Appendix E.

By introducing the dimensionless variable,

[TABLE]

the above Boltzmann moment equation for the temperature $T_{S}$ can be recasted into an alternative form that will used in the analysis:

[TABLE]

where $\widetilde{\langle\sigma v\rangle}_{SS\to\sum_{ij}{\rm SM}_{i}{\rm SM}_{j}}(x_{k})\equiv\langle\sigma v\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle_{SS\to\sum_{ij}\text{SM}_{i}\text{SM}_{j}}(x_{k})/T_{k}$ and $a(m_{X}/m_{S})^{4}/x^{4}\equiv(2-\delta_{H})\gamma/H(T)$. In Eq. (57), we do not distinguish $T_{X}$ from $T_{S}$. An unphysical result may occur when $y_{S}(x_{S})\gg y_{S}^{\rm eq}(x)$ and other interaction terms become much smaller than that involving $1-\delta_{\Gamma}$ which originates from the interaction, $S\leftrightarrow\text{SM~{}SM}$. The $(1-\delta_{\Gamma})$ term may result in $y$ to be monotonically increasing with $x$, even though it is highly small and $1-\delta_{\Gamma}$, which is less than 0.05 for $x>20$, vanishes in nonrelativistic limit. Actually, when we consider completely coupled temperature evolutions for $X$ and $S$, such an unphysical result will be suppressed by including the sizable $XS\leftrightarrow XS$ and $XX\leftrightarrow SS$ interactions. For instance, the terms, which are rewritten from Eq. (55) and not shown in Eq. (57), are given by

[TABLE]

with $y_{X}^{T}\equiv T_{X}/T$. When the DM freezes out, the effect of the $(1-\delta_{\Gamma})$ term is completely washed out due to a much larger DM energy density compared with the hidden scalar one. A more detailed treatment is described as follows.

For most of cases (the exception one will be discussed below) that the hidden sector is decoupled from the bath at temperature $T\approx m_{S}$ (i.e., $y_{S}(x_{S})>y_{S}^{\rm eq}(x)$ for $x\gtrsim 1$), the coefficient given in the square bracket of Eq. (58) is always much larger than that of the $S\leftrightarrow\text{SM~{}SM}$ interaction term involving $1-\delta_{\Gamma}$. In other words, the evolution of $T_{S}$ will closely follow $T_{X}$ for $x\gtrsim 1$. Moreover, after the time that the hidden sector is decoupled from the SM bath, the DM plays as an effective reservoir with respect to the hidden scalar. Because the dof ratio $g_{X}/g_{S}=3$ and $m_{X}>m_{S}$, the $S$ temperature change rate due to the $(1-\delta)$ term can thus be reduced by about $90\%$ and $75\%$ for the cases $m_{S}=0.8\,m_{X}$ and $m_{S}=0.99\,m_{X}$, respectively. As such, when the cannibal interaction becomes inactive, both $T_{S}$ and $T_{X}$, following $dT_{S,X}/dt+(2-\delta_{H})HT_{S,X}\approx 0$ and $\propto a^{-2}$, will approximately evolve with the same temperature. Based on the above reasons, we will neglect $(1-\delta_{\Gamma})$ term, i.e., simply take $\delta_{\Gamma}=1$, for the case that the hidden sector is decoupled from the bath at $T\approx m_{S}$.

Three remarks are in order. First, we have neglected the reheating of the bath due to the out-of-equilibrium decay of $S$. A thorough treatment of the reheating is beyond the scope of the present work, since the bath temperature is not a suitable variable to take into account the thermal evolution and the total comoving entropy, which is no longer conserved, increases after the out-of-equilibrium $S$ decay occurs. Second, the uncertainty due to reheating of the bath can be realized as follows. After decoupling, compared with the SM radiation energy density $\rho_{\rm SM}\propto a^{-4}$, the energy density of $S$ evolves as $\rho_{S}\propto a^{-3}$. Thus, the value of $\rho_{S}/\rho_{\rm SM}$ is about $1/120$ at temperature $\sim m_{S}$, but is increased to be $\sim(1/120)\times m_{S}/T$ at a later time with temperature $T$. Adopting the sudden-decay approximation, the bath temperature change due to reheating is less than 3% for the case with $\alpha<5\times 10^{-7}$, and $\sim 10\%$ for $\alpha=1\times 10^{-7}$, where we have used the temperature result ($T_{\rm de}^{\rm out}$) of the out-of-equilibrium $S$ decay from the next section. Third, unlike the other ones, for the case with $\alpha=1\times 10^{-7}$ and $m_{S}=0.8m_{X}$, because the down scattering results in $y_{S}>y_{X}$, the term involving $1-\delta_{\Gamma}$ can be comparable with the coefficient given in the square bracket of Eq. (58) in a very short period of time just before DM freeze-out. It may result in a small temperature difference between $X$ and $S$ in such a short period of time. However, after that the DM and hidden scalar will still evolve with the almost same temperature, $\propto a^{-2}$ until and after their kinetic decoupling, as the other cases. Here, for simplicity, we will not consider such effect for this case.

On the other hand, as for the case that the $\text{SM SM}\leftrightarrow S$ interaction is strong enough to maintain the hidden scalar in thermal equilibrium with the bath until a temperature, which may be below the DM freeze-out temperature, the cooling rate of the hidden scalar via $SS\to XX$ will be soon larger than the $\text{SM SM}\to S$ heating rate (described by the $(1-\delta_{\Gamma})$ term) in magnitude after the dark matter is kinetically decoupled from the hidden scalar. To have a more precise estimate for the temperature ($T_{S,{\rm end}}=m_{X}/x_{S,{\rm end}}$) below that the hidden scalar cannot be in thermal equilibrium with the bath, we will model a term as given in Eq. (59) to show the possible thermal flow due to the temperature difference $T_{S}>T_{X}$, which will occur after $X$ is kinetically decoupled from the hidden scalar. In the analysis, we improve the numerical result recursively. The approach is described as follows. We will first use the approximation $\delta_{\Gamma}=1$ to get the numerical result. From that we then extract the DM kinetic decoupling temperature, $T_{X}^{\rm kd}$, which will be defined and discussed in the next section. Adopting the obtained $T_{X}^{\rm kd}$, we use the exact value of $\delta_{\Gamma}$ and include the $XX\leftrightarrow SS$ term (rewritten from the second term of Eq. (55)),

[TABLE]

in Eq. (57) to have the improved solution for $x_{S,{\rm end}}$, where $x_{X}^{\rm kd}\equiv m_{X}/T_{X}^{\rm kd}$. Here the $XS\leftrightarrow XS$ term is negligible for this case, and the relation $T_{X}\simeq T_{X}^{\rm kd}\cdot(T/T_{X}^{\rm kd})^{2}$ is used when the DM and hidden scalar are kinetically decoupled from each other.

V Numerical results for the thermal evolution of the hidden sector

We present the numerical results for thermodynamic evolutions of the normalized yields (proportional to co-moving number densities) and hidden sector temperatures using two sets of masses for the hidden sector: (i) $m_{X}=80$ GeV, $m_{S}=0.8m_{X}=64$ GeV, and (ii) $m_{X}=80$ GeV, $m_{S}=0.99m_{X}=79.2$ GeV, where the latter one is the nearly degenerate case. These two sets of the hidden masses are capable of generating one-step cascade DM annihilation spectra that provide a good fit to the observed GC gamma-ray excess which will be further discussed in the next section.

To illustrate how this secluded DM model could be highly decoupled from the SM bath, for these two mass sets, we take small mixing angle $\alpha$, i.e. to be (1) $1\times 10^{-5}$, (2) $1\times 10^{-6}$, (3) $5\times 10^{-7}$, and (4) $1\times 10^{-7}$, respectively. Meanwhile, in the analysis, we set the parameter $g_{\rm dm}$, of which the value is relevant to $XX\to SS$ annihilation cross section, to account for the observed DM relic abundance (see Eq. (33)). Our results are shown in Figs. 3, 4, and 5.

The equilibrium number densities, described by the Boltzmann equations given in Eqs. (29) and (30), can be maintained by interactions, including $XX\leftrightarrow SS$, $SS\leftrightarrow\text{ SM SM}$, $S\leftrightarrow\text{ SM SM}$, and $3\leftrightarrow 2$ hidden sector cannibalization. The temperature evolution of the hidden scalar is affected by the kinetic energy transfer by interacting with the bath and with the DM. Such a kinetic energy transfer can be generated from the hidden scalar number changing interactions and from the elastic scattering.

The main results, categorized in terms of temperature scales relevant to the transition phases during the evolution of the hidden sector, are summarized as follows.

1. (i)

$x_{f}\equiv m_{X}/T_{\rm f}$ is the usual freeze-out temperature variable. For $T<T_{\rm f}$, the comoving DM number density tend to be conserved. From Eq. (22), we can estimate the freeze-out temperature, below which the DM production rate from $SS\to XX$ is overtaken by the dilution rate, giving the relation

[TABLE]

where $T_{X}=T_{S}$ before freeze out are functions of $T$, and $T_{S}$ can be determined by solving numerically the Boltzmann equations, Eqs. (29), (30) and (57). In Figs. 3, 4, and 5, $x_{f}$ is denoted by the magenta dot in each plot of the left panel, while the corresponding asymptotic yield $y_{X}^{\infty}$, described by Eq. (34), is depicted by the horizontal dotted line on the right panel.

For $T>T_{\rm f}$, the rate on the left hand side (LHS) of Eq. (60) is larger than the expansion rate, resulting in the detailed balance $\langle\sigma v\rangle_{XX\to SS}(T_{X})\,(n_{X}(T_{X}))^{2}=\langle\sigma v\rangle_{SS\to XX}(T_{S})\,(n_{S}(T_{S}))^{2}$. This implies that

[TABLE]

so that $X$ and $S$ (with $T_{X}=T_{S}$) have the same chemical potential, $\mu_{X}=\mu_{S}$, i.e., the hidden sector is in chemical equilibrium, for which the chemical potential can be non-zero if the hidden scalar undergoes an out-of-equilibrium decay before the DM freezes out (see (vi) for related discussions).

In the right panel of Figs. 3, 4, and 5, the dashed curves exhibit the results with the Boltzmann suppression (i.e., with $\mu_{X}=\mu_{S}=0$) for the hidden sector if $T_{X,S}=T$, while dot-dashed curves follow the Boltzmann suppression at their true temperatures, where curves in magenta and blue colors are for the DM and hidden scalar, respectively. 2. (ii)

$x_{\rm ann,2}\equiv m_{X}/T_{\rm ann,2}$, plotted as the black dot in the left panel of Fig. 3, corresponds to the bath temperature $T=T_{\rm ann,2}$ set by

[TABLE]

where $T_{S}$, the temperature of $S$, is a function of $T$ determined by Eq. (53), and the second term of the RHS is described by Eq. (37). Here, considering the relevant channels $SS\leftrightarrow\sum_{ij}{\rm SM}_{i}~{}{\rm SM}_{j}$, the value of $(T_{S}/a^{3})\,d(n_{S}a^{3})/dt=T_{S}(dn_{S}/dt+3Hn_{S})\propto dy_{S}/dt$ is about zero for $x\lesssim x_{{\rm ann},2}$, if the thermal equilibrium between the hidden scalar and bath can be maintained by this annihilation reaction. When the bath temperature falls below $T_{\rm ann,2}$, the kinetic energy changing rate of the hidden scalar due to variations of its temperature (the first term of RHS) and density (the second term of RHS) during the cosmic expansion becomes larger in magnitude than the heating rate transported from the bath via annihilations $\sum_{ij}{\rm SM}_{i}~{}{\rm SM}_{j}\to SS$, such that this kind of interactions cannot play the role to keep thermal equilibrium with the SM bath for $T<T_{\rm ann,2}$. In Figs. 3(b-1), 4, and 5, $x_{\rm ann,2}$ is too small to be visible in plots. 3. (iii)

$x_{\rm el}\equiv m_{X}/T_{\rm el}$, sketched as the orange dot in the left panel of Figs. 3, 4, and 5, denotes the temperature below which the hidden sector is elastically decoupled from the bath. For $T>T_{\rm el}$, the heating rate of the hidden sector, which gains energy by the elastic scattering $S~{}{\rm SM}\to S~{}{\rm SM}$, is larger in magnitude than the cooling rate due to the Hubble expansion,

[TABLE]

where $(2-\delta_{H})\gamma T$ is the kinetic energy-transfer rate from the relativistic SM particles to a hidden scalar particle via elastic scattering. Unlike the inverse hidden scalar annihilation into relativistic SM particles, of which the rate is reduced by the Boltzmann suppression of the $S$ number density (with $\mu_{S}=0$) for $T\lesssim m_{S}$, the kinetic energy-transfer rate of $S~{}{\rm SM}\to S~{}{\rm SM}$ is proportional to the relativistic SM density which is not suppressed. Therefore, in most cases, we have $x_{\rm el}>x_{\rm ann,2}$. Two remarks are in order. First, $x_{\rm el}<x_{\rm ann,2}$ may occur, if the annihilation cross section is largely enhanced by the resonant $s$-channel SM-Higgs exchange with $s\sim(2m_{S})^{2}\sim m_{h}^{2}$. Second, compared with the elastic scattering, because the $2\to 2$ annihilation is much more sensitive to $S$-$f$-$\bar{f}$ coupling, which is proportional to $s_{\alpha}$, a small enough mixing angle $\alpha$ may result in the $2\to 2$ annihilation decoupling to occur significantly before elastic decoupling. 4. (iv)

$x_{\rm de}\equiv m_{X}/T_{\rm de}$ is respectively denoted by a purple “X” and by the vertical dotted (red) line in the left and right panels of Figs. 3(a), 3(b) and 5(a). $T_{\rm de}$ denotes the temperature below which not only the requirements, $n_{S}(T_{S})=n_{S}^{\rm eq}(T)$ and $T_{S}=T$, need to be satisfied, but also the heating rate, generated from the inverse decay: $\text{SM~{}SM}\to S$, is larger in magnitude than the rate needed to keep the $S$ particles in kinetic equilibrium with the bath during the cosmic expansion,

[TABLE]

where, on the RHS, the first term is the rate for temperature variation with an unchanged comoving number density of $S$, while the second term, describing the thermal energy rate due to a change of the $S$ comoving number density, results from $S\leftrightarrow\text{SM~{}SM}$ interactions, for which its contribution to (T_{S}/a^{3})\,d(n_{S}a^{3})/dt=T_{S}\big{(}dn_{S}(T_{S})/dt+3Hn_{S}(T_{S})\big{)} is relatively negligible when $T\approx T_{\rm de}$ for the case shown in Figs. 3(a) and 5(a). As such, if $T_{\rm de}>\text{min}(T_{\rm el},T_{\rm ann,2})$ is satisfied, the hidden scalar (as well as the dark matter) can be still maintained in thermal equilibrium with the bath (i.e., $n_{S}(T_{S})=n_{S}^{\rm eq}(T)$ and $T_{S}=T$) until a later time, $x_{\rm S,end}(\equiv m_{X}/T_{\rm S,end})$, which could be larger than $x_{f}$. For this case, to have a more precise estimate, we have included a term given in Eq. (59) to show the possible thermal flow due to the temperature difference $T_{S}>T_{X}$, which occurs after $X$ is kinetically decoupled from the hidden scalar. $x_{\rm S,end}$, denoted as the blue dot in Figs. 3(a-1) and 5(a-1), is determined by

[TABLE]

where the second and third terms of the left hand side are respectively the rates originating from $SS\to XX$ and $XS\to XS$ collision terms of the Boltzmann moment equation, while the second term of the RHS is the contribution from \big{(}dn_{S}(T_{S})/dt+3Hn_{S}(T_{S})\big{)}T_{S}, resulting mainly from $XX\to SS$ due to the fact that $n_{X}(T_{X})>n_{X}^{\rm eq}(T_{X})$ after the DM freezes out. The DM temperature after kinetic decoupling follows $T_{X}(a)\simeq T_{X}^{\rm kd}\cdot(a_{X}^{\rm kd}/a)^{2}$, with $a$ being the cosmic scale factor and $a_{X}^{\rm kd}$ being the corresponding value at $T_{X}=T_{X}^{\rm kd}$. For the cases shown in Figs. 3(a) and 5(a), this relation can be rewritten as $T_{X}=T_{X}^{\rm kd}(T/T_{X}^{\rm kd})^{2}$. After DM kinetically decouples from the hidden scalar, the evolution of $T_{X}/T$ is sketched as the red line in the left panel of Figs. 3, 4, and 5, where $T_{X}^{\rm kd}\equiv m_{X}/x_{X}^{\rm kd}$, depicted as the red dot, is the DM kinetic decoupling temperature, featuring $T_{X}^{\rm kd}\leq T_{\rm f}$. A detailed discussion for $T_{X}^{\rm kd}$ will be given in Appendix E. Numerically, we obtain $x_{\rm S,end}\approx x_{X}^{\rm kd}$, as seen from Figs. 3(a-1) and 5(a-1).

On the other hand, provided that the relation given in Eq. (64) is satisfied but with $T_{\rm f}\ll T_{\rm de}<\text{min}(T_{\rm el},T_{\rm ann,2})$, i.e., $x_{f}\gg x_{\rm de}>\text{max}(x_{\rm el},x_{\rm ann,2})$, the hidden sector may be kinetically decoupled from the bath at $T\lesssim\text{min}(T_{\rm el},T_{\rm ann,2})$, such that at $x\gtrsim x_{\rm de}^{\rm out}$ (see (vi) for the definition), the $S$ first undergoes an out-of-equilibrium decay with a rate much larger than its inverse production rate $\text{SM~{}SM}\to S$ due to the fact that $n_{S}\gg n_{S}^{\rm eq}$, resulting in the RHS of Eq. (64) to be less than zero,

[TABLE]

It is interesting to note that, as the time evolves, the $n_{S}$ Boltzmann equation gives $n_{S}(T_{S})\to n_{S}^{\rm eq}(T)$, which is displaced in Figs. 3(b-2) and 5(b-2), so that as shown in Fig. 3(b-1) the requirement of Eq. (64) is possible to be met and $x_{\rm end}$ can thus exist. 5. (v)

$x_{\rm c}\equiv m_{X}/T_{c}$ corresponds to the bath temperature $T=T_{c}$, illustrated by the green dot in the left panel of Figs. 3, 4, and 5, and described by

[TABLE]

where $K_{3\to 2}$ shown in Eq. (53) is the kinetic energy released in a relevant $3\to 2$ process involving nonrelativistic $S$ and $X$. For $T>T_{c}$, this number changing interaction maintains the hidden sector, which is undergoing cannibalism, in kinetic equilibrium and in chemical equilibrium with $\mu_{X,S}=0$: $n_{X}(T_{S})\to n_{X}^{\rm eq}(T_{S})$, $n_{S}(T_{S})\to n_{S}^{\rm eq}(T_{S})$. We are interested in the cases, as given in Figs. 3(b), 4(a), 4(b) and 5(b), 5(c), 5(d), that the hidden sector is decoupled from the thermal bath and evolves with different temperature independently, before it becomes nonrelativistic. For these cases with $T\lesssim m_{S,X}$, the total comoving entropy density of the hidden sector tends to be conserved before the $S$ decay occurs. Moreover, during the cannibal process, the entropy density ratio for the SM, $s_{\rm SM}=(2\pi^{2}/45)h_{\rm SM}^{\rm eff}(T)T^{3}$, to the hidden sector, $s_{h}=(2\pi^{2}/45)h_{\rm h}^{\rm eff}(T_{S})T_{S}^{3}$, is constant, where $h_{\rm SM}^{\rm eff}$ and $h_{\rm h}^{\rm eff}$ are the effectively relativistic degrees of freedom of the SM and hidden sector, respectively. Thus, we find

[TABLE]

where

[TABLE]

and $s_{\rm SM}/s_{h}\simeq 30$ for decoupling at $T\sim m_{X}$. From this scenario of entropy conservation, the temperature ratio $T_{h}/T$ increases due to cannibalization and follows the dotted gray curve, illustrated on the left panel of Figs. 3, 4, and 5. The hidden sector temperature will deviate from the dotted curve earlier if the out-of-equilibrium decay of $S$ takes place before the end of cannibalization. As time evolves such that $T<T_{c}$, the cannibal process is inactive, and the out-of-equilibrium number densities of the hidden sector starts to be exponentially depleted (see Figs. 4(a-2), 4(b-2) and 5(c-2), 5(d-2)).

For $S$ with a lifetime longer than the inverse Hubble rate and during its epoch of cannibalization, the conservation of the total comoving entropy for the hidden sector gives $s_{h}a^{3}\simeq(\rho_{h}/T_{S})a^{3}\simeq$ constant. Therefore, the comoving number density of hidden sector as well as its temperature decreases logarithmically with the scale factor, i.e. logarithmically with the bath temperature parameter $x$,

[TABLE]

where $s=s_{\rm SM}+s_{h}$, and $a_{\rm out,h}$ (the cosmic scale factor) and $x_{\rm out,h}$ (the bath temperature parameter) correspond to the values at which the hidden sector starts to be out of equilibrium with the bath. The logarithmic dependence of the comoving number densities for $X$ and $S$ can be seen from Figs. 3(b-2), 4(a-2), 4(b-2) and 5(b-2), 5(c-2), 5(d-2). 6. (vi)

$x_{\rm de}^{\rm out}\equiv m_{X}/T_{\rm de}^{\rm out}$, denoted by the purple square in the left panel and by the vertical dashed (red) line in the right panel of Figs. 3, 4, and 5, is the temperature below which the $S$ undergoes an out-of-equilibrium decay, i.e., the second term in the RHS of Eq. (64) (see also Eq. (66)) is much larger than the term in the LHS of Eq. (64) in magnitude due to the fact that $n_{S}(T_{S})\gg n_{S}^{\rm eq}(T)$ at $T=T_{\rm de}^{\rm out}$. For this case, corresponding to a much smaller mixing angle as that given with $\alpha=5\times 10^{-7}$ or $1\times 10^{-7}$ in this paper, because $n_{S}(T_{S})\,T_{S}\gg n_{S}^{\rm eq}(T)\,T$ at a later time with $x_{f}\gg x>\text{max}(x_{\rm el},x_{\rm ann,2})$, we thus approximately define $x_{\rm de}^{\rm out}$ from Eq. (64) to satisfy

[TABLE]

where $x_{S,{\rm de}}^{\rm out}$ is the value of $x_{S}$, corresponding to $x=x_{\rm de}^{\rm out}$. When $x=x_{\text{de}}^{\rm out}$, the hidden scalar thus starts to undergo out-of-equilibrium decay at the cosmological time, (2H)^{-1}\approx\Gamma_{S}^{-1}\big{[}\frac{K_{1}(x_{S,{\rm de}}^{\rm out}\cdot m_{S}/m_{X})}{K_{2}(x_{S,{\rm de}}^{\rm out}\cdot m_{S}/m_{X})}\cdot 2/(2-\delta_{H})\big{]}^{-1}\approx\Gamma_{S}^{-1}. See also the related discussion in Sec II.

Because the number changing interactions between $S$ and $X$ affect the $S$ number density for a longer time interval, thus the estimation of the value of $x_{\rm de}^{\rm out}$ needs to be further improved. For simplicity, here we neglect the cannibal interaction between $S$ and $X$. Such an interaction results in a logarithmic dependence of the hidden sector comoving number densities on the temperature variable $x$. In the plots, we will use the definition for $x=x_{\rm de}^{\rm out}$ which satisfies $y_{S}(x_{\rm de}^{\rm out})/y_{S}(1)\simeq e^{-1}$ Yang:2018fje with the initial value $y_{S}(1)=y_{S}^{\rm eq}(1)$. The value of $x=x_{\rm de}^{\rm out}$ is estimated as follows. For the case with $m_{X}-m_{S}$ sizable enough (e.g. $m_{S}=0.8m_{X}$), the down-scattering rate, $XX\to SS$, can be significantly larger than the up-scattering rate, $SS\to XX$, such that after a sufficient time at $x=x_{\text{de}}^{\rm out}\gg 1$, we have $y_{X}^{\rm eq}/y_{S}^{\rm eq}\ll 1$ and $y_{X}\ll y_{S}$. Therefore we set the effective initial $S$ yield to be $y_{S}^{\rm in}(1)\approx y_{X}(1)+y_{S}(1)\approx y_{X}^{\rm eq}(1)+y_{S}^{\rm eq}(1)=\kappa y_{S}^{\rm eq}(1)$ with $\kappa\equiv[3(m_{X}/m_{S})^{2}K_{2}(1)/K_{2}(m_{S}/m_{X})+1]$, and approximate Eq. (30) as

[TABLE]

with

[TABLE]

Here the approximation in the last step of Eq. (74) is reasonable because $g_{*}^{1/2}/h_{\text{eff}}$, which is 0.109 for $T=m_{X}/20$ and 0.125 for $T=m_{X}/150$, weakly depends on $T_{\rm f}$ in the present study. Solving this equation, we obtain

[TABLE]

and $x_{\rm de}^{\rm out}\simeq\sqrt{1+2(1+\ln\kappa)/C}$. In terms of the cosmic time variable of the radiation dominated epoch,

[TABLE]

the solution of the normalized yield can be rewritten as $y_{S}\simeq\kappa y_{S}^{\rm eq}(1)e^{-\Gamma_{S}t}$. If $\kappa=1$, we have $t(\approx(2H)^{-1})=\Gamma_{S}^{-1}$ at $x=x_{\text{de}}^{\rm out}$, consistent with the that given in Eq. (72). For the case with $m_{X}\approx m_{S}$ (e.g. $m_{S}=0.99m_{X}$), summing Eqs. (29) and (30), we have

[TABLE]

Using the initial conditions: $y_{X}(1)=y_{X}^{\rm eq}(1)$, $y_{S}(1)=y_{S}^{\rm eq}(1)$, and $y_{X}^{\rm eq}(1)/y_{S}^{\rm eq}(1)=g_{X}/g_{S}=3$, and approximating $y_{X}+y_{S}\approx 4y_{S}$, we have

[TABLE]

with $x_{\rm de}^{\rm out}\simeq\sqrt{1+8/C}$. The solution can be given by $y_{S}\simeq y_{S}^{\rm eq}(1)e^{-\Gamma_{S}t/4}$. Note that a longer-lived $S$ will result in a larger $x_{f}$, i.e., a larger $y_{X}^{\infty}$, so that, to have a correct relic density, the dark matter annihilation cross section is generally boosted above the conventional WIMP value. This point will be further discussed in the next section.

For the case that the hidden sector is kinetically decoupled from the bath at $T\sim m_{X,S}$, the ending value ($x_{c}$) of cannibalization depends on the magnitude of $x_{\rm de}^{\rm out}$ since the number density of $S$ is exponentially depleted during decay. In Figs. 4(b) and 5(d), we show the cases with $x_{\rm de}^{\rm out}\sim x_{c}$, for which, when $S$ decays out of equilibrium, the $X$ and $S$ densities are exponentially depleted, instead of following the Boltzmann suppression with a zero chemical potential (see the dot-dashed curves in Figs. 4(b-2) and 5(d-2)). Note that, for this case, $X$ and $S$ are still in chemical equilibrium but with non-zero chemical potential before freeze-out. Moreover, it is also interesting to note that for $T<T_{c}$, we have $T_{X}=T_{S}$ and $T_{S}=T_{S}^{\rm c}\cdot(a^{\rm c}/a)^{2}$ even after thermal decoupling, where $a$ is the cosmic scale factor and $a^{\rm c}$ is its corresponding value at $T_{S}^{\rm c}=T_{S}(T=T_{c})$.

VI Discussions

Since we have considered the secluded vector dark matter model with DM mass $\sim{\cal O}\text{(80 GeV)}$ as an example to exhibit the thermodynamic evolution of the hidden sector, the related parameters in this model should be very likely constrained by the astrophysical and cosmological measurements. Therefore, before making conclusion, let us discuss the parameter space that can fit to the excess of GeV-scale gamma-rays emitted from the GC region and evade constraints from dwarf spheroidal observation, cosmic microwave background, direct detection, and big bang nucleosynthesis.

The differential gamma-ray flux from the one-step cascade DM annihilations is described by

[TABLE]

where $\langle\sigma v\rangle_{\rm LV}$ is the DM annihilation cross section in the low-velocity limit (consistent with $T\to 0$), $(dN_{\gamma}^{f}/dE)_{X}$ is the prompt gamma-ray spectrum produced per annihilation with final state $f$ in the DM rest frame, and the J-factor is the integral along the line of sight (l.o.s.) and over the region of interest (ROI) denoted by the solid angle $\Delta\Omega$. We use a Galactic DM density distribution which is a function of $r$, the distance to the GC, and parametrized by a generalized Navarro-Frenk-White (gNFW) profile Navarro:1995iw ; Navarro:1996gj ,

[TABLE]

where we adopt $r_{s}=20$ kpc, $r_{\odot}=8.5$ kpc, $\gamma=1.2$ and $\rho_{\odot}=0.4$ GeV/cm3 as the canonical inputs. Here “$-\gamma$” is the inner log slope of the halo density near the GC, and $\rho_{\odot}$ is the local DM density at a distance of $r_{\odot}$ from the GC. The gamma-ray spectrum in the DM rest frame can be expressed in terms of that given in the rest frame of the metastable mediator ($S$) by means of one-step Lorentz boost Elor:2015tva (see also Eq. (10) in Ref. Yang:2017zor ), where we use PPPC4DMID result Cirelli:2010xx ; Ciafaloni:2010ti to described the gamma spectra that are generated from the final state SM particle pair in the $S$ decay at rest. As for the parameter region with $m_{V}<m_{S}<2m_{V}$ (with $V\equiv W$ or $Z$), the three-body decay channel $S\to VV^{*}\to Vf_{1}\bar{f}_{2}$ is kinematically open and becomes much more important when $m_{S}$ is close to $2m_{V}$ (see also Fig. 1). In the $S$ rest frame, the gamma-ray spectrum generated from three-body decay channels can be obtained by boosting the gamma-ray spectra produced from $V$ at rest and from $V^{*}$ at rest, respectively. Because the description for this part, relevant to the parameter region of the gamma-ray line emission, is sophisticated and does not affect the conclusion of this paper, we will thus defer the details in a future study.

Using the GC excess result extracted by Calore, Cholis, and Weniger (CCW) Calore:2014xka from the study of Fermi-LAT Pass 7 data, where the gamma-ray spectrum covers the energy range between 300 MeV to 500 GeV in the $40^{\circ}\times 40^{\circ}$ square ROI around the Galactic center with latitude $|b|\leq 2^{\circ}$ masked out, we do the goodness-of-fit with a $\chi^{2}$ test statistic for the DM mass and annihilation cross section $\langle\sigma v\rangle_{\rm LV}$. In Fig. 6, two ratio values of $m_{S}=0.8m_{X}$ and $m_{S}=0.99m_{X}$ are used to show the GC excess result, where, taking $\rho_{\odot}=\text{0.4 GeV/cm}^{3}$ and $\gamma=1.2$, the blue regions with solid, dashed and dotted boundaries respectively satisfy p-value $\geq$ 0.3, 0.15, and 0.05, corresponding to $\chi^{2}/{\rm dof}=$ 24.9/22, 28.8/22, and 33.9/22. The best fit is denoted as the blue dot with p-value = 0.46 or 0.42, corresponding to $\chi^{2}_{\rm min}/{\rm dof}=$ 22.0/22 or 22.7/22, for the case of $m_{S}=0.8m_{X}$ or $m_{S}=0.99m_{X}$.

Further allowing variation of $\rho_{\odot}\in$ [0.25, 0.85] GeV/cm3 and $\gamma\in[1.15,1.25]$, the value of $\langle\sigma v\rangle_{\rm LV}$ can be raised (or lowered) extremely by a factor of 2.94 (or 0.194). For illustration, in Fig. 6, we also show the GC allowed region in orange color corresponding to $\rho_{\odot}=\text{0.25 GeV/cm}^{3}$ and $\gamma=1.15$, while that in pink color corresponding to $\rho_{\odot}=\text{0.85 GeV/cm}^{3}$ and $\gamma=1.25$. In Fig. 7, the GC fit together with other constraints is redrawn on the ($m_{X},g_{\rm dm}$) plane, where a larger $g_{\rm dm}$ is needed to account for the data for the nearly degenerate case because $\langle\sigma v\rangle_{\rm LV}$ vanishes in the limit $m_{S}\to m_{X}$.

We remark that a newer Pass 8 Fermi data set was analyzed in Ref. Linden:2016rcf , in which the authors showed that the considerable difference between Fermi Pass 7 and Pass 8 data appears only at low energies which might be due to the modeling for the point sources in various datasets TheFermi-LAT:2017vmf ; Linden:2016rcf .

In Fig. 6, the relic density of the conventional WIMP dark matter is accounted for by the narrow gray range, while above the gray range the nonconventional WIMP scenario, showing a boosted annihilation cross section, can be satisfied. The result can be also easily read from Fig. 8, where the nonconventional WIMP scenario corresponds to a small mixing angle $\alpha\lesssim 2\times 10^{-6}$, for which the hidden sector has kinetically decoupled from the thermal bath before it becomes nonrelativistic. The resulting nonconventional WIMP DM annihilation cross section that can account for the correct relic density is significantly boosted above the conventionally thermal WIMP value for $\alpha\lesssim 6\times 10^{-7}$ (or $\alpha\lesssim 1\times 10^{-6}$) if $m_{X}=80$ GeV, $m_{S}=64$ GeV (or $m_{X}=80$ GeV, $m_{S}=79.2$ GeV).

Fig. 6 shows constraints from the Fermi gamma-ray observations of dwarf spheroidal galaxies (dSphs) and the measurement of the cosmic microwave background (CMB). For the dSphs constraint, we have performed a combined likelihood analysis using the 6-year Fermi-LAT data of 28 confirmed and 17 candidate dSphs for gamma-ray energies within 500 MeV to 500 GeV Fermi-LAT:2016uux ; FermiLatDesData . In the likelihood analysis, we adopt the spectroscopically determined nominal J-factor for the individual target along with its error when possible, or use a predicted value from the distance scaling relationship with an uncertainty of 0.6 dex, otherwise Fermi-LAT:2016uux . See the detailed description in Ref. Yang:2018fje for the likelihood analysis. We also show the dSphs projection sensitivity denoted by the dashed red line by assuming that the 15-year data can be collected from 60 dSphs. For the CMB constraint which is complementary to that determined from dSphs observations, Planck sets a bound from temperature and polarization data (TT, TE, EE+lowP) at recombination to be Ade:2015xua

[TABLE]

where $\langle\sigma v\rangle_{\rm CMB}\simeq\langle\sigma v\rangle_{\rm LV}$ for s-wave DM annihilation, and the efficiency factor is

[TABLE]

Here, we use $f_{\rm eff}^{\gamma,e^{-}}(E)$ curve results suited for the “3 keV” baseline prescription shown in Ref. Slatyer:2015jla . Moreover, as the previous study for the GC excess, $(dN_{\gamma,e^{-}}/dE)_{X}$ generated from the one-step cascade DM annihilation is the photon/electron energy spectrum that can be obtained by boosting the spectra provided in PPPC4DMID. The current bound obtained from the CMB analysis seems to be much weaker than that from the Fermi-LAT dSphs data (see also Fig. 8).

In the present model, compared with the SM, we have 4 additional parameters, $m_{X},m_{S},\alpha$ and $g_{\rm dm}$. As shown in Figs. 3, 4, and 5, having the chosen masses for the DM and hidden scalar, and giving the magnitude of $\alpha$, we can fine-tune the value of $g_{\rm dm}$ in the numerical analysis of Boltzmann equations to obtain $y_{X}^{\infty}$, which matches the observed relic abundance determined by Eq. (27). Note that, in Eq. (27), $\langle\sigma v\rangle^{(0)}_{XX\to SS}$ is a function of $g_{\rm dm}$, and its $\alpha$-dependence is negligible in our study.

In order to have a more comprehensive understanding of the phenomenological constraints on the secluded DM that could exhibit a boosted DM annihilation cross section, as discussed in the previous section, using the two ratio values of $m_{S}=0.8m_{X}$ and $m_{S}=0.99m_{X}$ with $m_{X}=80$ GeV, we display the correct relic abundance as the black curve on the ($\alpha,y_{X}^{\infty}$) plane in Fig. 8. Since the $S$ decay width $\Gamma_{S}$ is a function of $\alpha$, we label its corresponding values on the top of the plots. On the other hand, for a obtained $y_{X}^{\infty}$, we can get $x_{f}$ from the relation given by Eq. $(\ref{eq:freeze-out})$, and further have the corresponding value $\langle\sigma v\rangle_{\rm LV}(\equiv\langle\sigma v\rangle^{(0)}_{XX\to SS})$ from Eq. (33) or from the value of $g_{\rm dm}$. The dependence of $\langle\sigma v\rangle_{\rm LV}$ on $\alpha$ is weak and thus neglected in the plots. All of the corresponding quantities are labeled in the plots. The range favored by observed features of the GC excess with variation of $\rho_{\odot}\in$ [0.25, 0.85] GeV/cm3 and $\gamma\in[1.15,1.25]$ is given in between the two horizontal dot-dashed (purple) lines.

As shown in Fig. 8, the LZ projected sensitivity can only reach the RHS of the dotted (blue) curve corresponding to the thermal WIMP region with $\alpha\gtrsim 3.7\times 10^{-3}$ for $(m_{X},m_{S})=$ (80 GeV, 64 GeV), or for $\alpha\gtrsim 5.5\times 10^{-3}$ for $(m_{X},m_{S})=$ (80 GeV, 79.2 GeV). If $\alpha$ is larger than the value denoted by the vertical dashed (magenta) line which corresponds to $T_{\rm de}=\text{min}(T_{\rm el},T_{\rm ann,2})$, the hidden sector particles can be well in the chemical and thermal equilibrium with the bath before freeze out. Nevertheless, the case of $(m_{X},m_{S})=$ (80 GeV, 64 GeV) with $\alpha\lesssim 6\times 10^{-7}$, or $(m_{X},m_{S})=$ (80 GeV, 79.2 GeV) with $\alpha\lesssim 1\times 10^{-6}$, clearly exhibits the boosted annihilation cross section capable of accounting for the correct relic density. The 95% C.L. limit from CMB is denoted as the dot-dashed (brown) line, while the 95% C.L. limit and project sensitivity for dSphs observations are shown as horizontal solid (red) and long-dashed (red) lines, respectively. The boosted annihilation cross section is thus stringently constrained by the current dSphs observations. The secluded DM scenario can be further tested by the dSphs projection.

Finally, we discuss the bound on the lifetime of the hidden scalar from the big bang nucleosynthesis constraint, i.e., a lower bound on the mixing coupling $\alpha$. It has been shown and discussed in Refs. Kawasaki:1999na ; Kawasaki:2000en ; deSalas:2015glj ; Hasegawa:2019jsa that the late-time entropy production by the massive particle decay could induce cosmological effects at the time $\sim\Gamma_{S}^{-1}\sim 1$ sec. We consider the case with a long-lived massive $S$ which starts to decouple from the SM bath at the temperature below $m_{S}$. After decoupling, the energy density of the nonrelativistic hidden scalar then decreases as $\rho_{S}\propto a^{-3}$, while the SM radiation energy density $\rho_{\rm SM}$ scales as $a^{-4}$. As a result, $\rho_{S}/\rho_{\rm SM}$, scaling linearly with $a$, is about $1/120$ at $T\sim m_{S}$, but becomes $\sim(1/120)\times m_{S}/T$ at a later time with temperature $T$. In other words, the universe can be rapidly dominated by the nonrelativistic hidden sector particles if $S$ is long-lived. When the hidden scalar decays out-of-equilibrium into SM particles, the universe becomes radiation-dominated again and the SM bath experiences the reheating due to the large entropy injection. While photons and charged leptons are quickly thermalized during reheating, the weakly interacting neutrinos are slowly produced in the bath. Since neutrinos decouple from the thermal bath after $T\lesssim 1.5$ MeV, they would not be well thermalized if the reheating temperature $T_{\rm RH}\sim{\cal O}(1)$ MeV. See the example shown in Fig. 3 of Ref. Kawasaki:2000en , where each neutrino follows the perfect Fermi-Dirac distribution very well for $T_{\rm RH}=10$ MeV, while the distributions are not in thermal equilibrium form for $T_{\rm RH}=2$ MeV (See also Fig. 4 in Ref. Hasegawa:2019jsa , where neutrino self-interaction and oscillation are included). If $T_{\rm RH}\lesssim 7$ MeV, the effective number of neutrino species $N_{\rm eff}$ becomes smaller than three (see Fig. 4 in Ref. Kawasaki:2000en or Fig. 1 in Ref. Hasegawa:2019jsa for reference). Note that in our model there are no additional relativistic particles present before or after BBN, although such particles can large the value of $N_{\rm eff}$. The deficit of the neutrino distribution functions due to the insufficient thermalization will decrease the interaction rate between proton and neutron, so that the helium nucleon fraction $Y_{p}\equiv 4n_{\rm He}/n_{b}$, and the deuterium ratio $D/H$ are thus enhanced.

Following Refs. Kawasaki:2000en ; Hasegawa:2019jsa , we define the reheating temperature of the SM bath to be $\Gamma_{S}=3H(T_{\rm RH})$. Using the approximation,

[TABLE]

$T_{\rm RH}$ can be related to the decay width of $S$ as

[TABLE]

where we have used $g_{\rm eff}=43/4$. Here, we quote a lower bound $T_{\rm RH}\gtrsim 4.1$ MeV at 95% C.L., corresponding to $m_{S}=$10 GeV$-$100 TeV, from the $Y_{p}+D/H$ analysis in the case of 100% hadronic decay (see Figs. 12 and 13 in Ref. Hasegawa:2019jsa ), which is suitable for our model. Further considering the neutrino self-interaction and oscillation, the same reheating bound is also required by $N_{\rm eff}=3.15\pm 0.23$ from Planck report Ade:2015xua (see Fig. 1 in Ref. Hasegawa:2019jsa ). Therefore, from Eq. (84), we can obtain the BBN constraint on the $S$ width to be $\Gamma^{-1}\lesssim 0.03$ sec. As such, we have $\alpha\gtrsim 1.13\times 10^{-10}$ for $(m_{X},m_{S})=$ (80 GeV, 64 GeV), or $\alpha\gtrsim 1.03\times 10^{-10}$ for $(m_{X},m_{S})=$ (80 GeV, 79.2 GeV).

VII Conclusions

Using the secluded vector dark matter model, we have presented a comprehensive study on thermodynamic evolutions of the hidden sector particles from the first principle. We have solved numerically the coupled Boltzmann moment equations for number densities and temperature evolutions of the hidden sector particles. Our formalism can be easily extended to a general secluded dark matter model.

Taking two mass sets: (i) $m_{X}=80$ GeV, $m_{S}=0.8m_{X}=64$ GeV, and (ii) $m_{X}=80$ GeV, $m_{S}=0.99m_{X}=79.2$ GeV, we have shown the detailed thermodynamics for which, while the dark matter in thermal equilibrium with the hidden scalar is secluded from the visible sector with small interaction rates in agreement with the limit from the direct detection and collider experiments, the hidden sector can be either in thermal equilibrium or out of equilibrium with the bath before the DM freezes out. The results are briefly summarized as below. More details about the thermodynamics of the hidden sector have been given in Sec. V.

For the case satisfying $T_{\rm de}\geq\text{min}(T_{\rm el},T_{\rm ann,2})$, the kinetic decoupling of elastic scattering $S\,\text{SM}\leftrightarrow S\,\text{SM}$ and/or annihilation $SS\leftrightarrow\text{SM SM}$ occurs only when the bath temperature is below $T_{\rm de}$ at which the heating rate of the hidden sector generated from the inverse decay $\text{SM~{}SM}\to S$ starts to overcome the dilution rate due to the cosmic expansion. As such, the nonrelativistic hidden sector can keep thermal equilibrium with the bath until freeze-out. Therefore, the DM is consistent with the conventional WIMP scenario, but can easily evade the searches from the colliders and direct detections (e.g. projected LZ measurement) for a small mixing angle $2\times 10^{-6}\lesssim\alpha\lesssim 4\times 10^{-3}$ as in the present model.

On the other hand, for the case that the hidden sector starts to be kinetically decoupled from the thermal bath at $T\sim m_{X,S}$ due to its weak couplings to the SM particles, the nonrelativistic hidden sector will first undergo a cannibal epoch, during which the total comoving entropy density of the hidden sector is approximately conserved before $S$ decays out of equilibrium. When out-of-equilibrium $S$ decay occurs, the hidden sector particles $X$ and $S$ are still in chemical equilibrium, but their densities, instead of following Boltzmann suppression with zero chemical potential, are exponentially depleted with non-zero chemical potential until freeze-out. We have shown that having a small mixing angle $\alpha\lesssim 6\times 10^{-7}$ (or $\alpha\lesssim 1\times 10^{-6}$) which corresponds to $m_{X}=80$ GeV, $m_{S}=64$ GeV (or $m_{X}=80$ GeV, $m_{S}=79.2$ GeV), the secluded DM annihilates into “long-lived” hidden mediators which later decay out of equilibrium with the bath, such that the resulting nonconventional WIMP-like DM annihilation cross section accounting for the observed relic density is boosted above the conventionally thermal WIMP value.

For the experimental constraints, we have shown the parameter space which yields a good fit to the GC excess data and is compatible with the LZ projected sensitivity, BBN bound, Planck CMB measurement and Fermi dSphs observation. Moreover, we expect that Fermi-LAT 15-yr dSph observations can explore the parameter region of the correct relic density described not only by the nonconventional WIMP scenario but also, if the DM and hidden scalar are not well degenerate, by the conventional WIMP one.

Acknowledgements.

This work was supported in part by the Ministry of Science and Technology, Taiwan, under Grant Nos. 105-2112-M-033-005 and 108-2112-M-033-002.

Appendix A The partial decay widths of the hidden mediator $S$

The main partial decay widths of the hidden scalar $S$ with mass $\lesssim$ 130 GeV are given by

[TABLE]

where $N_{c}^{q(\ell)}\equiv 3\,(1)$ for the quark (lepton), $G_{F}=(\sqrt{2}v_{H}^{2})^{-1}$, the coupling $g_{Sff}=-s_{\alpha}m_{f}/v_{H}$, $K_{q}=1+5.67\alpha_{s}(\mu)/\pi$ and $K_{g}=1+(215/12)\alpha_{s}(\mu)/\pi$ are the NLO QCD corrections Djouadi:2005gj , $\delta^{\prime}_{W}=1$, $\delta_{Z}^{\prime}=\frac{7}{12}-\frac{10}{9}\sin^{2}\theta_{W}+\frac{40}{27}\sin^{4}\theta_{W}$,

[TABLE]

with $x\equiv m_{V}^{2}/m_{S}^{2}$ Keung:1984hn ; Djouadi:2005gi , and $f_{S}(\tau)=[1+(1-\tau)f(\tau)]$ with

[TABLE]

Here, we will take the scale $\mu=m_{S}/2$.

Appendix B $2\to 2$ annihilation cross sections

B.1 The annihilation process for $XX\to SS$

In this secluded DM case, the relic density is determined by the thermally averaged annihilation cross section $\langle\sigma v_{\text{M\o l}}\rangle_{XX\to SS}$ which is also relevant to the indirect detection searches, where $v_{\text{M\o l}}$ is the Møller velocity. In the text, we have used $\sigma v\equiv\sigma v_{\text{M\o l}}$ for simplicity. The value of $\langle\sigma v_{\text{M\o l}}\rangle_{XX\to SS}$ equals to $\langle\sigma v_{\rm lab}\rangle_{XX\to SS}$ which is the result calculated in the rest frame of one of the incoming particles.

The diagrams for the $XX\to SS$ process are depicted in Fig. 9, where the $s$-channel annihilation via $h$ is negligible and does not shown. The resulting cross section is given by

[TABLE]

where $s$ is the center-of-mass energy squared, and

[TABLE]

Using the above result, the thermally averaged annihilation cross section for $T\lesssim 3m_{X}$ can be obtained by calculating Gondolo:1990dk ,

[TABLE]

with $K_{1,2}$ being the modified Bessel functions. At the indirect detection, we can take the approximation in the low-velocity limit, i.e, $\langle\sigma v_{\text{lab}}\rangle_{XX\to SS}=(\sigma v_{\text{lab}})_{XX\to SS}$ with the replacement $s=4m_{X}^{2}$.

B.2 The annihilation process for $SS\to\text{SM SM}$

The diagrams for the hidden scalar $S$ annihilation into the SM particles are depicted in Fig. 10. The resulting annihilation cross section is given by

[TABLE]

where $N_{c}^{f}\equiv 3\,(1)$ for $f\equiv$ quarks (leptons), $g_{SSS}$ is shown in Eq. (90), and the remaining couplings are

[TABLE]

On can further apply Eq. (91) to obtain thermally averaged value of the annihilation cross section. This result is relevant to the chemical equilibrium between the hidden sector and thermal bath in the early Universe.

Appendix C The thermal average $\langle\sigma v\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle_{SS\to\text{SM}_{1}\text{SM}_{2}}$

In this paper, for the process $S(p_{S})\,S(p_{S^{\prime}})\to\text{SM}_{1}\,\text{SM}_{2}$ that follows the Maxwell-Boltzmann distribution, its thermal average, $\langle\sigma v_{\text{M\o l}}\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle$, at temperature $T_{i}$ in the cosmic comoving frame is defined by

[TABLE]

where the Møller velocity is given by

[TABLE]

For a typical case that $\sigma v_{\text{M\o l}}=(\sigma v)_{0}$ is constant, because

[TABLE]

we thus have $\langle\sigma v_{\text{M\o l}}\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle=T_{i}\langle\sigma v_{\text{M\o l}}\rangle=T_{i}(\sigma v)_{0}$. For a general case, we can first rewrite the momentum-space volume element to be

[TABLE]

with $\theta$ being the angle between ${\bf p}_{S}$ and ${\bf p}_{S^{\prime}}$. As seen from Eq. (97) that $v_{\text{M\o l}}E_{S}E_{S^{\prime}}$ is Lorentz invariant, we can relate the Møller velocities in two different frames with and without a prime to be

[TABLE]

where the last step uses the fact that $p_{S}\cdot p_{S^{\prime}}=E_{S}E_{S^{\prime}}(1-{\bf v}_{S}\cdot{\bf v}_{S^{\prime}})$ which is Lorentz invariant. Thus, we can use the Møller velocity given in the laboratory frame, which is equivalent to the rest frame of one of the incoming particles, to obtain the Møller velocity defined in the cosmic comoving frame; for simplicity, in the following we will use $v_{\text{lab}}$ for the former and $v_{\text{M\o l}}$ for the latter. The resulting relation is

[TABLE]

As such, we get

[TABLE]

where we have used the fact that $\int_{-1}^{+1}\cos\theta\,d\cos\theta=0$. In order to calculate $\langle\sigma v_{\text{lab}}\cdot\frac{{\bf p}_{S}^{2}}{3E_{S}}\rangle(T_{i})$, we further change integration variables from $E_{S},E_{S^{\prime}},\theta$ to $E_{-},E_{+},s$, given by

[TABLE]

with the integration region

[TABLE]

Using and new variables and calculating the thermal average in terms of modified Bessel functions of the second kind, we obtain

[TABLE]

where we have used the recursive relation in the last step,

[TABLE]

Appendix D $3\to 2$ annihilations

We consider a thermally averaged cannibal annihilation cross section for $a(p_{a})\,b(p_{b})\,c(p_{c})\to d(p_{d})\,e(p_{e})$, where all particles resides in a hidden sector and keep the same temperature, $T_{X}=T_{S}=T_{h}$, during the interaction. The generic form defined through this paper is given by

[TABLE]

where $m\equiv 2$ if the final state particles, $d$ and $e$, are the same particle species, otherwise $m\equiv 1$. Here, the sum for the amplitude squared, $|M|_{abc\to de}^{2}$, has been taken over all internal degrees of freedom of the initial and final states. In the nonrelativistic limit, $E_{a}\approx m_{a},E_{b}\approx m_{b},E_{c}\approx m_{c}$, the cross section is approximately given by Berlin:2016gtr

[TABLE]

where $\overline{|M|^{2}}_{abc\to de}$ is the amplitude squared but with the initial state spin-averaged. To calculate $3\to 2$ thermally averaged annihilation cross sections for the nonrelativistic hidden sector particles with a temperature below their masses, i.e., $T_{h}<m_{S,X}$, we take the low-velocity approximation, $\langle\sigma v^{2}\rangle\simeq\sigma v^{2}$ and neglect its subleading corrections of order $T_{h}/m_{X,S}$. We show the diagrams in Figs. 11, 12, 13, 14, and 15, and summarize all the relevant results as below,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix E The kinetic decoupling of the DM from the hidden scalar

In this Appendix, we discuss the epoch that the DM and hidden scalar starts to be kinetically decoupled from each other. For this period of time, which is well after $3\to 2$ annihilation decoupling, the DM temperature follows the Boltzmann equation,

[TABLE]

where $\delta_{H}^{X}\approx 1$ is the same form as $\delta_{H}$ (defined in Eq. (38)) but with $T_{S}$ replaced by $T_{X}$, and the $XS\to XS$ momentum relaxation rate is given by

[TABLE]

Based on the fact that $-t\ll 4|{\bf p}|^{2}\sim T\ll m_{X}^{2},m_{S}^{2}$, we neglect $t$ in the calculation of the amplitude squared to obtain the approximate form of $\gamma_{X}$. Thus, the elastic scattering amplitude squared and summed over all internal degrees of freedom of initial and final spin states is given by

[TABLE]

Here the elastic scattering process $XS\to XS$ contains the amplitudes with a hidden scalar $S$ mediated in the $t$-channel and with $X$ in the $s/u$-channel. The contribution is dominated by the former one. We define the DM kinetic decoupling temperature $T_{X}^{\rm kd}$ below which the kinetic energy injection rate transferred by the elastic scattering $XS\to XS$ and/or by the annihilation $SS\to XX$ to the DM falls below the diluting rate arising from the Hubble expansion. $T_{X}^{\rm kd}$ approximately satisfies

[TABLE]

where the approximation $\langle\sigma v\cdot\frac{{\bf p}_{X}^{2}}{3E_{X}}\rangle_{XX\to SS}(T_{S})\simeq T_{S}\,\langle\sigma v\rangle_{XX\to SS}(T_{S})$ with an error less than 10% is used for the present nonrelativistic s-wave annihilation, and the second term of the RHS arises from the contribution of $XX\to SS$ to \big{(}dn_{X}(T_{X})/dt+3Hn_{X}(T_{X})\big{)}T_{X} for $T<T_{f}$. Here, before kinetic decoupling, $T_{X}=T_{S}$ are functions of $T$ and can be determined from solving the Boltzmann moment equations. The DM temperature after kinetic decoupling satisfies $T_{X}(a)\simeq T_{X}^{\rm kd}\cdot(a_{X}^{\rm kd}/a)^{2}$, where $a$ is the cosmic scale factor and $a_{X}^{\rm kd}$ is its corresponding value at $T_{X}=T_{X}^{\rm kd}$. The $T_{X}/T$ evolution after kinetic decoupling is sketched as the red line in the left panel of Figs. 3, 4, and 5, where $T_{X}^{\rm kd}\equiv m_{X}/x_{X}^{\rm kd}$ is depicted as the red dot.

Similarly, if we consider the temperature below that the DM and hidden scalar are decoupled from each other, the following terms need to be included in the RHS of Eq. (53), which is the temperature evolution equation of $S$,

[TABLE]

where the $XS\to XS$ momentum relaxation rate is given by

[TABLE]

and the following approximations have been used in the last step:

[TABLE]