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

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.


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) [1][2][3][4][5][6][7][8][9][10][11], for which, however, the allowed WIMP dark matter models have been also severely constrained by the current null results of the direct detection [12][13][14][15] 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 [16][17][18][19][20][21][22][23][24][25].
The mechanism for the secluded WIMP dark matter was discussed by Pospelov, Ritz, and Voloshin [16]. 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 [26].
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 [17][18][19][20][21][22][23][24]. 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 [27]. Moreover, during the period of time in which the dark sector is out of thermal equilibrium with the bath, if the 3 → 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. [28][29][30][31]. Alternatively, 3 → 2 DM annihilation mechanism is also relevant to the strongly-interacting massive particles (or called SIMP) [32] and elastically decoupling relic (or called ELDER) [33,34] 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 ∼ 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 [17,19]. 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 [15], 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 µ 's, which get the mass from the vacuum expectation value (VEV) of the hidden complex scalar field Φ S due to the spontaneously dark gauge symmetry U X (1) breaking, where the Z 2 symmetry, X µ → −X µ and Φ S → Φ * S , is imposed to stabilize the dark matter [17]. The relevant kinetic Lagrangian (L kinetic ) and the scalar potential (L scalar ) are given by where Φ H = (H + , H 0 ) T is the SM Higgs doublet, X µν = ∂ µ X ν −∂ ν X µ , and the covariant derivative is defined as D µ Φ S = (∂ µ + ig dm Q Φ S X µ )Φ S . Here g dm is the gauge coupling and Q Φ S is the U dm (1) charge of Φ S . After spontaneous symmetry breaking, the Higgs fields develop non-zero VEV's, where the CP-odd states, σ h and σ s , respectively becomes the longitudinal components of the Z boson and X µ ; the dark matter thus obtain a mass, m X = g dm Q Φ S v S . In the present paper, we will simply use Q Φ S = 1.
The scalar fields (φ h , φ s ) can be expressed in terms of mass eigenstates of physical Higgses (h, S) and the mass squared matrix in the former basis can be parametrized in terms of masses of the latter and the mixing angle α, The S decay rate (total width), Γ S , and twice of the Hubble rate, 2H, as functions of m X /T , where T is the temperature of the bath. The horizontal lines from up to down with colors blue, red, brown, and green are the S decay rates corresponding to α = 1 × 10 −5 , 1 × 10 −6 , 5 × 10 −7 , and 1 × 10 −7 , respectively, while the black dotdashed line stands for the Hubble rate. In this plot, using m X =80 GeV, the solid and dashed lines correspond to m S = 0.8 m X and 0.99 m X , respectively.
Here and throughout the paper, we adopt the abbreviations: s α ≡ sin α and c α ≡ cos α.
Using v H 246 GeV and m h = 125.18 GeV [35], we will take m X , m S , g dm and α as the independent parameters in the following analysis.
The branching ratios of the hidden scalar, S, with a mass of m S 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 m X 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 with the total energy density being Here, g eff,SM and g 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. (31) and (58), and discussed in Sec. V(iv), if the nonrelativistic scalar S is kept in kinetic equilibrium with the thermal bath via its inverse decay SM SM → S, then this kinetic energy injection rate to S will be larger than the Hubble cooling rate, i.e., roughly Γ S 2H.
In the right panel of Fig. 1, the S particles with the width corresponding to α = 1 × 10 −5 or 1 × 10 −6 can be in thermal equilibrium with the bath when m X /T 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 α = 5 × 10 −7 or 1 × 10 −7 , because the ratio of the Hubble cooling rate to heating rate, ∼ 2Hn S (T S )T S / Γ S n eq S (T ) T , is much larger than 1 due to the fact that n S (T S )T S n eq S (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 eq S (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 the S lifetime, (2H) −1 Γ −1 S . As shown in the right panel of Fig. 1, the corresponding out-of-equilibrium temperature is about m X /T ∼ 6 or 30 for α = 5×10 −7 or 1 × 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 gammaray excess data. For the first set, when the hidden sector with a sizable mass gap undergoes the cannibal process, the down-scattering rate, XX → SS, can be significantly larger than the upscattering one, SS → 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 → 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 where µ XN = m X m N /(m X + m N ) is the reduced mass of the dark matter (X) and nucleon (N ), Left panel: The direct search limit on the (m X , m S ) parameter space with the requirement m X > m S which is the region on the right hand side of the magenta solid line. For a given value of s α (≡ sin α) ≤ 1/ √ 2, the region above the corresponding dashed (brown) curve is allowed by XENON1T [12]. Here we use g dm = 0.173. Right panel: Upper bounds of sin α and g dm from XENON1T and LZ projected sensitivity [15] denoted by the solid and dashed lines, respectively, where the red line is for m X = 80 GeV, m S = 0.8m X , and the blue for m X = 80 GeV, m S = 0.99m X . and f N = q N |qq|N m q /m N 0.3 [36]. The parameter space constrained by XENON1T [12] is shown in Fig. 2, where in the right panel the bound by the LUX-ZEPLIN (LZ) projected sensitivity [15] 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 α is above the corresponding dashed line, where we have limited α ≤ π/4 which is suitable for the case with a small α. As for π/4 ≤ α ≤ π/2, the bound is same as that with a mixing angle = π/2 − α, because sin 2α = sin 2(π/2 − α).
The paramter constraint from the invisible Higgs decay, which is less than 25% at the 95% CL [35], is much weaker than that from direct detection. Meanwhile, for the present case, h → S S is kinematically forbidden.

IV. THERMAL EVOLUTION OF THE NONRELATIVISTIC HIDDEN PARTICLES
The evolution of the phase space distribution f h (with h ≡ X or S) of the hidden sector particles in the homogeneous isotropic Friedmann-Robertson-Walker Universe is described by the Boltzmann equation, where H is the Hubble expansion parameter, p = E 2 h − m 2 h 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, with µ h the chemical potential of the particle species h. In the present case, we consider m X m S ∼ O(10 − 100 GeV) and the thermal evolution that the elastic scattering XS ↔ XS can keep the X and S particles in thermal equilibrium (T X = T S ) until kinetic decoupling temperature where a is the cosmic scale factor and a kd X is its corresponding value at T kd X . For a hidden sector particle, h 1 , the generic form of one of the collision terms described by the interaction "h 1 h 2 · · · b 1 b 2 · · · ↔ h 1 h 2 · · · b 1 b 2 · · · " can be written as where h ( ) i is the particle of the hidden sector with temperature T X for X or T S for S, b ( ) i 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 are the symmetric factors in the initial and final states, respectively, ∆ and ∆ 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 + SM ↔ S + SM, and elastic self-scattering SS ↔ SS), the moment result can be non-vanishing (see Eq. (33) for instance), and an additional factor "1/2" needs to be added in C[f h 1 ] to avoid double-counting. Taking SSS ↔ XX as an example, we have S y ≡ 3!, and S y ≡ 2!. Moreover, we have ∆ ≡ 3, ∆ ≡ 0 for considering the Boltzmann equation of the S particles, while ∆ ≡ 0, ∆ ≡ 2 for X. Here the 1 ± 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 ± f i 1, and thus approximate the distributions as Basically, this is a good approximation even for high or low temperature. The collision term can be then given as the following form, where 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.
A. 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. (9) with "1" and integrating over the momentum space, 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 where K i is the modified Bessel function of the second kind with x S ≡ m X /T S and x ≡ m X /T , n eq i is the equilibrium number density with vanishing chemical potential, Γ S is the total decay width of S into SM final states, and σv i and σv 2 i are respectively the thermally averaged cross sections for 2 → 2 and 3 → 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 eq S and, meanwhile, relevant to S → SM i SM i and SS → SM i SM i on the right hand side (RHS) of Eq. (17) are functions of the bath temperature, "T ", which is explicitly indicated, whereas the remaining ones appearing in Eqs. (16) and (17) 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, with x i ≡ m X /T i and x ≡ m X /T being the temperature variables of the hidden sector particles and thermal bath, respectively, the yields Y i ≡ 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 σv XX→SS is the leading approximation of σv XX→SS which is s-wave. In the following discussion, we will simply use y i (x) ≡ y i (x i ; x). We use x ≡ m X /T as the evolution variable and trade the time derivative in the Boltzmann equations to be where the relativistic degrees of freedom, h eff and g 1/2 * are defined via Thus, we can further rewrite the Boltzmann equations to be where M pl ≡ (8πG) −1/2 = 2.44 × 10 18 GeV is the reduced Planck mass, and the equilibrium value of y i is given by Again, it should be noted that in Eq. (24) the terms involving y eq S 's and relevant to S → SM i SM i and SS → SM i SM i are functions only of "x", as shown explicitly, while other quantities appearing Eqs. (23) and (24) 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 which can be determined by matching the present-day DM relic abundance Ω DM = (0.1198 ± 0.0026)/h 2 [35,37], where Y ∞ X is related to y ∞ X = y X (x → ∞) (see also Eq. (18)), s 0 = 2891 cm −3 is the visible entropy density today, ρ c = 3H 2 0 /(8πG) is the critical energy density, and h 0.678 is the Hubble constant H 0 of the present day in units of 100 km 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. (23) can be approximated Solving the equation, we get where we use the fact that the value of y X at x = x f is significantly larger than y ∞ X , and we can approximate T X ≈ T S in the calculation (see Figs. 3, 4, and 5 for the temperature dependence in the next section).

B. 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: We adopt the definition of the temperature, which is suitable not only for the nonrelativistic case at low temperatures, T < m S , but also for relativistic case at high temperature, T m S /0.01. The Boltzmann moment equation of the hidden scalar's temperature can be formed by multiplying Eq. (9) with p 2 S /(3E S ) and then integrating over the momentum space. Thus, we arrive at the form of the temperature evolution equation, where which is approximately to be "0" for nonrelativistic particles or "1" for ultra-relativistic ones. Here corresponding to a process i(initial state) ↔ f (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 = constant from Eqs. (16) and (17)  The collision term resulting from S(p S,1 )S(p S,2 ) ↔ SM 1 (p 1 ) SM 2 (p 2 ) is given by where m! = 2 for identical final state particles or 1 otherwise, and we have used the energy conservation E 1 + E 2 = E S,1 + E S,2 , and the relation, Here the thermal average, for which the detailed description is provided in Appendix C, is given where v 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 σv· GeV, and becomes unity in the limit x → ∞. For a typical case with σv being constant, the ratio is equal to one. See also discussions in Appendix C.

The collision term due to
The collision term resulting from S(p S ) ↔ SM 1 (p 1 ) SM 2 (p 2 ) is given by where m! = 2 for identical final state particles or 1 otherwise, Γ S→SM 1 SM 2 is S → SM 1 SM 2 decay width, and which approaches to n eq S (T ) in a nonrelativistic limit. The result of Eq. (37) is also correct if replacing Γ S→SM 1 SM 2 with the relevant three-(or more-) body decay mode. Here, the result for XSS ↔ 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 m X,S , its temperature will decrease logarithmically with the cosmic scale factor during cannibalization (see Eq. (65) for discussion). The collision term resulting from X(p X )S(p S,1 )S(p S,2 ) ↔ X(p X )S(p S ) is given by where we have used p µ for the X, and have approximated three initial hidden particles (XSS) that annihilate or are produced in the 4. The collision term due to elastic scattering: S + SM ↔ S + SM Here we consider the elastic scattering, S(p S ) + SM(k) ↔ S(p S ) + SM(k ), where p µ S ( ) = (E S ( ) , p S ( ) ), k ( )µ = (ω ( ) , k ( ) ) 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, where the hidden scalar scattering with all relativistic SM fermions is taken into account. Under the typical condition m S T ∼ ω, this term can further reduce to a semi-relativistic Fokker-Planck-type equation [38][39][40][41][42][43], where the momentum relaxation rate is given by for which the sum runs over all relevant relativistic SM species, and the differential elastic scattering cross section is with |M Sf →Sf | 2 the square of the scattering amplitude summed over initial and final spin states.
Note that in Eq. (43), we have followed the approach given in Ref. [40] to adopt the t-average matrix (8k 4 ) −1 0 −4k 2 dt(−t)dσ/dt due to the fact that the scattering amplitude squared in our case vanishes around t = 0 in the relativistic limit m f → 0; therefore, it is unsuitable to take the result at zero momentum transfer of the t-channel as done in Ref. [38].
Taking into account the elastic scattering Sf → 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 with N f c ≡ 3 (1) for quarks (leptons), and the couplings shown in Eqs. (B2), (B5), (B6), and (B7).
Averaging over t for the scattering amplitude squared, we get the momentum relaxation rate to be Here the transferred momentum t in the denominator of the amplitude squared is neglected in the calculation consistent with the requirement −t < 4|k| 2 ∼ T 2 m 2 S . Therefore, for the case with the resulting elastic decoupling temperature ∼ 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 el (≡ m X /T el ) should be less than what is shown (see (iii) in Sec. V for the definition of x el ). However, such overestimation does not affect our conclusions.

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, where 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 = T S , we need to further include the following two terms to the RHS of Eq. (47), where the first and second terms are the kinetic energy-transfer rates by elastic scattering (XS ↔ XS) and by annihilation (XX ↔ 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, the above Boltzmann moment equation for the temperature T S can be recasted into an alternative form that will used in the analysis: In Eq. (51), we do not distinguish T X from T S . An unphysical result may occur when y S (x S ) y eq S (x) and other interaction terms become much smaller than that involving 1 − δ Γ which originates from the interaction, S ↔ SM SM. The (1 − δ Γ ) term may result in y to be monotonically increasing with x, even though it is highly small and 1 − δ Γ , 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 ↔ XS and XX ↔ SS interactions. For instance, the terms, which are rewritten from Eq. (49) and not shown in Eq. (51), are given by with y T X ≡ T X /T . When the DM freezes out, the effect of the (1 − δ Γ ) 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 ≈ m S (i.e., y S (x S ) > y eq S (x) for x 1), the coefficient given in the square bracket of Eq. (52) is always much larger than that of the S ↔ SM SM interaction term involving 1 − δ Γ . In other words, the evolution of T S will closely follow T X for x 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 − δ) 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 − δ H )HT S,X ≈ 0 and ∝ a −2 , will approximately evolve with the same temperature. Based on the above reasons, we will neglect (1 − δ Γ ) term, i.e., simply take δ Γ = 1, for the case that the hidden sector is decoupled from the bath at T ≈ 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-ofequilibrium 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 ρ SM ∝ a −4 , the energy density of S evolves as ρ S ∝ a −3 . Thus, the value of ρ S /ρ SM is about 1/120 at temperature ∼ m S , but is increased to be ∼ (1/120) × 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 α < 5 × 10 −7 , and ∼ 10% for α = 1 × 10 −7 , where we have used the temperature result (T out de ) of the out-of-equilibrium S decay from the next section. Third, unlike the other ones, for the case with α = 1 × 10 −7 and m S = 0.8m X , because the down scattering results in y S > y X , the term involving 1 − δ Γ can be comparable with the coefficient given in the square bracket of Eq. (52) 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, ∝ 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 SM SM ↔ 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 → XX will be soon larger than the SM SM → S heating rate (described by the (1 − δ Γ ) 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,end = m X /x S,end ) below that the hidden scalar cannot be in thermal equilibrium with the bath, we will model a term as given in Eq. (53) 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 δ Γ = 1 to get the numerical result.
From that we then extract the DM kinetic decoupling temperature, T kd X , which will be defined and discussed in the next section. Adopting the obtained T kd X , we use the exact value of δ Γ and include the XX ↔ SS term (rewritten from the second term of Eq. (49)), in Eq. (51) to have the improved solution for x S,end , where x kd X ≡ m X /T kd X . Here the XS ↔ XS term is negligible for this case, and the relation T X T kd X · (T /T kd X ) 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 α, i.e. to be (1)  while for T below that denoted by the blue dot, S is completely thermally decoupled from the bath. When T is below that denoted by the purple "X" (corresponding to the vertical dotted (red) line in the right panel), the heating rate due to the inverse decay of S can be larger than that needed to keep the hidden sector in thermal equilibrium with the bath. The purple square (corresponding to the vertical dashed (red) line in the right panel) stands for the temperature below which the S undergoes an out-of-equilibrium decay. When T falls below that denoted by the green dot, the heating rate of the cannibal process in the hidden sector is less than H. The conservation of the comoving hidden sector entropy is described by the dotted curve (see (i) x f ≡ m X /T f is the usual freeze-out temperature variable. For T < T f , the comoving DM number density tend to be conserved. From Eq. (16), we can estimate the freeze-out temperature, below which the DM production rate from SS → XX is overtaken by the dilution rate, giving the relation σv XX→SS (n eq X (T S )) 2 (n eq S (T S )) 2 (n S (T S )) 2 3H(T )n X (T X ) , 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. For T > T f , the rate on the left hand side (LHS) of Eq. (54) is larger than the expansion rate, resulting in the detailed balance σv XX→SS (T X ) (n X (T X )) 2 = σv SS→XX (T S ) (n S (T S )) 2 .
This implies that so that X and S (with T X = T S ) have the same chemical potential, µ X = µ 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 (ii) x ann,2 ≡ m X /T ann,2 , plotted as the black dot in the left panel of Fig. 3, corresponds to the bath temperature T = T ann,2 set by where T S , the temperature of S, is a function of T determined by Eq. (47), and the second term of the RHS is described by Eq. (31). Here, considering the relevant channels SS ↔ ij SM i SM j , the value of (T S /a 3 ) d(n S a 3 )/dt = T S (dn S /dt + 3Hn S ) ∝ dy S /dt is about zero for x x 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 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 ij SM i SM j → SS, such that this kind of interactions cannot play the role to keep thermal equilibrium with the SM bath for T < T ann,2 . In Figs. 3(b-1), 4, and 5, x ann,2 is too small to be visible in plots.
(iii) x el ≡ m X /T 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 el , the heating rate of the hidden sector, which gains energy by the elastic scattering S SM → S SM, is larger in magnitude than the cooling rate due to the Hubble expansion, where (2 − δ H )γ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 µ S = 0) for T m S , the kinetic energy-transfer rate of S SM → S SM is proportional to the relativistic SM density which is not suppressed. Therefore, in most cases, we have x el > x ann,2 . Two remarks are in order. First, x el < x ann,2 may occur, if the annihilation cross section is largely enhanced by the resonant s-channel SM-Higgs exchange with s ∼ (2m S ) 2 ∼ m 2 h . Second, compared with the elastic scattering, because the 2 → 2 annihilation is much more sensitive to S-f -f coupling, which is proportional to s α , a small enough mixing angle α may result in the 2 → 2 annihilation decoupling to occur significantly before elastic decoupling.
(iv) x de ≡ m X /T 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 de denotes the temperature below which not only the requirements, n S (T S ) = n eq S (T ) and T S = T , need to be satisfied, but also the heating rate, generated from the inverse decay: SM SM → S, is larger in magnitude than the rate needed to keep the S particles in kinetic equilibrium with the bath during the cosmic expansion, 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 ↔ SM SM interactions, for which its contribution to (T S /a 3 ) d(n S a 3 )/dt = T S dn S (T S )/dt + 3Hn S (T S ) is relatively negligible when T ≈ T de for the case shown in Figs.

3(a) and 5(a). As such, if
T de > min(T el , T 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 eq S (T ) and T S = T ) until a later time, x S,end (≡ m X /T 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. (53) 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 S,end , denoted as the blue dot in Figs. 3(a-1) and 5(a-1), is determined by where the second and third terms of the left hand side are respectively the rates originating from SS → XX and XS → XS collision terms of the Boltzmann moment equation, while the second term of the RHS is the contribution from dn S (T S )/dt + 3Hn S (T S ) T S , resulting mainly from XX → SS due to the fact that n X (T X ) > n eq X (T X ) after the DM freezes out. The DM temperature after kinetic decoupling follows T X (a) T kd X · (a kd X /a) 2 , with a being the cosmic scale factor and a kd X being the corresponding value at T X = T kd X . For the cases shown in Figs. 3(a) and 5(a), this relation can be rewritten as T X = T kd X (T /T kd X ) 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 kd X ≡ m X /x kd X , depicted as the red dot, is the DM kinetic decoupling temperature, featuring T kd X ≤ T f . A detailed discussion for T kd X will be given in Appendix E. Numerically, we obtain x S,end ≈ x kd X , as seen from Figs. 3(a-1) and 5(a-1).
On the other hand, provided that the relation given in Eq. (58) is satisfied but with T f T de < min(T el , T ann,2 ), i.e., x f x de > max(x el , x ann,2 ), the hidden sector may be kinetically decoupled from the bath at T min(T el , T ann,2 ), such that at x x out de (see (vi) for the definition), the S first undergoes an out-of-equilibrium decay with a rate much larger than its inverse production rate SM SM → S due to the fact that n S n eq S , resulting in the RHS of Eq. (58) to be less than zero, It is interesting to note that, as the time evolves, the n S Boltzmann equation gives n S (T S ) → n eq S (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. (58) is possible to be met and x end can thus exist.
(v) x c ≡ 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 i,j,k≡S,X K 3→2 σv 2 3→2 n eq i (T S )n eq j (T S )n eq k (T S ) 2H(T )n eq where K 3→2 shown in Eq. (47) is the kinetic energy released in a relevant 3 → 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 µ X,S = 0: n X (T S ) → n eq X (T S ), n S (T S ) → n eq S (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 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 SM = (2π 2 /45)h eff SM (T )T 3 , to the hidden sector, s h = (2π 2 /45)h eff h (T S )T 3 S , is constant, where h eff SM and h eff h are the effectively relativistic degrees of freedom of the SM and hidden sector, respectively. Thus, we find where 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 (ρ h /T S )a 3 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, where s = s SM + s h , and a out,h (the cosmic scale factor) and x 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).
(vi) x out de ≡ m X /T out de , 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. (58) (see also Eq. (60)) is much larger than the term in the LHS of Eq. (58) in magnitude due to the fact that n S (T S ) n eq S (T ) at T = T out de . For this case, corresponding to a much smaller mixing angle as that given with α = 5 × 10 −7 or 1 × 10 −7 in this paper, because n S (T S ) T S n eq S (T ) T at a later time with x f x > max(x el , x ann,2 ), we thus approximately define x out de from Eq. (58) to satisfy where x out S,de is the value of x S , corresponding to x = x out de . When x = x out de , the hidden scalar thus starts to undergo out-of-equilibrium decay at the cosmological time, (2H) −1 ≈ S . 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 out de 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 out de which satisfies y S (x out de )/y S (1) e −1 [31] with the initial value y S (1) = y eq S (1). The value of x = x out de 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 → SS, can be significantly larger than the up-scattering rate, SS → XX, such that after a sufficient time at x = x out de 1, we have y eq X /y eq S 1 and y X y S . Therefore we set the effective initial S yield to be y in S (1) ≈ y X (1) + y S (1) ≈ y eq X (1) + y eq S (1) = κy eq S (1) with κ ≡ [3(m X /m S ) 2 K 2 (1)/K 2 (m S /m X ) + 1], and approximate Eq. (24) as Here the approximation in the last step of Eq. (68) is reasonable because g 1/2 * /h eff , which is 0.109 for T = m X /20 and 0.125 for T = m X /150, weakly depends on T f in the present study.
Solving this equation, we obtain and x out de 1 + 2(1 + ln κ)/C. In terms of the cosmic time variable of the radiation dominated epoch, the solution of the normalized yield can be rewritten as y S κy eq S (1)e −Γ S t . If κ = 1, we have t(≈ (2H) −1 ) = Γ −1 S at x = x out de , consistent with the that given in Eq. (66). For the case with m X ≈ m S (e.g. m S = 0.99m X ), summing Eqs. (23) and (24), we have Using the initial conditions: y X (1) = y eq X (1), y S (1) = y eq S (1), and y eq X (1)/y eq S (1) = g X /g S = 3, and approximating y X + y S ≈ 4y S , we have with x out de 1 + 8/C. The solution can be given by y S y eq S (1)e −Γ S t/4 . Note that a longerlived S will result in a larger x f , i.e., a larger y ∞ X , 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 ∼ m X,S , the ending value (x c ) of cannibalization depends on the magnitude of x out de since the number density of S is exponentially depleted during decay. In Figs. 4(b) and 5(d), we show the cases with x out de ∼ 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 c S · (a c /a) 2 even after thermal decoupling, where a is the cosmic scale factor and a c is its corresponding value at T c S = T S (T = T c ).

VI. DISCUSSIONS
Since we have considered the secluded vector dark matter model with DM mass ∼ O(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 where σv LV is the DM annihilation cross section in the low-velocity limit (consistent with T → 0), (dN 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 ∆Ω. 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 [44,45], where we adopt r s = 20 kpc, r = 8.5 kpc, γ = 1.2 and ρ = 0.4 GeV/cm 3 as the canonical inputs. Here "−γ" is the inner log slope of the halo density near the GC, and ρ is the local DM density at a distance of r 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 [46] (see also Eq. (10) in Ref. [24]), where we use PPPC4DMID result [47,48] 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 ≡ W or Z), the three-body decay channel S → V V * → V f 1f2 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.

18.
m X (GeV)  while that in pink color corresponding to ρ = 0.85 GeV/cm 3 and γ = 1.25. In Fig. 7, the GC fit together with other constraints is redrawn on the (m X , g dm ) plane, where a larger g dm is needed to account for the data for the nearly degenerate case because σv LV vanishes in the limit m S → m X .
We remark that a newer Pass 8 Fermi data set was analyzed in Ref. [49], 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 [11,49].
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 α 2 × 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 α 6×10 −7 (or α 1×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 [50,51].
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 [50]. See the detailed description in Ref. [31] 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 [26] where σv CMB σv LV for s-wave DM annihilation, and the efficiency factor is Here, we use f γ,e − eff (E) curve results suited for the "3 keV" baseline prescription shown in Ref. [52]. Moreover, as the previous study for the GC excess, (dN γ,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 , α and g dm . As shown in Figs. 3, 4, and 5, having the chosen masses for the DM and hidden scalar, and giving the magnitude of α, we can fine-tune the value of g dm in the numerical analysis of Boltzmann equations to obtain y ∞ X , which matches the observed relic abundance determined by Eq. (27). Note that, in Eq. (27), σv (0) XX→SS is a function of g dm , and its α-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 (α, y ∞ X ) plane in Fig. 8. Since the S decay width Γ S is a function of α, we label its corresponding values on the top of the plots. On the other hand, for a obtained y ∞ X , we can get x f from the relation given by Eq. (54), and further have the corresponding value σv LV (≡ σv   8. Contour of having the correct relic abundance (black curve) on the (α − y ∞ X ) plane. On the top, the value of the S decay width Γ S , which is a function of α, is labeled. Here for a given α, we fine-tune the value of g dm from the numerical analysis of Boltzmann equations to obtain y ∞ X , which matches the observed relic abundance determined by Eq. (27). From each set of allowed parameters, we get x f from the relation given by Eq. (54), and further have the corresponding value σv LV (≡ σv (0) XX→SS ) from Eq. (27) or from the value of g dm . The dependence of σv LV on α is negligible. The 95% C.L. upper limit from CMB is denoted as the dot-dashed (brown) line. The 95% C.L. upper limit and project sensitivity for dSphs observations are shown as horizontal solid (red) and long-dashed (red) lines, respectively. The RHS of the dotted (blue) curve corresponds to projected reach by LZ. When α is larger than that denoted by the vertical dashed (magenta) line, the dark sector is well in the chemical and thermal equilibrium with the bath before freeze out. The region between the two horizontal short-dashed (green) lines provides a good fit to the GC gamma-ray excess. or (m X , m S ) = (80 GeV, 79.2 GeV) with α 1 × 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 α. It has been shown and discussed in Refs. [53][54][55][56] that the late-time entropy production by the massive particle decay could induce cosmological effects at the time ∼ Γ −1 S ∼ 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 ρ S ∝ a −3 , while the SM radiation energy density ρ SM scales as a −4 . As a result, ρ S /ρ SM , scaling linearly with a, is about 1/120 at T ∼ m S , but becomes ∼ (1/120) × 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  Fig. 4 in Ref. [56], where neutrino self-interaction and oscillation are included). If T RH 7 MeV, the effective number of neutrino species N eff becomes smaller than three (see Fig. 4 in Ref. [54] or Fig. 1 in Ref. [56] 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 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 ≡ 4n He /n b , and the deuterium ratio D/H are thus enhanced.
Following Refs. [54,56], we define the reheating temperature of the SM bath to be Γ S = 3H(T RH ). Using the approximation, T RH can be related to the decay width of S as where we have used g eff = 43/4. Here, we quote a lower bound T RH 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. [56]), which is suitable for our model. Further considering the neutrino self-interaction and oscillation, the same reheating bound is also required by N eff = 3.15±0.23 from Planck report [26] (see Fig. 1 in Ref. [56]). Therefore, from Eq. (78), we can obtain the BBN constraint on the S width to be Γ −1 0.03 sec. As such, we have α 1. For the case satisfying T de ≥ min(T el , T ann,2 ), the kinetic decoupling of elastic scattering S SM ↔ S SM and/or annihilation SS ↔ SM SM occurs only when the bath temperature is below T de at which the heating rate of the hidden sector generated from the inverse decay SM SM → 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 × 10 −6 α 4 × 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 ∼ 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-ofequilibrium 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 α 6 × 10 −7 (or α 1 × 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 "longlived" 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. via h is negligible and does not shown. The resulting cross section is given by where s is the center-of-mass energy squared, and Using the above result, the thermally averaged annihilation cross section for T 3m X can be obtained by calculating [60], where the Møller velocity is given by For a typical case that σv Møl = (σv) 0 is constant, because we thus have σv Møl · p 2 S 3E S = T i σv Møl = T i (σv) 0 . For a general case, we can first rewrite the momentum-space volume element to be with θ being the angle between p S and p S . As seen from Eq. (C2) that v Møl E S E S is Lorentz invariant, we can relate the Møller velocities in two different frames with and without a prime to be σv Møl = σv Møl where the last step uses the fact that p S · p S = E S E S (1 − v S · v S ) 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 lab for the former and v Møl for the latter. The resulting relation is As such, we get where we have used the fact that +1 −1 cos θ d cos θ = 0. In order to calculate σv lab · p 2 S 3E S (T i ), we further change integration variables from E S , E S , θ to E − , E + , s, given by with the integration region Using and new variables and calculating the thermal average in terms of modified Bessel functions of the second kind, we obtain where we have used the recursive relation in the last step, K n+1 (x) = K n−1 (x) + 2n x K n (x) . where δ X H ≈ 1 is the same form as δ H (defined in Eq. (32)) but with T S replaced by T X , and the XS → XS momentum relaxation rate is given by Based on the fact that −t 4|p| 2 ∼ T m 2 X , m 2 S , we neglect t in the calculation of the amplitude squared to obtain the approximate form of γ X . Thus, the elastic scattering amplitude squared and summed over all internal degrees of freedom of initial and final spin states is given by Here the elastic scattering process XS → 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 kd X below which the kinetic energy injection rate transferred by the elastic scattering XS → XS and/or by the annihilation