Constraining Non-thermal and Thermal properties of Dark Matter

We describe the evolution of Dark Matter (DM) abundance from the very onset of its creation from inflaton decay under the assumption of an instantaneous reheating. Based on the initial conditions such as the inflaton mass and its decay branching ratio to DM, reheating temperature, and the DM mass and interaction rate with the thermal bath, the DM particles can either thermalize (fully/partially) with the primordial bath or remain non-thermal throughout their evolution history. In the thermal case, the final abundance is set by the standard freeze-out mechanism for large annihilation rates, irrespective of the initial conditions. For smaller annihilation rates, it can be set by the freeze-in mechanism, also independent of the initial abundance, provided it is small to begin with. For even smaller interaction rates, the DM decouples while being non-thermal, and the relic abundance will be essentially set by the initial conditions. We put model-independent constraints on the DM mass and annihilation rate from over-abundance by exactly solving the relevant Boltzmann equations, and identify the thermal freeze-out, freeze-in and non-thermal regions of the allowed parameter space. We highlight a generic fact that inflaton decay to DM inevitably leads to an overclosure of the Universe for a large range of DM parameter space, and thus poses a stringent constraint that must be taken into account while constructing models of DM. For the thermal DM region, we also show the complementary constraints from indirect DM search experiments, Big Bang Nucleosynthesis, Cosmic Microwave Background, Planck measurements, and theoretical limits due to the unitarity of S-matrix. For the non-thermal DM scenario, we show the allowed parameter space in terms of the inflaton and DM masses for a given reheating temperature, and compute the comoving free-streaming length to identify the hot, warm and cold DM regimes.


INTRODUCTION
There is overwhelming astrophysical and cosmological evidence for the existence of Dark Matter (DM) in our Universe (for a review, see 1). Assuming the standard CDM (cosmological constant+Cold Dark Matter) picture of the Universe, the recent measurements from the Planck mission yield the current matter density in the Universe to be 4.9% in the form of baryonic matter and 26.6% as non-baryonic, non-luminous DM, while the remaining 68.5% is in the form of Dark Energy [2]. Despite all the compelling evidence from its gravitational interaction, the origin and nature of DM are still unknown, and resolving these issues is one of the main goals of modern cosmology as well as particle physics.
On the other hand, cosmological observations such as Planck [3] are strongly pointing toward an epoch of primordial inflation (for a review, see [4]), which is considered to be one of the best paradigms to create the seed perturbations for the DM particles to form the observed large-scale structures [5]. Inflation not only stretches the primordial perturbations on large scales but also dilutes all matter, and therefore, it is important that the inflaton must excite the Standard Model (SM) degrees of freedom (d.o.f) after the end of inflation for the success of Big Bang Nucleosynthesis (BBN) [6].
Irrespective of the origin of the inflaton field, whose potential leads to inflation, the inflaton could decay into the SM d.o.f, and also directly to the DM particles. The process of creating the entropy happens after the end of inflation during reheating or preheating (for reviews, see 4,7). If the inflaton is a SM gaugesinglet field φ, it can in principle couple and decay to DM particles χ with some unknown branching ratio 1 . As a concrete example, one can envisage a visible-sector inflation scenario within the Minimal Supersymmetric Standard Model (MSSM), where inflation can be driven by the superpartners of quarks and leptons [9][10][11]. Thus, if the DM can have a significant coupling to the inflaton or the moduli field, it can be created rather efficiently from the direct decay of the inflaton or the moduli [12,13] 2 . The initial abundance of such DM particles could be large enough to overclose the Universe, unless their interaction rate is sufficiently larger than the Hubble expansion rate in the early Universe, i.e., χ H(t), to make them annihilate efficiently into the SM d.o.f. The requirement of not to overproduce the DM poses a stringent constraint on its mass and interaction rates. In fact, a small branching fraction of the inflaton energy density to DM particles can be sufficient to overclose the Universe [18].
In this paper, we seek a model-independent way to analyze the thermal and non-thermal properties of the DM directly produced from the inflaton decay, in terms of their masses, the initial inflaton branching ratio and the strength of DM coupling to the thermal bath. For this purpose, we will make the following minimal assumptions: 1. The inflaton decays into the SM d.o.f and the DM in a perturbative scenario; hence on kinematic grounds, m χ < m φ /2. We do not consider non-perturbative DM production processes during the coherent oscillations of the inflaton, e.g., superheavy DM with m χ m φ for large enough amplitude of the inflaton field [8], or fragmentation of the inflaton condensate associated with a global symmetry [19,20]. 2. The process of reheating is instantaneous, and the SM d.o.f produced in the inflaton decay quickly thermalize to achieve a local thermodynamic equilibrium (LTE), i.e., both chemical and kinetic equilibrium. Thus we can define a unique reheat temperature T R [21] at which the Universe is dominated by relativistic species. Since we assume the DM to be part of the decay products of the inflaton, its initial number density (n χ ) will be determined in terms of the number density of the inflaton field (n φ ) and the branching ratio of the inflaton decay to DM (B χ ) which should be small in order to have the standard radiation-domination epoch immediately after reheating.
We should note here that the analysis presented in this paper based on the simplified picture mentioned above may not be completely valid if the reheating process is not instantaneous and a significant amount of DM is produced during reheating through inelastic scatterings between high energy inflaton decay products and the thermal plasma. However, a proper treatment of this issue must include the details of the thermalization process which involve some model-dependent subtleties. In particular, it is important to know when the inflaton decay products thermalize during the inflaton oscillations. This depends on how the DM interacts with the ambient thermal plasma, such as the SM d.o.f.
These thermal interactions lead to some finite temperature effects on DM production rate which will be model-dependent and has to be studied separately. In addition, there exist certain physical situations where the thermalization of the inflaton decay products can happen very late, even beyond the complete decay of the inflaton. This issue has to be dealt again separately, because the DM can still be created before the epoch of thermalization, but instead of thermal corrections, there will be finite momentum effects due to the hard-hard, hard-soft and soft-soft scatterings between the inflaton decay products, all of which must be accounted for in a field-theoretically consistent manner. Besides these issues, there could also be an epoch of preheating during which the DM can be created abundantly, but again this has to be studied in a modeldependent setup. Our main goal in this paper is to study the DM abundance from the inflaton decay in a model-independent perspective, and in this regard, we make the simple assumption that the thermalization has happened instantly right at the time of reheating. A more exhaustive analysis of DM creation from the inflaton decay, addressing all the subtleties mentioned above, will be presented elsewhere.
With the assumptions made above, the evolution of the DM particle number density, can be completely described in a modelindependent way by its thermally averaged interaction rate σ v with the thermal bath. Depending on the size of σ v , we consider the following three possible scenarios: 3 1. For large enough σ v , the DM particles will quickly reach LTE with the bath, thus losing their initial abundance, and follow the equilibrium distribution until their reaction rate eventually drops below the Hubble rate, after which they will freeze out as a 'thermal relic' with a constant comoving number density. This is the standard WIMP scenario [26] in which the final relic abundance is independent of the initial conditions or the details of the production mechanism. Depending on their mass and interaction rate, they could freeze out as a cold, warm, or hot relic [26]. It is well-known that σ v ∼ 10 −26 cm 3 s −1 naturally gives the observed cold DM relic density [2], almost independent of the DM mass. 2. If the interaction of DM particles with the thermal bath is too small to bring them into full LTE, they will decouple from the bath soon after being produced. Hence, if they are produced abundantly, the final number density will remain large, thus leading to overclosure of the Universe. However, if their initial abundance is negligibly small, the interactions with the bath, although feeble, could still produce some DM particles whose 3 Again a concrete example is MSSM inflation in which case the lightest supersymmetric particle, e.g., gravitino or neutralino, could be excited directly from the inflaton decay or its decay products [22], besides the SM d.o.f. Since gravitinos mostly interact via Planck-suppressed interactions, their abundance will freeze out soon after their production and will be mainly determined by the reheat temperature, while neutralinos have weak interactions and can be quickly brought into kinetic equilibrium (though not necessarily chemical equilibrium) with the bath. Therefore irrespective of how the neutralinos were initially created, their final abundance is always set by the thermal decoupling temperature, as long as T R ≥ m χ [23]. On the other other, for low reheating temperatures below the standard freeze-out temperature T F ∼ m χ /20, neutralinos could be a non-thermal DM candidate [24,25].
final abundance freezes in at some point as the interaction rate eventually becomes smaller than the expansion rate. This is the FIMP (Feebly Interacting Massive Particle or Frozen-In Massive Particle) scenario [27]. 3. For extremely small annihilation rates, the DM particles are never in thermal contact with the bath, and are practically produced decoupled in an out-of-equilibrium condition, and remain non-thermal throughout their evolution. This leads to a super-WIMP (SWIMP)-like scenario [28], where the final abundance is primarily determined by the initial conditions which, in our case, are set by the inflaton mass, reheat temperature and branching ratio [29,30] 4 . Note that it is also possible to have a non-thermal DM with chemical equilibrium, provided the reheat temperature is smaller than the usual freeze-out temperature so that the DM decouples during reheating [25]. We do not consider this case here since we have assumed instant reheating.
For each of the above three scenarios, we study the evolution of the DM number density by numerically solving the Boltzmann equation, and obtain their final relic abundance as a function of their mass and interaction rate for cold, warm as well as hot relics. We highlight a generic fact that inflaton decay to DM inevitably leads to an overclosure of the Universe for a large range of parameter space, and provides a generic constraint for models of DM with an arbitrary coupling to the inflaton field. For a given reheat temperature and initial abundance, we show the overclosure region as a function of the DM mass and annihilation rate in a model-independent way. For the thermal WIMP case, we show the complementary constraints on the (m χ , σ v ) parameter space by taking into account various theoretical as well as experimental limits from unitarity, dark radiation, indirect detection, BBN, Cosmic Microwave Background (CMB) and Planck. For the non-thermal production of DM from inflaton decay, we show that a large fraction of the (m φ , m χ ) parameter space leads to an overclosure for a generic class of hidden sector models of inflation. This an important result in pinning down the nature of DM from particle physics point of view and on the allowed region of the inflaton-DM coupling and the branching ratio. The rest of the paper is organized as follows: In section 2, we briefly review the evolution of DM as governed by the Boltzmann equation. In section 3, we discuss the production of DM from inflaton decay: thermal production (both freeze-out and freeze-in scenarios) in section 3.1, and non-thermal production in section 3.2. In section 4, we discuss various experimental/observational constraints on DM. In section 5, we present our numerical results for both thermal and non-thermal scenarios. Our conclusions are given in section 6.

EVOLUTION OF DM: A BRIEF REVIEW
The microscopic evolution of the number density n χ for any species χ , and the departure from its thermal equilibrium value n χ,eq , can be computed exactly by solving the Boltzmann Equation [26] dn χ dt where σ v is the thermally averaged total annihilation rate, σ being the total (unpolarized) annihilation cross section, and v being the relativistic relative velocity between the two annihilating particles 5 . In the absence of Bose-Einstein condensation or Fermi degeneracy, one can neglect the quantum statistical factors, and write the number density as where g χ is the number of internal (e.g., spin or color) degrees of freedom of χ , T is the temperature, and μ χ is the chemical potential of species χ (energy associated with change in particle number) which we assume to be zero for the equilibrium number density n χ,eq . It is useful to express Equation (1) in terms of the dimensionless quantities Y χ = n χ /s and Y χ,eq = n χ,eq /s to scale out the redshift effect due to the expansion of the Universe. Here, is the entropy density and g s is the effective number of relativistic degrees of freedom contributing to the total entropy density.
Recall that in the early Universe with radiation domination, g s is same as the relativistic degrees of freedom g ρ contributing to the energy density, and also appearing in the Hubble expansion rate: where m Pl = 1.22 × 10 19 GeV is the Planck mass. Henceforth, we will not distinguish the two, and will take g ρ = g s ≡ g which is valid for most of the thermal history of the Universe 6 . Assuming an adiabatic and isentropic (constant entropy per comoving volume) expansion of the Universe, Equation (1) can be rewritten as Srednicki al. [33] and Gondolo and Gelmini [34] dY 5 For the non-relativistic case, v is approximated by the relative velocity v r = |v 1 − v 2 |, while in the general case, it is usually taken to be the Møller velocitȳ Here we use the manifestly Lorentz-invariant , where p 1 , p 2 are the four-momenta of the annihilating particles with mass m 1 and m 2 respectively [32]. 6 g s and g ρ differ only when there are relativistic species not in equilibrium with photons which happens in the SM for temperatures below the electron mass when the neutrinos have already decoupled from the thermal bath, and e ± pair-annihilation transfers entropy only to the photons, thus making g s slightly higher than g ρ today.
with the introduction of a new independent variable x = m χ /T. The current number density Y χ (x 0 ) of the species χ is obtained by integrating Equation (5) from x = 0 to x = x 0 ≡ m χ /T 0 , where T 0 = 2.7255(6) K is the present temperature of the CMB photons [35]. Knowing Y χ (x 0 ), we can compute the relic density of χ , conventionally defined as the ratio of its current mass density, ρ χ (x 0 ) = m χ s 0 Y χ (x 0 ), and the critical density of the Universe, ρ c = 3H 2 0 /8π . Using the current values for the entropy density s 0 = 2889.2 cm −3 (T 0 /2.725 K) 3 , and the critical mass density ρ c = 1.05375(13) × 10 −5 h 2 GeVcm −3 [36], we obtain Equation (5) is a form of the Riccati equation for which there is no general, closed-form analytic solution. Therefore, the current density Y χ (x 0 ) in Equation (6) has to be obtained either by numerically solving Equation (5) or by approximating it with an analytic solution in some special cases. In the standard CDM cosmology, the thermal relics decouple from the thermal plasma in the radiation-dominated era after inflation, and the decoupling occurs at some freeze-out temperature T F when the annihilation rate χ = n χ σ v drops below the Hubble expansion rate H. Depending on the exact value of x F = m χ /T F , one can have the following three scenarios:

NON-RELATIVISTIC CASE
For x F > ∼ 3, the DM particles are mostly non-relativistic when they decouple from the thermal plasma. This leads to the usual cold DM scenario with free streaming lengths of sub-pc scale [1], as favored by the standard theory of large-scale structure formation [37,38]. Analytic approximate formulas for their relic abundance have been derived in the non-relativistic limit x F 1 [34,[39][40][41][42]. The key point is that the actual abundance Y χ tracks the equilibrium abundance Y χ,eq during early stages of evolution (for x < ∼ x * ), while at late stages (x > ∼ x * ), Y χ,eq is exponentially suppressed and has essentially no effect on the final abundance Y χ (x 0 ). Here x * is some intermediate matching point (not the freeze-out point x F , as commonly assumed) where the deviation from equilibrium starts to grow exponentially. After solving for x * iteratively as a function of m χ , σ v and g * (the relativistic degrees of freedom at x = x * ), Equation (5) can be integrated from x = x * to x = x F dropping the Y χ,eq term, to finally obtain an improved analytic solution for the relic density (in the s-wave limit) [42]: where the subscript * means the values evaluated at x = x * , and The analytic result in Equation (7) agrees with the exact numerical result within ∼ 3%, almost independent of the DM mass. Note that for an arbitrary l-wave annihilation, the above formalism can be repeated by Taylor-expanding σ v in powers of v 2 r ∼ 1/x.

RELATIVISTIC CASE
In the other extreme limit, where the freeze-out occurs when the χ particles are still relativistic (x F 1), their current relic abundance Y χ (x 0 ) is approximated by the equilibrium abundance at freeze-out Y χ,eq (x F ) [26,43]. In this case, Equation (2) gives n χ,eq = (ζ (3)/π 2 )g eff T 3 , where g eff = g χ (3g χ /4) for bosonic (fermionic) χ , and ζ (x) is the Riemann zeta function. Using the entropy density in the relativistic limit as given by Equation (3), we obtain Y χ,eq (x F ) = 0.28g eff /g(x F ) which is insensitive to the details of freeze-out. From Equation (6), the present relic density is then given by Relativistic DM particles in our Universe will lead to large damping scales > ∼ 10 Mpc (roughly the size of typical galaxy clusters), thereby suppressing the growth of small-scale structures. They would predict a top-down hierarchy in the structure formation [44,45], with small structures forming by fragmentation of larger ones, while observations have shown no convincing evidence of such effects, thereby imposing stringent upper limits on these "hot" DM species. For instance, the SM neutrino contribution to the non-baryonic DM relic density is currently constrained to be ν h 2 ≤ 0.0062 at 95% confidence level (CL) [36]. Thus, hot DM cannot yield the total observed DM density in our Universe [2], and if it exists, 7 must coexist with other cold/warm components. For example scenarios of such multi-component DM, (see 18, 47-50).

SEMI-RELATIVISTIC CASE
In the intermediate regime x F ∼ 1, the χ particles are semirelativistic when they decouple from the thermal bath. The improved analytic treatment of Steigman et al. [42], as discussed in section 2.1, is not applicable in this case, since the thermally averaged cross section σ v involves multiple integrals, and cannot be expanded in a Taylor series of the velocity-squared. One way is to approximate the cross section by interpolating between its relativistic and non-relativistic expressions. Following this approach, it was shown [51] that the Maxwell-Boltzmann distribution can still be used to compute σ v , and the more appropriate Fermi-Dirac or Bose-Einstein distributions are only needed for the calculation of the freeze-out abundance Y χ,eq (x F ). Note that although the current observations do not rule out the possibility of the whole DM density being comprised of warm DM species (see e.g., [52][53][54]), there exist strong constraints from observations of early structure, in particular from Lyman-α forest data [55,56].
On the other hand, if the interaction of the DM particles with the thermal bath is not large enough, they may not come into full LTE before they decouple from the plasma. In such cases, their current relic density Y χ (x 0 ) in Equation (6) will also depend on the initial abundance, and hence, on the production mechanism. This is discussed in the following section with a simple production mechanism.

DM FROM INFLATON DECAY
As discussed in section 1, we assume that the DM particles χ directly couple to the inflaton field φ so that it can be produced in the perturbative inflaton decay for m χ < m φ /2. 8 The initial energy density stored in the inflaton field is ρ φ ≈ n φ m φ which is transferred to the decay products at the end of inflation, thereby (re)heating the Universe with a temperature T R . Assuming the Universe to be radiation dominated immediately after inflation, the total energy density is given by ρ r = (π 2 /30)gT 4 R . 9 Hence, the initial DM number density is given by where B χ is the branching ratio of the inflaton decay to DM. Since we are interested in model-independent constraints on the DM parameter space, we keep our discussion general in terms of the branching ratio, without specifying its exact formula in terms of the DM-inflaton couplings, their masses, and the n-body decay kinematics (for n ≥ 2, depending on the specific DM candidates). Once produced, depending on the strength of their interaction with the thermal plasma, they could either thermalize fully/partially with the thermal bath or could remain non-thermal throughout their evolution. In the former case, their current relic density will be determined by their freeze-out abundance, independent of the initial abundance set by the inflation parameters. This is also true for the freeze-in scenario, provided the initial abundance is small compared to the thermal abundance. In the non-thermal case, however, their final number density is essentially the same as their initial abundance, only redshifted by the Hubble expansion rate. These two different scenarios are discussed below in somewhat details, with some numerical examples.

THERMAL DM
In this case, depending on their thermal annihilation rate σ v , the χ particles can either reach full LTE (i.e., both kinetic and chemical equilibrium) with the plasma before decoupling or decouple from the plasma before the full equilibrium could be established. The former case occurs for large annihilation rates, which enable the χ particles to attain equilibrium soon after their 8 This is the necessary condition required solely due to kinematic reasons, and could be sufficient, for instance, for fermionic DM coupling to the inflaton through a φχχ term in the Lagrangian. For more complicated inflaton decay chains involving many particles, a more stringent kinematic condition may be required. Also we do not consider derivative couplings of the inflaton with DM, which could be the case for axion DM, for example. 9 Note that the assumption ρ φ = ρ r = (π 2 /30)gT 4 R which determines the reheat temperature may not be correct for all possible inflaton or moduli coupling to the SM d.o.f. This definition is correct for large inflaton coupling to matter, i.e., α φ ≥ 10 −7 [21]. Typically the moduli coupling to the SM d.o.f and DM will be very small, i.e., α φ ∼ (m φ /m pl ) 2 . Therefore, a more rigorous treatment of the reheating scenario is required for the case of moduli decay. production. In this case, the χ particles follow the equilibrium distribution until they freeze out at a certain stage, depending on the exact value of the interaction rate. Thus, the initial abundance is irrelevant for their final relic density. In the latter case, the annihilation rates are not large enough to bring the χ particles into full LTE, and hence, their final abundance is determined by the annihilation rate as well as the initial abundance given by Equation (10). For given inflaton and DM masses, the final relic abundance χ h 2 will exceed the observed value for a large reheat temperature and/or large branching ratio of the inflaton to DM, thus overclosing the Universe. If the initial abundance is small, the DM particles can still be produced from the thermal plasma unless the interaction rate is utterly negligible. The dominant production in this case occurs at temperatures T > ∼ m χ when the interaction rate is still larger than the Hubble rate, and as the interaction rate drops below the Hubble rate, the relic abundance will freeze-in. We discuss below both freeze-out and freeze-in scenarios for the DM produced from inflaton decay, and give a numerical example for each case to illustrate the magnitudes of the interaction cross section, as compared to the well-known thermal WIMP scenario.

Freeze-out
In this case, the final relic abundance of the DM species is set by the freeze-out abundance which is determined by the freezeout temperature. This is obtained by solving the Boltzmann Equation (5) for Y χ . To calculate the freeze-out abundance more precisely, we track the evolution of the quantity χ = (Y χ − Y χ,eq )/Y χ,eq which represents the departure from equilibrium. From Equation (5), the evolution equation for is obtained to be of the form where χ,eq = n χ,eq σ v = Y χ,eq s σ v is the equilibrium annihilation rate, and H(x) can be readily obtained from Equation (4). For a Maxwell-Boltzmann distribution, the equilibrium number density and thermally averaged annihilation cross section are respectively given by Gondolo and Gelmini [34] Y χ,eq (x) = 45 where K n (x) is the n-th order modified Bessel functions of the second kind, and √s is the center-of-mass energy. Strictly speaking, Equation (13) is only applicable for the non-relativistic case with x 1. However, as noted in Scherrer and Turner [39] and Drees et al. [51], this is a good approximation (within 3% accuracy) even for the semi-relativistic case with x ∼ 1. For the relativistic case x 1, the final abundance is simply the equilibrium abundance, as given in Equation (9).
Following the strategy developed in Steigman et al. [42] to solve Equation (11) for , we note that in the early stages of evolution, Y χ tracks Y χ,eq closely, and hence, , d /dx 1. In this case, the left-hand side of Equation (11) can be safely dropped, thus leading to As the χ particles start freezing out with increasing x, increases exponentially, eventually becoming much larger than 1. Thus for some intermediate value of x = x * , ∼ O(1), and for x > x * , it grows exponentially. We define x * when (x * ) ≡ * = 1/2, 10 and solve Equation (14) iteratively for x * as a function of m χ , σ v and g * . For the logarithmic derivative of of g(T), we use the calculations of Laine and Schroeder [57] for the SM relativistic d.o.f. For the cases with no phase transition around is almost constant, and hence, this term can be ignored in Equation (14). Once the value of x * is found, we can determine T * = m χ /x * and Y χ (x * ) = (3/2)Y χ,eq (x * ) (corresponding to * = 1/2). The actual freeze-out temperature T F is somewhere below T * , since at T = T * , ( χ /H) * is still larger than 1 [42]. For x > x * , Y χ Y χ,eq , and hence, the Y 2 χ,eq term in Equation (5) can be dropped. Integrating from x = x * to x = x 0 , we obtain the present relic abundance: which can be used in Equation (6) to compute χ h 2 . To perform the integration in Equation (15), we need to know the x-dependence of σ v using Equation (13) which is one of the key quantities that determine the current relic density. In general, one can find an ansatz for σ v which smoothly interpolates between the non-relativistic and relativistic regimes. For simplicity, we will use the ansatz for an s-wave annihilation of two Dirac fermions [51]: where α χ denotes the coupling constant of the four-fermion interaction, which we will treat as a free parameter. This approach works well for DM species that freeze out between 0.5 x F 15. For x F > 15 (roughly corresponding to m χ > 10 MeV), the particles are already non-relativistic at decoupling, and hence, one can expand σ v in a Taylor series in terms of the averaged relative velocity: 10 As verified in [42], other alternative choices of * change the final result only by about 0.1%.
For s-wave annihilation, only the first term is considered, and in this case, Equation (15) simplifies further to finally yield the relic density given by Equation (7). We note that this approximation of using a constant value for σ v also works well in the semi-relativistic case, and induces an error of only about 6%, as compared to using the ansatz given by Equation (16). From Equation (15), it is clear that the final abundance is inversely proportional to the thermal annihilation rate. Thus, the larger the cross section, the longer the DM particles stay in equilibrium with the thermal bath, and hence, the lower the final abundance. This is true for both cold and warm DM cases, while for the hot DM case, the freeze-out is insensitive to the interaction cross section, as discussed in section 2.2.
The dependence of the current relic abundance on the annihilation rate for the thermal DM which has frozen out is illustrated in Figure 1. Here we have chosen m χ = 100 GeV. The solid black line shows the equilibrium distribution which is constant in the extreme relativistic regime (x 3), and exponentially suppressed in the non-relativistic regime (x 3), as can also be seen from Equation (12) by taking the asymptotic limits of the Bessel function. For large enough annihilation rates, the DM particles quickly thermalize, thereafter following the equilibrium evolution until their freeze-out, and the final relic abundance is independent of the initial abundance. The observed relic density as measured by Planck, shown as the horizontal band, is obtained for the thermal annihilation rate of σ v = 2 × 10 −26 cm 3 s −1 , as shown by the solid red line. As the annihilation rate decreases, the DM freezes out earlier (with smaller x F ), thus giving a larger relic density.

Freeze-in
In this scenario, the DM particles are very weakly coupled to the bath, and hence, cannot reach full thermal equilibrium with the bath before decoupling. However, the feeble interactions with the thermal bath (either directly [27] or mediated by a portal [58]) could still populate the DM, until the interaction rate drops below the Hubble rate when the DM abundance will freeze in. In this case, the final abundance is directly proportional to the interaction strength; the larger the interaction cross section is, the more DM particles are produced. In this sense, freeze-in can be viewed as the opposite process to freeze-out.
The final relic density in the freeze-in scenario will in general be determined by both the interaction cross section and the initial abundance which in turn depends on the reheat temperature and the branching ratio of the inflaton in our case. Note that the decoupling in this case occurs for small values of x F , where the equilibrium abundance Y χ,eq is independent of x, as can be seen from Equation (12): Also since the DM particles decouple very soon after being produced, the annihilation cross section as well as the number of relativistic degrees of freedom can be treated as constant with respect to x during this short period of time. Hence, the general Boltzmann Equation (5) can be approximated in this case to the following simple form: where A and B are constants in x. Equation (19) has a simple analytic solution in terms of the initial values x i = m χ /T R and Y χ,in , where the latter can be obtained from Equations (10) and (3): In the limit x → ∞, the expression for Y χ (x) simplifies further, and the final relic density can then be obtained using Equation (6). This has two contributions: where the first term represents the non-thermal contribution which only depends on the initial abundance, and the other two terms represent the thermal contribution which also depend on the interaction rate. Note that the analytic expression (21) is valid as long as m χ T R otherwise the thermal production will be delayed to lower values of temperature (or higher values of x) when the equilibrium distribution in Equation (19) may no longer be flat, but exponentially decaying. For the freeze-in scenario, it is usually assumed that the initial abundance is negligible, so that the final abundance is solely determined by the interaction strength in Equation (21), as in the freeze-out scenario. This is illustrated in Figure 1 for a typical choice of parameters: m χ = 100 GeV, m φ = 10 13 GeV, T R = 10 TeV, and B χ = 10 −15 so that the initial abundance given by Equation (20) is negligible. The different dashed lines in Figure 1 correspond to the freeze-in scenario with various interaction rates, and hence, different final abundances. Note that the final abundance increases with increasing interaction rate, in contrast with the freeze-out scenario (the solid lines) where the final abundance decreases with increasing interaction rate. As shown here, the observed relic abundance shown by the gray horizontal band can be obtained in the freezein scenario for an interaction rate of 10 −47 cm 3 s −1 , which is much smaller than the typical value of 2 × 10 −26 cm 3 s −1 , as in the freeze-out scenario.
We should mention here that there could be other thermal production mechanisms for the DM in specific models, depending on its interaction with the SM particles and/or the model construction for the beyond SM sector. For instance, a keV-scale sterile neutrino DM can be produced by the Dodelson-Widrow mechanism [59], which is very similar to the freeze-in mechanism discussed above.

NON-THERMAL DM
For very small cross sections, the DM particles are produced already decoupled from the thermal bath, and hence, the thermal production in Equation (21) is negligible compared to the initial abundance, which could be sizable for large branching ratios. In this case, the annihilation rate, and hence, the right-hand side of Equation (5) can be neglected, thus leading to dY χ /dx 0. Hence, the final relic abundance is completely determined by the initial one given by Equation (20). Using the general expression (6), this yields the non-thermal relic DM density which can also be identified with the first term on the righthand side of Equation (21). Thus for super-weak interaction rates, the final abundance only depends on the reheat temperature and inflaton branching fraction for given DM and inflaton masses 11 . Some illustrative cases for the non-thermal DM are shown in Figure 2 for two typical values of the branching ratio B χ = 10 −5 and 10 −15 . The choice of small values of B χ will be justified below. The various contours show the reheat temperature values required to obtain the correct relic density χ h 2 = 0.12 for given values of the inflaton and DM masses. These plots were obtained by numerically solving the Boltzmann Equation (5) for a typical annihilation rate σ v = 10 −60 cm 3 s −1 (see section 5 for details) following the procedure mentioned above, but the results agree quite well with the approximate analytic formula given in Equation (22). From Figure 2 it is clear that as the inflaton branching fraction increases, the allowed

FIGURE 2 | The colored contours show the reheat temperature values required to give the correct relic density for non-thermal DM as a function of the inflaton and DM masses, for a given inflaton branching ratio.
range of the DM mass shifts to lower values in order to satisfy the observed relic density, in accordance with Equation (22). We have shown the results for the inflaton mass m φ in the range 10 3 -10 13 GeV, the reheat temperature T R between 1 and 10 9 GeV and for the DM mass m χ ≤ m φ /2. Note that a late-time entropy production would induce various cosmological effects, leading to a lower limit on the reheat temperature of about 1 MeV from BBN constraints [61,62], which, when combined with the CMB and large scale structure data, improves to about 4 MeV at 95% CL [63]. However, for T R below the QCD phase transition scale, QCD = O(100) MeV, the estimation of the number density of particles from inflaton decay suffers from large QCD uncertainties due to the non-perturbative hadronization effects. Hence, we have used a conservative value, which is one order of magnitude above the hadronization scale, as the minimum value of T R in Figure 2 (left panel) 12 .
For m χ m φ , the non-thermal DM directly pair-produced from the inflaton decay will have a large velocity at the time of matter-radiation equality (t eq ), unless the reheat temperature is sufficiently high to make the velocity small due to redshift. The comoving free-streaming length of a non-thermal DM at matterradiation equality is given by 12 One should also note that for T R below the electroweak phase transition scale, EW = O(100) GeV, one cannot use the standard electroweak sphaleron processes to explain the observed baryon asymmetry, and must invoke a post-sphaleron baryogenesis mechanism [64][65][66][67][68].
where a(t) is the scale factor, t D is the time at inflaton decay, and is the magnitude of the velocity of the DM particle. Integrating Equation (23), and requiring that λ fs < ∼ 1 Mpc, from Lyman-α constraints (for warm/cold DM), one obtains a lower limit on the reheat temperature [69] T R > ∼ 5 × 10 4 GeV g 200 Combining this with Equation (22), and requiring that χ h 2 ≤ 0.12 to satisfy the observed relic density for cold/warm DM relics, we derive an upper limit on the branching ratio of the inflaton decay to DM: B χ < ∼ 0.01(g/200) 1/4 . This is complementary to what is already expected from the fact that for a standard Cosmology, B χ must be small in order to have a radiation-domination epoch immediately after reheating, followed by matter domination only at a late stage.

EXPERIMENTAL CONSTRAINTS
In this section, we summarize the various experimental constraints on the DM properties relevant for our analysis.

OVERCLOSURE
For any DM candidate, we must ensure that it does not lead to an overclosure of the Universe. Thus, we set the upper limit on the relic density of our χ particles coming from inflaton decay using the observed value CDM h 2 = 0.1199 ± 0.0027 (68% CL; Planck + WP) [2], which combines the Planck temperature data with WMAP polarization data at low multipoles. We do not set a lower limit on χ since for the cases in which the χ particles do not account for the total observed abundance, the remaining fraction can be obtained by invoking a hidden-sector/multi-component DM scenario (see e.g., 18, 47-50).

UNITARITY
The partial-wave unitarity of the scattering matrix, together with the conservation of total energy and momentum, impose a generic upper bound on the cross section of thermal DM annihilation into the j-th partial wave [70]: where v r = 2 1 − 4m 2 χ /s is the relative velocity between the two annihilating particles in the center-of-mass frame with total energy √s . Assuming that the s-wave piece with j = 0 dominates in the partial-wave expansion, we obtain an upper bound on the thermally averaged annihilation rate σ v as a function of the DM mass from Equation (13), where σ is replaced with (σ 0 ) max from Equation (26). Since the current abundance of a non-relativistic thermal relic scales as χ ∝ 1/ σ v , the observed DM relic density constrains the mass of the thermal relic to be m χ < ∼ 130 TeV to satisfy the unitarity bound. Note however that this bound may not be applicable when the higher partial-waves are not suppressed, as is the case when the DM particles decouple from the thermal bath while still being relativistic.

PLANCK
Precision measurements of the CMB angular power spectrum by Planck put stringent constraints on the number of effective neutrino species (N eff ), which parametrizes the total radiation energy density of the Universe: where ρ γ = (π 2 /15)T 4 is the energy density of photons, and the neutrino-to-photon temperature ratio T ν /T γ = (4/11) 1/3 assuming exactly three neutrino flavors and their instantaneous decoupling. In the standard cosmological model, T ν /T γ is slightly higher than (4/11) 1/3 due to partial reheating of neutrinos when electron-positron pairs annihilate transferring their entropy to photons, thus giving N eff = 3.046 [71]. Now if the DM species remains in thermal equilibrium with the neutrinos or electrons and photons after neutrino decoupling, and transfers its entropy to them during its annihilation after it decouples at a later stage, it can increase or decrease the value of N eff as we decrease the DM mass. Using the constraints on N eff from Planck [2], together with the helium abundance Y p , [72] derived a robust lower bound of 2-10 MeV on the thermal DM mass, depending on whether it is a fermion (Dirac/Majorana) or scalar (real/complex) and whether it was in equilibrium with neutrinos or with electrons and photons.
Another generic lower bound on the cold DM mass can be obtained using the CMB and matter power spectrum observations which place an upper bound on the DM temperatureto-mass ratio: T/m χ ≤ 1.07 × 10 −14 (1 + z) 2 [73]. Evaluating this bound at matter-radiation equality with a redshift of z eq = 3391 [2] and T γ,eq 0.77 eV [26], we obtain a lower limit of m χ > ∼ 6.5 keV, which is much weaker than the limit derived in Boehm et al. [72] using N eff .

DARK RADIATION
The Planck constraints on N eff can also be used to set an upper limit on the amount of dark radiation at decoupling. From Equation (27), the radiation energy density apart from the photon and SM neutrino contribution is given by where γ h 2 = 2.471 × 10 −5 (T/2.725) 4 is the CMB radiation density [36], and N eff = N eff − 3.046. Using the 95% CL measured value of N eff = 3.30 +0.54 −0.51 from Planck+WMAP-9 polarization data+SPT high-multipole measurement+Baryon Acoustic Oscillation measurements from large scale structure surveys (Planck+WP+highL+BAO) [2], we obtain an upper limit on the amount of dark radiation from Equation (28): dark h 2 ≤ 4.46 × 10 −6 . This also sets the upper limit on the relic density of hot DM species. In order to obtain the mass range in which the thermal DM species decouple while being relativistic, we calculate their free-streaming length [74]: where z is the redshift, and I 0 = [ ∞ 0 f (0) (y)dy]/[ ∞ 0 y 2 f (0) (y)dy] is a dimensionless ratio, given in terms of the comoving energy distribution function f (y) = 1/ exp [ y 2 + x 2 + 1] with x = m χ /T and the superscript (0) refers to the current value of the distribution. From the Ly-α constraints, we require λ fs (0) < ∼ 1 Mpc for cold/warm DM candidates. For concreteness, we will impose the dark radiation upper limit from Equation (28) for the parameter space corresponding to λ fs (0) > 2 Mpc.

BBN AND CMB
The late annihilation of DM particles (after freeze out) can deposit hadronic and/or electromagnetic energy in the primordial plasma, thereby altering the history of nucleosynthesis (BBN) [75][76][77] and recombination (CMB) [78][79][80][81][82][83][84][85]. These effects depend only on the type and rate of energy injection into the thermal bath, thus allowing to set rather model-independent bounds on the annihilation rate, especially for DM masses in the MeV-GeV range. During nucleosynthesis, the injection of hadronic and/or electromagnetic energy can affect the abundance of nuclei via (i) raising the neutron-to-proton ratio and therefore the primordial 4 He abundance, and (ii) high energy nucleons and photons disassociating nuclei. During recombination, the injected electromagnetic energy ionizes hydrogen atoms, which results in an increased number of free electrons, causing the broadening of the surface of last scattering, and results in scaledependent changes to the CMB temperature and polarization power spectra, especially in the low multipole modes. The precision measurements of BBN and CMB from WMAP and Planck data have been used to set upper bounds on the DM annihilation cross section σ v , as a function of the DM mass [75][76][77][81][82][83][84][85].

INDIRECT DETECTION
After freeze-out, the relic annihilations of WIMP DM may be indirectly observed by searching for their annihilation products such as charged particles, photons and neutrinos (for a review, see 86). In fact, a number of indirect detection experiments have observed an excess of electrons and positrons in the charged cosmic ray flux, and this was recently confirmed with the precision measurements by AMS-02 [87]. Assuming a possible DM contribution to this positron excess and using the high quality of AMS-02 data, [88] has performed a spectral analysis to put stringent constraints on the DM annihilation cross section for various leptonic final states 13 . Similar constraints were obtained for the DM annihilation into hadronic final states [93,94] in order to explain the absence of a corresponding excess in the cosmic-ray antiproton flux in the PAMELA data [95,96].
The DM annihilation to various SM final states can also lead to an observable photon flux which can be produced either by direct DM annihilation ("prompt" gamma-rays) or by inverse Compton scattering and synchrotron emission of the electrons and positrons created in the DM annihilation. These photon signals are preferentially searched for in regions with high DM densities and/or regions with reduced astrophysical background. The Fermi-LAT, with its unprecedented sensitivity to gamma rays in the MeV-TeV energy range, has performed deep searches for line spectrum (mono-energetic gamma-rays due to direct DM annihilation) [97] as well as continuum spectrum (through DM annihilation into intermediate states) [98] 14 . They have derived additional constraints on the DM annihilation cross section from the isotropic diffuse gammaray emission in the galactic halo [100], nearby galaxy clusters [101], and nearby dwarf spheroidal galaxies [102]. Similar constraints were also derived from the galactic center region for various DM density profiles [103]. Complementing the Fermi-LAT range toward higher energies, the HESS collaboration has performed a number of DM searches up to multi-TeV DM masses [104][105][106].
The DM annihilation can also produce neutrinos which, like gamma-rays, can travel essentially unabsorbed through the galaxy, and can be observed at large neutrino detectors on Earth. Constraints on the DM annihilation rate were derived by the 13 A DM interpretation of the AMS-02 positron excess is still viable, if the DM annihilates to four-lepton final states [89][90][91][92]. 14 There exists yet another class of spectral signature, namely, box-shaped gamma-ray spectrum, which arises if the DM annihilates/decays into intermediate particles which further decay into photons [99]. The cross-section limits derived using this feature are currently comparable to those obtained using the line-like spectral feature.
IceCube experiment from the upper limits on the high-energy neutrino fluxes from the galactic halo [107], galactic center [108], dwarf galaxies and clusters of galaxies [109]. These limits are currently somewhat weaker than the gamma-ray limits for low DM masses, but become competitive at larger DM masses. Combining the Fermi-LAT data on the diffuse gamma-ray and the IceCube data on diffuse neutrino flux, robust constraints were derived on the DM annihilation rate for heavy DM masses (1 TeV -10 10 GeV) [110].

RESULTS AND DISCUSSION
Using the model-independent approach outlined in section 3, we solve the Boltzmann equation (5) numerically for the evolution of DM produced from inflaton decay. Here we assume an s-wave annihilation, and take the annihilation rate σ v to be a free parameter. 15 Both thermal and non-thermal regions are identified in the (m χ , σ v ) parameter space. Our results are shown in Figures 3-6 for a fixed inflaton mass m φ = 10 13 GeV. We consider two typical values of the reheat temperature T R = 10 9 GeV and 10 4 GeV, and branching ratios of the inflaton decay B χ = 10 −5 and 10 −15 for our illustration purposes. We have considered the DM masses only below the reheating temperature, and do not analyze scenarios in which DM could be produced during preheating or reheating (e.g., the WIMPzilla scenario 114).
For each case shown in Figures 3-6, we calculate the current relic density of the thermal DM to show the overclosure region (blue-shaded) which rules out a wide range of the parameter space, irrespective of the initial choice of parameters. In the remaining allowed parameter space, we identify the region with very low annihilation rates belonging to the non-thermal DM scenario (green-shaded region) since for such extremely small interaction rates, the DM particles cannot achieve LTE, and decouple soon after being produced essentially with their initial abundance, as discussed in section 3.2. So the overclosure condition for non-thermal DM will be determined solely by the initial conditions, and this will be discussed in section 5.2.

THERMAL CASE
In the thermal DM regime, the region above the overclosure region with large annihilation rate belongs to the freeze-out scenario, while in the white region below the overclosure one with small annihilation rate belongs to either freeze-out or freezein scenario, depending on the initial conditions. The observed value of the relic density is obtained at the boundary between these regions with the overclosure region (shown by the solid and dotted blue lines). The thermal freeze-out region with large annihilation rate is severely constrained by many experimental searches, as discussed in section 4, some of which are shown by the shaded regions 1-15 in Figures 3-6, and also summarized below: FIGURE 3 | Model-independent constraints on the DM annihilation rate as a function of the DM mass for both thermal and non-thermal production mechanisms. Here we have chosen m φ = 10 13 GeV, T R = 10 9 GeV, and B χ = 10 −5 as the initial parameters for the DM evolution. The blue-shaded region is excluded from relic density constraints, and the observed relic density is obtained at its boundary (shown by the solid and dotted blue lines). The green-shaded region at the bottom represents the non-thermal DM scenario, while in the rest of the parameter space, the DM can be fully/partly thermalized with the primordial bath. The various colored-shaded regions in the thermal region are excluded (under certain assumptions) by the constraints given in section 4; see text for details. • Region 1 is excluded by the dark radiation constraint, as discussed in section 4.4. In this region, the comoving freestreaming length is greater than 2 Mpc, thus corresponding to a hot DM regime, while the relic density χ h 2 exceeds the upper limit of 4.46 × 10 −6 derived from the Planck data [2]. • Region 2 is excluded by the recent Planck measurements of the effective number of neutrino species, as discussed in section 4.3, assuming that the DM interacts with neutrinos or electrons and photons after the neutrino decoupling. This sets a robust lower bound of order of MeV on the thermal DM mass with large interaction rates. The precise value of the lower bound depends on whether the DM is a scalar or fermion and on whether it couples to neutrinos or to electrons and photons. The bound shown by region 2 assumes a Majorana fermion DM coupling to neutrinos. Note that [72] had originally derived this limit for a cold DM candidate, but this is generically applicable as long as the interaction rate is large enough to keep the DM in LTE after the neutrino decoupling, thus transferring entropy at a late stage and affecting N eff . • Regions 3 and 4 are excluded by the BBN data, as discussed in section 4.5, and assuming DM annihilation into electronpositron and quark-antiquark pairs respectively [77]. Similarly, the region 9 is excluded by constraints derived from a combination of the CMB power spectrum measurements from Planck, WMAP9, ACT, and SPT, and low-redshift datasets from BAO, HST and supernovae [85]. • Region 5 is excluded by the Fermi-LAT limit at 95% CL derived using the diffuse gamma-ray flux from a combined analysis of 15 dwarf spheroidal galaxies, for an NFW DM density profile and assuming DM annihilation into tau-antitau final states [102]. Region 7 is excluded by the 3σ Fermi-LAT limit obtained using the diffuse gamma-ray emission in the Milky Way halo, assuming an NFW DM distribution and for annihilation into bottom-antibottom quark pairs [100]. Region 8 is excluded by a similar analysis using the Fermi-LAT data from galactic center [103]. The corresponding limits for other SM final states are weaker, and are not shown here for clarity purposes. • Region 6 is excluded by the Fermi-LAT 95% CL upper limit on the cross section of DM annihilation to two photons from a dedicated search for the gamma-ray line spectrum [97]. Region 13 is excluded from a complementary search for the line spectrum by HESS [106]. Note that these limits, although very stringent, can be evaded in most of the popular WIMP DM models, since the direct annihilation to photon final states is suppressed due to loop effects. • Region 10 is excluded by the measurements of the antiproton flux from PAMELA, and assuming the DM annihilation to bb final states [94]. These limits are applicable only for hadronic final states. Similarly, region 15 is excluded by the 95% CL upper limits, derived from the AMS-02 data, on the DM annihilation cross section for e + e − final state [88]. The corresponding limits for other leptonic final states are somewhat weaker, and hence, are not shown here. • Region 11 is excluded by the IceCube upper limit on the DM annihilation cross section for neutrino final states for the Virgo galaxy cluster including subhalos [109]. The corresponding limits for other final states as well as from searches in galactic halo [107] and galactic center [108] are somewhat weaker. • Region 12 is excluded by the cascade gamma-ray constraints obtained using the Fermi-LAT diffuse gamma-ray background data up to very high energies [110]. The corresponding limits derived using the IceCube high-energy neutrino data are stronger at higher DM masses, but weaker than the unitarity constraint (see section 4.2). • Region 14 is excluded by the unitarity constraints [70], as discussed in section 4.2, which sets an upper limit on the CDM mass of about 130 TeV for the allowed region, and rules out heavy thermal DM, even with annihilation rates many orders of magnitude below the thermal annihilation rate. This theoretical constraint is the most stringent one for very heavy DM, and is applicable as long as the DM is produced thermally.
Note that for the indirect detection constraints, we have shown only a few of them (typically the most stringent ones) for illustration purposes. Most of these limits have limited applicability, as they were derived assuming DM annihilation into a particular final state, and could be evaded in specific models where some of these annihilation channels might be suppressed due to various reasons. Also note that additional constraints on the annihilation cross section for a given DM mass might be derived using possible correlations with the DM direct detection cross section limits [115] and collider search limits from mono-jet [116,117] and mono-photon [118,119] final states with large missing energy.
In the absence of a collider signal for DM, model-independent constraints can be derived on the mass and interactions of a generic WIMP DM candidate from direct and indirect detection searches [120]. The other allowed thermal DM parameter space, namely, the region with very low interaction rates such as the FIMP scenario, is hard to constrain from the existing experimental limits. Various experimental tests of the freeze-in mechanism by measurements at colliders or by cosmological observations were outlined in Hall and Jedamzik [27]. However, these signals depend very much on the particular freeze-in scenario under consideration, and hence, it is difficult to derive model-independent constraints in the (m χ , σ v ) parameter space, except for the generic dark radiation constraint as shown in Figures 3-6. Just to give an example of additional model-specific bounds, a keV-scale sterile neutrino DM, which has a small interaction rate due to its mixing with the active neutrinos, could radiatively decay to an active neutrino and a photon which will lead to a mono-energetic X-ray line [121], the absence of which puts severe constraints on such keV-scale DM models, including their production mechanisms [122].

NON-THERMAL CASE
Now we move on to discuss the non-thermal DM region (greenshaded) in Figures 3-6. As already discussed at length in section 3.2, the final relic density of these DM particles is solely determined by the initial conditions, which in our case, are set by the inflaton and DM masses, the reheat temperature and the branching ratio of the inflaton decay to DM. For fixed reheat temperature and inflaton branching ratio, we show in Figures 7,  8 the contours for relic density computed using Equation (22) in the (m φ , m χ ) plane. We also calculate the comoving freestreaming length using Equation (23), and identify the regions with λ fs < 10 kpc as cold DM (CDM), with λ fs > 2 Mpc as hot DM (HDM), and the rest as warm DM (WDM). Note that there is no well-defined boundary between these regions, and we have just chosen some typical values derived from various astrophysical data [56,74] for our illustration purposes. We find that the observed DM relic density can be satisfied for a narrow parameter space in the CDM region (the boundary between the blue and orange regions), and the region above this is excluded due to overclosure constraints. For the HDM case with T R = 10 9 GeV and B χ = 10 −5 (Figure 7, left panel), an additional portion of the parameter space (blue-shaded region at bottom-left corner) is ruled out due to the dark radiation constraint, as discussed in section 4.4.
Similar to the thermal DM case, additional constraints can be derived for specific non-thermal DM candidates. For instance, a popular class of such candidates, known as the Weakly Interacting sub-eV particles (WISPs) such as axions and axion-like particles which often arise as the Nambu-Goldstone bosons associated with some global symmetry breaking, can be constrained from various low-energy experiments involving lasers, microwave and optical cavities, strong electromagnetic fields or torsion balances [123].

CONCLUSION
In this paper, we have investigated the thermal and non-thermal properties of DM from inflaton decay in a model-independent manner, assuming that the reheating and thermalization of the FIGURE 7 | The color-coded contours show the relic density of non-thermal DM produced from inflaton decay as a function of the inflaton and DM masses for a fixed reheated temperature T R = 10 9 GeV and fixed inflaton branching ratios B χ = 10 −5 and 10 −15 . We identify the cold, warm and hot DM regions in each case by assuming that the corresponding free-streaming length given by Equation (23) should be < 10 kpc, between 10 kpc -2 Mpc, and above 2 Mpc respectively. The blue-shaded region for the CDM case is excluded by the overclosure constraint, as discussed in section 4.1. The additional blue-shaded region in the HDM case is ruled out by the dark radiation limit, as discussed in section 4.4. ambient plasma have happened instantly at a given unique temperature. In the thermal DM scenario, the relic abundance of the DM species is determined by the freeze-out abundance, irrespective of the initial conditions or production mechanism, provided its interaction with the thermal bath is large enough to bring it into LTE soon after its production. For smaller interaction rates when the DM does not attain full LTE, but can still be produced from the thermal bath, one can also obtain the correct relic density through freeze-in mechanism. On the other hand, if the interaction rate is negligibly small so that the DM remains decoupled from the thermal bath from the beginning, the relic density is essentially determined by the initial conditions. Assuming that the DM has a non-zero coupling to the inflaton so that it can be directly produced from the inflaton decay, we have investigated all the above scenarios by tracking the evolution of the DM species from the very onset of its production. We have numerically solved the Boltzmann equation for DM number density, and have shown that the inflaton decay to DM inevitably leads to an overclosure of the Universe for a large range of parameter space, especially for non-thermal DM scenarios. This is an important constraint for hidden sector DM models with an arbitrary DM coupling to the inflaton. For the thermal DM scenario with large annihilation rates, we show the complementary constraints on the DM parameter space from various experimental searches. On the other hand, the other viable regions for both thermal and nonthermal DM candidates with very small interaction rate remain mostly unexplored.