Imprints of MeV Scale Hidden Dark Sector at Planck Data

New light species can contribute to the number of effective relativistic degrees of freedom ($N_{\rm eff}$) at Cosmic Microwave Background (CMB) which is precisely measured by Planck. In this work, we consider an MeV scale thermally decoupled non-minimal dark sector and study the imprint of the dark sector dynamics on the measurement of $N_{\rm eff}$ at the time of CMB formation. We have predicted the allowed region of model parameter space in the light of constraints arising from the measurements of both $N_{\rm eff}$ and dark matter relic density by Planck. It turns out that the impact of the dark sector dynamics on $N_{\rm eff}$ is significant in case of a non-hierarchical mass spectrum of the dark sector particles.

New light species can contribute to the number of effective relativistic degrees of freedom (N eff ) at Cosmic Microwave Background (CMB) which is precisely measured by Planck. In this work, we consider an MeV scale thermally decoupled non-minimal dark sector and study the imprint of the dark sector dynamics on the measurement of N eff at the time of CMB formation. We have predicted the allowed region of model parameter space in the light of constraints arising from the measurements of both N eff and dark matter relic density by Planck. It turns out that the impact of the dark sector dynamics on N eff is significant in case of a non-hierarchical mass spectrum of the dark sector particles.

I. Introduction
The increasing tension due to non-detection of dark matter at several direct and indirect detection experiments [1][2][3][4][5] motivate to explore alternative scenarios beyond weakly interacting massive particle (WIMP) paradigm [6][7][8][9][10]. One of the attractive proposals is freeze in mechanism [11,12] where the DM abundance results from decays or annihilations of visible sector particles. In this scenario, the DM never reaches thermal equilibrium with the Standard Model (SM) bath due to its feeble interaction strength. Such DM candidates are popularly dubbed as feebly interacting massive particle (FIMP). Interestingly, the final abundance of the FIMP dark matter is sensitive to the initial production mechanism in contrast to WIMP. A different kind of scenario has been discussed in [13][14][15][16] where the dark matter is a part of a secluded dark sector and the annihilation of the dark matter into the other dark sector particles set the relic abundance. This is known as the secluded dark sector freeze-out.
The presence of any beyond Standard Model (BSM) light species can leave detectable imprints at Planck satellite experiment. The Planck data provides the number of relativistic degrees of freedom at the time of Cosmic Microwave Background (CMB) to be N CMB eff = 2.99 +0. 34 −0.33 with 95% confidence limit (C.L) [17]. In SM alone, an accurate analysis of neutrino decoupling by including exact collision terms, effect of neutrino oscillations, and finite temperature QED estimates N SM eff = 3.04 at the time of (CMB) [18][19][20][21]. Earlier efforts of constraining light BSM degrees of freedom can be found in [22][23][24][25][26][27][28][29][30]. These studies are based on the assumptions that SM neutrino decouples instantaneously and the new BSM particle is in thermal equilibrium either with the neutrino bath or with the photon bath. A more detailed analysis in [31][32][33][34] by incorporating non-instantaneous SM neutrino decoupling imposes lower bound on masses of BSM degrees of freedom from the measurement of N CMB eff . In addition to that, during Big Bang Nucleosynthesis (BBN), non-negligible abundances of the additional light species increase the Hubble parameter which in turn may alter the abundances of Helium-4 and Deuterium. The measured abundances of Helium-4 and Deuterium give rise to an upper bound on number of relativistic degrees of freedom at the time of BBN (N BBN eff ) to be 2.878 ± 0.278 with 68.3% C.L [35]. Thus light species with mass O(1)MeV can also contribute to the N BBN eff [24,[36][37][38]. The measurement of effective relativistic degrees of freedom can be used as a probe of light dark matter as well. A concise discussion on other possible modes of probing a light dark sector can be found in [39][40][41]. In [31,[42][43][44][45][46][47], it has been argued that a light dark matter can contribute sizeably to N CMB eff if it remains in thermal contact with either neutrino or electron bath till late time. This poses a lower bound on the MeV scale DM which differs depending on the spin of the dark matter candidate.
In this work, we consider an MeV scale non-minimal dark sector consisting of two real gauge singlet scalar fields. The lightest scalar is assumed to be neutrinophilic. The heavier scalar, being stable in the Universe lifetime can be identified as a viable dark matter candidate. Both the scalars interact very feebly with the visible sector and hence never attain thermal equilibrium with the SM bath, thus form a secluded dark sector. First, the non-thermal production of lighter neutrinophilic scalar takes place from the inverse decay of SM neutrinos and subsequently its annihilation yields the dark matter. The production rate of DM depends on two factors. One is the abundance of the mediator particle and secondly DM coupling strength to the mediator also determines the efficiency of DM yield. After the production, the dark matter annihilates back to the lighter scalar and in fact a stronger coupling leads to formation of internal dark sector equilibrium with a different temperature other than that of SM. Finally, when the ratio between dark sector interaction strength and the expansion rate of the Universe comes down below unity, the dark matter decouples and freezes out.
We anticipate that such a non-trivial dynamics of MeV scale dark sector can have two fold impacts on the observation of Planck measurement of N CMB eff . Firstly, if the lighter scalar decays to SM neutrinos during or after neutrino decoupling, it raises the neutrino bath temperature and in turn increases N CMB eff . Secondly, when mass hierarchy between dark matter and the lighter scalar is small, it may increase the N CMB eff further. In order to address such possibilities, we have performed a detailed study by solving system of coupled differential equations for the evolutions of number densities and the temperatures for the relevant species. We also include constraints on DM relic density in our present analysis. The bound on the measurement of ∆N CMB eff are expected to be more stringent with the next generation CMB experiments with improved sensitivity. In fact the proposed CMB Stage IV (CMB-S4) experiment [48] has the ability to strengthen the upper bound of ∆N CMB eff upto 0.06 at 95% C.L and can potentially probe our proposed setup.

II. Basic set up
In this section, we describe the particle contents of our proposed set up and their interaction pattern. We consider an extension of SM with two real gauge singlet scalar fields S and φ. The effective Lagrangian of our interest is given by, where we have considered the relevant mass scales (m S and m φ ) and the coupling coefficients (λ and κ) as real and positive. In our framework, the scalar field φ is purely neutrinophilic whereas S is the DM candidate. The stability of S can be ensured by imposing a Z 2 symmetry under which S has odd charge. The last term of Eq. 1 describes the interaction of the dark sector with the neutrino bath. Here the couplings of φ with the neutrinos are flavor universal and i represents the SM neutrino flavors. In our analysis, we work in the regime m φ < m S . In the opposite limit where the DM is lighter than the mediator particle, the dark matter S has negligible impact on N CMB eff . This is because in this case we need smaller values of the portal coupling to satisfy the relic density constraint, which in turn implies that the impact on N CMB eff is negligible. The effective vertex φνν is clearly not invariant under SU (2) L ⊗ U (1) Y . This can be originated from an SU (2) L ⊗ U (1) Y gauge invariant set up which is valid at higher energy scale. In appendix A, we have discussed two such realisations, guided by the introduction of a discrete symmetry Z 4 in both cases, under which φ transforms nontrivially. This choice automatically forbids the interaction terms involving φ and other SM particles at linear order in φ. Additionally, there exist other renormalisable terms such as φ 2 |Φ| 2 , S 2 |Φ| 2 in the Lagrangian, where Φ is the SM Higgs doublet. These may lead to new production channels for the dark sector (DS) particles. However, in the present work we are interested in a specific scenario where the DS gets populated dominantly from the SM neutrinos (for example in [49]). This can be realized by either choosing the Higgs portal coupling (φ 2 |Φ| 2 etc.) parameters to be tiny enough or assuming the reheating temperature of the Universe much lower than the Higgs boson mass scale [50,51], making the Higgs abundance Boltzmann suppressed in the early Universe. In both the cases, the DS does not equilibrate with the SM bath through Higgs portal.
In the present study, we aim to probe the secluded non-minimal dark sector (DS) with field contents S and φ at Planck experiment. The dark sector is secluded when S and φ never reach thermal equilibrium with the SM bath.
The non-thermalisation of dark sector with SM bath implies, where σv νiνi→φ is the thermally averaged cross section for the inverse decay of φ and H represents the Hubble parameter of the Universe. A conservative check for Eq. 2 can be performed at the freeze in temperature T ν ∼ m φ of φ yield. Then using the standard analytical expressions for σv νiνi→φ , n eq νi and Hubble parameter, the Eq. 2 can be converted to, where y 2 = i κ 2 . We have confirmed numerically that Eq. 3 is a reasonable check to keep the DS and SM neutrino out of thermal equilibrium for any temperature of the SM neutrino bath. We start with the assumption that both φ and S have zero initial abundance. The inverse decay of SM neutrino yields φ and conversion of φ takes place through φφ → SS process. At a later stage when number density of S is sufficient, the reverse process SS → φφ turns effective (see Fig. 1). If the rate of such forward and reverse interactions in dark sector are fast enough, a dark equilibrium can be formed with a new temperature T D = T S = T γ , T ν . Assuming the dark sector is internally thermalised, we have estimated the dark sector temperature T D by solving the Boltzmann equation of total dark sector energy density ρ D = ρ S + ρ φ . The equilibration of the dark sector can be ensured by the following condition.
where we have considered σv SS→φφ to be s-wave dominated. Simultaneous requirements of dark sector equilibrium and non-thermalisation of dark and visible sectors put constraints on the model parameters and the consequences of these constraints will be discussed in section III.
We are particularly interested in the dynamics of a MeV scale secluded dark sector in presence of dark S − φ equilibrium. The impact of such MeV scale dark sector on Planck can be parameterized by the effective number of relativistic species at CMB as defined by, where ν are the energy densities of photon and neutrino baths respectively. Before neutrino decoupling, the temperatures of the neutrino bath and the photon bath are equal. In the post neutrino decoupling era, T γ and T ν evolve separately. If all the BSM particles are having mass much larger than eV scale, they do not contribute to the energy budget during CMB and hence in Eq. 5, one can write, Similarly one can define the number of relativistic degrees of freedom at the time of BBN and it is given by ( considering all BSM particles are non-relativistic during BBN)

III. Boltzmann equations and Numerical analysis
To study the effect of the dark sector dynamics on ∆N eff , we need to solve the Boltzmann equations for comoving number densities Y i (i = φ, S), ξ = T ν /T γ , and ξ D = T D /T γ . In deriving the above Boltzmann equations, we have used the following assumptions.
• The temperatures of all the neutrino flavors are same i.e. T νe = T νµ = T ντ ≡ T ν .
• The entropy of the photon bath is conserved since at late time φ can only decay into a pair of neutrinos. Using the conservation of entropy in the photon bath, we have defined Y i = n i /s γ and the relic abundance of the DM is calculated accordingly 1 .
• The interaction rate between the dark sector particles such as S and φ is sufficient enough to keep them in thermal equilibrium and they share a common temperature T D (see discussion above Eq. 4).
Using the above assumptions, first we write the evolution of ξ as a function of x ≡ m S /T γ . where Here, C el (T γ , T ν ) and C inel (T γ , T ν ) are the relevant collision terms for the elastic and inelastic processes such as The matrix amplitude square of these processes and the generic form of collision terms for elastic, inelastic, and decay processes are given in decay processes are given in table III of appendix B and appendix C respectively. The energy densities of electron, photon, and neutrino are respectively denoted by ρ e , ρ γ , and ρ ν . p e denotes the pressure of electron. The evolution of the co-moving number densities of φ (Y φ ), and S (Y S ) are governed by the following two coupled Boltzmann equations.
Here, s γ is the entropy density of the photon bath, and Y eq φ(S) (T D ) denotes the equilibrium co-moving number density of φ(S) at hidden sector temperature T D . The Hubble parameter H and h eff are defined as follows.
where M Pl = 1.22 × 10 19 GeV and g * s,γ (T γ ) is the relativistic degrees of freedom contributing to the entropy density of the photon bath. The energy densities of photon and neutrino baths are indicated by ρ γ , ρ ν respectively. Note that, in the expression of the Hubble parameter, we have included the energy densities of dark sector particles in addition to the contributions from SM fields. On the right hand side (RHS) of Eq. 10 and Eq. 11, σv rel for the thermally averaged annihilation cross section of SS → φφ at temperature T D . The analytical form of σ SS→φφ is given by, where √ s is the centre of mass energy of the collision. Γ φ T D and Γ φ Tν on the RHS of Eq. 10 are the thermally averaged total decay width of φ at temperature T D and T ν respectively. A general expression of the thermally averaged decay width at temperature T i is given by 16π is the total decay width of φ. In deriving the Boltzmann equations for Y φ and Y χ , we have also taken into account the effect of non-zero chemical potential for dark sector particles as evident from the presence of square of the ratio (Y eq S (T D )/Y eq φ (T D )) in the RHS of Eq. 10.
As we discussed earlier, the dark sector is connected to SM by the portal coupling y and our choice of y keeps dark and the visible sector out of thermal equilibrium. Initially inverse decay of φ populates the dark sector and subsequently φφ → SS process occurs. The Boltzmann equation that governs the evolution of total dark sector energy density (ρ D = ρ φ + ρ S ) is read as where p D stands for the pressure and C νν↔φ includes the collision integrals for both decay and inverse decay processes of φ. We use the assumption of Maxwell-Boltzmann statistics and write ρ D and p D as function of T D as given by In the RHS of the above equation, the quantities, ρ eq φ and p eq φ represent the energy density and pressure at equilibrium as given by where once again, the ratio in Eqs.15-16 takes care of the non-zero chemical potential during the temperature evolution for the i th species. Utilizing these, we obtain, where,  In Fig. 2, we demonstrate the evolution patterns of T D , Y φ , Y S , and ∆N CMB eff as functions of the SM temperature T γ . We have fixed m φ and y at 5MeV and 2 × 10 −10 respectively and considered two specific DM masses which are 5.1MeV and 12MeV. The corresponding λ value for each of the DM mass is fixed by the relic density requirement (Ω S h 2 ∼ 0.12). The dark sector temperature starts from zero initial value and increases at the early stage of its evolution due to the production of φ fromν i ν i → φ process. The temperature T D continues to grow until it reaches a maximum value, after which it starts to redshift like radiation upto T D m φ due to the expansion of the Universe. Around T γ ∼ m φ , we notice that the redshift for T D gets slower. This occurs since the production rate of νν → φ The symbol ' ' stands for the λ value that yields correct relic abundance [17]. In all three figures, the grey region is disallowed from the measurement of ∆N CMB eff by Planck within 2σ limit. reaches maximum around this temperature. For m S = 5.1MeV we also observe the impact of dark matter freeze out in the evolution of T D , where late time conversion process SS → φφ slows down the redshift of T D for the second time. This pattern is not visible for m S = 12MeV as in this case the dark matter freeze out occurs at a higher temperature when the νν → φ process is still active. Finally after the freeze out of the dark matter and decay of φ, the dark sector temperature redishifts like a non-relativistic matter field with usual scale factor (a) dependence T D ∝ a −2 .
In Fig. 2(b), we show the evolutions of Y S (solid lines) and Y φ (dashed lines) as functions of T γ considering two different choices of (m S , λ) as indicated by black and blue curves for m S = 5.1MeV and 12MeV, respectively. The corresponding values of λ are chosen accordingly to yield correct relic abundance with a larger DM mass requiring lower value of λ to obey the relic density bound. On the other hand, in Fig. 2(c), we show the variation of ∆N CMB eff as a function of T γ for the same choice of parameters as considered in Fig. 2(b). It is observed from Fig. 2 that the dark matter S of mass 12MeV freezes-out much earlier and at the time of DM freeze-out, the production of φ is still in progress. As a result, when φ starts decaying to a pair of SM neutrinos, the DM has already decoupled from φ bath and the final value of ∆N eff is governed mostly by m φ and y. However the situation turns a bit different as we decrease the mass hierarchy of S and φ. In case of m S = 5.1MeV, decoupling of S is delayed in comparison to the earlier scenario and at the time of decay of φ, SS → φφ is still active with S having equivalent comoving abundance of φ around T D ∼ m φ . As a result, the depletion of φ has been counterbalanced by the production SS → φφ, which in turn enhances ∆N CMB eff . Thus, ∆N CMB eff increases as we reduce the mass gap between S and φ. Next, the dependence of the parameter λ on ∆N CMB eff is shown in Fig. 3 for y = 2 × 10 −10 and three different values of m φ . The lowest chosen value of m φ is 3MeV in Fig. 3 and for this choice, φ remains out of equilibrium from SM neutrino bath as can be confirmed from Eq. 3. For larger m φ with same y, the Eq. 3 is easier to satisfy. We consider two different hierarchical patterns of dark sector particles which are (m S − m φ )= 0.5MeV and 2MeV. Let us define δ = (m S − m φ ) for convenience. The dark sector reaches internal equilibrium for the chosen ranges of λ in Figs. 3(a)-3(c) and it has been numerically checked. We find that ∆N CMB eff is insensitive to λ. This is because when the dark sector is internally thermalized, the abundance of φ at the onset of its decay is mostly governed by y, m φ and m S . We also notice ∆N eff increases with the decrease in mass gap between S and φ as earlier observed in Fig. 2(c). In all three subfigures of Fig. 3, the point highlighted by ' ' indicates the required λ to yield correct relic density for DM [17]. Another important outcome of Fig. 3 is related to the dependence of m φ on ∆N CMB eff . When m φ is fixed at 3MeV, the predictions for ∆N CMB eff considering both the mass differences δ = 0.5MeV and 2MeV are below the Planck limit. With the increase of m φ ∼ 7MeV with δ = 0.5MeV, the prediction for ∆N CMB eff enters into the disallowed region. This happens since larger m φ increase the decay width of φ and causes enhanced rate of entropy injection into the SM neutrino bath. Further increase of m φ to 10MeV shows different pattern and δ = 0.5MeV again makes a comeback to the allowed region of ∆N CMB eff . This is due to the fact that a sufficient large m φ would make φ to decay early with less impact on evolution of T ν . On the other hand for δ = 2MeV, the predictions for ∆N CMB eff always remain inside the Planck favored region.
As we find from Fig. 2(c), the predictions on ∆N CMB eff for m φ = 5MeV and m S = 5.1MeV are disfavored when y = 2 × 10 −10 . Now, let us examine the effects of varying y for the same set of m φ and m S . In Fig. 4, we show the variation of ∆N CMB eff as a function of T γ for three different choices of the portal coupling y by fixing m φ and m S at 5MeV and 5.1MeV respectively. The choices for λ are made to obtain the best fit value of observed relic abundance Ωh 2 = 0.1200 ± 0.0012 [17] after formation of internal dark equilibrium. Earlier we found that y = 2 × 10 −10 is ruled out by present CMB data which is observed once again in Fig. 4. However reducing y suppresses the energy injection rate and thus decreases ∆N CMB eff and hence would be allowed by the Planck data. In Fig. 5, we have shown the evolutions of N BBN eff following the definition of Eq. 6 considering three different values of m φ with the m S is fixed by m S = m φ + 0.5MeV and y = 2 × 10 −10 . The magnitude of λ is determined by the requirement of obeying the relic density bound. It turns out that N BBN eff at T BBN ∼ 1MeV increases with m φ however stays inside the allowed 2σ limit [35]. Importantly, for lighter m φ , N BBN eff is smaller than its SM value which is approximately three. This owes to the fact that for the lightest m φ = 3MeV, as considered here, the production of φ continues even during BBN from the inverse decay of φ. This reduces the temperature of the SM neutrino bath at BBN and we observe a dip in the value of N BBN eff at T BBN ∼ 1MeV for m φ = 3MeV. Overall, we find that BBN requirements do not impose any serious constraints to our model parameter space. A detailed analysis in order to study the impact of a secluded non-minimal dark sector on the abundances of the primordial light elements (e.g. Helium-4 and Deuterium [55,56]) is beyond the scope of the present article and will be pursued elsewhere.
Finally in Fig. 6, we show the predictions for ∆N CMB eff in m S − λ plane. We have fixed m φ and y at 5MeV and 2 × 10 −10 respectively. These particular choices ensure that φ does not equilibrate with SM neutrinos as followed from Eq. 3. The overabundant region of DM is highlighted with light grey color. We have numerically found that the condition for dark sector equilibrium following Eq. 4 with m S 20MeV poses much weaker constraint compared to the one from observed DM relic in the m S −λ plane. Hence we do not show it in Fig. 6. The ∆N CMB eff lines for different 2σ lower limit of N eff at BBN mφ= 3MeV, mS= 3.5MeV, λ= 3.49× 10 -6 mφ= 7MeV, mS= 7.5MeV, λ= 6.58× 10 -6 mφ= 10MeV, mS= 10.5MeV, λ= 9. with Tγ is shown for three different sets of (m φ , ms) with y = 2 × 10 −10 . The corresponding λ value for each set is fixed by the relic density requirement. The magenta and blue horizontal solid lines denote the upper and lower limit of ∆N eff at the time of BBN at 2σ [35] respectively. The perpendicular dot-dashed line represents TBBN = 1MeV. m S are also depicted in Fig. 6. As explained earlier the predictions for ∆N CMB eff remain insensitive to variation in λ. Clearly, one can see that the value of ∆N CMB eff is maximum when m S turns closer to m φ . We also observed that as m S turns larger, its impact on ∆N CMB eff gets reduced. The latest Planck 2σ bound on ∆N CMB eff restricts m S to be larger than around 7 MeV for m φ = 5MeV and y = 2 × 10 −10 .

IV. Summary and Conclusion
In this work, we examine the scope of probing a MeV scale secluded dark sector from the observation of ∆N CMB eff by Planck. For the purpose, we have considered a simple structure of secluded dark sector comprising of two SM gauge singlet scalars S and φ. The dark sector never reaches thermal equilibrium with the SM sector. The lighter scalar φ is neutrinophilic and we identify the heavier one S as our DM candidate. Initially the dark sector particle φ gets populated from the SM sector by the ν i ν i → φ process and afterwards φφ → SS process produces S. Depending on the interaction strength between S and φ, dark thermalisation can occur which causes freeze-out of S via SS → φφ process. We constrain our parameter space from the observations of both dark matter relic density and measurement of N CMB eff by Planck ensuring internal dark thermal equilibrium. We particularly emphasize the impact of the dark sector parameters on the prediction of N CMB eff . Below we summarize few important observations that came out from our analysis.
• We obtain unique predictions for the ∆N CMB eff for a particular DM mass that satisfies the relic density bound with suitable value of λ, provided mass of the other scalar and its coupling strength with the SM neutrinos are fixed.
• The dark matter mass has a nontrivial role on the prediction of ∆N CMB eff . The predictions for ∆N CMB eff turn out to be maximum in the nearly degenerate spectrum of dark sector particles and consideration of an increased hierarchy between the dark sector particles reduces the impact of dark matter mass on ∆N CMB eff . We also find that prediction for ∆N CMB eff is insensitive to the dark sector interaction rate, provided the dark thermal equilibrium is reached.
• Increasing m φ upto a certain value keeping other parameters (y, λ and δ) fixed, enhances ∆N CMB eff hence restricted by the present Planck limit. When m φ is too large, the effect on ∆N CMB eff of φ diminishes.
• For a particular m φ with a constant y, we are able to impose a lower bound on the DM mass from the Planck 2σ data on ∆N CMB eff . As an example we find m S 7 MeV when m φ = 5MeV and y = 2 × 10 −10 . Additionally, a lower bound on λ can be derived such that the DM relic abundance remain below 0.12 as experimentally favored.
In summary, an MeV scale secluded non-minimal dark sector is difficult to probe at dedicated dark matter search experiments. In view of this, we propose that the measurement of ∆N CMB eff can be used as a tool to test such a decoupled dark sector. Our work indeed shows that some of the parameter space yield observable ∆N CMB eff and are already within the reach of the Planck sensitivity. Upcoming CMB stage IV experiments with improved sensitivity should be able to probe/refute the allowed region of our model parameter space further.

V. Acknowledgements
SG would like to thank Anirban Biswas for some insightful discussions. SG would also like to thank the University Grants Commission (UGC), Government of India, for providing financial support as a senior research fellow. AKS is supported by NPDF grant PDF/2020/000797 from Science and Engineering Research Board (SERB), Government of India. In this appendix, we discuss two SU (2) L ⊗ U (1) Y invariant frameworks that give rise to φνν vertex effectively at MeV scale.
Case I: We consider a variant of type-I seesaw framework with three gauge singlet right handed (RH) neutrinos. We impose Z 4 discrete symmetry to write the desired Lagrangian. The charge assignments of the relevant fields under the new Z 4 are presented in table I. We propose the following Lagrangian for generating the φνν vertex.
where c φ1 , c φ2 are dimensionless O(1) coupling coefficients and Φ = iσ 2 Φ * . After electroweak symmetry breaking, the SM neutrino mass matrix is given by, where The active-sterile neutrino mixing is θ αi , quantified as θ αi m D M −1 R αi [57]. Due to this active-sterile mixing, the first two terms in Eq. A1 can be translated to φνν vertex with the coupling coefficient being proportional to ∼ (θ 2 µ1 + θ 2 e2 ). This set up is able to provide correct order of active neutrino mass (∼ 0.1eV) for M 33 , M 12 ∼ O(1)GeV and Y ∼ O(10 −7 ) with the desired value of φνν coupling coefficient ∼ O(10 −10 ). A detailed study on neutrino masses and mixing is beyond the scope of the present paper. To prevent the production of DS particles from tree level RH neutrino decay (N → φν), we have assumed that the reheating temperature (T RH ) of the Universe is less than 1 GeV which is consistent with the lower limit on the T RH ( 4MeV) from big bang nucleosynthesis [50].
Case II: We may also start with the following SU (2) L ⊗ U (1) Y invariant Lagrangian at dimension six level to effectively generate the φνν vertex [58], where c 2 is the dimensionless O(1) coupling coefficient, Λ is the cut off scale andL l = iσ 2 L C l with l = {e, µ.τ }. In  this case also, we consider an additional discrete symmetry Z 4 and tabulate the charge assignments of the relevant fields in table II. After electroweak symmetry breaking one obtains the vertex φνν with the coupling coefficient being proportional to c2v 2 Λ 2 where v = 246GeV.

C. Collision Terms
• The collision term for the evolution of the energy density of A for a generic inelastic process such as A(p 1 )+B(p 2 ) → C(p 3 ) + D(p 4 ) is given by where Here the distribution function of the i th species is denoted by f i . p 1 , p 2 , p 3 , p 4 are four momenta of A, B, C, D respectively. |M 2 | AB→CD is the matrix amplitude square of AB → CD process and dΠ i = d 3 p i 2E i (2π) 3 .
• The collision term for the evolution of the energy density of A for a generic elastic process such as A(p 1 ) + B(p 2 ) → A(p 3 ) + B(p 4 ) is given by where Here |M 2 | AB→AB is the matrix amplitude square of AB → AB elastic process.
• For A(p 1 ) → B(p 2 ) + C(p 3 ) process, the collision term for the evolution of the energy density of A is as follows: where |M 2 | A→BC is the matrix amplitude square for A → BC process and