When Freeze-out occurs due to a non-Boltzmann suppression: A study of degenerate dark sector

Exponential suppression or commonly known as the Boltzmann suppression in the number density of dark matter is the key ingredient for creating chemical imbalance prior to the usual thermal freeze-out. A degenerate/quasi-degenerate dark sector can experience a different exponential suppression in the number density analogous to the radioactive decay law leading to a delayed freeze-out mechanism of dark matter known as the co-decaying dark matter. In this work, we study the dynamics of a multicomponent dark matter from thermally decoupled degenerate dark sector in a hidden U$(1)_{X}$ extension of the Standard Model. We compute the relic density of dark matter frozen-out through the co-decaying mechanism by solving four coupled Boltzmann equations. We demonstrate how temperature $T^\prime $ of the dark sector changes due to all types of $3\rightarrow 2$ and $2\rightarrow 2$ interactions along with the eternal expansion of the Universe. We find that $3\rightarrow 2$ interactions enhance $T^\prime$ by producing energetic particles in the dark sector while the excess heat is transferred by $2\rightarrow 2$ interactions to the entire dark sector. As the direct detection is possible only through the feeble portal couplings, we investigate the neutrino and $\gamma$-ray signals from dark matter annihilation via one step cascade processes and compare our results with the measured fluxes of atmospheric neutrinos by Super-Kamiokande and diffuse $\gamma$-rays by Fermi-LAT, EGRET, INTEGRAL collaborations. We find that the present scenario easily evades all the existing bounds from atmospheric neutrino and diffuse $\gamma$-ray observations for degenerate dark sector. However, the constraints are significant for quasi degenerate scenario.


Introduction
The gravitational effects established the fact that almost 25% energy budget of our Universe is made of an unseen matter commonly known as the dark matter (DM). In particular, the satellite borne experiments like WMAP [1] and Planck [2] while measuring the temperature anisotropy in the Cosmic Microwave Background radiation (CMB) have established the present value of dark matter relic density Ω DM h 2 = 0.120 ± 0.001 [2]. Other indirect evidences like the rotation curves of spiral galaxies [3], the gravitational lensing of distant objects [4], the bullet cluster observation [5] etc. also strongly indicate the presence of more matter than only the visible matter. In spite of having clear evidences, the nature and properties of dark matter is an open problem to date. There exist many particle physics models in the literature, which have one or more candidates of dark matter in their particle spectrum. Most of these beyond Standard Model (BSM) theories have focused on weakly interacting massive particle (WIMP) dark matter [6][7][8][9], a popular and well studied class of dark matter candidates. In the WIMP paradigm, it is assumed that dark matter was in thermal equilibrium (both chemical and kinetic equilibrium) with the Standard Model (SM) particles and the chemical equilibrium was lost as dark matter became non-relativistic. When expansion rate of the Universe becomes dominant over the DM DM → SM SM interaction rate, the dark matter number density gets frozen-out and changes only by expansion of the Universe afterwards. The most important property of a WIMP dark matter is that it has sizeable interaction strength (in the weak scale) with the SM particles that predicts observable number of events in the experiments over the possible backgrounds. However, the parameter space of this theoretically well motivated scenario are now acutely constrained due to non-observations of any signature at various direct detection experiments over the last two decades [10][11][12][13][14][15][16][17][18]. Several near future experiments like DARWIN [19], XENONnT [20], LUX-ZEPLIN [21] will have the required sensitivities to probe the remaining WIMP parameter space above the neutrino floor [22], a region dominated by coherent neutrino-nucleon scatterings.
In the light of these observations, many well motivated proposals have been suggested which can explain the null results of direct detection experiments and also remain compatible with the current Planck result. The "secluded sector dark matter" [23][24][25][26][27][28][29][30][31][32] is one of such alternatives. In the secluded sector scenario, dark matter is connected with SM through a metastable mediator and the coupling between the mediator and SM is small to evade the constraints from the direct searches and also from collider experiments. However, the coupling is sufficient enough to establish the kinetic equilibrium with the SM sector. Since the metastable mediator can decay into the SM species, it should decay well before the onset of big bang nucleosynthesis (BBN) to ensure that the observables (e.g. effective number of neutrino species, abundance of deuterium and helium) during the BBN era remains unaltered. In this scenario, if the dark matter mass is greater than the mediator mass, annihilation of dark matter into the metastable mediators is possible and that determines the dark matter relic abundance at the present era while the dark matter is still a thermal relic. This type of dark matter though secluded from the SM sector due to the tiny portal coupling with the SM, can be detected via indirect searches. Most of the work on secluded sector dark matter have been done on the basis of the assumption that the dark sector is in kinetic equilibrium with the SM bath. In general, this condition can be relaxed by assuming that the dark sector is decoupled from the SM sector in the early Universe while it was relativistic. Since, the dark sector is not in kinetic equilibrium with the SM, the entropy is conserved separately in each sector and the temperature evolution in the dark sector is in general different from the SM bath. Now, if dark sector has sufficiently strong interactions among its species then beyond 2 → 2 scatterings e.g. inelastic scatterings like 3 → 2, 4 → 2 processes etc. can be active. As a result, the dark sector enters into the cannibal phase [33][34][35][36] and during this phase, the dynamics of temperature evolution of the dark sector is drastically different from the evolution of the SM temperature.
In this work, our principal objective is to carry out a detailed phenomenological analysis of the secluded sector dark matter involving both 2 → 2 and 3 → 2 scatterings. In order to implement this, we have considered a dark sector which has a U(1) X gauge invariance. In this dark sector, two left chiral Weyl fermions ξ 1 L and ξ 2 L singlet under SU(3) c ⊗ SU(2) L ⊗ U(1) Y gauge group and having U(1) X charges +1 and −1 respectively are added. Moreover, we have included a complex scalar η which is also singlet under the SM gauge group but has a nonzero U(1) X charge +2 to break this additional U(1) symmetry spontaneously. After symmetry breaking in the dark sector, we have two Majorana fermions χ 1 and χ 2 , a dark scalar h d and a neutral gauge boson Z . Among them, the Majorana fermions are absolutely stable and will play the role of dark matter. However, other two species h d and Z are feebly coupled to the SM sector through the kinetic mixing between two U(1) groups (one of them is U(1) Y of the SM) and the scalar mixing α between the SM Higgs boson and η. However, these portal couplings are such that two sectors are not in kinetic equilibrium and Z , h d can decay out of equilibrium to the SM particles. Since we do not have any prior knowledge about the species of dark sector (except dark matter) and the scale of interactions among the members of dark sector, it is natural to assume all the dark sector couplings are of same order. This simplifying assumption results in a degenerate/quasi-degenerate dark sector which can drastically modify the entire freeze-out mechanism of the dark matter. It has been observed that due to degeneracy among the species there is no Boltzmann suppression in number densities of nonrelativistic dark species which is a key ingredient for a dark matter proceeding towards freeze-out. Instead, here each species develops a nonzero chemical potential which discards the Boltzmann suppression. However, the necessary suppression in number density arises when either Z or h d (or both) starts to decay into the SM particles and thereby introduces a new kind of exponential suppression e − Γ 4H to the number densities of dark species, where Γ is the decay width and H is the Hubble parameter. This results in the violation of chemical equilibrium in the dark sector and consequently triggers the freeze-out process of dark matter. As the freeze-out is initiated only after a sufficient time interval from the beginning of decay of species other than dark matter, which are feebly coupled to the SM fields, the process naturally leads to a delayed freeze-out of dark matter.
This novel mechanism was proposed in [37] under the name of "co-decaying dark matter".
Thereafter a few studies have been performed on this topic concentrating on astrophysical implications like small scale structure formation [38] and formation of primordial black hole [39]. In this co-decaying framework, we have considered a temperature T , different from the SM temperature T , in the dark sector, which is a natural choice as the two sectors are not in thermal contact. We have done a detailed analysis of temperature evolution of the dark sector taking into account all types of 2 → 2 and 3 → 2 scattering processes among the dark species. In this process, we have numerically solved four coupled Boltzmann equations, among which three equations are for evolution of number densities of χ 1 + χ 2 , Z and h d while the last one is for T and have computed the relic density of dark matter frozen-out via the co-decaying mechanism. We have found that 3 → 2 processes have important implications in the evolution of T as they increase the temperature by introducing highly energetic particles in the dark sector. The excess heat thus generated via 3 → 2 scatterings spreads across all the species through the 2 → 2 scatterings so that the entire dark sector has a common temperature. Here, we have chosen the portal couplings and α which connect the dark sector with the visible sector in such a way that the observations during the BBN era remain unaltered. Moreover, in this framework, we have also studied extensively the prospects of indirect signatures of our dark matter candidates (χ 1 and χ 2 ) through neutrinos and γ-rays. For that we have considered the most dominant processes which are one step cascade processes and have the following generic structure like χ i χ i → A(A → X Y ). The intermediate states are Z for neutrinos and h d for γ-rays. This type of cascade process generally produces box-shaped spectrum. However, as we are considering degenerate/quasi degenerate dark species, the resulting spectrum takes a line like shape with energy m χ /2 of the outgoings particles. Moreover, γ ray flux from the Final State Radiation (FSR) and the Inverse Compton Scattering (ICS) has also been discussed.
We have compared our results with the observed neutrino flux by the Super-Kamiokande [40] detector while the gamma-ray flux has been compared with diffuse γ-ray background measured by the Fermi-LAT [41], the EGRET [42], and INTEGRAL [43] collaborations. We have found that in both cases the present scenario easily satisfies all the existing bounds arising from diffuse background γ-rays and atmospheric neutrinos for degenerate dark sector. However for quasi degenerate dark sector, the parameter space is constrained from CMB observations and measurement of positron flux at AMS-02 experiment. The rest of the paper has been organised in the following manner. In Section 2 we have discussed the present model describing our dark sector. A detailed discussion on the dynamics of the dark sector following co-decaying mechanism is given in the first part of Section 3. The necessary Boltzmann equations and related discussions are given in Subsection 3.1. Various constraints on the portal couplings have been shown in Subsection 3.2. Numerical results that we have obtained by solving four coupled Boltzmann equations are presented in great detail in Subsection 3.3. The indirect detection constraints on this framework from astrophysical neutrinos and diffuse background γ-rays are discussed in Section 4. Finally, we summarise in Section 5. A comprehensive discussions on the Boltzmann equations, all 2→2 and 3 → 2 inelastic scatterings, necessary vertex factors, an approximate analytical form of the dark matter relic density and a short note on the thermal averaged cross section in the degenerate limit have been given in Appendices A-F.

Model
In this section, we discuss an anomaly free U(1) extension of the SM which has two Majorana dark matter candidates. In order to proceed further, let us discuss the model in detail. We extend the fermionic sector of the SM by two left-chiral Weyl fermions (ξ 1 L , ξ 2 L ) which are singlet under the SM gauge group but are charged under the new U(1) X gauge group. Besides the two fermions, we need at least one complex scalar (η) with nonzero U(1) X charge for breaking of new gauge symmetry spontaneously and this results in a massive neutral gauge boson in the particle spectrum. The complete field contents and their individual charges 1 Table 1. The U(1) X charges of ξ 1 L and ξ 2 L are equal to +1 nd −1 respectively and this is an economic choice leading to the cancellations of both [U(1) X ] 3 and [Gravity] 2 U(1) X anomalies. This can easily be understood as follows. Suppose, q 1 and q 2 are the charges of ξ 1 L and ξ 2 L respectively under U(1) X . Then the anomaly cancellation conditions for axial vector anomaly and mixed gauge-gravitational anomaly are given by The general solution (real) for the above set of equations is q 1 = −q 2 and we choose q 2 = −1 and q 1 = −q 2 = 1. On the other hand, the gauge invariance of the Majorana mass terms for both ξ 1 L and ξ 2 L demands the charge of η to be -2. Now, we shall write the Lagrangian of our model which is invariant under the entire symmetry groups i.e. SU(3) c ⊗ SU(2) L ⊗ U(1) Y ⊗U(1) X . The total Lagrangian L, composed of the SM Lagrangian (L SM ) and the additional part involving newly added gauge, Yukawa and scalar sectors respectively, is given by where, Where, Φ is the SM Higgs doublet, ξ i L = C ξ i L T , (i = 1, 2) and C is the charge conjugation operator. The field strength tensor of the extra U(1) X gauge symmetry is X µν = ∂ µ X ν − ∂ ν X µ while the corresponding tensor for the U(1) part of the SM is B µν . The Covariant derivatives in the above Lagrangian, needed to restore the gauge invariance under the local U(1) X symmetry, have the usual definition as D µ = ∂ µ + ig X q X X µ where, q X is the U(1) X charge of that particular field on which D µ is acting and g X is the new gauge coupling. In Eq. 2.3, the first two terms are responsible for the Majorana masses for both ξ 1 L and ξ 2 L respectively when η gets a nonzero VEV. There could also be a Dirac mass term like m 12 ξ 1 L c ξ 2 L involving both ξ 1 L and ξ 2 L , which conserves U(1) X charge. However, since the Dirac mass term triggers mixing between ξ 1 L and ξ 2 L , for simplicity, one can easily avoid such term by assuming one of the two chiral fermions as Z 2 odd. The electroweak symmetry and U(1) X symmetry are broken when both Φ and η acquire VEVs v and v X respectively and in the unitary gauge they can be represented as The second term in Eq. 2.1 with coefficient /2, represents the kinetic mixing between U(1) Y and U(1) X and is not forbidden by any symmetry of the Lagrangian. Therefore, before proceeding further, we shall write Eq. 2.1 in the canonical form. To do this, we perform a basis transformation from the "un-hatted" basis to the "hatted" basis by a non-orthogonal transformation, which is given by Assuming << 1 (supported by various experimental observations [44][45][46][47]), one can write B µ ≈ B µ − X µ and X µ ≈X µ . In spite of restoring the canonical form in the Lagrangian forB µ andX µ , the mixing among the neutral gauge bosons W 3 µ ,B µ andX µ will again reappear when we substitute the transformation relation for B µ in the covariant derivative of Φ. After spontaneous symmetry breaking these mixing terms are solely responsible for all four off-diagonal elements (except the mixing between W 3 µ andB µ ) of the 3 × 3 neutral gauge boson mass matrix which has the following form in the basis B µ W 3 µX µ : Here, g 1 and g 2 are the gauge couplings of U(1) Y and SU(2) L respectively while as mentioned earlier is the co-efficient of the kinetic mixing term in Eq. 2.1. To diagonalize M 2 GB , we first rotate the basis vector B µ W 3 µX µ T at the Weinberg angle (θ W ) in theB µ − W 3 µ plane. After this rotation, both W 3 µ andB µ have changed to Z µ = cosθ W W 3 µ − sinθ WBµ and the orthogonal state A µ = sinθ W W 3 µ + cosθ WBµ with m 2 A = 0 while the third stateX µ remains unaffected. Finally, we apply another rotation in the Z µ −X µ plane at an angle θ 1 and get the diagonal matrix in the basis (A µ Z µ Z µ ): where, masses of the new physical states are where, Z µ is the usual SM Z-boson having mass M Z = 91.1876 ± 0.0021 GeV [48] while the new gauge boson corresponding to U (1) X is represented by Z µ . Unlike the SM, here we have one more mixing angle θ 1 along with the Weinberg angle θ W . The two mixing angles θ W and θ 1 are given by For completeness, we have written below the transformation relation between the physical basis and the gauge basis in a matrix form as Let us now look at the fermionic sector. After the breaking of U(1) X one can easily rewrite Eq. 2.3 in terms of two Majorana fermions χ 1 and χ 2 as where, each Majorana field is defined as χ i = ξ iL + ξ c iL , and the corresponding Majorana mass In the scalar sector, after the EWSB and U(1) X breaking, the masses of two physical CP-even scalars, which are the admixtures of H and ζ, are generated. On the other hand, both the CP-odd scalars become would-be Goldstone bosons corresponding to Z and Z bosons respectively. The CP-even scalar mass matrix in the basis (H ζ), using Eq. 2.4 and Eq. 2.5, is given by The scalar mass matrix M 2 scalar can be diagonalised using an orthogonal transformation The two physical CP-even scalars h and h d have masses m h and m h d respectively and are given by The new scalar sector mixing angle α is given by (2.14) Here, h is our SM-like Higgs boson having mass m h = 125.10 ± 0.14 GeV [48][49][50]. Using the transformation relations in Eq. 2.9 and Eq. 2.12 along with the definition of two Majorana fermions

Dynamics of the dark sector
As discussed in the previous section, the present model has a very rich dark sector having a scalar h d , two Majorana fermions χ 1 , χ 2 and a dark gauge boson Z . Among these particles, both χ 1 and χ 2 are absolutely stable and thus they are automatic choice for our dark matter candidates. The remaining two species h d and Z have tiny couplings with the SM particles through mixings (parametrised by the mixing angles θ 1 and α respectively) which are assumed to be extremely small θ 1 , α << 1. In this situation, if we consider a "democratic choice" that all the dark sector couplings are of the same order 2 , a degenerate dark sector can easily be achieved. However, such a simplifying assumption has deep impact on the cosmological evolution of dark matter and associated particles. In this case, instead of the usual freeze-out process with the Boltzmann suppressed dark matter number density in the non-relativistic regime, we are encountered with a different type of exponential suppression ∝ e −Γt in the number density followed by late freeze-out of dark matter. Here, Γ is the decay width of associated dark sector particles into the SM particles. This is known as co-decaying dark matter scenario as first proposed in [37]. Before going to the detailed phenomenological analysis of the dark sector, we would first like to discuss on the co-decaying dark matter scenario briefly.
Here we need at least two species in the dark sector which are either degenerate or quasidegenerate in mass and the dark sector has decoupled from the SM while it was relativistic. Now, let us consider two species A and B in the dark sector and furthermore let us assume B is feebly connected to the SM, so that later on it can slowly decay (Γ B H) into the SM particles. In our model, both χ 1 and χ 2 are playing the role of A while h d and Z are mimicking the species B. Initially, after decoupling from the SM sector at a temperature T d , the dark sector particles A and B remain in thermal equilibrium among themselves due to the scattering AA BB and strong individual self-interactions respectively. Moreover, they maintain the thermal equilibrium even after the dark sector temperature T drops below m, the common mass scale of A and B. This is particularly due to the reason that as A and B are degenerate in mass, the scattering AA BB remains kinematically viable even after T becomes less than m. Besides, there can be 3 → 2 interactions as well involving both A and B species, which are also in chemical equilibrium 3 . During this phase, due to these cannibalism, the chemical potential of each species in the dark sector is zero. However, 3 → 2 processes decouple before 2 → 2 scatterings as interaction rates of the former are quadratically suppressed by both number density and velocity of the initial state particles. Let us assume that the temperature at which all 3 → 2 interactions are frozen-out is T c and the corresponding SM temperature is T c . Therefore for T < T c , both A and B develop a nonzero chemical potential which helps them to get rid of the Boltzmann suppression in the non-relativistic regime. Now, the chemical equilibrium of AA BB interactions will be lost if the species B starts decaying into the SM particles and thus reducing the number density of B by an exponential factor e −Γ B t , where Γ B is the decay width of B. Once the chemical equilibrium is lost (AA / BB) at T Γ , the freeze-out of dark matter (A) occurs when the interaction rate of AA → BB at a particular temperature T f goes below the corresponding expansion rate of the Universe controlled by the Hubble parameter H. Therefore, the final abundance of A in this co-decaying scenario depends on both the annihilation cross section of AA → BB as well as the decay width of B into the SM particles. Moreover, as we will see later the freeze-out of A, for certain values of model parameters, occurs well after it enters into the non-relativistic regime (i.e. x f = m/T >> 20) if B has longer lifetime. This requires large annihilation cross section of A compared to that of the usual freeze-out cases (∼ 3 × 10 −26 cm 3 /s) and this opens the possibility of indirect detection prospects of A in the present and upcoming neutrino and γ-ray detectors, which has been discussed elaborately in Section 4. 2 In an unknown dark sector, the assumption of equal order for all couplings is a realistic one. Typically here we need yi/ √ 2 2 gX √ 2λX . 3 Although, 2 → 3 scatterings are kinematically suppressed for T < m, the interaction rate for forward and backward processes still could-be of same order. Now, we will try to understand the dynamics of the dark sector in more details by computing certain thermodynamic quantities analytically. We will compare our analytical predictions with those obtained from numerical solutions of the Boltzmann equations in the next section. As we mentioned earlier, the dark sector has decoupled from the visible sector when they are relativistic. As a result, the entropies of the dark sector and the visible sector are separately conserved until either Z or h d or both start to decay out-of equilibrium into the SM species. Suppose, t d , t c and t Γ are the time scales for decoupling of the dark sector from the SM, freeze-out of 3 → 2 processes and beginning of decay of metastable mediators (h d , Z ). The corresponding hidden sector temperatures 4 are already defined above as T d , T c and T Γ respectively with T d = T d . Before the freeze-out of the 3 → 2 processes i.e. during the cannibal phase with temperature T > T c , chemical potential of each dark sector species is zero. Now, using the entropy conservation and the conservation of total number of dark sector species in a co-moving volume for a temperature T lying in the range T c > T > T Γ , one can have, where s is the total entropy density of the dark sector and n = n χ 1 + n χ 2 + n Z + n h d is the total number density of all dark species. Therefore, from Eq. 3.1, we can write Here a(T ) is the cosmic scale factor at a temperature T . Now, using the second law of thermodynamics for a non-relativistic species, we can easily find the ratio of entropy density to number density at an arbitrary temperature T as where, ρ (T ) and P (T ) are total energy density and total pressure of all the dark sector particles at T . Moreover, in the above we have considered that chemical potential for each species is same 5 and equal to µ (T ). Substituting the expressions of ρ (T ), P (T ) for a non-relativistic dark sector in Eq. 3.3, we get Here m is the common mass scale between all degenerate dark sector species. Finally, substituting Eq. 3.4 in Eq. 3.2 and considering µ (T c ) = 0, we get the following expression of chemical potential for each species at any arbitrary temperature T (T c > T > T Γ ), Eq. 3.5 clearly shows that each species in the dark sector develops a temperature dependent chemical potential which pauses the Boltzmann suppression in number density for a non-relativistic species between T c and T Γ . Now we want to calculate how T changes with respect to the SM temperature T between T d and T Γ . The nature of T between T d and T Γ is vastly different in two regimes separated by T c (this is the temperature where all 3 → 2 interactions are frozen-out). In these two regimes, the chemical potential µ (T ) for each species behaves differently with T . In the first domain bounded between T d and T c , µ = 0 while the temperature dependence in the second domain is given in Eq. 3.5. Let us first consider T lying between T d and T c . At the time of decoupling of the dark sector from the SM, both the sectors are relativistic and they have a common temperature T d . Now, using the entropy conservation separately for both the sectors at a temperature T < T d , we have where g s (T ) is the number of relativistic degrees of freedom in the visible sector contributing to the entropy density at T while the corresponding quantity for the dark sector has been denoted by g s (T ). The SM at T is radiation dominated however, the dark sector is non-relativistic at T (T < m). Using s (T ) given in Eq. 3.3 for µ (T ) = 0, we get where, x = m/T and x = m/T . After a few mathematical simplifications we obtain the following expression of x valid for x 1, Therefore, due to the logarithmic dependence, the hidden sector temperature before T c varies slowly with respect to the SM temperature. This can be understood in the following way. In this regime, due to the expansion of the Universe there is a rapid decrease of T with respect to T as the former is proportional to T 2 . However, the decrement in T due to expansion has been compensated partially by heat generation in the dark sector via all the 3 → 2 scatterings active at temperature T > T c . The resultant effect is a logarithmic dependence on T . Finally, we will show the behaviour of T with respect to T in the other domain where T lies between T c and T Γ . Following the similar procedure as we have done in the previous case, i.e. using the entropy conservation at T c and T respectively for both the sectors and substituting the expression of µ (T ) as given in Eq. 3.5, we find the following expression of x as given bellow As there is no other source active for heating the dark sector, the variation of x follows the usual redshift of temperature of a non-relativistic species due to expansion of the Universe, which scales as a −2 (T ).

The Boltzmann equations
Now, we will formulate the system of Boltzmann equations for our dark sector described in Section 2. The dark sector particles Z and h d have important roles on the thermal evolution of the dark matter candidates χ 1 and χ 2 as the freeze-out of the latter is initiated only after the beginning of decay of the former into the SM particles. Moreover, since our dark sector has decoupled from the SM at a temperature T d , it is not in kinetic equilibrium with the SM thereafter and its temperature T , beyond T d , is different from that of the SM denoted by T . Thus, we have four coupled Boltzmann equations, where first three equations are corresponding to the evolution of number densities of χ 1 + χ 2 , Z and h d respectively while the last equation describes the evolution of T for T < T d . The four coupled Boltzmann equations which we need to solve to extract the physics of the dark sector are given below, where, in the last equation 6 we have substituted the right hand side of Eq. 3.10. Here n χ = n χ 1 +n χ 2 is the total dark matter number density while that of Z and h d are denoted by n Z , n h d respectively.
The corresponding equilibrium number density of a species i at temperature T is shown by n eq i (T ) and H is the Hubble parameter. The quantity σv T a b→c d is the thermal averaged annihilation cross section for the process a b → c d and it depends on the temperature of the species appearing in the initial state which in our case is T . Here v is the magnitude of the relative velocity between 6 In non-relativistic limit, this equation can also be derived from the more commonly used equation of energy density (see Eq. 12 of [37]) by putting ρχ mχnχ + 3nχT 2 .
a and b. The 2 → 2 annihilation processes in the dark sector are The Feynman diagrams and the corresponding cross sections for these processes are given in Fig. 1 and Appendix E respectively. Here We have taken contributions from both the dark matter candidates for each annihilation channels and these have been incorporated by the quantity σv χχ→jj On the other hand, Γ i (i = Z , h d ) is the total decay width of the species i and the thermal average of Γ i has been indicated by Γ i T which like the thermal average of annihilation cross section depends on the temperature of the parent particle. However, in Eq. 3.11 and Eq. 3.12, we also have the quantities like Γ i T for Z and h d respectively. These terms actually represent the contributions coming from inverse decay. In that case the initial state particles are the SM particles whose temperature is different from T , as both Z and h d can decay only into the SM particles. Moreover, δ(T ) in Eq. 3.13 has the following expression where m χ is the mass of any of the dark matter species (i.e. m χ = m 1 = m 2 ) and its distribution function is denoted by f χ (p, T ) with p and E p being the magnitude of three momentum and energy respectively. g χ is the internal degrees of freedom which in our case is equal to 2. In Fig. 2, we have shown how δ(T ) varies with the temperature T for two different values of m χ such as m χ = 1 GeV and 100 GeV respectively. From this figure, it is clearly evident that in both cases δ(T ) tends to a very small value when T << m χ (i.e. in the non-relativistic regime). However, δ(T ) increases as T increases and finally saturates to a value of unity in the ultra-relativistic limit (i.e. T >> m χ ).
In the present work, as we have studied the dynamics of a dark sector in the non-relativistic regime, we have set δ(T ) to be equal to zero while solving the above set of Boltzmann equations. A detailed derivation of all the four Boltzmann equations are given in Appendix A. Finally, the last term in the right hand side of Eq. 3.13 contains all the 2 → 2 and 3 → 2 processes having significant impact on the evolution of dark sector temperature. The contribution from relevant 2 → 2 scatterings is denoted by the function F 2→2 (T ). In Appendix B.1, we have derived the analytical expression of F 2→2 (T ) for a generic process like χ i χ i → jj. More importantly, the dark sector temperature T has significant effect from inelastic 3 → 2 scatterings which are shown in Figs Let us define three dimensionless variables namely, co-moving number density Y i = n i (T )/s(T ), x = m 0 /T and ξ = T /T , where m 0 is a reference mass scale which can also be the mass of any dark sector particle. As expected, the final solutions are independent of the choice of m 0 . In terms of these dimensionless variables, the above set of Boltzmann equations take the following form.
In the above, g ρ and g s are effective number of degrees of freedoms associated with the energy and entropy densities of the Universe while g 1/2 * is defined as In Eq. 3.15-Eq.3.18, we have replaced the dark sector temperature T by the newly defined variables ξ and x. The right hand side of Eq. 3.18 has been obtained after substituting the expression of F(T ) 2→2 given in Appendix B.1 into Eq. 3.13 and defining a new quantity σv χχ→jj for an interaction process Thermal average ( Thermal average ( with x concurrently. Here, χ i can be any of the dark matter candidates χ 1 or χ 2 . From this plot, we can see that for large x (i.e for smaller temperature T ), the thermal averaged annihilation cross section becomes independent of x and it coincides with the value of σv χ i χ i →Z Z in the non-relativistic limit denoted by ( σv χ i χ i →Z Z ) NR . It is mainly due to the fact that for larger values of x , the velocity independent term (i.e the s-wave term) of σv χ i χ i →Z Z dominates over the other terms which, at that regime, are heavily velocity suppressed (e.g. p-wave, d-wave terms). However, unlike σv χ i χ i →Z Z , the quantity σv χ i χ i →Z Z keeps on decreasing with increasing x after attaining a peak at T ∼ m χ . The continuous decreasing nature of σv χ i χ i →Z Z can be easily understood if we notice the dotted line indicating the behaviour of 1 T Here, in the large x limit both dotted and dash-dotted curves are overlapping with each other which implies that in the non-relativistic regime (i.e. for x >> 1), the quantity Therefore, from Eq. 3.20 it becomes very straightforward to figure out the decreasing nature of σv χ i χ i →Z Z . The similar thing for a p-wave annihilation cross section χ i χ i → h d h d has been depicted in the right panel of Fig. 3. Here, unlike the previous case of a s-wave process, the thermal averaged annihilation cross section σv χ i χ i →h d h d decreases in the large x limit as the former is proportional to v 2 (neglecting O(v 4 ) term and beyond) which is getting lower with x . Nevertheless, here also 1 T non-relativistic regime and the decreasing nature of σv remains. Now, if we look at Eq. 3.18 once again, we can understand the effects of various terms on ξ (and hence on T ). Let us start from the right hand side. The second term in the right hand side is the only source term which raises the temperature and this is due to 3 → 2 inelastic scatterings. These scatterings convert non-relativistic species into a fewer number of relativistic species having more kinetic energy compared to the mean kinetic energy which is ∼ T . Thus, the term F Y 3→2 (T ) heats up the dark sector and shuts off when all inelastic scatterings are frozen-out. Here, the function F Y 3→2 (T ) represents the contribution from all relevant 3 → 2 scatterings as encountered previously in Eq. 3.13, except all the number densities in F 3→2 (T ) should now be replaced by the corresponding comoving number densities. The excess heat generated due to these inelastic scatterings is smeared over the entire dark sector by the 2 → 2 scatterings so that all the species have a common temperature. This has been incorporated in the first term on the right hand side of Eq. 3.18. Therefore, as the 2 → 2 scatterings reduce the rate of increase of temperature, the first term appears with a -ve sign. Moreover, it is already evident from the discussions given above along with Fig. 3 that the effect of 2 → 2 scatterings become sub-leading at large x. Finally, there is another source of reduction of dark sector temperature which has always been present in the background. It is the expansion of the Universe and the second term in the left hand side confirms this effect. In the large x limit when all the terms in the right hand side become insignificant, ξ varies as 1/x (or T ∝ a −2 (T )) due to the effect of expansion of the Universe.

Constraints on the portal couplings
The portal couplings ( , α) and the mass of the mediators (m Z , m h d ) can be constrained from various experimental observations (See [51] for example). In this section we discuss the relevant constraints in m Z -and m h d -sin α planes.

BBN constraint
In our model, Z and h d can decay into the SM fields through the kinetic mixing parameter and h − h d mixing angle α respectively. Therefore, the lifetime of Z and h d are controlled by these portal couplings. To ensure the observations during the BBN epoch remain unaltered, the lifetime of these fields must be less than one second i.e. they must decay before the BBN substantially. This sets a lower limit on each portal coupling. The constraints from the BBN on and α are shown in Fig. 4 by the light blue coloured regions.

Thermalisation condition
On the other hand, we cannot choose very large value for the portal couplings since the upper limits of these couplings are set from the requirement of early kinetic decoupling of the dark sector from the visible sector. Now, at the time of decay of Z and h d , However, in order to ensure that the dark sector will not re-equilibrate with the SM bath, we have adopted a more conservative lower limit x Γ > 5m 0 /m j , 7 as mentioned in [37]. Therefore, the requirement x Γ > 5 sets an upper limit on the portal couplings and α. The disallowed region of the parameter spaces are indicated by light green color in Fig. 4.

Direct Detection
Due to the presence of kinetic mixing portal, our DM candidates χ i can interact with the nucleon exchanging Z and Z bosons. Since χ i is a Majorana fermion, therefore only the axial vector couplings are present in our model. Consequently, both Z and Z mediated DM-nucleon scattering 7 As in this work we are considering a degenerate dark sector, the condition in the present case is effectively xΓ > 5.
contributes to the spin dependent (SD) scattering cross section. Moreover, the DM-nucleon scattering can be mediated through h and h d also as there is h-h d mixing in our model and it contributes to the spin independent (SI) elastic scattering cross section. The spin dependent DM-nucleon scattering cross section mediated by Z and Z is given by the following expression [52,53]. (3.21) In the above, are the axial vector couplings of χ i (q) with Z and Z bosons respectively. The quantity µ χn is the reduced mass of the DM-nucleon system. The contribution of proton (neutron) to the nuclear spin is denoted by S p (n) . The values of ∆ n q and ∆ p q are: The spin independent DM-nucleon scattering cross section mediated by h and h d is given by [52] σ In the above expression Z(A) is the atomic (mass) number of a nucleus while m p and m n are the mass of proton and neutron respectively. Moreover, [54].
For both the cases discussed above an effective DM-nucleon scattering cross section is defined to compare with the current bounds from XENON1T [18,55], CDMSlite [56,57], CRESST-III [58]. In Fig. 4, we show the spin dependent (left panel) and spin independent (right panel) direct detection constraint for g X = 1 by two purple coloured regions.

Supernovae cooling
Dark sector particles with mass 200 MeV can be constrained from the observation of supernova SN 1987A [59,60]. The constraints are derived based on the following two conditions.
i) The emissivity of the hidden sector particles must be smaller than that of the observed neutrino signal. This is known as "Raffelt Criterion" [61] and the constraint derived from this condition is known as "cooling constraint". This condition puts an upper bound on the portal couplings.
ii) The mean free path of the hidden sector particle must be smaller than the radius of the supernova (SN). The constraints derived from this condition put a lower bound on the portal coupling and this is known as "Trapping condition".
In this work, we have considered the bound from SN 1987A cooling [59,62,63] relevant for the dark gauge boson Z in our model. Similarly h d can also contribute to the SN 1987A cooling because of the presences of h − h d mixing angle α [64]. In both panels of Fig. 4, the red coloured regions indicate the excluded parameter space from the SN 1987A cooling.

Beam dump experiments
The dark gauge boson Z with mass 1GeV and 10 −8 −10 −2 can be constrained from the result of beam dump experiments. The excluded parameter space [65] is shown in the left panel of Fig. 4 using magenta color. The dark scalar h d of mass 1 GeV and portal coupling sin α 10 −4 − 10 −2 can be constrained from proton beam dump experiments by CHARM collaboration [66]. Following [67,68], we have derived the bound in m h d − sin α plane and the excluded parameter spaces is shown in the right panel of Fig. 4 with magenta color.

Electroweak Precision Observables
The presence of kinetic mixing parameter can alter the decay width of Z boson and also all the couplings involving Z boson. To study the constraints from Electroweak Precision Observables (EWPO), we calculate the S and T parameters for our model in terms of the measured Z boson mass and the Weinberg angle by comparing our Lagrangian with the effective Lagrangian formulation of Zf f current [69,70]. The best fit values of S = 0.06 ± 0.09 and T = 0.10 ± +0.07 [71] are used in deriving the constraint and it has been shown in the left panel of Fig. 4 by the orange coloured contour.

Dark matter Self Interaction
We have calculated the momentum transfer cross section σ T for DM self interaction processes mediated by Z and h d . For multicomponent DM, we define an effective momentum transfer cross section σ eff 3) and from the Bullet cluster observation it is bounded as σ eff T < 1.25 cm 2 (m χ /g) [72]. Using this relation, we have found g X < 1.87 × 10 4 (m χ /1 GeV) 3/4 and our parameter space of interest is well below the exclusion limit.

Numerical results
We have numerically solved four coupled Boltzmann equations given in Eq. 3.15-3.18 to find the evolution of comoving number density of each dark sector species with temperature T . The initial condition 8 for solving these equations is at x = 0.1, ξ = 0.1 (i.e. x = 1) and Y i = Y eq i . We have checked that with g X ∼ 1 both sectors are kinetically decoupled from each other at x = 0.1 for the considered ranges of portal couplings. Using the solutions of coupled Boltzmann equations, one can easily estimate the observable parameter namely the dark matter relic density from the following well known relationship Ω χ h 2 = 2.755 × 10 8 (m χ /GeV) Y 0 χ , where Y 0 χ is the asymptotic value of Y χ at the present era. The numerical results are presented in Fig. 5. As we have discussed in the Section 3 that the degeneracy among species of the dark sector leads to a completely different freeze-out dynamics compared to the usual thermal freeze-out of dark matter and this has been known as the co-decaying framework. Therefore, in this numerical analysis we have considered mass degeneracy in the dark sector by appropriately tuning (see footnote 3) the relevant couplings.
In the upper panel of  number density Y eq χ is also shown by the black dash-dotted line. Initially, the number density of dark matter follows the equilibrium number density. Thereafter it deviates from equilibrium, starts depleting and finally freezes out to a particular density. The time of departure from the equilibrium depends on when the decay modes of Z and/or h d into the SM particles open up. This creates chemical imbalance in the 2 → 2 interactions between χ i and other species (Z , h d ) leading to the freeze-out of χ 1 and χ 2 . For example, in the left and middle plot for m = 0.1 GeV and 1 GeV respectively, the freeze-out occurs at larger values of x (x > 100) and much larger values of x (depending on the parameter ξ which is at least 0.1 or even less for x > 100 as seen from the plots in the lower panel) while in the right plot dark matter freezes-out relatively early. This is because, in the plot for m = 30 GeV we have chosen higher values of g X and and hence decay of Z starts much earlier compared to the other two cases. The relevant model parameters such as mass (m), gauge coupling (g X ), kinetic mixing parameter ( ) and scalar sector mixing angle (α) in all three plots are adopted in such a way so that we achieve the right dark matter relic abundance.
The evolution of dark sector temperature with x for the same choices of masses have been depicted in the lower panel of Fig. 5. The trivial variation of the SM temperature T with x (inverse of T ) is also shown by a black dashed line in each plot. From these plots it is seen that at first T decreases with x and thereafter the variation of T with x becomes almost constant between 1 x 100 (in the left and middle plot) and 0.3 x 10 (in the right plot) respectively. This is the "phase of cannibalism" where all the inelastic 3 → 2 scatterings among the dark sector species are active and generating excess heat in the dark sector. A detailed discussion about these inelastic scatterings are given in Appendix B.2 along with all possible Feynman diagrams in Figs. 16-25. As mentioned earlier, the heat generated due to 3 → 2 scatterings spreads equally over the entire dark sector by the 2 → 2 scatterings so that all species share a common temperature. However, expansion of the Universe red-shifts the temperature and it is always present in the background. Due to these two opposite effects, the dark sector temperature becomes a slowly varying function of x. Actually, in this regime T varies logarithmically with T as shown in Eq. 3.8. Moreover, it is also seen from these plots in the lower panel of Fig. 5 that after this plateau region, T sharply decreases with x as a result of both 2 → 2 scatterings as well as the Universe's expansion. Finally, for large x when all the scatterings are frozen-out, T redshifts as non-relativistic matter due to the effect of expansion of the Universe. We would like to spend a few sentences on the chemical potential of species in the dark sector. From the earlier discussion in Section 3, we know that during the number changing interactions by 3 → 2 processes, the chemical potential of each species is zero. After freeze-out of these inelastic scatterings at a temperature T c , each species develops a nonzero chemical potential which helps them to get rid of the exponential suppression in number density even if they become nonrelativistic. We have derived the expression of chemical potential (µ(T )) for T c > T > T Γ (Eq 3.5) using the conservation of entropy and number density in a co-moving volume. It shows that µ(T ) of a species first increases sharply as T departs from T c and thereafter saturates to the mass m of that particular species when T T c . The exactly similar nature for µ has been obtained from the numerical simulations by solving the Boltzmann equations and it has been shown in Fig. 6. Here, we have shown the variation of µ with x (inverse of T ) for three different values of m = 0.1 GeV, 1 GeV and 30 GeV respectively. We have computed µ from the well known relationship between µ and T for a species i as µ = T ln n i n eq i where, the number density n i is obtained by solving the Boltzmann equations. We have noticed that the analytical expression given in Eq. 3.5 remains valid between T c and T Γ . However, the results we have found from the numerical simulations show that µ continues to maintain its saturation value well beyond the temperature T Γ at which decays of Z and h d into the SM particles begin. In Fig. 7, we have demonstrated the dependence of important model parameters on the comoving number density. In the first plot at the top-left position, we have shown the effect of gauge coupling g X on Y χ for three different values of g X = 0.1 (green solid line), 0.477 (red solid line) and 1.0 (cyan solid line) respectively. The corresponding equilibrium values of Y χ have been indicated by the dashed lines of same colors. Note that g X = 0.477 corresponds to correct relic abundance. This plot clearly shows that increasing g X leads to much suppressed final abundance and late freezeout of dark matter. For example, changing g X by one order of magnitude (g X = 0.1 to g X = 1.0) reduces Y χ by more that three orders i.e from Y χ = 10 −7 to Y χ < 10 −10 . In the second plot at the top-right position shows the effect of kinetic mixing parameter on Y χ . Here also we have chosen three different values of namely = 10 −11 (green), 10 −10 (red) and 10 −9 (cyan) respectively. From this plot one can notice that the impact of on Y χ depends on the value of . For example, altering from 10 −11 to 10 −10 does not reflect a significant change in Y χ . However, increasing further by one order of magnitude (i.e. from 10 −10 to 10 −9 ) results in a noticeable difference in the final abundance. Moreover, in the latter case much earlier freeze-out of dark matter has occurred. This is mainly due to the fact that increasing enhances the decay width Γ Z which reduces the lifetime of Z . Hence, Z decays much earlier into the SM particles and creates enough chemical imbalance in the dark sector which thereafter leads to earlier freeze-out of the dark matter. The dependence of another portal coupling α on Y χ is exactly similar to the earlier case of . It is shown in the bottom-left plot. Analogous to the previous cases, we have adopted three different values of scalar mixing angle α = 10 −8 rad (green), 10 −7 rad (red) and 10 −6 rad (cyan) respectively. The equilibrium values of Y χ in each case has been shown by the dashed line. Here also we see much earlier dark matter freeze-out and more suppressed final abundance when we increase α from 10 −7 rad to 10 −6 rad. Finally, in the bottom-right plot we show dependence of mass on Y χ . In this plot, we have considered four different values of m as indicated in the plot legend and portal coupling parameters are kept fixed to a particular value as = 10 −8 and α = 10 −6 rad. Moreover, we have also varied the gauge coupling g X so that the final value of Y χ for each value of m reproduces the correct relic abundance. We observe that for a fixed portal couplings as we increase m, we require higher values of g X to achieve Ω χ h 2 in the right ballpark value. This parametric dependence of Y χ on m and g X will be more clear when we analyse Fig. 9.
In both the panels of Fig. 8 we present the allowed parameter space that we have obtained by solving four coupled Boltzmann equations and comparing our results with the observed value of relic density in 2σ range i.e. Ω DM h 2 = 0.120 ± 0.002 [2]. In the upper panel, we have shown the − g X parameter space for two distinct dark matter masses 1 GeV and 30 GeV respectively. The corresponding variation of the third parameter α has been indicated by the respective colour bar in each plot. The constraints coming from BBN, prevention of thermalisation of the dark sector with the SM bath and measurement of diffuse γ-ray fluxes by the INTEGRAL space telescope have also been indicated. From the plot in the top-left position, we can notice that allowed range of decreases sharply as g X increases from 0.1 to 1.0. Additionally, the production of correct relic abundance at the present era demands smaller α for higher g X . For example, we require α 10 −9 rad when g X 0.8 while higher mixing angles α ∼ 10 −6 rad are necessary for smaller couplings (g X 0.2). More or less similar feature has also been observed for the other dark matter mass (m = 30 GeV) in the top-right plot. In this case, reduction of parameter space of for increasing g X is not so sharp and the variation of α in the scattered region mostly concentrated around 10 −7 rad α 10 −6 rad. Moreover, unlike the plot for m = 1 GeV, we have not depicted any constraint in the − g X plane for m = 30 GeV from INTEGRAL as in this case the corresponding bound is on the higher values of gauge coupling (g X >> 1). The allowed regions in the − α parameter space have been shown in the lower panel of Fig. 8 for m = 1 GeV and 30 GeV respectively. We observe a definite pattern in the allowed values of and α in order to satisfy the correct relic abundance of dark matter.
In the bottom-left plot, almost entire range of α allowed from BBN and thermalisation condition is also allowed from the relic density criterion for 10 −10 . It implies that for these values of the chemical equilibrium in the dark sector is disturbed mostly by the decay of Z . In contrast, an opposite situation is also observed when the scalar mixing angle α 10 −7 rad. In this case, the entire range of satisfies relic density in 2σ range expressing Γ h d dominance in the freeze-out dynamics of dark matter. In the intermediate region, decays of both Z and h d are responsible for the departure of chemical equilibrium in the dark sector. The similar feature is seen for m = 30 GeV (bottom-right plot), however, in this case only very small portion of − α plane is allowed. In both these plots, possible variation of gauge coupling g X is indicated in the respective colour bar. Similar to the plots in the upper panel, here we have also shown the other constraints in the − α plane from BBN and thermalisation using two different colours. Moreover, we have tried to find the allowed values of g X when mass of the dark sector species varies between 100 MeV to 100 GeV. The result has been shown in Fig. 9 where varies between 10 −12 to 10 −8 and is indicated by the colour bar. The green shaded region is the disallowed as x Γ is less than five here while the light yellow and the light blue shaded regions are ruled-out from the measurements of diffuse γ-ray background by EGRET and INTEGRAL respectively. A detailed discussion on the indirect detection prospects of the co-decaying dark matter scenario has been given in Section 4. We have obtained this m − g X parameter space for a fixed value of scalar mixing angle α = 10 −6 rad. From this plot we can see that generally as m increases we need higher values of g X (i.e. higher annihilation cross section) to satisfy the observed relic abundance (within 2σ range) except for a couple of values of m near 200 MeV and 2 GeV respectively where decay width 9 of Z and h d suddenly enhance due to the crossing of kinematic thresholds for new decay modes. In this particular case, the decay width Γ h d increases substantially due to opening of ss and cc channels 9 For 350 MeV m 5 GeV, the hadronic decay modes of Z and h d are also present. Following [73,74], we have checked the effect of these hadronic decay modes on the relic density and we find that the impact of inclusion of these hadronic decay modes is not very significant in m − gX plane since relic density depends on g 4 X .
and it dominates the freeze-out process of dark matter over Γ Z for < 10 −10 and α = 10 −6 rad (see Fig. 10 and related discussions). The sudden enhancements in Γ h d are compensated by the curtailments in dark matter annihilation cross sections through g X as Ω χ ∼ m σv Γ Z (h d ) (see Appendix D for an approximate analytical expression of relic density). Let us note that in Fig. 9, the threshold effects are prominent only for ss and cc channels. This can be understood clearly from Fig. 10. In this figure we have plotted the total decay width of h d (Γ h d ) for α = 10 −6 rad and the total decay width of Z (Γ Z ) for = 10 −8 , 10 −10 as a function of mass m. It is evident from Fig. 10 that for m 5 GeV, the decay width Γ h d is dominant over Γ Z for all values of we have considered and in this region the relic density is controlled by Γ h d only.

���
Since the threshold effect for the bb channel is mild in Γ h d , the corresponding sharp feature is absent in Fig. 9. The situation becomes different for the case of lighter dark matter where m < 5 GeV. In this case, the freeze-out process of dark matter is not entirely dominated by Γ h d . Instead, the chemical imbalance prior to the freeze-out can be created either by the decay of h d or by Z or both depending on the portal couplings α and . As a result, we have two distinct regimes for m 3 GeV in the m − g X plane where the allowed points forming a line corresponds to the small values of < 10 −10 and for these set of allowed g X and m, the decay width of h d controls the relic density. Consequently, the sharp features are clearly visible for < 10 −10 due to the threshold effects of cc and ss channels. Moreover, as the threshold effects are mild in Γ h d for the other decay modes like µ + µ − and τ + τ − , we have not observed any peak for m = 2 m µ and 2 m τ respectively in Fig. 9. On the other hand, the scattered points in the low mass regions are for 10 −10 where from Fig. 10 it is seen that Γ Z Γ h d and hence the relic density has been controlled by the dynamics of Z instead of h d . Additionally, we would also like to mention that the threshold effect in Γ h d is more prominent in comparison with that in Γ Z due to the presence of the fermion mass in the Yukawa coupling of h d with the SM fermions. when there is a finite mass splitting between the χ (species A) and Z , h d (species B). Here, we have kept m Z and m h d fixed at 1 GeV and consider the m χ = m Z /(1 − δ m ). In Fig. 11 we have shown the variation of Y χ with x for four different values of δ m = 0, 10 −4 , 10 −3 and 10 −2 respectively. We find that the features of the co-decaying scenario (i.e. freeze-out at large x, x f >> 20) remains still valid for mass splitting δ m ≤ 10 −2 . Beyond which we do not find any stable numerical solution of Y χ . The plots shown in Fig. 11 are drawn for g X = 0.477, = 10 −10 and α = 10 −7 rad so that we get the observed relic abundance for δ m = 0 (red line). As we gradually increase δ m from zero, freeze-out of χ occurs earlier with a reduced final abundance. From the above figure it is clearly seen that for increasing δ m from 0 to 10 −2 , the comoving number density of χ reduces nearly two orders of magnitude while the corresponding x f (inverse of the freeze-out temperature) changes from ∼ 400 to ∼ 250.

Indirect signature: γ-rays and Neutrinos from dark matter annihilation
In this work, we also investigate the possibility of detecting γ-ray and neutrino flux from cascade processes 10 having a general structure like For the γ rays, A, X, and Y are h d , γ and γ respectively whereas for the neutrinos in the final state, the intermediate particle, instead of h d , is the U (1) X gauge boson Z . This kind of cascade process generates a polynomial box-shaped spectrum of outgoing particles where shape of the polynomial depends on the polarisation of the intermediate particle A. In the rest frame of A, the energy of X (or Y ) is m A /2. In the Laboratory frame, assuming non-relativistic nature of dark matter, the energy of outgoing particle is given by 1) 10 We are not considering the direct production ofνν and γ γ from dark matter annihilations because those processes are suppressed by the portal couplings such as , λ .
where θ Lab is the angle between the direction of motion of A and the direction of emission of X(Y ). It is clear from Eq. 4.1 that the resulting spectrum has a sharp cut-off at θ Lab = 0 and θ Lab = π and width of the spectrum depends on mass splitting between the incoming dark matter particle and intermediate state A, i.e. ∆E X = E max X − E min X = m 2 χ − m 2 A . Thus in the exact degenerate limit of χ i and A, this box-spectrum will be transformed into a line. In the present case, as we have mentioned in the previous sections, all the dark sector species are degenerate. Therefore, line like signal of neutrinos and γ rays can be originated from the annihilation of dark matter which is frozen-out through co-decaying scenario. In the present scenario, the final state radiation and the inverse Compton scattering can also result in γ-ray flux from DM annihilation. In Section 4.1, indirect detection prospect of our DM candidate from annihilation is discussed for mass splitting parameter δ m = (m χ − m Z (h d ) )/m χ = 0 whereas the same is discussed in Section 4.2 in the light of non zero δ m .
The differential energy spectrum of X, produced from the decay of an arbitrarily polarised intermediate particle A, has the following general form Here, m is the polarisation of A and n is the number of X particle produced in the final state. Particular forms of the function f m (E X /m χ ) depends on the spin of the intermediate state. The differential flux of X coming from annihilation of χ i in the Milky Way halo and measured at the earth is given by In the above σv χ i χ i →AmAm is the thermal averaged annihilation cross-section of χ i into a pair of A with polarisation m, the quantity Br(A m → X Y ) is the branching ratio of A m → X Y while the definition of Br m is Br m = σv χ i χ i →AmAm σv χ i χ i →AA . Since we have two dark matter candidates of identical properties, they will contribute equal amount to the local dark matter density ρ at the solar location. Hence, an extra 1/4 factor has appeared in the above expression of differential flux. The total differential flux arising from annihilations of both the dark matter candidates is noting but the above expression multiplied by an extra factor of 2. The J-factor averaged over solid angle ∆Ω for dark matter annihilation is defined in terms of the galactic co-ordinate system (b, l) as where, R = 8.5 kpc is the distance of the solar system from the centre of our Milky way galaxy and s is the line of sight (l.o.s) distance which is integrated over 0 to s max . The expression of s max is given by Moreover, Eq. 4.4 can be framed in terms of angle θ GC where θ GC being the angle between the galactic centre and the line of sight distance. In this co-ordinate systemJ ann takes the following form,J The upper limit of the l.o.s distance in this co-ordinate system is s max = R 2 MW − sin 2 θ GC R 2 + R cos θ GC . In Eqs. 4.4 and 4.6, ρ = 0.3 GeV/cm 3 is the dark matter density at the solar neighbourhood, R MW = 40 kpc is the size of the Milky Way (MW) galaxy halo and the solid angle ∆Ω represents the field of view of the detector. Moreover, throughout our analysis, we have assumed the DM density follows the standard Navarro-Frenk-White (NFW) [75] density profile.

γ-ray signal from dark matter annihilation
In our model, γ-ray flux is composed of three different components such as γ-ray line signal, final state radiation (FSR) from charged particle final states and inverse Compton scattering (ICS) of the charged particles with the CMB photons. We have discussed all the three components in the following sections.
• γ-ray line signal: The monochromatic photons coming from dark matter annihilation/decay are considered as one of the important smoking gun signatures of dark matter. Motivated by this, we have discussed the prospect of γ-ray line in the present scenario due to dark matter annihilation at the galactic centre. In our model, the dark scalar h d mixes with the SM Higgs. Thus, a line shaped γ-ray spectrum can be generated form χ i χ i annihilation through one step cascade process 11 like χ i χ i → h d (h d → γγ). Since the intermediate state in this case is a spin 0 boson h d , the energy spectrum of the emitted photons is a Dirac delta function for m χ = m h d . The differential photon spectrum is given by Here, the prefactor 2 in dN Line γ dE γ arises because for each decay of h d , two photons are emitted. Now, using Eq. 4.7, we can calculate the γ-ray flux from Eq. 4.3. The differential photon flux is given by where Br(h d → γγ) is the branching ratio of h d into a pair of γ and σv χχ→h d h d is the dark matter annihilation cross section into a pair of h d which is velocity suppressed in the 11 The lifetime of h d is much smaller than tBBN in the parameter space in which we are interested in. For example, for m h d = 10 GeV, α = 10 −6 , τ h d 10 −5 s << tBBN.
non-relativistic limit (see right panel of Fig 3). The analytical form of σv χχ→h d h d in the non-relativistic limit for exactly degenerate dark sector is given by Here v is the magnitude of relative velocity of the initial state particles and v is the average thermal velocity of DM. This is to be noted that the leading order term is proportional to v 3 instead of v at threshold (δ m = 0). This is because of the absence ofâ 0 term in the cross section (see Eq. F.5 of Appendix F for details).
• Final State Radiation: As discussed in [76,77], the FSR contribution to γ-ray flux can constrain our parameter space in m χ − g X plane. To study the γ-ray flux from FSR, we have considered the cascade process χχ → Z (h d )Z (h d ) and the subsequent decay of Z (h d ) → e + e − + FSR. The energy spectrum of photons in the center of mass (CoM) frame of the DM is given by 12 [79] dN FSR . (4.10) Here f denotes the final state fermion having mass m f , α EM is the fine structure constant and x = E γ /m χ where E γ is the energy of the emitted photon in the CoM frame of DM annihilation. Therefore, using Eq. 4.3 and Eq. 4.10 we can write down the differential photon flux from FSR as In the non-relativistic limit σv χ i χ i →h d h d is given in Eq. 4.9 and whereas σv χ i χ i →Z Z is given by the following relation.
Here also the leading order term of σv χ i χ i →Z Z is proportional to v 3 . It is due to the fact that for δ m = 0, the coefficientâ 0 in Eq. F.5 vanishes as it is proportional to δ 3/2 m (see Appendix F for details). 12 Here, we have not considered the angular dependence in the differential FSR spectrum. However the effect of the angular dependence is not very significant as discussed in [78].
• Inverse Compton Scattering: In the present scenario, high energy electron positron pair can be produced from DM annihilation via one step cascade processes through mixing parameters and α. The Inverse Compton scattering between these charged particles and ubiquitous CMB photon can produce high energy gamma ray fluxes which may be within the reach of current γ-ray detectors. To study this prospect, we have calculated the γ-ray flux following [76,80].
Finally, the total photon flux will have all three contributions from different sources as In Fig. 12, we have shown the E 2 γ weighted total differential photon flux as a function of the energy of emitted photons for three benchmark points which satisfy the relic density constraint. Here the fluxes are calculated for the Region of Interest (ROI) having (330 The computed flux has been compared with the available data for diffuse gamma-ray from EGRET (Fig. 4 of [42]), Fermi-LAT (Fig. 6 of [41]) and INTEGRAL (Fig. 7 of [43]) collaborations. The ROIs considered in above are adopted from the respective experimental collaboration to be as close as possible with the experimental data.
It is clearly seen from Fig. 12 that the γ-ray flux due to dark matter annihilation calculated from the ROI for EGRET is larger than that for Fermi-LAT and INTEGRAL. We know that for a cuspy density profile like the NFW profile the value ofJ ann increases sharply as we move towards low latitude. Therefore, the value ofJ ann for the ROI of the EGRET collaboration is much larger compared toJ ann for the Fermi-LAT and INTEGRAL. In the left panel, we have considered δ m = 0 whereas the results in the right panel are for δ m = 10 −3 . The enhancement in the photon flux as seen in the right panel compared to that in the left panel has been discussed in the light of nonzero mass splitting in Section 4.2.

Neutrino line from DM annihilation
In our model, the kinetic mixing portal can produce an interesting neutrino signal from dark matter annihilation. For the exact degenerate case (i.e. δ m = 0), the neutrino spectrum is a line spectrum. Here we have considered the following cascade process χ i χ i → Z (Z →ν µ ν µ ) and the emitted neutrinos are monochromatic. The neutrino energy spectrum is given by where m is the polarisation of the decaying particle i.e. Z . Note that the polarisation of Z has no effect on the spectrum of monochromatic neutrinos. However, the polarization of Z has non-trivial effect for the quasi degenerate case as discussed later. The 2 factor in Eq. 4.14 is for ν µ andν µ . Using Eq. 4.3, we can easily write the differential neutrino flux for ν µ where the energy spectrum andJ ann are obtained from Eq. 4.6 and Eq. 4.14 respectively.
As the dark matter candidates are non-relativistic, we have found Br 0 = 0. Therefore, in Eq. 4.15 we have used [81,82] Br −1 = Br +1 = 1/2. Since the annihilation cross section σv χ i χ i →Z Z is velocity suppressed for δ m = 0 and the branching ratio of Z intoν µ ν µ is extremely small, the resulting neutrino flux is well below the observed flux by Super-Kamiokande detector.

Effect of nonzero mass splitting (δ m = 0) on the Indirect detection prospect
If we consider the mass splitting parameter δ m = 0, the γ-ray and ν µ fluxes change significantly. In the following, we have discussed the γ-ray and ν µ fluxes for δ m = 0.
• γ-ray signal for δ m = 0: As discussed earlier, the parameter δ m = 0 breaks the degeneracy between the DM χ i and the mediators Z and h d . As a result, the γ-ray energy spectrum becomes a box like spectrum instead of a line like one for the δ m = 0 case. Thus, the energy spectrum of the photons from the annihilation of DM into two h d s and the subsequent decay of h d into two photons is given by Here, r h d = m h d /m χ and z = E γ /m χ . The cut-off energies for the box-spectrum are E γ min Lab and E γ max Lab and these can be calculated from Eq. 4.1. Origin of the factor 2 is mentioned in Eq. 4.7. Similar to the case of δ m = 0, the total photon flux for δ m = 0 also contains three components such as γ-ray box spectrum and γ-rays from FSR and ICS. Therefore, using Eq. 4.3 and Eq. 4.13, we calculate total photon flux for nonzero mass splitting. Note that the line contribution in Eq. 4.13 should be replaced by the box contribution for non zero δ m . The main difference in this case compared to the previous one (δ m = 0) is that now the annihilation cross section σv rel χ i χ i →Z Z is not velocity suppressed and it scales as 13 δ 3/2 m . Therefore, the final γ-ray flux also scales as δ 3/2 m . In the right panel of Fig. 12, the total photon flux from DM annihilation is plotted for δ m = 10 −3 . As we can see from the plot that the computed flux is several order of magnitude larger than the case for δ m = 0 and at some energy range with certain choices of model parameters it overshoots the observed flux. This is primarily because of the presence of the s-wave term in σv rel χ i χ i →Z Z .
• Neutrino line from DM annihilation for δ m = 0: For δ m = 0, the neutrino spectrum due to one step cascade process is a polynomial spectra with sharp cut-off at E min ν Lab and E max ν Lab and the shape of the polynomial depends on polarisation of Z . Therefore, following [83], we can calculate the differential energy spectrum of neutrinos f m (E ν /m χ ) for an arbitrarily polarised Z having polarisation m. The functional forms of the spectrum for different polarisation states of Z are given by m . In Fig. 13, we have shown the variation of E 2 ν weighted differential neutrino flux for ν µ as a function of E ν for δ m = 10 −3 and compared it with the observed ν µ flux by the Super-Kamiokande (SK) detector. We have found that the neutrino signal from dark matter annihilation is several order of magnitude lesser than that of the observed flux by the SK detector at low energy. This is primarily due to the fact that branching ratio of Z toν µ ν µ becomes heavily suppressed at low value of m Z (m Z 10 GeV). Therefore, although the annihilation cross-section σv χ i χ i →Z Z is not velocity suppressed in the non-relativistic limit for δ m = 0 (see left panel of Fig. 3), the resulting neutrino flux lies well below the observed flux due to additional suppression in the branching ratio of Z →νν channel. We have divided the full range of cos θ GC into ten bins as shown by the SK collaboration [84] and combined fluxes from all such bins to show the variation with E ν . As we have assumed the dark matter density profile to be the cuspy NFW profile, therefore, we would expect the neutrino flux from dark matter annihilation has an angular variation. For the cuspy profile, the dark matter density and hence the J-factor increases as we go towards the galactic centre. Therefore, the number of annihilations are also very large near the galactic centre because of the huge amount of dark matter. Thus, a large amount of neutrino flux due to dark matter annihilation from the nearby regions of the galactic center is expected. In the left panel of Fig. 14 we show the variation of differential neutrino flux as a function cos θ GC for δ m = 10 −3 , m χ = 5 GeV and m χ = 50 GeV. This plot reveals that the differential flux dΦ ν dE νµ is higher for the neutrinos arising from annihilation of 5 GeV dark matter. The energy E νµ is always larger for neutrinos coming from m χ = 50 GeV dark matter than those are from m χ = 5 GeV dark matter. However, the differential flux proportional to σv χ i χ i →Z Z /m 2 χ becomes heavily suppressed as m χ increases.
Finally, in the right panel of Fig. 14, we have presented the variation of dark matter annihilation cross section into neutrinos as a function of dark matter mass. Since we have considered an one step cascaded process to produce a pair ofν µ ν µ , therefore we have plotted σv χ i χ i →Z Z × Br(Z →ν µ ν µ ) as a function of dark matter mass (m χ ) for two different values of the dark sector gauge coupling (g X ) and compared it with the current upper limit on DM DM →νν channel as given in [85]. Here also we have considered the mass splitting parameter δ m = 10 −3 . From this plot we can notice that the effective annihilation cross section of dark matter into a pair of neutrinos in the present scenario is several orders of magnitude smaller than the current upper limit. Thus, our model easily avoids the indirect detection constraint onνν channel by the SK collaboration.
It is to be noted that for δ m = 0, rate of production of SM charged particles from DM annihilation is sufficient enough to modify the ionisation history of Hydrogen and Helium at the time of CMB. The observation from the Planck [86] experiment puts a strong constraint on DM mass as discussed in [87]. Moreover, another constraint is coming from the measurement of positron flux by the AMS-02 detector [88] as discussed in [89]. In Fig. 15, we show the constraints in the m − g X plane from CMB measurements and AMS-02 experiment for δ m = 10 −3 . We have also depicted two relic density satisfied lines in the m − g X plane for two fixed values of the kinetic mixing parameter = 10 −8 (dashed line) and = 3 × 10 −9 (dotted line) respectively. From this figure one can clearly notice that the relic density satisfied parameter space for = 10 −8 and m 1 GeV is allowed by the AMS-02 measurements of positron flux in cosmic rays as well as CMB measurements by Planck. However, for = 10 −8 , the relic density satisfied parameter space do not satisfy the criterion x Γ > 5 required for preventing the thermalisation of the dark sector with the SM bath (see left panel of Fig. 4). The thermalisation criterion is more relaxed for the portal coupling < 10 −8 . As a result, we have found a parameter space for = 3 × 10 −9 (dotted line) and 1 GeV m 10 GeV which produces correct dark matter relic density and at the same time is allowed from CMB measurements by Planck and positron flux measured by AMS-02. Additionally, we would like to note that here we have chosen α = 10 −7 rad to avoid the constraint arising from the thermalisation criterion (see right panel of Fig. 4).

Summary and Conclusion
In this work, we have considered an anomaly free minimal U(1) X extension of the SM by considering two SM gauge singlet left chiral Weyl fermion of opposite U(1) X charge. The U(1) X symmetry of the dark sector is broken by the vacuum expectation value of a complex scalar η. Thus, the dark sector contains a spin zero scalar h d , two spin half Majorana fermions χ i (i = 1 , 2) acting as the dark matter candidates and a spin one massive gauge boson Z . Apart from that, the dark sector can talk to the visible sector through the kinetic mixing portal and also through the h − h d mixing angle α. Since, the couplings within the dark sector particles are unknown, we choose all the couplings in the dark sector are of the same order. This "democratic" choice would lead to degenerate/quasi-degenerate dark sector which has very rich phenomenology. Under the assumption of the feeble portal couplings, the dark sector is kinetically decoupled and it evolves with a different temperature compared to the visible sector. Due to the degeneracy within the dark sector, the number density of dark species are exponentially suppressed only when the mediator particles (Z , h d ) start to decay into the SM fields. Since, the decays are out-of-equilibrium, the exponential suppressions in number densities arise much later compared to the standard freezeout mechanism leading to a delayed freeze-out of dark matter and this is known as co-decaying dark matter. In this scenario, we solved three coupled Boltzmann equations (for Y χ 1 +χ 2 , Y Z , Y h d ) along with the temperature evolution equation to study the dynamics of kinetically decoupled dark sector taking into account all the 2 → 2 and 3 → 2 processes within the dark sector. We show that the presence of 3 → 2 processes introduces the cannibal phase and it changes the dynamics of the dark sector temperature evolution significantly. The dark sector dynamics essentially depends on four parameters i.e. g X , , α, m χ and we have presented the parameter spaces allowed from the relic density constraint. Apart from that, we show parameter space for the portal couplings allowed from the BBN observations, criterion for kinetically decoupled dark sector, beam-dump experiments, SN1987A cooling, direct detection, and Electroweak Precision Observables.
Furthermore, we have investigated the prospect of detecting our dark matter candidates due to annihilation around the galactic centre of the Milky-Way galaxy. Due to the presence of the portal couplings, both χ 1 and χ 2 can produce neutrino and γ-ray line spectrum from one step cascade processes. We have calculated the neutrino flux from a cascade process like χ i χ i → Z (Z → ν µ ν µ ) and compared our result with the observed atmospheric neutrino flux by Super-Kamiokande neutrino detector. We have also calculated the γ-ray flux from DM annihilation. In our scenario, γ-ray flux composed of three components such as γ-ray line signal from χ i χ i → h d (h d → γγ), final state radiation from χ i χ i → X(X → e + + e − + FSR, X = h d , Z ), and ICS of final state charged particles with CMB photons. The calculated total flux is then compared with the available diffuse background γ-ray data from Fermi-LAT, EGRET, and INTEGRAL. We have found that the allowed parameter space obtained after satisfying the relic density bound easily satisfies all the constraints coming from the indirect detection and CMB anisotropy measurement by Planck for the degenerate dark sector. However, a certain region of parameter space for the quasi degenerate dark sector has already been ruled-out from the measurements of diffuse γ-ray flux by INTEGRAL, the CMB anisotropy, and positron flux from AMS-02.

Acknowledgements
Authors would like to acknowledge Alejandro Ibarra for a very helpful and informative email communication regarding indirect detection. One of the authors SG would like to acknowledge University Grants Commission (UGC) for financial support as a senior research fellowship.

A A brief discussion on Boltzmann equation
In this section we have derived the necessary Boltzmann equations for the dark sector. The evolution of phase space distribution function f A (p, T ) of a species A having p and T as magnitude of three momentum (p ≡ | p|) and temperature respectively is governed by the Liouville equation which has the following form Here H(t) is the Hubble parameter and C[f A (p, T )] is the usual collision term for the evolution of f A (p, T ). Due to homogeneity and isotropy of the Universe, the distribution function depends on p only and not on the individual components. Now, for a generic process like A 1 + A 2 + ..... + A n → B 1 + B 2 + B 3 + .....B m , the moment for any quantity (say O(p κ )) for A κ is given by where p κ and g κ are three momentum and internal degrees of freedom of A κ respectively. Assuming the phase space distribution function vanishes at the boundary, Eq. A.2 can be written as In the above, we have ignored the Pauli-blocking and the stimulated emission terms for fermions and bosons as these terms will not be significant in our case where we have a non-relativistic dark sector obeying the Maxwell-Boltzmann distribution. Here, p α s and k β s are the three momenta of the initial and the final state particles and the corresponding four momenta are denoted by P α , K β respectively. The Lorentz invariant phase space measure dΠ α = g α d 3 p α 2E pα where α stands for the initial state particles. For the final state particles the index α and p α in dΠ α are replaced by β and k β respectively. Moreover, the Lorentz invariant matrix amplitude square |M| 2 for the process .B m is averaged over spins of both initial and final state species. Furthermore, the quantity O A i is a function of T only and it is defined as Let us consider a particular case when there are x and y number of identical A species in the initial and the final state respectively i.e. the process we are considering is as follows In this situation we need to keep in mind that under exchange of p 1 , p 2 , ... p x in the initial state and k 1 , k 2 , ... k y in the final state the process remains unaltered. Hence, in order to take into account that effect the collision term must be divided by x! and y! to avoid overcounting. Now, if we want to find the time evolution of O A (T ) for the species A, then we need to write the Boltzmann equations (Eq. A.3) for all the As (both for initial as well as final state) and then add them. After a few algebraic simplifications, we have arrived at the following form In the right hand side of the above equation, we have used the freedom of permutation symmetry of p α s (for α = 1 to x), k β s (for β = 1 to y) and 1 ≤ r ≤ x, 1 ≤ t ≤ y. Let us note in passing, if particle composition in the initial and the final states are exactly identical (e.g. elastic scattering) then the collision term must be divided by an extra 2! since under the exchange of {initial state} ↔ {final state}, the process remains unchanged. Now, we will derive the Boltzmann equations for our dark sector species χ 1 , χ 2 , Z and h d . Let us first consider the species χ 1 and a particular annihilation channel χ 1 (P 1 )χ 1 (P 2 ) → Z (K 1 )Z (K 2 ). For this process, the Boltzmann equation for the number density of χ 1 follows from Eq. A.5 using O = 1 and O = n χ 1 (T ) as Here we consider χ 1 is a Majorana fermion. Using f Z (k β , T ) = n Z n eq Z Exp(−E k β /T ) and f χ 1 (p α , T ) = n χ 1 n eq χ 1 Exp(−E pα /T ) and also the conservation of energy in Eq. A.6 we get, Now, using the definition of σv for an annihilation process as given in [7], we can write the collision term of the above equation in a more compact form as with K 1 √ s T being the modified Bessel function of second kind and order one while s = (P 1 + P 2 ) 2 is one of the Mandelstam variables. The cross section σ χ 1 χ 1 →Z Z appearing above has the usual definition, and n eq χ 1 (T ) can be obtained from Eq. A.4 by using f Aκ = f eq χ 1 (E p 1 , T ) = Exp (−E p 1 /T ). In our model, χ 1 has two annihilation channels and hence the complete Boltzmann equation for χ 1 is Since both the dark matter candidates χ 1 and χ 2 have identical interactions, thus the Boltzmann equation for χ 2 is similar to that of χ 1 . Now, if they are degenerate in mass then their number densities will also be identical. In that case, it is needless to consider number densities of χ 1 and χ 2 separately and instead we will consider the total dark matter number density n χ = n χ 1 + n χ 2 with n χ 1 = n χ 2 . Therefore, in terms of total number density n χ , the Boltzmann equation for the dark matter can easily be obtained from Eq. A.11 as The other two dark sector species Z and h d have both decay and annihilation modes. Therefore, the collision terms of both Z and h d have contributions from annihilations and from decays as well. Therefore, following the similar procedure as we have discussed above, one can write the Boltzmann equations for Z and h d which are given below, where, Γ j T is the thermal averaged total decay width of the species j (j = Z , h d ) and it depends on dark sector temperature T . The expression of Γ j T in terms of total decay width Γ j is given by The last terms in Eqs. A.13 and A.14 represent the contribution from the inverse decay. In our model, both Z and h d have decay modes only into the SM particles which have temperature T . Hence, these inverse decay terms depend on the SM temperature instead of the dark sector temperature T like the others terms. The equilibrium number density n eq j (T ) of the species j at T can be found from Eq. A.4 using the Maxwell-Boltzmann distribution at temperature T .

A.1 Evolution of dark sector temperature (T )
Our next task is to derive the fourth Boltzmann equation describing evolution of the dark sector temperature T with respect to that of the SM. All the dark sector species are in kinetic equilibrium and they share a common temperature T before the freeze-out of dark matter candidates at T f . Therefore, in order to compute the temperature of the dark sector, we need to find the temperature of a particular species. Here we consider the temperature of our dark matter candidates. As we have seen in the previous section that for two identical dark matter candidates, we do not need to consider them separately in the Boltzmann equations. Similarly, in the temperature Boltzmann equation also, we will use the total dark matter density n χ rather than the individual ones. The temperature T of any species is defined as [90] T = g χ n χ (T ) where, m χ and g χ are the mass and internal degrees of freedom of dark matter respectively. The energy E p = p 2 + m 2 χ and n χ is the total number density of dark matter. Now considering O(p) = p 2 3E p , we can easily derive the temperature evolution equation from Eq. A.5 which is written below and The functions F(T ) 2→2 and F(T ) 3→2 in the right hand side of Eq. A.17 contain contributions from all 2 → 2 and 3 → 2 processes which can affect the temperature of the dark sector. The general form of these functions is shown in Eq. A.5. Their actual expression depends on the particular interaction process. We have presented a detail discussion about these functions in Appendix B.

B Calculation of F(T ) 2→2 and F(T ) 3→2
In this section, we have derived the specific form of the functions F(T ) 2→2 and F(T ) 3→2 encoding all the information about the relevant 2 → 2 and 3 → 2 processes which can change the temperature of χ i . Since all the species of the dark sector share a common temperature, we have considered only those interaction processes which have impact on the temperature of χ i .

B.1 Calculation of F(T ) 2→2
In our model, the 2 → 2 processes within the dark sector that can affect the evolution of T are Here we derive the function F 2→2 (T ) for a generic process like χ i χ i → jj where j can be either Z or h d . The master equation for writing the expression of F 2→2 (T ) can be found from the right hand side of Eq. A.5 as where, as defined earlier P α s are the four momentum for the initial state while that of the final state are denoted by K β s. The other quantities are already defined in the previous section. Eq. B.1 is symmetric under the exchange of P 1 ↔ P 2 and K 1 ↔ K 2 and using this exchange symmetry between P 1 ↔ P 2 one can further simplify as Now, following the similar procedure as we have done below Eq. A.6, a more compact form of the function F 2→2 (T ) is given below, Now, following the prescription given in [7], the expression of p 2 1 3E p 1 σv χ i χ i →jj T can be simplified further by using the new variables instead of E p 1 , E p 2 and θ. The transformation relations are given below where θ is the angle between p 1 and p 2 . The region of integration for these new variables are In terms of the new variable s the compact form of p 2 1 3E p 1 σv χ i χ i →jj T , after a few mathematical simplifications, can be obtained as This integral form of p 2 1 3E p 1 σv χ i χ i →jj T agrees well with the expression given in [36]. Finally, while incorporating this term in the Boltzmann equation (i.e. in F 2→2 (T )) we have to take contributions from both the dark matter candidates (sum over the index i) and replace individual densities n χ i by the total density n χ . Therefore, the final expression of F 2→2 (T )| χχ→jj for a annihilation process χχ → jj is given by

B.2 Calculation of F(T ) 3→2
In our case, there are several 3 → 2 processes which can change the dark sector temperature (T ) and we have considered all of them.
We first derive the contribution to F(T ) 3→2 from one such process namely χ i χ i h d → χ i χ i . Thereafter, we will present the final expressions of other 3 → 2 processes. The Feynman diagrams for the scattering χ i χ i h d → χ i χ i are shown in Fig. 16. The contribution to F(T ) 3→2 from χ i χ i h d → χ i χ i is given by As mentioned in [35], the quantity σv 2 a+b+c→d+e for a 3 → 2 scattering process a + b + c → d + e is defined as where |M| 2 a+b+c→d+e is the Lorentz invariant matrix amplitude square of a + b + c → d + e averaged over both initial and final states and the prefactor 1/m! is due to m number of identical particles in the final state. In the non-relativistic limit (i.e. E p 1 m a , E p 2 m b , E p 3 m c and p 1 = p 2 = p 3 0), Eq. B.8 can be further simplified as The expression of σv 2 has the following form in the non-relativistic and quasi degenerate limit σv 2 Here by quasi degenerate limit we want to mention that the above cross section has been calculated for m 1 = m 2 = M and m h d = m Z = M (1 − δ m ) with δ m 1. This kind of quasi degenerate dark sector is necessary for the co-decaying mechanism to work successfully [37].
Therefore, in the non-relativistic limit, the quantity F(T ) 3→2 | χ i χ i h d →χ i χ i can be written as where, 2 k 2 1 3E k 1 NR indicates the value of the quantity within the brackets in non-relativistic limit and it can be expressed as a function of masses of the initial state particles as Finally, the net contribution to F(T ) 3→2 from both the dark matter candidates χ 1 and χ 2 for the scattering process χ i χ i h d → χ i χ i is given by The Feynman diagrams for the scattering χ i Z Z → χ i Z are shown in Fig. 17 The contribution to F(T ) 3→2 is given by where, λ is the Kallen function which has the following definition and M in = 2 m Z + m i is the total mass of all the initial state particles.
The Feynman diagrams for the 3 → 2 scattering χ i Z Z → χ i h d are shown in Fig 18. The corresponding expression of σv 2 in non-relativistic and quasi degenerate limit is given by σv 2 This 3 → 2 scattering process has the following contribution to F(T ) 3→2 , with M in = 2 m Z + m i , the total mass of all the initial state particles.
The Feynman diagrams for the scattering χ i h d h d → χ i h d are shown in Fig. 19. The expression of σv 2  Figure 19. Feynman diagrams for This 3 → 2 scattering has the following contribution to F(T ) 3→2 , with M in = 2 m h d + m i is the total mass of all the initial state particles.
The Feynman diagrams for the scattering χ i Z h d → χ i Z are shown in Fig. 20 and the expression of σv 2 χ i Z h d →χ i Z in non-relativistic and quasi degenerate limit is given by The contribution to the function F(T ) 3→2 due to this inelastic scattering is given by where, M in = m i + m Z + m h d and λ is the Kallen function defined in Eq. B.16.
All the Feynman diagrams for this inelastic scattering is shown in Fig. 21 and the corresponding σv 2 is given by σv 2 Contribution to the function F(T ) 3→2 due to the 3 → 2 scattering χ i Z h d → χ i h d is given below Fig. 22. In the non-relativistic and quasi degenerate limit σv 2 for this scattering has the following expression  Figure 22.
All the relevant Feynman diagrams for the 3 → 2 scattering h d h d Z → χ i χ i are shown in Fig. 23.
In the non-relativistic and quasi degenerate limit, σv 2 for h d h d Z → χ i χ i can be expressed as In the function F 3→2 (T ), this particular 3 → 2 scattering has the following contribution Feynman diagrams contributing to the scattering Z Z h d → χ i χ i are depicted in Fig. 24 and corresponding σv 2 in non-relativistic and quasi degenerate limit is given below The expression of F 3→2 (T )| Z Z h d →χχ is given by where, m in = 2m Z + m h d .

B.2.10 Z Z Z → χ i χ i
Feynman diagrams for this inelastic scattering are depicted in Fig. 25. The expression of σv 2 for Z Z Z → χ i χ i is given by It is needless to mention here that σv 2 for this 3 → 2 scattering is also computed in non-relativistic and quasi degenerate limit. Z Z Z → χ i χ i inelastic scattering has following contribution to where λ is the Kallen function (Eq. B.16) and M in = 3m Z . Moreover, there is another 3 → 2 scattering χ i χ i Z → χ i χ i involving dark matter candidates χ i s however, these scatterings have null effect to the function F 3→2 (T ) in the non-relativistic limit.

C Vertex factors
The relevant vertex factors are listed in table 2.

Vertex
Vertex factors

D Approximate analytical form of Ω A
In this section, we have derived an approximate analytical expression for relic density of dark matter in co-decaying scenario. We have considered a general scenario where the dark sector has two degenerate species A and B as discussed in Section 3. The species A is our dark matter candidate while B can decay into the SM particles in out of equilibrium. As mentioned in Section 3, T d , T Γ and T f are the visible sector temperatures corresponding to decoupling of the dark sector from the SM, beginning of decay of the species B and the freeze-out of A respectively. The corresponding temperatures in the dark sector are denoted by T d , T Γ and T f with T d = T d . The relic density of A is defined as where, ρ c = 3H(T 0 ) 2 /(8πG N ) = 3.714×10 −47 GeV 4 is the critical density of the Universe and G N is the Newton's gravitational constant. In the second step, we have used the conservation of total number of A species per comoving volume after freeze-out to the present era (T 0 ) while the conservation of entropy per comoving volume between T f to T 0 has been utilised in the last step. We would like to note that all the dark sector thermodynamic variables are denoted with a prime while those for the visible sector have no prime. Using the freeze-out condition i.e. n A (T f ) σv AA→BB (T f ) H(T f ) in Eq. D.1 we get, f . In order to understand the parametric dependence of Ω A analytically, we have neglected the variation of degrees of freedoms (g ρ , g s ) with temperature. We have taken into account these effects while performing the numerical analyses shown in Section 3.3. Now, we need to know x f to calculate the relic density of A using Eq. D.2. For that, we require the expression of n A (T f ) which is related to n A (T Γ ) as This is the unique feature of the co-decaying scenario where only after the departure from chemical equilibrium an exponential suppression in number density, different from the Boltzmann suppression for a non-relativistic species, arises. In the last step we have used H(t) ∼ 1 2 t in the radiation dominated era and t f >> t Γ . Therefore, the number density n A (T f ) is given by The expression n A (T Γ ) can be found using the second law of thermodynamics (Eqs Here we have considered that the number densities and chemical potentials of species A and B are equal at T = T Γ . After T Γ , the species B has started to decay and thus the number densities of A and B deviates from each other. In the last step, we have used the fact that entropy per comoving volume is separately conserved in both the sectors between T d and T Γ . Moreover, ξ d is the ratio of degrees of freedom for the dark sector to the visible sector at the temperature T d . Now, substituting the expression of n A (T Γ ) in Eq. D.4 and considering the freeze-out for the species A once again, we get an equation for T f as .
(D. 6) In terms of x f , the above equation can be written as It is a transcendental equation of x f which can be solved numerically using iterative method and substituting the exact expression of σv AA→BB . However, an analytic solution for x f can be approximately given if we assume that σv AA→BB is independent (or does not depend significantly) of temperature. For example, the s-wave scattering at small T (see the plot in left panel of Fig. 3). Therefore, the approximate expression of x f is given by Here W 0 is the principal value of Lambert W function. Finally, substituting the expression of x f in Eq. D.2 we get an analytical expression of relic density (Ω A ) for the co-decaying scenario as E 2 → 2 scattering cross sections We have listed the expressions of relevant 2 → 2 scattering cross sections of the dark sector species below. The corresponding Feynman diagrams are shown in Fig. 1.

F Thermal average for degenerate initial and final state
In the non-relativistic limit, the matrix amplitude square for the process χ i (P 1 )χ i (P 2 ) → X(P 3 )X(P 4 ) (X = Z , h d ) can be written as, Where, the coefficients a i = a i (δ m , θ, m χ i ) and δ m = m χ i − m X m χ i , cos θ = p 1 . p 3 | p 1 || p 3 | . In the above we have used κ = s − 4m 2 << 1 (definition of κ is same as in [7]) for the non-relativistic regime.
Therefore the cross section of χ i χ i → XX is given by In the above, the dependence of κ in σv χ i χ i →XX is different for two scenarios. This is because the prefactor in σ i.e. s − 4m 2 X s − 4m 2 χ i ∝ κ −1/2 for δ m = 0 whereas it becomes independent of κ for δ m = 0.
Let us note that we have absorbed the phase space factors into a i s and define a new quantitŷ a i . Thus following the prescription given in [7], we have arrived at the following expression of the thermal average for χ i χ i → XX in the non-relativistic limit.
Now using x = m χ i /T 2/v 2 where v is the average DM thermal velocity, it is clearly seen that the thermal averaged annihilation cross section depends on either the even powers of v for δ m = 0 (corresponds to the away from threshold condition) or the odd powers of v for δ m = 0 (corresponds to the threshold condition).
In our model, the coefficientâ 0 for the process χ i χ i → Z Z vanishes when δ m = 0 (sincê a 0 ∝ δ 3/2 m ). Therefore, σv χ i χ i →Z Z is proportional to v 3 at threshold whereas it has a velocity independent leading order term away from threshold (δ m = 0). In contrast, the other annihilation channel of DM χ i χ i → h d h d has the coefficientâ 0 = 0 for all values of δ m , which implies that the thermal averaged annihilation cross section σv χ i χ i →h d h d is proportional to v 2 and v 3 for δ m = 0 and δ m = 0 respectively.