WIMPs, FIMPs, and Inflaton phenomenology via reheating, CMB and $\Delta N_{eff}$

In this paper, we extensively analyzed the reheating dynamics after inflation and looked into its possible implication on dark matter (DM) and inflaton phenomenology. We studied the reheating through various possible channels of inflaton going into massless scalars (bosonic reheating) and fermions (fermionic reheating) via non-gravitational and gravity-mediated decay processes. We further include the finite temperature effect on the decay process. Along with their precise roles in governing the dynamics, we compared the relative importance of different temperature-corrected decay channels in the gradual process of reheating depending on the reheating equation of state (EoS), which is directly related to inflaton potential. Particularly, the universal gravitational decay of inflaton is observed to play a very crucial role in the reheating process for a large range of inflaton decay parameters. For our study, we consider typical $\alpha$-attractor inflationary models. We further establish the intriguing connection among those different inflaton decay channels and the CMB power spectrum that can have profound implications in building up a unified model of inflation, reheating, and DM. We analyze both fermion and scalar DM with different physical processes being involved, such as gravitational scattering, thermal bath scattering, and direct inflaton decay. Gravitational decay can again be observed to play a crucial role in setting the maximum limit on DM mass that has already been observed earlier in the literature [52]. Depending on the coupling strength, we have analyzed in detail the production of both FIMP and WIMP-like DM during reheating and their detailed phenomenological implications from the perspective of various cosmological and laboratory experiments.


I. INTRODUCTION
Reheating is a phenomenon that has been studied quite extensively over the years. It is the phase that bridges the two paradigms of cosmology, namely, inflation [1][2][3] and the standard Big Bang. While inflation sets the uniform initial condition for all the causally disconnected patches of exponentially large homogeneous space of the size of our present universe, the subsequent phase fills the spaces with visible hot matters homogeneously distributed through the process called reheating [4][5][6][7]. In the standard scenario, reheating is the physical process through which the inflaton field decays into the visible matter fields. The endpoint of the reheating is when the standard radiation-dominated universe begins, which sets the proper initial conditions for Big Bang nucleosynthesis (BBN). From these chronological cosmological events, it is obvious that the state of our present universe must be non-trivially dependent on the process of reheating. Depending on the nature of inflaton and its decay, reheating dynamics can be effectively described by parameters such as the equation of state of inflaton (EoS) and its coupling with other fields. Investigation of this phase is still an ongoing effort since its theoretical inception proposed in [4][5][6][7]. Since there is no way to directly probe this phase with the present experimental techniques, it is important instead to understand this phase through various possible physical processes and look for direct/indirect process-dependent observables which can be probed directly/indirectly in future experiments. Our present study comes with two main objectives: firstly, study the reheating through various possible decay channels; namely, i) inflaton (ϕ) decaying into bosonic radiation through g r 1 ϕs 2 , g r 2 ϕ 2 s 2 interaction, ii) inflaton decaying into fermionic radiation through h r ϕff interaction, iii) inflaton decaying into all fundamental fields through minimal s-channel graviton exchange interaction, (1/M P )(h µν T µν ϕ + h µν T µν s/f ) [8,9], and identify the region of their dominance in the parameter space of bosonic coupling g r i and fermionic coupling h r . T µν is the energy-momentum tensor and h µν is the tensor metric perturbation. To this end, we would like to point out that in the context of gravity-mediated decay, the effect of non-minimal Ricci curvature (R) coupling ξϕ 2 R has been considered in the reference [10]. However, the contribution of such a term has been shown to be negligible for dimensionless coupling ξ < 1 and hence will be ignored in this paper. On top of those couplings, we further included the finite temperature effect on the decay process and compared it with the zero-temperature case for all the cases. Our second objective is to study the dark matter (DM) production during reheating considering different physical processes via gravitational scattering, thermal bath scattering, and direct inflaton decay to DM via g x 1 ϕS 2 , g x 1 ϕ 2 S 2 interactions for scalar DM S, and h x ϕF F interaction for fermionic DM F . From the phenomenological perspective, DM is assumed to be an integral part of the visible standard model components in quantum field theoretic framework . In this framework it is just the DM mass and the cross-section which are shown to be sufficient to explain the current abundance of DM. However, a large number of attempts over the years to detect [33][34][35][36][37] such particle appeared to go in vain. Therefore, going beyond the existing framework of both experimental and theoretical approaches to understanding dark matter may seem to be important [38][39][40][41]. Towards this endeavor recent proposal of graviton mediated DM production [42][43][44][45][46][47][48][49][50][51] has been shown to have some promising universal features, and it cannot be ignored in any DM studies. Interestingly such gravitational production has been observed to set a limit on the maximum possible value of the DM mass [52,53]. In the standard DM literature, two distinct DM production mechanisms exist in the early universe. For the standard WIMP (Weakly Interacting Massive Particle) scenario, the DM is assumed to be in thermal equilibrium with the radiation bath. During the course of background evaluation, DM particles became thermally decoupled from the bath (known as the Freeze-out mechanism), and the present value of abundance is achieved [54][55][56][57][58][59][60][61]. In the second type, known as FIMP (Feebly Interacting Massive Particle) scenario, the DM is assumed to remain out of equilibrium with the radiation bath and produced due to decay of other fields throughout. During the course of background evaluation, the decay channel ceases to produce DM at some point (known as the Freeze-in mechanism), and the present value of abundance is achieved [62][63][64][65][66][67][68][69][70]. In this paper, we will explore those mechanisms in the context of the early universe with a non-trivial reheating phase. Apart from understanding the very nature of DM, such studies actually open up the possibilities of looking for the signature of reheating and, most importantly, the nature of inflaton through the physics of DM.
Most of the DM phenomenological studies were confined to the early radiation-dominated universe [54][55][56][57][58][59][60][61]. DM physics during reheating has gained significant interest only recently [25,26,62,[72][73][74][75][76][77][78][79][80][81]. Primary motivation of this attempt has two-folds: a) to construct a unified framework where inflaton is assumed to be an integral part of DM model building, and b) explore the physics of reheating and its impact on DM physics. Finally, analyze and constrain the inflaton, reheating and DM parameters through the constraint on extra relativistic degrees of freedom in terms of ∆N eff at the time of Big-Bang Nucleosynthesis (BBN) [82][83][84][85], and different DM searches in the astrophysical/cosmological/laboratory experiments [86][87][88][89][90]. Keeping those two-fold motivations in mind, we study in detail the parameter space wherein both WIMP and FIMP-type mechanisms can be realized. In this study, we will further see how universal gravitational DM production during the initial stage of the reheating phase plays an important role in constraining the DM parameters. WIMP production during reheating has been considered very recently in [91,92]. However, detailed studies taking into account the physics of inflation and the subsequent reheating processes are still lacking. In this paper, we fill this gap and study in detail its constraints and significance in the context of present and future experiments.
We organize our paper as follows: In section II, we will discuss the basic setup for the perturbative reheating processes for different decay channels and their connection with the inflationary parameters. As mentioned earlier, we include the effect of finite temperature corrections in the decay widths for various decay channels. In section IV, we discuss in detail the reheating dynamics due to two different bosonic decay channels. We identify the parameter regions with respect to reheating EoS, where reheating can be successfully achieved depending upon the strength of the different decay channels under consideration. We further discuss the possible constraints on those inflationary coupling parameters with respect to the CMB observations. In section V, we discuss in detail the fermionic reheating scenario where the reheating occurs due to the inflaton decaying into fermionic radiation. Similar to the bosonic reheating case, we analyze and constrain the inflaton-fermion coupling parameter through the perspective of inflationary observable and CMB constraints. In section VI, We consider various possible scenarios corresponding to DM production. We elaborately discuss DM production from direct inflaton decay and thermal radiation separately. In section VII and VIII, we discuss their possible constraints from the perspective of various theoretical and experimental bounds. Finally, we conclude with some future directions.

II. PERTURBATIVE REHEATING: GENERAL SET UP
In the first half of our present paper, we will discuss in detail the single-sector reheating, which has been studied earlier for some special cases [25,26,75,[93][94][95]. However, we will perform a comprehensive study on this with many significantly new results. In the second half, we include dark matter production and explore different possible production mechanisms and their observable possibilities in a unified framework.
After the end of inflation, the inflaton field starts to oscillate with a decaying amplitude. During this phase, the inflaton field transfers its energy into massless fields termed as radiation (R) through various decay channels, ϕ → ss/f f and ϕϕ → ss, where we assume explicit coupling between inflation and daughter fields. In addition, there exists universal gravitational coupling between inflaton and daughter field through s-channel graviton (h µν ) exchange interaction, (1/M 2 P )h µν T µν , which has recently been shown to influence the reheating dynamics [53,96]. In our present analysis, we shall include such universal contributions as well. When we talk about gravitational scattering, we mainly consider the process ϕϕ → h µν → ss/f f . In Fig.1, we have shown the Feynman diagram for all possible interactions between inflation (ϕ), and radiation (R). We consider inflaton coupled with scalar (s) and fermion (f ) as massless radiation for the case of bosonic and fermionic reheating, respectively. Compared to the direct decay with free coupling parameters, decay due to gravitational interaction is universal and hence will always be present. Setting all the direct inflaton coupling to zero implies purely gravitational scattering, named as gravitational reheating (GR), which has been studied in detail by two of the present authors in [53,96]. Therefore, we consider following general interaction Lagrangian, h µν M P (T µν ϕ + T µν s/f ) + g r 1 ϕs 2 + g r 2 ϕ 2 s 2 + h r ϕf f, where g r 1 is the coupling parameter for the trilinear interaction with mass dimension unity, g r 2 is the dimensionless coupling parameter for quartic interaction and h r is the Yukawa coupling.

A. Finite temperature decay widths and Boltzmann equation
The standard assumption of any reheating studies is that radiation is always thermalized among its constituents throughout the entire process. Hence, inflaton decay products must encounter a finite temperature thermal bath which will modify the decay width [97][98][99][100]. Such finite temperature effect has already been discussed in [101][102][103] for the special cases with a matter-dominated reheating phase. Our goal is to generalize those studies for the arbitrary inflaton equation of state (EoS) w ϕ , which has not been studied earlier. For the bosonic decay of inflaton, the thermal effect introduces Bose enhancement factor into the decay width, which makes bosonic reheating efficient compared to the zero temperature one. On the other hand, for fermionic decay of inflaton, the thermal bath induces an additional Pauli blocking factor into the decay width, which makes fermionic reheating less efficient compared to the zero-temperature one. Including the finite temperature effect, we list up the following decay widths 1 for various decay/scattering channels as [75,95,[103][104][105] Γ s/f = Γ gr ϕϕ→ss = ρ ϕ m ϕ 1024πM 4 p (1 + 2f B (m ϕ /T )) , [50,95] (3) [50,95] (4) where f B/F (z) = 1 e z ∓1 are the equilibrium Bose-Einstein (B)(−) and Fermi-Dirac (F)(+) distribution function. The last two decay width expressions are the gravity-mediated inflaton decay to other fields. In these expressions, we ignore the effect of thermal mass correction due to self-interaction. m ϕ (t) corresponds to time-dependent inflaton mass defined as m 2 ϕ (t) = ∂ 2 V (ϕ)/∂ϕ 2 for a generic inflaton potential V (ϕ). It is obvious from the expression that when the temperature of the radiation bath is greater than the inflaton mass, i.e. (T rad > m ϕ (t)), the Bose-enhancement or Pauli-blocking is effective. Due to the Bose enhancement, the decay rate is enhanced for the Bosons, and due to Pauli blocking, the decay rate is suppressed for the Fermions. Since the radiation particles are massless, the total decay width for the gravitational sector is mainly associated with the scalar particles. For fermionic particles, the decay width for gravitational interaction ∝ m 2 f , and that is the reason for not taking the contribution from the fermionic particles in the thermal bath for gravitational production. In addition, in the case of vector boson production, we need to introduce the mass term to break the conformal invariance; thus, massless vector boson production is impossible. Therefore, the total gravitational scattering rate to radiation production will be simply Γ gr ϕϕ→RR = Γ gr ϕϕ→ss . One important point is to note that we consider only one type of scalar particle in the radiation bath, and all other particles are assumed to be fermions or vector bosons, so their gravitational production is unimportant. It has already been established that for inflaton equation of states w ϕ > 0.65, the gravitational scattering can alone reheat the universe (see, for instance, reference [96]). In this paper, we will include the explicit inflaton coupling with radiation and do a comprehensive combined analysis to figure out the complete parameter space for successful single-sector reheating. The general set of Boltzmann equations for single sector reheating are [74,106] ρ ϕ + 3H(1 + w ϕ )ρ ϕ + (Γ s/f + Γ gr ϕϕ→RR )(1 + w ϕ )ρ ϕ = 0, ρ r s/f + 4Hρ r s/f − Γ s/f (1 + w ϕ )ρ ϕ = 0 (6) ρ r gr + 4Hρ r gr − Γ gr ϕϕ→RR (1 + w ϕ )ρ ϕ = 0.
Where ρ ϕ is the inflaton energy density. ρ r s/f corresponds to radiation energy density from direct inflaton decay to either scalar (s) or fermion (f ) via the coupling parameters (g r i , h r ) respectively. ρ r gr is the radiation energy density produced through universal gravitational scattering. For solving this set of equations numerically, we define dimensionless comoving variables Φ = ρ ϕ A 3(1+w ϕ ) /(m end ϕ ) 4 , and R rad = ρ rad A 4 /(m end ϕ ) 4 . Where R rad is the radiation produced from either direct decay or gravitational scattering from inflaton. The derivative is taken with respect to the cosmic time defined through the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric ds 2 = −dt 2 + a(t) 2 (dx 2 + dy 2 + dz 2 ). The normalized scale factor is defined as A = a/a end , with a end as the scale factor at the end of inflation, and the Hubble expansion parameter H =Ȧ/A. At the end of inflation i.e., A = 1, ρ ϕ = 3V (ϕ end )/2 and the energy densities of all the other components are set to be, ρ r s/f /gr = 0. Hence the appropriate initial condition for the above set of Boltzmann equations is, Where (ϕ end , m end ϕ ) are the values of the inflaton field and its mass at the end of inflation, which we define later.

B. Relating reheating and inflationary parameters through CMB
The connection between inflation and reheating parameters are established through the initial conditions Eq.8. We consider the well known α-attractor E-model [108][109][110] inflaton potential, where Λ is the mass-scale fixed by the CMB power spectrum, which is typically of the order 8×10 15 GeV, and the parameter (α, n) controls the shape of the potential. Through out our analysis we have taken α = 1, although our analysis is not much sensitive within the allowed range of α values from Planck and BICEP/Keck combined results [130,131]. Using the slow roll parameters ϵ(ϕ) , one obtains physically measurable quantities namely scalar spectral index (n s ) and the tensor to scalar ratio (r) defined at a particular pivot scale k, Another important quantities are inflationary e-folding number (N k ) and the Hubble constant H k defined for a particular pivot scale as, where A s represents the amplitude of the inflation fluctuation. Throughout our analysis, we assume the central value of A s ∼ 2.1 × 10 −9 from Planck [111]. Here ϕ k and ϕ end are the values of the inflaton field at the point of horizon crossing for a particular pivot scale and at the end of the inflation, respectively. ϕ end is obtained from the ending condition of the inflation ϵ(ϕ end ) = 1. The field and the potential value at the end of inflation take the following form Combining Eqs.10-12, one can find the field value ϕ k and the mass scale Λ as, After the end of the inflation, inflaton field oscillates around its minima, and the reheating phase begins. At the minimum of the potential, we expand the inflaton potential (Eq.9) in the limit ϕ ≪ M p as where β = 2/3αM 2 p . The field-dependent mass becomes, Using the envelope approximation at any instant of time, the envelope value of the field ϕ 0 must satisfy V (ϕ 0 ) ≃ ρ ϕ (t) [75]. Under this approximation, the inflaton mass can be written as Using the virial theorem, one can further calculate the equation of states (EoS) w ϕ as a function of the power of the inflation potential [75] w ϕ = n − 1 n + 1 .
Through the background dynamics and entropy conservation, we can connect the inflationary parameters defined above with the reheating parameters. The reheating period is effectively described by very few parameters viz. reheating temperature (T re ), reheating e-folding number (N re ), and the equation of state of inflation ω ϕ . In general, the end of reheating is marked at the point when the condition ρ ϕ (A re ) = ρ r (A re ) is satisfied. Where ρ r is the total radiation density constituted of all massless daughter fields. The reheating temperature T re can be expressed in terms of the radiation temperature T rad as where g ⋆ (T re ) is the effective number of relativistic degrees of freedom at the end of reheating, and we take g ⋆ (T re ) = 100 though out the paper. Combining the above two equations, we can get the one-to-one correspondence between the coupling constant (g r i ,h r ) and reheating temperature T re , where i = 1, 2 corresponding two different inflaton-Boson couplings described in the introduction. One can further obtain a constraint relation between reheating temperature T re and the present CMB temperature (T 0 ) by considering a physical assumption that from the end of reheating to the present time, the co-moving entropy density of the universe is conserved. This condition leads to [112], Where k/a 0 = 0.05 Mpc −1 is the CMB pivot scale, the present CMB temperature T 0 = 2.725 K and a 0 is the present day scale factor. Combining the above Eq.20 and 19, we can essentially obtain an indirect connection between the coupling constant (g r i ,h r ) and the inflationary spectral index n s , which parameterizes the anisotropies in the CMB fluctuations. In the next section, we moved our discussion to the possible constraints on the perturbative reheating scenario from the fragmentation of the inflaton that leads to self-resonance.

III. THE EFFECT OF SELF-RESONANCE IN REHEATING PARAMETER SPACE
Homogeneous oscillation of the inflation condensate can be unstable, leading to self-resonance. In Figs. 2 and 7, we have shown the region in the w ϕ versus coupling parameter space where self-resonance may be important. It has been pointed out that self-resonance can be sufficient to reheat the universe (except w ϕ = 0) without any coupling to the other fields with inflaton [128,129], however, that strictly depends on the values of α. When α << 1/6, self-resonance is efficient, and the RD universe is established within less than an e-fold of expansion after inflation end. However, for α > 1/6, self-resonance is not very efficient, and it takes many e-folds to give rise to radiation dominated universe. In the reference [129], the authors provide an estimation for the number of e-folds calculated from the end of inflation to the beginning of the radiation domination for the α-attractor E-model, where ∆k/k is the fractional width of the resonance band, d is its dimensionless "strength" and δ = 0.126 is a dimensionless number which is independent of n. This bound can be used as an upper bound on the transition duration between the inflation end and the radiation-dominated states of the universe. Self-resonance can significantly modify the dynamics for which a comprehensive analysis is required. Here, we have shown the region where it can play its role. If the perturbative reheating is completed with e-folding number < ∆N sr , our analysis will not be affected. Otherwise, self-resonance will significantly affect the parameter space, which needs to be taken into account, and we defer this for our future study. Using the bound from Eq. 21, in Figs. (2), (6), (7), and (10) we have shown regions where the self-resonance effect is important. We have found that self-resonance is effective when 0 < w ϕ ≤ 0.75 for α = 1.

A. Dynamics: probing different decay channels
During reheating, if the dominant decay channel of the inflaton is through massless Bosons, we call it bosonic reheating. For our phenomenological purpose, we discuss two different non-gravitational bosonic channels ϕ → ss and ϕϕ → ss of inflaton along with gravitational scattering process ϕϕ → h µν → ss. The total radiation component will be composed of all these thermalized decay products produced from different decay channels. Most of the studies on this single sector reheating were done at zero temperature, except for a few very special cases [101,103]. Therefore, a more realistic consideration would be to take into account finite temperature correction due to radiation baths. We compare the results with the zero temperature case in all the figures to quantify such effects. The detailed analytic calculation for both zero and finite temperature cases and the description of the dynamics are given in Appendix-A.
We solve the Boltzmann equations for the general inflaton equation of state w ϕ and scan the entire inflatonscalar coupling (g r i ) parameter space and figure out one of our most significant results shown in Fig.2. The parameter space (w ϕ , g r i ) can be generically divided into five regions marked in different colors: i) Light-cyan region is where reheating is entirely controlled by the gravity-mediated decay channel (gravitational reheating), ii) Yellow region is controlled by mostly inflaton-Boson coupling, iii) Light-red region is where successful reheating cannot be achieved as reheating temperature T re < T BBN , iv) Pink region where initial parametric resonance will be the important and v) Light-blue region where the self-resonance of the inflaton will be important. In our present paper, we ignore both effects and defer for our future study. The system has two competing effects due to direct and gravity-mediated inflaton decay. Based on their relative dominance, we observe three distinct regions of coupling g r i where reheating evolution will be different: 1) Case-I: Entire reheating dynamics will be dominated by direct inflaton decay, 2) Case-II: Both the decay processes will partially dominate the reheating dynamics, 3) Case-III: Entire reheating dynamics will be dominated by gravity mediated decay (gravitational reheating [96]). These three cases immediately suggest the existence of two critical coupling values for every individual inflaton-scalar radiation coupling g r i , which set the phase boundaries among those three cases in (w ϕ , g r i ) plane. If the coupling g r i > G 1, th ci ≃ G 1 ci , the reheating evolution will be according to case-I. Where, G 1 ci is computed without thermal effect, for it turned to be same for both with and without thermal effect (see detailed derivation in the Appendix-A). The critical coupling strength G 1 ci is found to be, If we lower the couplings g r i below G 1, th ci , the gravitational scattering starts to reveal its presence in the early phase of the reheating process, and the complete takeover happens if the non-gravitational coupling strength is lower than a new critical coupling which we denoted as G 2, th ci . Therefore, if coupling strength in between G 2, th ci < g r i < G 1, th ci , the reheating evolution will be according to case-II. The expressions of G 2, th ci for different interaction are calculated as, where, A gr re is the scale factor defined at reheating end for the gravitational reheating scenario (see, for this instance, Eqn.A17).
Finally, as pointed out just above, if g r i < G 2, th ci , the reheating evolution will be according to case-III, which we call gravitational reheating. Detailed analysis on this possibility have been discussed in [96]. We will now dwell on these three cases and discuss their thermal histories in detail: Case-I: Coupling strength g r i > G 1, th ci : In this regime, direct decay of inflaton into radiation controls the entire reheating process. In the left panel of Fig.4, we showed the evolution of the different energy components with the coupling parameter. Since the radiation bath of temperature T rad is produced from the decay products of homogeneous inflaton background, the typical energy of the bath particles will be of the order of inflaton mass m ϕ . And hence for the condition T rad > m ϕ (t), the thermal effect will be dominant. For any reheating dynamics, there exist two important energy scales of importance, and those are maximum radiation temperature (T max rad ) and the reheating temperature (T re ). Given the inflation model under consideration, we have two free parameters namely, the inflaton equation of state (ω ϕ ) and the inflaton-Boson coupling (g r i ), where "i" stands for two different bosonic decay channels mentioned earlier. Depending upon the evolution of (T rad , m ϕ ), and consequently the behavior of thermal effect, we have observed rich reheating histories. In the following section, we lay bare the detailed discussions on those for different cases in different temperature regimes.
When T max rad > m end ϕ : This is the situation which typically occurs for large value of inflaton-scalar coupling mostly in the pink region of Fig.2. Since the maximum radiation temperature T max rad is greater than the inflaton mass m end ϕ defined at the end of inflation, the thermal effect influences the reheating dynamics significantly. Details of this finite temperature effect will be further controlled by the parameters (w ϕ , g r i ) and the timedependent inflaton mass m ϕ (t). We will first discuss the situation when T rad > m ϕ (t) throughout the entire period of reheating. However, important to remember that such a condition does not satisfy though out the entire reheating parameter range, and this can be observed from Fig.2.
It is observed that in this region the ratio T rad /m ϕ (t) varies as A (9w ϕ −1)/2 (A (11w ϕ −3)/2 ) for ϕ → ss(ϕϕ → ss). Such variation of the ratio indicates that there exists a critical value ω c ϕ = 1/9 for the decay channel ϕ → ss and ω c ϕ = 3/11 for the decay channel ϕϕ → ss, above which T rad > m ϕ (t) condition is always maintained. With this condition, the finite temperature decay widths can be approximated as, T BB N = 10 -2 G e V  figure) as a function of w ϕ . Dashed and solid lines correspond to without and with the thermal effect being taken into account in the decay process respectively. The yellow and pink shaded regions indicate the explicitly coupling-dominated region where the decay channel controls the reheating temperature. The light-cyan region corresponds to gravitational reheating. The light-red region corresponds to Tre < TBBN ≃ 10 MeV. The light-gray region corresponds to the no reheating region where inflaton energy density falls faster than radiation energy density, and successful reheating is not impossible. The pink region corresponds to the non-perturbative regime where bounds on are obtained from resonance condition of Mathieu equation for scalar field [113][114][115]. The light-blue region corresponds to the self-resonance domination region.
and with this the reheating temperatures are estimated for two different decay channels as, Where, A re is the normalized scale factor at the end of reheating (detailed derivations and expressions can be found in Appendix-A). To this end we would like to point out an interesting fact associated with the maximum radiation temperature T max rad , which generically satisfies T max rad > T re . We observed the existence of a critical value of w t ϕ = (1/3, 3/5) for two different bosonic decay channel ϕ → ss and ϕϕ → ss (see, for instance, Eqs.A23 and A24). If w ϕ < w t , the maximum radiation temperature (T max rad ) satisfies the usual condition mentioned above, Where T r, max s indicates the maximum radiation temperature. Surprising result emerges, however, for w ϕ > w t ϕ case, for which the evolution of radiation and background conspire in such a manner that at the end of reheating maximum radiation temperature becomes equal to the reheating temperature, T max rad ≃ T re (for better visualization, see the left most plot of Fig.4). Such behavior has been observed before considering the phenomenological expression of the decay rate as a function of temperature [116]. The implication of this specific case could be interesting to study. Again, when EoS stays within 0 ≤ w ϕ < w c ϕ , due to initial high radiation temperature thermal correction will have a significant effect. As the reheating progresses, such effect diminishes with the complete takeover by the zero temperature dynamics at a certain value of scale factor A c , which depends on the inflaton equation of state as follows, where A max and ρ r, max s are defined in Eqs.A26,A25. After this crossover happens, the radiation energy density simply follows Eqs.A10. In the DM studies, we will see such an intermediate scale will have non-trivial dependence on its abundance. We find the associated reheating temperature as where In reheating model building, gravitational contribution to the radiation sector is universally present along with the non-gravitational one. Therefore, it is natural to expect that the T max rad for a given reheating model can not assume an arbitrarily low value. In fact, due to universal gravitational contribution, there exists a lower limit on T max rad , which is set by the gravitation reheating T r, max gr ≃ 10 11 → 10 12 GeV [53,96]. The small variations are due to different values of ω ϕ . Therefore, the minimum possible value of T max rad simply turns out as T r,max gr . In the following discussion, we now consider the regime where T max rad > T r, max gr but less than m end ϕ . When T r, max gr < T max rad < m end ϕ : The parameter region (see Fig.2) wherein this condition is satisfied belong to the yellow region. However, for this case, initially T max rad < m end ϕ , and hence there is no thermal effect initially, and the ratio T rad /m ϕ ∝ A −(3−27w ϕ )/8 (A −(9−33w ϕ )/8 ) for ϕ → ss(ϕϕ → ss) respectively. Thus, for w ϕ < w c ϕ (defined before), the finite temperature effect will never be significant. Consequently, the reheating dynamics will be the same as that of the zero temperature, and details of such dynamics are described in the Appendix-A (for example, the radiation energy density evolves following Eq.A10). However, if the inflaton equation of state satisfies the condition w ϕ > w c ϕ , the finite temperature effect (T rad > m ϕ ) is expected to occur at some intermediate radiation temperature T c = T rad (A c ) with the scale factor A c during reheating. We have where A max and ρ r, max s are defined in Eqs.A13, A14. After this crossover happens, the radiation energy density simply follow the Eqs.A23 and A24. In the DM studies, we will see such an intermediate scale will have non-trivial dependence on its abundance.
For case-I, above two temperature regimes will be possible. To this end let us point out an another important situation that deserves detailed discussion is a special case when T max rad = T r, max gr , which indicates that the initial  phase of reheating must be dominated by the graviton mediated inflaton decay. And such situation arises only when non-gravitational inflaton coupling g r i < G 1, th ci . This condition, therefore, belongs to the other two cases of coupling ranges mentioned before.
Case-II: Coupling strength in between G 2, th ci < g r i < G 1, th ci : In this coupling range, gravitational interaction drives the dynamics of reheating at the beginning. The gravitational production is nearly instantaneous and happens just at the beginning of reheating. In fact, this case is true for the entire range of coupling with g r i < G 1, th ci , and T max rad = T r, max gr condition always holds. However, since the non-gravitational bosonic coupling is non-zero, during the reheating process gravitational coupling, non-gravitational coupling, and thermal effect of the produced radiation undergo interesting interplay among themselves. Let us illustrate the following two different possibilities of thermal history in this context, • The thermal effect may start dominating the early phase when gravity-mediated decay controls the re-heating. During this phase the ratio behaves T rad /m ϕ (A) ∝ A 3w ϕ −1 . Hence, for w ϕ > 1/3, it is clear from this ratio that the thermal effect cannot be ignored. And we found a particular value of scale factor A g c , after which Bose enhancement starts affecting the dynamics, where, A max is the scale factor at which T rad = T max rad = T r, max gr , and maximum radiation energy density (ρ r, max gr ) is obtained from gravitational decay (see, for instance, the last expression of Eq.A13 and Eq.A15). After this point, the radiation energy density simply varies as A −4 . However, as reheating proceeds towards the end, there is another crossover from gravitational decay domination to non-gravitational decay domination, which will happen at the value of scale factor, defined as where C(w ϕ ) = . The reheating temperature assumes the same form as Eq.25. For better visualization of the evolution of different energy components see the middle panel of Fig.4. In summary, the reheating dynamics can be read off as follows: at the beginning gravitational sector dominates the porcess with radiation temperature varies as T rad ∝ A −1 → as reheating proceeds the thermal effect starts to play its role but with the same temperature evolution T rad ∝ A −1 → non-gravitational coupling takes over the process, and the radiation temperatures vary as for the decay process ϕ → ss (ϕϕ → ss)).
• There may be a situation where the thermal effect will be important only during non-gravitational production. For this case, reheating proceeds from gravity-mediated decay to explicit inflaton decay domination, and the transition occurs at the scale factor, where ρ r, max gr and ρ r, max s are defined in Eqs.A14, A15. Once reheating process starts to dominate by the explicit inflation coupling, the scale factor beyond which the thermal effect starts working is followed by Eq.30. Finally, the decay channel defines the reheating temperature (see, for instance, Eq.25). In summary, dynamics can be described as follows: reheating proceeds through gravity-mediated decay with no finite temperature effect (T rad ∝ A −1 ) → non-gravitational decay dominates the phase with negligible thermal effect with radiation temperature varies as for decay process ϕ → ss (ϕϕ → ss)) → non-gravitational coupling domination with significant thermal effect ( Case III: when g r i < G 2 ci : The bath temperature always falls as A −1 , and thermal effect is observed to play no role throughout. The gravity-mediated decay of inflaton controls the entire dynamics of reheating, and the scenario is termed as gravitational reheating, ant that will occur only for w ϕ > 0.65 (see the light-cyan region of Fig.2). The reheating temperature is defined when ρ ϕ = ρ r gr , and the condiction gives  where A gr re is the normalized scale factor at the end of gravitational reheating. In the right panel of Fig.4, we have shown the dynamical behaviour of different energy components, and in Table-I, showing the evolution of the bath temperature with A for non-gravitational reheating and gravitational reheating.
In the subsequent section, we will focus on the possible constraints on the inflaton coupling strengths depending on the CMB (n s ) and reheating parameter T re .

B. Inflaton phenomenology: Constraining reheating and bosonic decay parameters:
For illustration, we consider five different values of the inflaton equation of state w ϕ = (0, 0.2, 0.5, 0.82, 0.99). For each w ϕ , we have plotted :(i) T re vs n s , (ii) g r i vs n s , (iii) g r i vs T re . We compare the results with and without the thermal effect for all cases.
i) Reheating temperature (T re ) in terms inflationary (CMB) parameter (n s ): From the Fig.5 we observe that the evolution of reheating temperature in terms of the inflationary scalar spectral index (n s ) is insensitive to the finite temperature correction of inflaton decay width. Such behavior of the reheating temperature has already been reported in [101]. The generic feature is that for w ϕ < 1/3, (see Fig.5), T re increases with increasing n s , and as a consequence the reheating e-folding number N re decreases with n s . This indicates the existence of a maximum scalar spectral index n max s corresponding to the maximum reheating temperature T max re = 10 15 GeV and that is called instantaneous reheating. Similarly, the minimum reheating temperature T re = T BBN ∼ 10 MeV [119][120][121], corresponds to a minimum allowed spectral index n min s for a given the inflaton equation of state w ϕ . On the other hand, for w ϕ > 1/3, one finds the opposite feature: maximum T re corresponds to the minimum spectral index n min s and vice versa. When the reheating phase is dominated by purely gravitational interaction, the minimum possible reheating temperature fixes the maximum possible value of n s for the equation of state w ϕ ≥ 0.65. For example, as shown in the Fig.5, for ω ϕ = (0.82, 0.99), we obtain T min re ≃ (10 3 , 10 6 ) GeV, respectively. From Fig.(5), it is clear that thermal feedback to the decay rate does not affect the variation of reheating temperature with n s . In Table II,      (ii) Constraining inflaton couplings with bosonic radiation (g r 1 , g r 2 ): One of the most important findings of our present analysis is illustrated in Fig.2. The figure clearly depicts different regions in the parameter space of (w ϕ , g r i ), where the effect of different inflaton decay channels on the reheating process can be understood. At this point, let us reiterate different regions again: i) the light-cyan region is where reheating is entirely controlled by the gravity-mediated decay channel (gravitational reheating), ii) the yellow region is controlled by mostly inflaton-scalar coupling, iii) the light-red region is where successful reheating cannot be achieved as reheating temperature T re < T BBN , and iv) pink region where initial parametric resonance will be important which we ignored in this paper. Furthermore, in the right panel of Fig.2, there is a light-gray region where reheating is not possible for the decay process ϕϕ → ss. For this process, if the thermal effect is subdominant at the beginning (m end ϕ < T max rad limit), the ratio between inflaton and radiation energy density varies as (1−5w ϕ ) (see, for instance, Eqs.A1 and A10), and hence if w ϕ < 1/5, the universe will always be inflaton dominated irrespective of the value of inflaton-scalar coupling g r 2 .
In addition to that, if the thermal effect starts dominating from the beginning (m end ϕ > T max rad limit), the ratio varies as for instance, Eqs.A1 and A10) which implies that if w ϕ < 3/13 ∼ 0.23 achieving radiation domination is not possible. However, for extremely large coupling, parametric resonance may have some effect.
From Fig.2, a generic feature can be observed, and that is related to the monotonic decrease of g r i with w ϕ for a fixed reheating temperature T re . The reason behind this behavior can be understood as follows: with increasing w ϕ , inflaton energy density dilutes faster, and hence to achieve the reheating condition ρ ϕ = ρ r s , one needs to lower the coupling. Furthermore, m ϕ (t) decays faster with increasing w ϕ , and for both types of bosonic decay channels (ϕ → ss and ϕϕ → ss), the production rate goes as ∝ 1/m ϕ (t), which will boost up the production. As a result, to keep reheating temperature fixed, one needs to lower the value of g r i again. Due to very nature of the bosonic particles the finite temperature correction in the decay width enhances the particle production rate from the inflaton condensate. As discussed, this physical fact is imprinted in the reheating dynamics and is further reflected in the parameter plot shown in Fig.6. Finite temperature correction naturally increases the effective decay width of the inflaton to scalar radiation, and consequently, one needs to lower the values of the dimensionless coupling parameters g r 1 =g r 1 m end ϕ and g r 2 as compared to their zero temperature case to have successful reheating. This can be observed in Fig.6 both with respect to reheating temperature (T re ) (lower two plots) and CMB spectral index n s (upper two plots).
The maximum limiting value of the coupling parameters will naturally be set by the maximum possible reheating temperature T max re ≃ 10 15 GeV, where all the lines converge (see Fig.6). If the reheating dynamics are controlled directly by the inflaton-radiation coupling, the minimum possible value of the coupling will be set by the minimum reheating temperature. However, such a limit on the inflaton coupling is observed to be dependent on the finite temperature correct, which will be discussed in detail. When the radiation temperature T rad ≫ m ϕ (t), the thermal effect significantly influences the radiation dynamics and consequently affects on the possible constraints on the coupling parameter as compared to the zero temperature case. It can be observed that higher the value of ω ϕ , more will be the effect of finite temperature correction on the thermal bath. For w ϕ = 0, the effective mass of the inflaton m ϕ (∼ 10 13 ) remains constant; as a result, the thermal effect manifests (see left two plots of Fig.6) itself only very near and above the reheating temperature ∼ 10 13 GeV (or for n s > 0.9645). On the other hand, for w ϕ > 0, the rate of decrease of effective inflaton mass m ϕ ∝ ∂ 2 ϕ V increases with increasing ω ϕ such that the condition T rad > m ϕ (t) becomes easier to satisfy even at a lower temperature. For example, for w ϕ = 0.2, the above condition begins to satisfy (see left Fig.6) when n s > 0.9639 (T ∼ 10 8 ) GeV, and accordingly, the finite temperature effect (solid line) manifest itself after T re ∼ 10 8 GeV. For w ϕ = 0.50, 0.82, 0.99, the condition T rad > m ϕ (t) can be observed to satisfy throughout the whole range of reheating temperature.
To this end, we would like to elaborate on the finite temperature effect for low reheating temperatures. The author in the reference [101] claims that the thermal effect will be insignificant at low reheating temperatures. However, generically such an effect depends on the evolution of the ratio m ϕ /T rad . And we found that the finite temperature effect can be significant at low reheating temperature for the higher equation of state w ϕ = (0.20, 0.50, 0.82, 0.99) (see solid and dotted lines in Fig.6), for which inflaton mass undergoes non-trivial evolution. As can be seen from the second row of Fig.6, for the higher equation of state finite temperature effect becomes more prominent at lower reheating temperature mainly because inflaton mass can become significantly smaller during the course of reheating. Moreover, for higher equation state (w = 0.82, 0.99), the gravitational reheating has been observed to give a lower limit of the reheating temperature (10 3 , 10 6 GeV), which is again found to be directly corresponding to specific inflationary scalar spectral index n s . The spectral indices associated with those temperatures are n s = 0.96855 and 0.9681 for w = 0.82 and 0.99, respectively. When n s reaches these values, the coupling parameter tends towards zero, i.e., gravitational scattering solely controls the reheating dynamics.

A. Dynamics: probing different decay channels
During reheating, if the dominant decay channel of the inflaton is through massless fermions, we call it fermionic reheating. For this purpose, we consider inflaton decaying only into massless Fermion though the standard Yukawa decay channel ϕ →f f along with gravitational scattering process ϕϕ → h µν → ss. Instantaneous thermalization of those different components are assumed throughout. To study the evolution of the radiation energy density, we took the finite temperature effect arising due to Pauli blocking (see, for instance, Appendix-B for details calculation). Similar to bosonic reheating, for the fermionic case, we identified distinct regions in (w ϕ , h) plane depending upon different physical processes involved in controlling the reheating dynamics (see Fig.7). For this case, we have plotted separately with and without the finite temperature effect and observed the quantitative change in parameter space due to the finite temperature effect where reheating would be successful. The parameter space (w ϕ , h) is again divided into four regions marked in different colors: i) Light-cyan region is where reheating is entirely controlled by the gravity-mediated decay channel (gravitational reheating), ii) Yellow region is controlled by mostly inflaton-Fermion coupling, iii) Light-red region is where successful reheating cannot be achieved as reheating temperature T re < T BBN , and iv) Pink region signifies initial parametric resonance domination which we ignored in this paper. As we discussed for bosonic reheating, based on whether gravitational or non-gravitational sectors dominate the dynamics, we have a fermionic critical coupling H c which sets a boundary for two distinct scenarios: 1) Case-I: The entire reheating dynamics is dominated by the Yukawa coupling. For this case, the fermionic coupling parameter is in the range h r > H c . The critical coupling H c is identified by equating the maximum energy densities from non-gravitational and gravitational sector ρ r, max f = ρ r, max gr . One can find the expression for the critical coupling for fermionic reheating as The maximum energy densities for both gravitational and non-gravitational sectors appear at the initial stage of reheating, and as the radiation bath associated with gravity mediated process dilutes faster than the nongravitational one, for h r > H c reheating dynamics always have explicit fermionic coupling domination.
2) Case-II: For this case coupling parameter satisfies h r < H c and both sectors partially dominates the reheating phase. However, when EoS w ϕ lies above 0.65, gravity-mediated decay controls the reheating phase termed as gravitational reheating. In our succeeding discussion, we will discuss these two cases in detail: Case-I: Coupling strength h r > H c : In this coupling regime, non-gravitational decay of inflaton into fermionic radiation controls the entire reheating process. In order to give analytical estimation and compare our result with the zero temperature scenario, we consider two separate regimes of the maximum radiation temperature (see Fig.8 for its depiction), and those are as follows: When T max rad > m end ϕ : If the maximum radiation temperature is greater than m end ϕ , the coupling parameter associated with this region entirely lies in the pink region of Fig.7 (non-perturbative resonance-dominated region). Due to strong inflaton-Fermion coupling, the gravitational channel is naturally subdominant throughout the reheating, and the thermal effect is non-negligible from early stage of reheating. However, its effectiveness through out the reheating process will depend on how T rad /m ϕ evolves. Since T rad > m ϕ at the initial stage, decay width can be approximated as, With the aforementioned decay width, one can estimate the behavior of radiation energy density, The above equation suggests that the evolution of the radiation component is entirely different in two different regimes • w ϕ > 7/15 : Most of the production occurs at the initial reheating stage, and temperature decreases with the scale factor as A −1 . Since the ratio T rad /m ϕ behaves as A 3 w ϕ −1 , T rad always remains greater than m ϕ , and hence the finite temperature effect will be significant till the end of reheating. The reheating temperature for this case assumes the following form, • 0 ≤ w ϕ < 7/15 : For this case, the ratio T rad /m ϕ ∝ A 3 10 (5w ϕ −1) induces two different evolution history with regard to the finite temperature effect. It turns out that when EoS stays within 0 ≤ w ϕ < 1/5, due to initial high radiation temperature thermal correction will have a significant effect. As the reheating progresses, such effect diminishes with the complete takeover by the zero temperature dynamics at a certain value of scale factor, which depends on the inflaton equation of state as follows, After this intermediate scale factor, the dynamics is governed by the zero temperature decay channel following the Eq.B3 (see without thermal effect section of Appendix-B for details calculation) till the end of reheating and eventually equating ρ ϕ = ρ r f , we find the associated reheating temperature as, On the other hand when w ϕ is in between 1/5 < w ϕ < 7/15, the thermal effect will be non-negligible throughout the entire reheating history, and we have the associated reheating temperature . When T r, max gr < T max rad < m end ϕ : For this case, the coupling parameter mostly lies in the pink region of Fig.7. Similar to the previous case, whole reheating dynamics are governed by non-gravitational coupling. Since T rad < m ϕ at the initial stage, the thermal effect is minimal. As reheating progresses, depending on the ratio T rad m ϕ , the thermal effect may state dominating the dynamics. Thus initially, the radiation component evolves as The aforementioned equation clearly suggests that the radiation component behaves differently for w ϕ < 5/9 and w ϕ > 5/9. Let us discuss these two cases, • w ϕ > 5/9 : Maximum production happens initially, and the bath temperature falls as A −1 . For this case, T rad /m ϕ ∝ A 3 w ϕ −1 and hence finite temperature effect will be important near the final stage of reheating, and the reheating temperature can be expressed as, • 0 ≤ w ϕ < 5/9 : For this case, T rad /m ϕ ∝ A −3/8(1−5w ϕ ) implies two different evolution histories depending on the value of equation state greater or less than 0.2. For 0 ≤ w ϕ < 0.2, the thermal correction will be subdominant, and the reheating temperature can be simply read off from Eq.40. Whereas for EoS w ϕ > 0.2, the thermal effect will start to dominate at some intermediate time within the reheating phase, which we call the crossover point, and reheating temperature is given by Eqn.41.
Case-II: Coupling strength h r < H c : In this coupling regime, the maximum temperature is always controlled by the gravitational sector T max rad = T r, max gr . This condition generally satisfies within the entire allowed region shown in Fig.7 except the parametric resonance dominated region shaded in pink. Evolution of the different energy densities in two different regimes h r > H c and h r < H c with two distinct values of inflaton equation of state w ϕ (1/3, 0.82) are shown in Fig.9. Depending on the inflaton equation of state, here also we have the following three different scenarios.
• w ϕ > 0.65 : For this case, the gravitational sector governs the entire reheating phase, and we termed this as gravitational reheating. The parameter space where this condition is met is shaded in light cyan in both the Figs.7. The reheating temperature can be followed from Eq.A17, which depends only on the reheating equation of state.
• 5/9 < w ϕ < 0.65 : This case turned out to be within the light red shaded region in the (h, w ϕ ) plane of Fig.7. As the figure suggests, reheating temperature evolved into below BBN temperature, which does not support the standard cosmological constraints.
• 0 ≤ w ϕ < 5/9 : In this case, the competition between two sectors of production, along with the finite temperature effect, leads to two different physically distinguishable reheating dynamics. In the (w ϕ , h) plane, the condition under consideration lies in the light yellow region of Fig.7. Here we have two different possibilities depending on how the thermal effect plays its role during the reheating history. As discussed for the bosonic reheating case 1) The thermal effect starts to influence the reheating dynamics in its early stage (T rad > m ϕ ) during the gravitational decay of inflaton. For this case, the behavior of the radiation component during nongravitational sector domination is simply followed by Eq.37, and we have reheating temperature as in Eq.41.
2) The thermal effect starts dominance during the later stage of reheating when it is governed by nongravitational inflaton decay. The scale factor associated with the point where the thermal effect starts to influence the dynamics can be the same as the Eq.44, and reheating temperature is given by Eq.41.

B. Inflaton phenomenology: Constraining reheating and fermionic decay parameters
Similar to the bosonic reheating case, to illustrate our results in terms of inflationary parameter n s and to see how the coupling strength behaves as a function of reheating reheating temperature, we have taken five different sample values of w ϕ = (0, 0.2, 0.5, 0.82, 0.99) and compare the results for with and without thermal effect.
i) Reheating temperature (T re ) in terms inflationary (CMB) parameter (n s ): The qualitative relation between reheating temperature and inflationary parameters remains similar to that of the bosonic case discussed earlier. Therefore, the possible bounds on inflationary parameters such as spectral index n s , the maximum inflationary e-folding number N max k do not depend on the details of the reheating dynamics but only the reheating temperature. As a result, bounds on the inflationary parameters remain the same as bosonic reheating (see, for instance, Table-II). This can easily be read off from the left plot of Fig.5. In addition, Fig.(5), further indicates that the thermal feedback on the decay process does not affect the reheating temperature variation with n s .
(ii) Constraining inflaton couplings with fermionic radiation (h r ): For the fermionic reheating, the parameter space in Fig.7 illustrates different regions in (w ϕ , h r ) plane, where the effect of different inflaton decay channels on the reheating process can be read off. We have given two plots with and without finite temperature effects for clear depiction.
An interesting distinction can be observed as compared to the scalar reheating case is that for fermionic reheating inflaton-Fermion coupling h r does not vary monotonically with respect to w ϕ given a fixed reheating temperature. There exists a critical value of w ϕ ≃ 7/15(5/9) for with (without) thermal effect, below which  6), the only difference is that here we have plotted for fermionic reheating and the self-resonance dominated region is same for both w ϕ = 0.5 and 0.2 shaded by orange-shaded region.
one requires a higher h r value for a higher equation state for a fixed reheating temperature. And this can be understood from the behavior of Fermion production rate ∝ Γ ϕ→f f ρ ϕ ∝ (h r ) 2 m ϕ (t)ρ ϕ . With increasing w ϕ , the effective mass of the inflaton decreases faster with time; hence, to achieve a fixed reheating temperature, h r needs to be enhanced. However, this simple physical argument is no longer tenable after w ϕ ≥ 7/15 (5/9) for with (without) thermal effect. For such cases, most of the production happens initially, and the radiation energy density simply dilutes as A −4 , which is slower than that of the inflaton energy density. In such a situation with increasing w ϕ and fixed reheating temperature, we need a lower the value of h r to satisfy the reheating condition ρ ϕ = ρ r f . Due to its intrinsic nature, the finite temperature Fermion bath diminishes its own production rate from the inflaton condensate. Consequently, for successful reheating one needs higher values of the dimensionless coupling parameters h r as compared to the zero temperature case. This can be clearly observed from Fig.10, and such behavior is opposite to that of the scalar reheating case. The qualitative behavior of the fermionic coupling in terms of the spectral index and reheating temperature are the same as that of the scalar reheating case. For example, for w ϕ = 0 the coupling parameter h r with thermal effect begins to affect only at very high temperatures at around ∼ 10 13 GeV, and in terms of the spectral index, the deviation is visible for n s > 0.9645. On the other hand, since the effective mass of the inflaton varies as A −3w ϕ (see, for instance, Eq.A7), for w ϕ > 0, m ϕ (t) decrease faster with increasing ω ϕ and the T rad > m ϕ (t) condition becomes easier to fulfill even at a lower radiation temperature. As a result, for w ϕ = 0.2, the thermal effect starts dominating even at smaller radiation temperature T rad ≥ 10 8 GeV when n s > 0.9639. The situation is entirely different though for w ϕ > 5/9. For EoS greater than 5/9, the maximum radiation production happens initially, and hence thermal effect will be dominant from the beginning, which indicates maximum radiation temperature T max rad > m end ϕ . From Fig.10, we can clearly see that for EoS w ϕ = (0.82, 0.99), the deviation between the results for with and without thermal effect start visible when T max rad ∼ m end ϕ .

VI. FIMPS AND WIMPS DURING REHEATING AND OBSERVATIONAL CONSTRAINTS
Discussion on DM will be considered in three parts. In the first part, we discuss the production of DM exclusively from the inflaton through non-gravitational and gravitational interaction. We mainly point out the constraints on the inflaton-DM coupling and DM mass from both theory and observation. Since it is produced solely from the inflton decay, we call these as FIMP like DM. In the second part, we discuss the production from the thermal bath assuming an effective radiation to DM annihilation cross-section ⟨σv⟩, added with the universal gravitational production discussed in the first section. For this case, we will have both the Freeze-in and Freeze-out production scenarios depending upon the strength of ⟨σv⟩. The DM produced due to thermal Freeze-out from the radiation bath will be generally called WIMPs. On the other hand, DM produced by the decay process from the radiation bath via the Freeze-in mechanism will be called FIMPs. In the third part, we discuss about the experimental constraints on various reheating and DM scenarios. ∆N eff = extra radiation energy density(ρ DM ) energy density of single SM neutrino species(ρ ν ) T =T BBN = 43 7 We will be discussing two production mechanisms. For FIMP like DM, we intend to separately discuss its production from the direct inflaton decay and radiation bath. For inflaton decaying into DM, production freezes during reheating or at the end of reheating, depending on the decay channels and DM mass. On the other hand, for DM from the thermal bath, its production rate crucially depends on the radiation production rate. Therefore, in this case, also freeze-in occurs mostly during or at the end of reheating, depending on DM mass. Overall, for the freeze-in mechanism, we, therefore, can always express If the DM is relativistic after Freeze-in, both radiation (ρ rad ) and relativistic DM energy density (ρ DM ) fall as a −4 . The ratio ρ DM /ρ rad , therefore, stays constant between reheating and BBN. Hence, the above equality Eq.46, holds true generically for any FIMP scenario. For WIMP, on the other hand, the situation becomes very different but simpler. For such cases, till it freezes out, DM remains in equilibrium with the thermal bath. Therefore, any relativistic DM being in the thermal bath at the time of BBN always behaves like an additional degree of freedom. Therefore, ∆N eff naturally transforms into [126,127] ∆N eff = 4 7 j x = 0.571 scalar DM 1.14 fermionic DM Where j x is the DM's intrinsic number of degrees of freedom, WIMPs mass lighter than T BBN ∼ 10 MeV always behaves like dark radiation at the time of BBN. Hence, Therefore, all the WIMPs of mass m x ≤ 10 MeV violate the BBN bound of ∆N eff (see Eq. (47)).

B. Freeze-in production of dark matter from inflaton decay and constraints
Similar to massless radiation production, we will now discuss DM production during reheating, considering various decay channels for both scalar and Fermion DM. As discussed in the introduction, we considered two categories of production channels: (i) DM production from inflaton through gravitational scattering and (ii) production through explicit coupling-dependent decay channel.
The governing equation for the DM number density (n x ) produced from direct inflaton decay takes the following formṅ where Γ ϕ x is the inflaton decay width to DM and ⟨E x ⟩ ϕ the average energy of the DM particles.

freeze-in production via Gravitational interaction
Gravitational freeze-in production of DM is universal in nature and hence will always be present in any inflationary scenario. In this work, we have considered scalar (S) and Fermion (F) DM (see Fig.11 for relevant Feynman diagrams for gravitational scatting from inflaton). The decay rates associated with the gravitational production are [50,51] Since, DM are feebly coupled with the radiation bath, the thermal correction to the decay width will be unimportant. Using the above equations in Eq.48, we have obtained the following solutions for the number density, Therefore, the gravitational contribution to the DM abundance is calculated as 2 Where the suffix "g" stands for production due to the gravitational scattering process. Ω r h 2 ≃ 5 × 10 −5 is the present value of the radiation abundance. It is clear from the expression above that gravitational production depends on H end , m end ϕ , and DM mass m x . Given an inflation model, inflationary parameters such as H end , m end ϕ are fixed by CMB observation. Therefore, w ϕ and DM mass m x are the only free-parameter. Therefore, the present DM abundance can completely fix the DM mass once a particular inflaton equation of states w ϕ is assumed. We will later observe that the aforesaid mass will set a maximum possible limit on the DM mass, which we symbolized as m g,max x , for a large range of coupling and the inflaton equation of state. To this end, let us reiterate that due to its universal nature, the gravitational contribution to DM must be added to all the additional production processes we take up in the following sections.

freeze-in production via inflaton direct decay
We introduce different inflation coupling to DM. We consider three types of possible interaction: g x 1 ϕS 2 , g x 2 ϕ 2 S 2 , h x ϕF F (see Fig.11 for relevant Feynman diagrams), and corresponding decay widths are, Using the above equation in Eq. (48), we have obtained the following solution of number density, Unlike the previous gravitational production case, we now have additional coupling parameters g x i , h x in the problem. Therefore, the present DM abundance will provide the constraint equation between (m x , g x i /h x ) once we fix a particular reheating history by fixing (w ϕ , T re ) and g r i /h r . Depending upon the DM mass, we will have two different expressions for the DM abundance. If m x > m ϕ (A re ), the DM freezes in before the end of reheating at some intermediate scale factor A = A re (m x /m ϕ (A re )) −1/3w ϕ for w ϕ ̸ = 0, and, if m x < m ϕ , the DM freezes in after the end of reheating. The contribution to the DM abundance for different direct decay channels are calculated as ( m x < m ϕ (A re )), We introduce a new symbol, M ϕ = 2n(2n − 1)βΛ 2/n ϵ w/1+w and m ϕ = M ϕ T is the inflaton mass defined at the end of reheating. Note that for ϕϕ → SS (with w ϕ < 0.2) DM production channel, most of the DM production happens at the initial phase of the reheating similar to the gravitational production. On the other hand if m x > m ϕ (A re ), we have Once we obtain products from the direct decay, the total DM abundance can be expressed as In order to constrain the coupling parameters, we consider the total DM abundance.
3. Constraining the inflaton-DM couplings (g x i , h x ) and DM mass (mx) Fig.(12) depicts detailed allowed parameter space for which present DM abundance is satisfied. Let us emphasize again that while plotting the DM abundance, we take into account the contribution from both the direct and the universal gravity-mediated decay of inflation. For the sake of presentation, we have considered two sample values of the inflation equation of states w ϕ = (0.0, 0.50) with reheating temperatures range between maximum and minimum values T re = (T BBN , 10, 10 6 , 10 10 , 10 15 ) GeV.
In each plot for a fixed (T re , w ϕ ) we see the maximum limit on DM mass (m g,max x ) [52,53] which is due to gravitational interaction as mentioned before. However, since DMs are produced from the inflaton decay, purely kinematic constraints can also set the upper limit to be m end ϕ for some specific cases when m g,max x > m end ϕ which is observed for w ϕ = 0 (see Fig.12). It is natural to expect that for m x < m g,max x , the DM from the decay channel solely controls the abundance. Upon m x approaching m g,max x value, the gravitational contribution starts dominating the abundance till m x = m g,max x . However, the lower bound on the DM mass will be fixed either from observation or theory. For example, the lower bound has been observed to be controlled by reheating temperature and the physical processes of reheating under consideration. In general, for this scenario lower the DM mass, the larger would be the inflaton-DM coupling to achieve the correct DM abundance. However, the inflaton-DM coupling should be bounded from above so that the universe should be radiation dominated from the end of reheating (T re ) to the matter-radiation equality (T eq ≃ 10 −9 GeV). Hence, there exists an upper limit of the DM coupling, above which we always get the DM-dominated universe after the inflaton domination, or in other words one never achieve the radiation-dominated universe. Since there is a one-to-one correspondence between the DM coupling and the DM mass, a lower limit of the DM mass corresponds to the upper limit of the DM coupling, and that can be obtained by equating the DM energy density and the radiation energy density at the time of matter-radiation equality. The upper limit of coupling are calculated as which does not depend on the details of the reheating history but on the nature of DM, reheating temperature, and equation of states. The corresponding lower limit on the DM mass can similarly be calculated as (see Fig.12) However, such theoretical lower bound will be further constrained by the observational BBN bound on ∆N eff ≤ 0.50 (yellow shaded region). Using this in Eq. (47) one readily infers that the relativistic DM energy density should be less 8% to the over all energy budget at the end of reheating. Upon incorporating such observational constraints we arrive at the following modified expression of lower limit of the DM mass as To this end, we would like to point out the fact that there exists Lyman−α bound on the DM mass m Lyman x > 5 × 10 −6 GeV. However, such bound on DM depends non-trivially on the details of its phase-space distribution and equation of state. Therefore, we defer this discussion in detail for our future studies. From our discussion so far, we have obtained two broad conditions on the DM mass, say m x > m ϕ (A re ) and m x < m ϕ (A re ). When the DM mass satisfies the condition m x < m ϕ (A re ), its abundance decrease with increasing reheating temperature, as shown in Eq.(54) (see also Fig.(12)). As a result, in order to achieve correct abundance for a fixed DM mass, the inflaton-DM coupling must be increased for larger reheating temperature with an exception (see Fig. 12) for ϕϕ → SS production channel with w ϕ < 1/3. The reason is that for such a situation, the ϕϕ → SS channel produces DM only during the initial stage of reheating. On the other hand, when the DM satisfies m x > m ϕ (A re ), the slope of the figure changes (see the bottom plot of Fig.12), and the co-moving DM freezes in at any point during reheating due to kinematics reason where m x ∼ m ϕ . As a consequence, the mass dependency of the abundance Ω x h 2 also changes (see Eq. (55)).

C. Freeze-in and Freeze-out production of DM from radiation bath
In this section, we will discuss DM production exclusively from the thermal bath. In addition, the universal gravitational production of DM will always be present, which can not be ignored. The associated Boltzmann equations (see, for instance, Eqn.7) for the freeze-in production scenario arė n r x + 3Hn r x + ⟨σv⟩ (n r x ) 2 − (n r x,eq ) 2 = 0.
And for the freeze-out scenario, since the cross-section is strong enough, the gravitationally produced DM from inflaton and radiation occurs at the initial reheating stage and reaches thermal equilibrium within a short period. Thus, the Boltzmann equations associated with DM take the following forṁ Where Γ r ϕ = Γ th s/f +Γ gr ϕϕ→RR , γ S = 1.9×10 −4 for scalar DM, γ F = 1.9×10 −3 for fermionic DM [50]. In the above expression, the fourth and last terms are associated with the DM production through gravitational scattering from inflaton, and radiation bath respectively (see Fig.11 for relevant Feynman diagrams for gravitational scattering from the thermal bath (R)). ⟨E x ⟩ r = m 2 x + 9 T 2 rad is the average energy per DM particle produced from the thermal bath [74]. ⟨σv⟩ be the thermally averaged cross section times velocity. n r x,eq be the equilibrium number density of the DM, which can be expressed as where, T rad is the temperature of the radiation bath and j x be the internal degrees of freedom of DM and K 2 (x) is the modified Bessel function of the second kind. The expression of the DM relic abundance in terms of radiation abundance [122,123] where T F be the temperature of the radiation bath at the very late times A F , when both the radiation and DM energy density became freezes, and N r x = n r x A 3 is the co-moving number density of DM. We constrain the DM parameter space (m x , ⟨σv⟩) in terms of (w ϕ , T re ). The population of the DM particles produced from the thermal bath strongly depends on the scattering cross-section ⟨σv⟩. If the scattering cross-section is large enough, the produced DM particles reach thermal equilibrium, and when the bath temperature falls below the DM mass, the number density of DM freezes out-this mechanism is known as the freeze-out mechanism [54][55][56][57][58][59][60][61]. On the other hand, if the scattering cross-section is small enough, the DM can never reach thermal equilibrium, and this mechanism is called the freeze-in mechanism [62,[64][65][66][67][68][69][70]. In this paper, we will discuss both production mechanisms and analyze the parameter space needed to satisfy the correct relic.

Freeze-in from radiation bath: bosonic and fermionic reheating
For the freeze-in from the thermal bath, the DMs will never be in thermal equilibrium, and hence DM number density generically satisfies n r x ≪ n r x,eq . Dynamical equation Eq.61 then transformed into simplified form in terms of co-moving DM number density N r x = n r x A 3 as, The above equation suggests that the production rate is simply proportional to the square of the equilibrium DM number density. In the m x < T rad limit, it behaves as ∝ T 6 rad . Thus freeze-in production from radiation bath naturally follows the way radiation temperature evolves up to the point m x ∼ T rad . The production for the masses m x > T rad will naturally be Boltzmann suppressed. As we have extensively discussed, the evolution of the bath temperature is conditioned non-trivially not only by the production process and its constituents, but also by the bath temperature itself. Therefore, details of the reheating history is expected to have interesting impact on DM evolution and its final abundance. Throughout our analysis, we will provide a detailed analysis of DM phenomenology and its dependence on the reheating history for different physical situations discussed before.
♣ : Freeze-in from the bosonic radiation bath Earlier, we discussed different possible bosonic reheating histories. In this section, we will quote our findings of the DM abundances for those different reheating histories. The detailed calculations are shown in the appendix-C. As discussed for the bosonic reheating case, depending on the range of inflaton-scalar coupling, we have three different cases, Case-I: Coupling strength g r i > G 1, th ci : In this regime, direct decay of inflaton into radiation controls the entire reheating process, and as discussed we have the two broadly classified thermal histories, When T max rad > m end ϕ : For this case, the direct inflaton decay channel controls the reheating dynamics and the thermal effect is effective throughout the reheating period for w ϕ > w c ϕ , so the temperature evolves according to Eq.A23,A24. With this reheating background, we now find the DM abundance for two different mass ranges.
When m x < T re , the present-day DM abundance can be obtained as, for g r 2 ϕ 2 s 2 with w ϕ > 3/11 When the DM mass is lower than the reheating temperature, kinematically, DM production is expected to continue even after reheating until the point when m x ∼ T rad . However, it is important to note that freeze-in production of DM from the radiation bath typically follows the evolution of the bath itself. In most cases, the comoving radiation energy density freezes at the end of reheating. Therefore, for analytical calculation, it is safe to take DM production up to the end of reheating even for m x < T re . On the other hand, there are some situations where radiation production happens initially, which is visible in most of the scenarios where w ϕ < 3/11 for ϕϕ → ss reheating process. For such case, DM production similarly happens instantaneously at the end of inflation, and its number density turned out to be independent of mass but depends on the maximum radiation temperature T rad rad . The resulting expression for the abundance is given in the last expression of Eq.65. Therefore, for this particular reheating, since maximum production happens initially, the above expression of the abundance will remain the same even for m x > T re . However, this should not be true in general.
Generically if one considers the mass m x > T re , the DM is naturally expected to be produced during reheating until m x ∼ T rad , and the abundance for different decay channels are obtained as At this point, we would like to elaborate the exceptional case for w ϕ > w t ϕ , for which the reheating temperature equals the maximum radiation temperature (see, for instance, Eq.A23). If the EoS satisfies w ϕ > w t ϕ and m x > T re , the DM production remains always suppressed and the dominating contribution comes from the initial stage of reheating. On the other hand, for w ϕ < w t ϕ where maximum radiation temperature appears at the initial phase of reheating and for m x > T re DM production occurs till the point m x ∼ T rad .
As discussed earlier, when w ϕ < w c ϕ , we found an intermediate temperature scale T c (at the point A c defined in Eq.27) above which bath temperature dynamics is controlled by the thermally corrected decay width. Hence, for m x > T re , since the freeze-in occurs during the reheating epoch itself, two different possibilities arise.
If m x > T c > T re , the DM will freeze in during the early phase of reheating, where the evolution of bath temperature is controlled by thermally corrected production rate, and the abundance will take the following form, Whereas for T c > m x > T re , the DM will freeze in during the later part of the reheating phase when finite temperature effect is diminished, and the abundance assumes different form as, And if m x < T re , the DM abundance at the present time can be written as When T r,max gr < T max rad < m end ϕ : For this case, direct inflaton decay channel controls the entire reheating dynamics. Similar to the discussion above for m x < T re , the abundance can be written as However, for m x > T re freeze-in will naturally occur during reheating, and for w < w c ϕ the abundance can be found to be the same as Eq.68. On the other hand for w > w c ϕ the abundance will be same as Eq. (66) for m x < T c , and for m x > T c is Case-II: Coupling strength in between G 2, th ci < g r i < G 1, th ci : As discussed earlier, for this coupling range, the gravitational interaction drives the dynamics of reheating at the initial stage. Therefore, the maximum temperature is always controlled by the gravitational sector T r, max gr = T max rad . In this coupling range, a crossover temperature scale T s (at the point A gr→ngr -see Eq.33) exists across which gravitational decay dominated to non-gravitational decay-dominated reheating occurs. When m x < T re , the DM will be produced up to the end of reheating (except for g r 2 ϕ 2 s 2 with w ϕ < 3/11), and hence the maximum production occurs at the end of reheating, and the abundance follows Eqs.65, 69. However, if T s > m x > T re and freeze-in occurs during the decay channel-dominated phase, and the final abundance follows the same form as expressed in Eqs. (66), and (68) with T c being replaced by T s . However, if m x > T s > T re and freeze-in happens during the universal gravitational decay-dominated phase, and we have and Case III: when g r i < G 2 ci : For this case the gravity-mediated decay of inflaton controls the entire dynamics of reheating (termed as gravitational reheating). This particular phase is realized for w ϕ > 0.65 (see the light-cyan region of Fig.2). Similar to the previous case, T r,max gr = T max rad will always holds. Since gravitational radiation production happens only at the beginning of reheating, most of the DM production is also expected to happen at the initial phase of the reheating, and the abundance for such a case is calculated to be, ♣ : Freeze-in from the Fermionic radiation bath Details of the fermionic reheating has been discussed. In this section, we will quote our findings of the DM abundances for those different reheating histories. The detailed calculations are shown in the appendix-C. The distinct behaviour of fermionic reheating, as opposed to bosonic one, mainly arises for the higher value of inflation equation of state, say w ϕ > 7/15 (5/9) for with (without) thermal effect. In addition to that, if w ϕ > 7/15 (5/9), most of the fermionic radiation production occurs during the initial stage of reheating due to its specific behaviour of production rate from the inflaton. Similar to the bosonic one, for fermionic reheating depending on the range of inflaton-Fermion coupling we have two different possibilities, Case-I: Coupling strength, h r > H c : As discussed before, for this coupling regime, non-gravitational decay of inflaton into radiation controls the entire reheating process. In the following subsection consider different temperature regimes When T max rad > m end ϕ : For this case, the thermal correction in the decay rate will be dominant from the beginning. Depending upon the equation of state we have two different possibilities of evolution of the radiation component: (i) 1 > w ϕ ≥ 7/15, the radiation production mainly takes place at the initial stage, (ii) 0 ≤ w ϕ < 7/15, radiation production happens throughout the reheating phase (see, for instance, Eqn.37).
• 1 > w ϕ ≥ 7/15 : As just pointed out, radiation production occurs at the initial stage and hence the comoving radiation density becomes constant early. Following the radiation evolution, the DMs are also produced at the beginning of reheating. As a result the abundance (see, Eqn.75), naturally be controlled by the T max rad , as follows • 0 ≤ w ϕ < 7/15 : For this range of EoS, the ratio T rad /m ϕ varies as A − 3 10 (1−5w ϕ ) which indicates that the thermal effect is dominant throughout the entire reheating phase when 7/15 > w ϕ > 1/5. However, for 0 ≤ w ϕ < 1/5, there exists an intermediate temperature scale T c with the scale factor A c (see Eq.V) above which the thermal effects drops down drastically. Let us discuss two possible scenarios in this context: 1.Dominant finite temperature effect during entire reheating period for EoS 1/5 < w ϕ < 7/15 : We found two different sub-possibilities depending on EoS. a) When EoS is in the range of 9/25 < w ϕ < 7/15, comoving DM freezes at the initial stage of reheating, and due to that, there is an explicit maximum temperature dependence in the DM abundance. Henceforth, for both m x > T re and m x < T re , we have the same DM abundance expression, and that is, b) When EoS lies between 1/5 < w ϕ < 9/25, the expression for DM abundance will be different for m x < T re (comoving DM freezes at the end of reheating) and m x > T re (comoving DM freezes at any intermediate point during reheating where m x ∼ T rad ). For m x < T re , DM abundance is found to be, and for m x > T re we have, 2. Finite temperature effect will not be dominant during the entire reheating period for EoS 0 ≤ w ϕ < 1/5 : As already discussed earlier, for this range of EoS, there exits an intermediate temperature scale T c across which thermal effect drops down ( see discussion around Eq.39). Therefore, for m x > T re (comoving DM freezes in during reheating), we have two different possibilities: i) When m x > T c > T re , the DM freezes in before thermal to non-thermal domination crossover, and we get ii) When T re < m x < T c , the comoving DM freezes after thermal to non-thermal domination crossover, and abundance assumes, when m x < T re , the Comoving DM freezes at the end of reheating, and the abundance can be expressed as T r, max gr < T max rad < m end ϕ : In this case, as described earlier in detail (see, for instance, sec-V) there is no thermal effect at the beginning of reheating. Depending on the reheating background, there are two different possibilities (see, for instance, Eq.42): • 1 > w ϕ ≥ 5/9 : The bath temperature always falls as A −1 since the radiation component is frozen at the beginning of reheating. For the freeze-in mechanism from the thermal bath, it is expected that the DM component follows the radiation and freezes at the beginning, irrespective of its mass. Thus the expression for the abundance will be exactly the same as defined earlier in Eq.75.
• 0 ≤ w ϕ < 5/9 : In this scenario, the thermal effect is subdominant at the beginning, but there is a chance of the thermal effect being dominant at any intermediate scale where T rad ∼ T c . However, depending upon the way the finite temperature effect made its presence on the DM evolution, the EoS range is divided into three sub-ranges, 1.Sub-dominant finite temperature effect during the reheating for EoS 0 ≤ w ϕ < 1/5 : If the thermal correction is not applicable throughout reheating, then the abundance follows the Eq. (81) for m x < T re and Eq.80 for m x > T re .
2. Dominant finite temperature effect at intermediate temperature T c for 1/5 < w ϕ < 3/7 : Here we have two different cases depending on different EoS regimes: i) In the presence of the thermal effect when EoS lies within 9/25 < w ϕ < 3/7, most of the DM production occurs before the intermediate temperature scale T c . Hence, for any value of m x < T c , the DM always freezes-in its production near around the T c , and DM abundance assumes the form, Whereas, for the same EoS range when m x > T c > T re , the comoving DM will freeze in during the initial phase when the thermal effect would be subdominant, and the abundance takes the following form ii) For the range of EoS 1/5 < w ϕ < 9/25, the DM production continues up to the end of reheating, and for m x < T re , the abundance follows the Eq.77. Whereas, when m x > T re , the DM production continue up to m x ∼ T rad . Therefore, for T c > m x > T re , the DM abundance will follow the Eq.78 and for m x > T c > T re , Eq.83 will be the abundance expression.
3. For EoS 3/7 < w ϕ < 5/9: For this case, it is observed that the comoving DM freezes immediately after the reheating begins irrespective of DM mass, and the abundance of the DM takes the following form, Case-II: Coupling strength, h r < H c : As discussed before, for this coupling regime, gravitational decay of inflaton into radiation controls the entire reheating process. Therefore, T max rad = T r, max gr will always be the case. Depending on the different EoS, there are the following possibilities, • 1 > w ϕ > 0.65: Since for this case h r < H c and w ϕ > 0.65, the gravitational sector is the governing reheating dynamics, termed as gravitational reheating. In gravitational reheating following the radiation component, the DM component will also freeze just at the beginning of reheating, irrespective of its mass, and the abundance will be of the same form as expressed in Eq.74.
• 0 ≤ w ϕ < 0.65: For this case, purely gravitational production will not be sufficient to reheat the universe. Hence, to have successful reheating, one needs to have non-gravitational production during the later stage of reheating, and the reheating temperature is defined by non-gravitational fermionic coupling. We found three interesting cases, which are as follows: 1.Reheating temperature T re < T BBN for EoS, 5/9 < w ϕ < 0.65 : It is observed that if the EoS lies in this range, both gravitational and non-gravitational production happen very early in reheating phase. Due to subsequent expansion, the reheating temperature turns out to be always < T BBN (see the red shaded region of Fig.7). However, from the Fig.7 it is clear that for EoS between 0 ≤ w ϕ < 5/9, we have non-trivial dynamics 2.Dominant finite temperature effect at intermediate temperature T s for EoS, 1/5 < w ϕ < 5/9 : In the range of EoS, the finite temperature effect start to dominate at some intermediate temperature scale T s during reheating. Now, if m x < T re , the DMs are expected to freeze after the end of reheating, and consequently, its abundance is calculated to be the same as given in Eq.77. If T s > m x > T re , the DMs freeze in during the later phase of reheating when non-gravitational decay dominates, and the abundance assumes the form of Eq.78. Finally, if m x > T s > T re , the DM will freeze during the initial part of the reheating phase when gravity-mediated decay controls the reheating, and for such case, the abundance has been calculated as 3.For EoS 0 < w ϕ < 1/5 : When m x < T re , the DM abundance is the same expression as Eq.81 and when T s < m x < T re , the abundance is same as defined in Eq.80. Again, if m x > T s , the DM will freeze in during the initial gravitational channel domination sector; in such a case, the DM abundance follows the below equation

Freeze-out from the bosonic and fermionic radiation bath
For freeze in mechanism, the interaction cross-section ⟨σv⟩ is so small that DM can never reach thermal equilibrium. However, if ⟨σv⟩ is large enough, the DM can strongly interact with the SM bath and be in thermal equilibrium. The background expansion eventually helps the DM freeze out from the bath at a certain temperature T f . Conventionally these DMs are called WIMP. The freeze-out temperature T f is defined as, During reheating, inflaton decays into radiation; hence, entropy is not conserved. Due to this physical situation, one can find two distinct situations: Freeze out after reheating: If the mass of the DM m x < T re , it will freeze out after reheating, i.e., during the radiation-dominated era. Moreover, after the freeze out, the co-moving number density N r x = n r x A 3 will be much larger than the comoving equilibrium number density N r x, eq . Thus after freeze-out happens, one can neglect N r x, eq in comparison with N r x and from the Eq.61, one can find 3 dN r 3 Near the DM thermal freezes-out temperature, the non-thermal gravitational production can become dominant as compared to thermal one. This is usually interpreted as the re-annihilation phase [71]. We thank the anonymous referee for pointing out this important fact. In our analysis, we ignored this effect. However, it is to be noted that if freeze-out occurs after reheating, such an effect turns out to be subdominant. The reason for this is that during radiation domination, the inflation energy density is negligible compared to the radiation energy density implying << ⟨σv⟩(n r x,eq ) 2 at T f . On the other hand, for T f > Tre, such an effect can become important, particularly for the higher DM mass range which we discussed in the appendix-F.

Integrating from freeze-out point (A = A f ) to the present time (A =
N 0 is identified as the present day comoving number density leading to Ω x h 2 ∝ 1/⟨σv⟩. Hence, the abundance decreases with increasing ⟨σv⟩. Freeze out during reheating : Alternatively, DM freeze-out could occur during reheating if m x > T re .
, after utilizing this in Eq.61 4 Freeze-out occurs during reheating at some intermediate scale factor A f < A re . Therefore, integrating the above equation for the number density from A f to A re , the comoving number density at the end of reheating is Therefore, the current abundance can be written as, ; We will evaluate the abundance for the aforementioned two cases for different reheating models discussed earlier in detail. Freeze out temperature : The freeze-out temperature T f in general can be computed from Eq.87 by assuming H(T f ) ∝ T k f as, The general solution of the above equation is expressed in terms of Lambert function W −1 (q) of branch −1 with argument q, We further assume the approximate ⟨σv⟩ being independent of temperature in the above expression and throughout our paper. If freeze-out happens after the reheating, one will simply have K = K 0 = √ ϵ √ 3Mp⟨σv⟩jx (2π/m x ) 3/2 , and k = 2. Consequently, the DM parameter space (m x vs⟨σv⟩) turns out to be the same for all reheating temperatures and inflaton equations of states. However, if freeze-out happens during reheating, the expression of (K, k) will differ. It is observed that generically we can express K(T re , T c/s ) = K 0 T µ re T ν c/s , where (µ, ν) will assume different values for different reheating history and that will be our subject of our subsequent discussion.
♣ : Freeze-out from the bosonic radiation bath Case-I: Coupling strength g r i > G 1, th ci : Similar to the freeze-in case, let us discuss two different regimes. When T max rad > m end ϕ : If the inflaton equation of state w ϕ > w c ϕ , the abundance due to freeze-out can be calculated as, for g r 2 ϕ 2 s 2 , and ν = 0, µ = 2 − k, k = 4 To find out the analytical expressions of the DM abundance, here we neglect the re-annihilation of DM. where T f is the freeze-out temperature which we already defined in Eq. (94).
However, if w ϕ < w c ϕ , as has already been discussed for the bosonic reheating, there exists an intermediate temperature scale T c at which the ratio T rad /m ϕ goes less than unity, and the thermal effect becomes subdominant at the later part of the reheating phase. Therefore, if T f > T c , DM freezes out during early reheating phase and the abundance is calculated to be However, if T c > T f , the freeze-out occurs during the later phase of reheating, and the abundance assumes a different form as T r,max gr < T max rad < m end ϕ : For this case, if the inflaton equation of state w ϕ < w c ϕ , the effective mass of inflaton remains to be greater than the radiation temperature, and the finite temperature effect will be subdominant throughout. The DM abundance for such case will be same as Eq.97, and the evolution of the comoving number density is depicted in Fig.(13). Whereas for w ϕ > w c ϕ , a new temperature scale T c emerges as before. If DM mass happens to satisfy the condition m x < T c , DM freezes out with radiation temperature T < T c , and the abundance assumes the same form as expressed in Eq. (95). On the other hand if m x > T c , freeze-out occurs during the initial non-thermal phase, and DM abundance becomes, Case-II: Coupling strength in between G 2, th ci < g r i < G 1, th ci :This reheating history is described in the bosonic reheating section. The gravity-mediated decay channel controls the initial reheating dynamic up to T = T s , and then the non-gravitational decay channel controls the reheating dynamics. As a result, initially, DM production is driven by the gravity-mediated decay channel (up to T = T s ) and then driven by the non-gravitational decay channel (see middle Fig.13). For w ϕ > w c ϕ , if the DM freezes out during the late nongravitational decay channel domination phase, i.e., T f < T s , then the DM abundance follows the Eq.95, and if the DM is frozen out during gravity mediated reheating phase T f > T s , the DM has following abundance, Tre Ts Tre Ts Again for w ϕ < w c ϕ , if the DM is frozen out during the late decay channel domination phase, the abundance has the Eq. (97), and if the DM is frozen out during the initial gravity-mediated reheating phase, the DM abundance has the following equation Tre Ts 4 for, g r 1 ϕs 2 , ν = 4(1+w ϕ ) Tre Ts Case-III : For this reheating scenario, the gravitational sector is the governing reheating dynamics termed as gravitational reheating (GR). The evolution of co-moving number density is shown in the left plot in Fig. (13). The DM abundance is ♣ : Freeze-out from the Fermionic radiation bath • w ϕ > 5/9 : For w ϕ > 5/9, the bath temperature always behaves as T rad = T re (A re /A), and using this equation into Eq.(92), the abundance is, • 0 < w ϕ < 5/9: For this range of equation of state, we discuss three different possibilities as follows: Case-I: Coupling strength h r > H c : T max rad > m end ϕ : Depending upon the evolution of radiation, we will have different behavior of the DM abundance in terms of reheating temperature. For 7/15 < w ϕ < 5/9 the radiation temperature behaves A −1 , and abundance follows the Eq.102. On the other hand if 1/5 < w ϕ < 7/15, radiation temperature behaves T rad = T re (A/A re ) − 3(1+5w ϕ ) 10 , and using this in Eq.92, we get Finally, if the 0 ≤ w ϕ < 1/5 similar to scalar reheating case, the intermediate temperature scale T c , leads to two different possibilities. If m x > T c , DM freezes out during the early reheating phase, and we have, where On the other hand, if m x < T c , DM freezes out during the late reheating phase when the finite temperature effect is subdominant, and we have T r,max gr < T max rad < m end ϕ : For this case, finite temperature correction is subdominant at the beginning. In fact, for lower EoS, namely 0 < w ϕ < 1/5, such finite temperature effect will be subdominant throughout the reheating, and hence the abundance turned out to be same as Eq. (105). On the other hand, for w ϕ > 1/5, the finite temperature effect will become dominant after an intermediate temperature scale T c , and if the DM freeze-out temperature satisfies T f < T c , the abundance obeys the Eq.103. Conversely, if the DM freezes out during the early non-thermal phase with T f > T c , one will have Coupling strength h r < H c : For this case, initially, the reheating phase is governed by the gravitymediated decay channel up to a point A gr→ngr , and the later part of the phase is dominated by direct inflaton decay. If the freeze-out happens during the later phase, the abundance will be the same as Eq.103 for w ϕ > 1/5 and Eq.104 for w ϕ < 1/5. If freeze-out occurs at the early phase of reheating, the abundance assumes, Tre Ts Tre Ts

VII. DM PARAMETER SPACE, (⟨σv⟩, mx) FOR FIMPS AND WIMPS FROM RADIATION BATH
In this section, we will discuss in detail the parameter region where DM production mechanisms are at play for different reheating histories. For the sake of completeness and understanding the DM parameter space, we consider three different EoS (w ϕ = 0.0, 0.20, 0.50) with five different reheating temperature (T BBN , 10, 10 6 , 10 10 , 10 15 ) GeV which includes both minimum and maximum reheating temperature. An important point we want to infer is that in the final contribution to DM abundance, we ignore any possible contribution from the non-gravitational inflaton-DM coupling but include the universal gravitational production from both inflaton and radiation scattering. And as we pointed out before, such universal production has been observed to set a maximum limit on the DM mass m g,max x (≤ m end ϕ of course) specifically for the freezein production. Interestingly for a given w ϕ , we found a universal feature of lowest possible DM mass within m min x ≃ 150 − 300 eV (see Eq. (D3)) irrespective of its nature and the reheating histories for which freeze-in and freeze-out mechanism coincides during the radiation-dominated era. Below this mass scale, all the DM has under abundance today. However, the critical cross-section ⟨σv⟩ crit (for analytical expressions, see Appendix-D) at which this coincidence occurs depends on the reheating temperature, which is represented by filled black circles in Figs.14-16. Moreover, another critical cross-section exists for higher DM mass (m max x ) where the freeze-in and freeze-out mechanism coincides during reheating. Unlike the lower DM mass bound, m max x has been found to have non-trivial dependence on the reheating temperature but also depends on the EOS w ϕ , reheating background. All these features are clearly observed in Figs.14, 15, 16. In some reheating temperatures, there is no meet point of freeze-in and freeze-out; it is due to the gravitational DM, which gives the overabundance for the freeze-in mechanism where we expect the meet point occurs. When m max x < m g,max x condition is satisfied, then the abundance condition Ω x h 2 = 0.12 gives rise to a closed contour in (m x , ⟨σv⟩) plane due to aforementioned two critical cross-sections for two different masses (m min x , m max x ), where the smooth transition happens from freeze-in to freeze-out or vice versa. For example, see Fig. (14) where we have found close contours for T re = (10 −2 , 10, 10 6 ) GeV and in Fig. (15) for T re = (10 −2 , 10) GeV with bosonic reheating (for fermionic reheating we only have closed contour for T re = 10 GeV). On the other hand, if m end ϕ > m max x > m g,max x , the critical cross-section disappears for higher mass values; for such cases, the maximum allowed DM mass is turned out to be m g,max x (m max x ) for freeze-in (freeze-out) mechanism. Again, when m max x > m end ϕ (it happens for higher reheating temperatures), one can find the freeze-in and freeze-out meet point, but the close contour is not formed. To this end, we would also like to point out that DM annihilation is dominant but under-abundant in the region outside the closed contours, whereas shaded regions are under-abundant for open contours. The minimum and the maximum DM masses up to which DM can give present abundance are 10 −7 and 7 × 10 16 GeV, respectively. But, further, some parameter space is ruled out by the ∆N eff at BBN (yellow shaded region) and by unitarity bound of ⟨σv⟩ ≤ 8π/m 2 x [74,124,125]. When WIMP DM has a mass scale in the order of T BBN or below, it violates the ∆N ef f bound at BBN. Again for WIMP DM, freeze-out happens during RD (m x < T re ) if the mass scale lies above 10 5 GeV, it violates the unitarity bound of ⟨σv⟩, which is shown by the light purple color region. Using two bounds, for the freeze-in mechanism, the maximum and the minimum allowed scalar (fermionic) DM masses are 10 −7 (10 −7 ) GeV, 4 × 10 16 (10 15 ) GeV respectively. FIG. 16: Here we have plotted the ⟨σv⟩ as a function of dark matter mass mx for w ϕ = 0.50 for five different reheating temperature Tre = BBN, 10, 10 6 , 10 10 , 10 15 GeV. The plot is on the left side for scalar dark matter and the right side for fermionic dark matter. The solid lines for the freeze-in mechanism and the dotted lines for the freeze-out mechanism. The yellow-shaded region is ruled out by the ∆N eff bound at BBN, and purple shaded region is ruled by the unitarity bound.
When m x < T re , we obtained some generic behavior of the DM abundance specifically for freeze-in production. Moreover, in the freeze-in production mechanism, if DM freezes during the radiation-dominated era, ⟨σv⟩ behaves as ∝ 1/m x (see Eqs. 65,69,70). On the other hand, when m x < T re , for the freeze-out scenario, even though we do not have such simple relation, it turned out that ⟨σv⟩ generically increases slowly with increasing m x .

VIII. WIMPS, EXPERIMENTAL BOUNDS AND CONSTRAINTS ON REHEATING
In this section, we would like to discuss our results from the perspective of some indirect experimental constraints with their future projected sensitivity limit on the DM cross-section and its mass plane. BBN constraints on the effective number of relativistic degrees of freedom ∆N ef f put direct constraints on the possible lower limit on the DM mass. In all the ⟨σv⟩ Vs m x plots, the yellow shaded regions depict the forbidden region coming from BBN bound. For WIMP-like DM, the condition of the ∆N ef f at BBN sets the approximate lowest possible DM mass m x ∼ T BBN . In addition, for 2 → 2 scattering processes, the unitarity constraints further provide a bound on the maximum possible mass m x ∼ 10 5 GeV for thermally produced DMs which freeze out, particularly during the RD era (T re > m x ). However, if freeze-out occurs during reheating era (T re < m x ), DM mass can be large, which is observed from the pink dotted curve of the third and fourth plot of Fig.17. Therefore, reheating dynamics plays a significant role in setting the possible values of maximum DM mass, and it depends non-trivially on the inflationary parameter such as inflaton equation of state. In Fig.(17), we have shown where our WIMP DM parameter space lies with the available sensitivity curve for indirect experiments within the mass range (0.01 − 10 5 ) GeV such as (i) PLANCK CMB measurement [86] labeled as CMB shown in the light yellow region in the mass range of 0.1 to 10 GeV. (ii) Alpha magnetic spectrometer (AMS)-02 experiment [87,88] provided cosmic ray positron (e + ) and antiproton (p) data with unprecedented precision in the mass range 5 → 10 3 GeV. (iii) Combined bounds (labeled as Combined) from Fermi-LAT, HAWC, HESS, MAGIC, and VERITAS experiment [89] provided a upper limit on the DM annihilation cross-section in the case of two annihilation channelsbb (purple dashed) and τ + τ − (red dashed). (iv) CTA experiment [90] is a ground-based gamma-ray detector that could be capable of searching the WIMP DM at the TeV scale with large sensitivity. The brown (orange) dashed line shows its expected sensitivity for WIMP annihilation tobb (W + W − ). In the near future, CTA might be able to probe a significant portion of TeV scale WIMP DM that freeze-out during either RD or reheating era. In order to illustrate our result, we choose w ϕ = (0.0, 0.20, 0.50) with the background of both bosonic (upper plot) and fermionic reheating (lower plot). The solid and dot-dashed lines are for T re = T BBN , 10 GeV, respectively. The black solid lines correspond to those DMs which freeze out during the RD era and for which one requires ⟨σv⟩ ∼ O(10 −9 ) GeV −2 to get correct present-day DM relic. Moreover, most of the parts of this black line are ruled out by the experimental bound (except CTA since it is a future experiment) as it lies inside the regime of the sensitivity curves. For bosonic reheating, if w ϕ = 0.5, for both T re = T BBN and T re = 10 GeV, the WIMPs follow the black line (see, for instance, first and the second plot of Fig. 17) to achieve the correct relic. For lower EOS, w ϕ = 0.0, 0.20 with low reheating temperature, WIMP-like DM freezes out during reheating, and one requires the lower value of ⟨σv⟩ to obtain present-day DM relic. Therefore, the DM parameter space for freeze-out during reheating is still very much allowed. Thus in order to detect those DMs in the near future, we need to increase the sensitivity of the experiments.

IX. CONCLUSIONS AND OUTLOOK
In this paper, we performed a detailed analysis of the physics of reheating after the end of inflation. The physics after the end of inflation is expected to play a significant role in every aspect of the late time universe. Apart from the well-established correspondence between the physics of CMB and the early inflationary era, the intriguing effects of reheating on our present universe can not be avoided. Unlike inflation, broadly reheating is similar to the usual phases of the standard big-bang model, except for the fact that it occurs right after the end of inflation. Therefore, experimentally it is challenging to look for its direct signatures. Further, it is generically argued that decoding any physics information is challenging because of non-linear thermalization processes during this phase. Over the years, however, endeavor toward understanding this phase has gained significant interest due to its rich new physics contents. Reheating is the phase that naturally encodes the physics of inflaton itself. The way inflaton decays into different fundamental fields is expected to be imprinted into different cosmological observables such as CMB anisotropy, DM, and gravitational waves in terms of new physics. Therefore, reheating could be an interesting playground for phenomenological studies of DM and inflation in a unified framework that has not been studied extensively in the literature, and this is the main objective of our present paper. We have addressed two main topics: Inflaton phenomenology: In the first part, we have studied in detail the dynamics of reheating separately through inflaton decaying into scalar field (bosonic reheating) and decaying into fermions (fermionic reheating). Apart from the direct decay term, we further include two important effects, namely, the universal gravitymediated decay of inflaton, and the finite temperature correction on different decay channels. The effects of gravitational decay have already been analyzed before [96]. It has been shown that purely gravitational reheating sets a lower bound on the reheating temperature when the inflaton equation of state satisfies w ϕ > 0.65. Such a lower limit on the reheating temperature further sets constraints on the inflaton direct coupling parameter below which gravity-mediated decay will be the dominant channel for reheating process (see the cyan shaded region in Fig.2, 7). The cyan zone of these plots is the important outcome of our present analysis. There is an equation of state-dependent critical values of the coupling constant below which reheating dynamics become completely independent of direct inflaton coupling. Combining this gravitational decay with the finite temperature correction in the direct decay, we observe intriguing interplay among those different physical effects during reheating. The thermal effect is expected to influence the dynamics when radiation temperature satisfies the condition T rad > m ϕ (t). Depending on the strength of the inflaton-radiation coupling, we observed several interesting reheating histories which have not been reported earlier. Let us outline the main interesting outcomes in the following discussion ♠ We particularly observed a new reheating history for which reheating temperature (T re ) and maximum radiation temperature (T max rad ) coincide when the equation of state w ϕ > (1/3, 3/5) for two different bosonic decay channels ϕ → ss and ϕϕ → ss (see, Fig.4 for depiction) respectively. ♠ Another interesting observation in the case of fermionic reheating is that when EoS lies above 5/9(7/15), maximum radiation production takes place at the initial stage of reheating. In this respect, therefore, fermionic reheating turns out to be qualitatively similar to gravitational reheating. The bosonic reheating through ϕϕ → ss decay process deserves special mention with regard to the fact that, if w ϕ < 3/13(1/5) with (without) thermal effect, successful reheating can not be achieved (see the grey shaded region of Fig.2). ♠ The phenomenological constraints on the inflaton coupling are shown in Fig.2 for bosonic reheating and in Fig.7 for fermionic reheating. The plots depict an interesting connection between the inflaton coupling parameters (g r i ,h r ) and inflationary parameter n (power of the inflaton potential at its minimum), which is translated into effective EoS w ϕ through the relation (n − 2)/(n + 2) during reheating. Contrary to the earlier claim [101], we observe that the thermal effect appeared to be most significant at low reheating temperatures for non-zero inflaton equation state, w ϕ = (0.20, 0.50, 0.82, 0.99) (see solid and dotted lines in Fig.6, 10), for which inflaton mass undergoes non-trivial evolution. ♠ Thermal correction in the production rate significantly modifies the production rate (enhances for bosonic channels and diminishes for fermionic channel), which is imprinted in the radiation temperature evolution. Better visualization of how the thermal correction affects the evolution of the radiation temperature with respect to scale factor is shown in Tables-I, III. ♠ For higher equation state (w = 0.82, 0.99), the gravitational reheating has been observed to give a lower limit on the reheating temperature (10 3 , 10 6 ) GeV, which is associated with the fixed scalar spectral index n s = 0.9685 and 0.9681, respectively. When n s reaches these values, the coupling parameter tends towards zero, i.e., gravitational scattering solely controls the reheating dynamics.
DM phenomenology: In the second half of the paper, we have extensively analyzed DM phenomenology in the context of production during reheating. We have discussed all known DM production mechanisms, namely: (i) gravitational production from both the inflaton and radiation scattering, (ii) production from direct inflaton decay through various decay channels, (iii) production from the thermal bath through effective average crosssection times velocity ⟨σv⟩ in both freeze-in and freeze-out mechanism with assuming ⟨σv⟩ as a free parameter. The universal gravitational production of DM is always incorporated throughout our analysis. In the following, we summarize some of the main results and important outcomes of our analysis: ♠ When DMs are produced from direct decay of inflaton (see Fig.12), the final abundance appears to depend only on the coupling strength, DM mass (m x ), and reheating temperature (T re ). We call this FIMP-like dark matter. For a fixed inflation EoS and reheating temperature, the minimum DM mass is fixed by the ∆N eff bound at BBN, and the maximum allowed DM mass is fixed by either gravitational decay or the inflaton mass. ♠ For the production of DM from the thermal bath, we have considered both freeze-in (FIMP) and freeze-out (WIMP) mechanisms. Since the DM is produced from the thermal bath, the DM production strictly depends on the evolution of the bath temperature. Due to the thermal effect, the bath temperature evolves differently during reheating, and the thermal effect directly influences DM production and its final abundance. ♠ In the context of both bosonic and fermionic reheating, we have discussed the DM production (both FIMP and WIMP) and derived the analytical expressions of the DM abundance. In Section-VII, we have shown in detail the DM parameter space (⟨σv⟩ Vs. m x ) for both FIMP and WIMP production mechanisms. Interestingly, the lowest possible DM mass for both the mechanisms assumes a universal theoretical value m x ∼ 10 −7 GeV, which is independent of T re , w ϕ and the nature of DM. However, constraints from different physical considerations such as lyman-α, ∆N ef f bounds may not satisfy this universal bound. For WIMPs, the minimum allowed value of DM mass is T BBN due to restriction from the ∆N ef f . Depending on the reheating parameters such as reheating temperature, inflation EOS w ϕ , and the background reheating dynamics, a maximum allowed DM mass also exists. In most of the scenarios, the maximum allowed mass is given by a particular mass scale where there is a transition happening from freeze-in to Freeze-out mechanism or vice versa (see, for instance, black circle of Figs. 14-16). ♠ For the DM production from the radiation bath, if one incorporates the unitarity bound, the DM mass lies above 10 5 GeV is ruled out in most of the cases. However, for both fermionic and bosonic reheating, we still have some parameter space even for m x > 10 5 GeV (see, for instance, the third and fourth plot of Fig. 17), which are consistent with the unitarity bound and satisfies the correct relic. The reason behind this is that if freeze-out happens during reheating (m x > T re ), the freeze-out temperature shifts towards the maximum radiation temperature with the increasing mass, and we need to lower the cross-section to meet the correct relic. Such lowering of the cross-section with higher DM mass (m x > 10 5 GeV), therefore, naturally evades the aforementioned unitarity bound. ♠ Important distinction of fermionic decay width from that of the bosonic one is its proportional behavior in terms of inflaton mass, which dilutes faster with increasing w ϕ . Due to this, one finds a critical w ϕ ≃ 7/15 considering the thermal effect from the beginning of reheating, below which the inflaton decay process occurs throughout and above which it occurs at the beginning of the reheating. These two different reheating histories have had their distinct effects on the DM production, which is reflected on ⟨σv⟩ V s m x parameter plane (see, for instance, third and forth plot of Figs.15, and fifth and sixth plots of 16). This phenomenon can also be explained from the Eqs.102-105. For w ϕ = 0.5 > 7/15, a fixed value of T re < m x and T f ∼ m f , the relic abundance behaves as Ω x h 2 ∝ m 1/4 x /⟨σv⟩ (see Eqn.102), which implies that with increasing m x , ⟨σv⟩ increases (see, for instance, fifth and sixth plots of Fig.16). However, for w ϕ = 0 < 7/15, for a fixed T re if freeze-out happens during reheating one has following relation m x ∝ ⟨σv⟩ −3 (see, for instance, Eqn.105). This clearly suggests that with increasing DM mass above reheating energy scale, ⟨σv⟩ decreases (see, for instance, Fig.14). ♠ The oscillating inflaton condensate may be fragmented due to self-resonance (see, for instance, Refs. [128,129]), which is expected to cause reheating dynamics non-trivial. It is observed that such an effect is not significant for n = 1 (w ϕ = 0). However, for higher n, self-resonance effects turned out to be important and may lead to radiation dominated (RD) universe depending on both α and n values. When α << 1/6, selfresonance is efficient, and the RD universe is established within less than an e-fold of expansion after inflation end. Whereas for α > 1/6 and n > 2, the self-resonance effect increasingly becomes inefficient, and the inflaton condensate remains intact for a long time without any other decay channel. It is in this parlance our present study is important. In our study, we considered α = 1 and further checked that the final results, namely the estimates of reheating temperature for different inflaton couplings and the whole DM parameters space, are not much sensitive with the increase of α within the limit set by recent Planck and BICEP/keck data [130,131]. For higher α values, even though the self-resonance effect would be inefficient (as pointed out in [128,129]), it may have a non-trivial effect on the reheating parameter space. Another important way forward of our present work would be to construct particle physics-motivated UV complete models that can be directly probed through laboratory and cosmological experiments. One may further study the gravitational wave dynamics in this context and explore the potential signature of reheating through the present-day gravitational wave spectrum. Finally, the analysis of cosmological perturbation will be interesting to work on, which will be a good probe of early universe physics. We defer all these topics for our future work.
Mp . The field-dependent mass becomes, Using the envelope approximation at any instant of time, the envelope value of the field ϕ 0 must satisfy V (ϕ 0 ) ≃ ρ ϕ (t). Under this approximation, the inflaton mass at the inflation end is given by After utilizing Eqns.(A5), (A6) and (A1), one can find the evolution of the inflaton mass as Inserting the above equation in Eqn.(A3), we can re-write decay width in terms of the scale factor as Upon substitution of the above decay width (Eqn.A8) in Eqn.A2, we have the evolution of the energy density associated with the bosonic non-gravitational channel as where, H end = ρ end ϕ /3M 2 p . Since during reheating, the background dynamics is mainly governed by the inflaton field, we approximate H ≃ ρ ϕ /3M 2 p = H end (a/a end ) −(3+3w ϕ )/2 , which is used to determine the above equation. After straightforward integration of Eq.(A9), we obtain, We will follow the similar procedure to find the expression of gravitational radiation energy density, and the production width associated with the gravitational sector can be expressed as Combining Eqns. (7) and (A11), we have the following solution for the radiation energy density associated with the gravitational sector ρ r gr (A) = The total energy density of the radiation bath is simply sum of the contribution from both gravitational and non-gravitational sector ρ r tot (A) = ρ r s (A) + ρ r gr (A) and the corresponding radiation temperature T rad = 30ρ r tot /(π 2 g ⋆ ) 1/4 . Where g ⋆ is the total relativistic degrees of freedom associated with the thermal bath. Now let us find the maximum radiation energy density for both gravitational and non-gravitational sectors separately.
At the beginning of the reheating, the individual components of radiation (both non-gravitational and gravitational) set to be ρ r s (A = 1) = ρ r gr (A = 1) = 0. Within a very short time, during the initial stage of reheating, both the components attain a maximum value when d ρ r s/gr d A = 0, which is associated with a specific value of the scale factor That, in turn, fixes the maximum value of the radiation component as ρ r, max Since the gravitational scattering process is a kind of irreducible background, which will be present regardless of the coupling-dependent decay channel. However, a critical coupling exists G 1 ci , above which the non-gravitational decay channel always controls the reheating dynamics. The maximum gravitational production happens initially, and the temperature associated with it falls off faster than the non-gravitational one. Thus, the critical coupling G 1 ci , above which value the non-gravitational sector dominates throughout reheating, is defined from the condition ρ r, max s = ρ r, max gr . Since in the gravitational sector domination, the maximum radiation temperature comes out to be 10 11 → 10 12 GeV, which is less than the inflaton mass m ϕ (A max ); thermal effect not in working situation. Therefore the expression for critical coupling for with and without thermal effect be the same G 1, th ci = G 1 ci (see, for instance, Eqn.22 ). In the limit of g ci > G 1 ci , the maximum radiation temperature is followed by Eqn.A14, whereas, for g ci < G 1 ci gravitational sector determines the maximum radiation temperature (see, Eqn.A15). If the bosonic coupling strength is less than G 1 ci , the gravitational sector can not reheat the universe before BBN energy scale for the equation of state w ϕ < 0.65. Therefore those (g r i , w ϕ ) parameter space are forbidden from BBN constraints (see, for instance, light red region of Fig.2). However, for w ϕ > 0.65, the gravitational sector can successfully reheat the universe and sets the lower bound of reheating temperature that is defined in Eq.(A17). There exists another critical coupling parameter G 2 ci , indicating below which only the gravitational sector defines the reheating temperature (see, for instance, light cyan region of Fig.2). We have the expression of G 2 ci as follows for g r 1 ϕs 2 and w ϕ ≥ 0.65 , for g r 2 ϕ 2 s 2 and w ϕ ≥ 0.65 .

(A16)
and when the gravitational sector controls the reheating temperature, the reheating temperature can be expressed as where A gr re is the normalized scale factor at the end of gravitational reheating. When the coupling parameter g r i > G 2 ci , the decay channel controls the reheating temperature and that can be defined as where for g r 2 ϕ 2 s 2 . (A19)

With thermal effect
If the radiation temperature is much smaller than the inflaton mass scale, we can safely use the zerotemperature decay width. But, for the cases where radiation temperature is much higher than the inflaton mass scale, the thermal effect becomes very efficient. In the limit T rad >> m ϕ (t), the inflaton decay rate for different decay channels can be expressed as , . (A20) When the bosonic non-gravitational channel dominates over the gravitational one, we can approximate radiation temperature as T rad = These aforementioned equations represent the modified expression of the bosonic radiation energy density when the thermal effect is functional (T rad >> m ϕ (t)). The maximum value of this radiation component associated with the scale factor A max at which the dρ r s /dA = 0. We have the maximum radiation energy density and the scale factor A max associated with the maximum radiation energy density can be expressed as for g r 1 ϕs 2 , 2 3−5w ϕ 2 3(5w ϕ −1) for g r 2 ϕ 2 s 2 . (A26) Now let us define reheating temperature in this context. For the coupling g r 1 ϕs 2 , we have Whereas for the coupling g r 2 ϕ 2 s 2 T re = 9M 4 p (1 + w ϕ )H 3 end (g r 2 ) 8πϵ(5w ϕ − 1)(m end ϕ ) 4 .
(A28) and the radiation temperature T rad = (ρ r tot /ϵ) 1/4 . After the end of the inflation, the energy budget associated with the thermal bath begins to increase from zero initial value and attain a maximum value when dρ r f /gr /dA = 0, then starts to fall off. We obtain the maximum radiation energy density corresponding to both gravitational and non-gravitational sector as for h r ϕf f , for gravitational scattering . (B5) And the scale factor correlated with these maximum energy densities turns out as Since the gravitational scattering process is always present, there exists a critical value of coupling below which the gravitational sector always defines maximum radiation temperature T max rad . Equating ρ r, max f = ρ r, max gr , one can find the expression for the critical coupling for fermionic reheating as (B7) Therefore, in the limit of h r > H c , T max rad can be approximated as T max rad ≃ ρ r, max f /ϵ 1/4 , whereas for h r < H c , T max rad ≃ ρ r, max gr /ϵ 1/4 . Another interesting fact is that the behavior of the radiation component associated with the non-gravitational sector is different for w ϕ < 5/9 and w ϕ > 5/9 (see, for instance, Eq.37). For w ϕ > 5/9, the radiation component decays as A −4 as the gravitational sector. Thus, if your Yukawa coupling strength h r < H c and w ϕ > 0.65, then reheating dynamics is governed by the gravitational sector, which is termed as gravitational reheating. Moreover, the reheating temperature is followed by Eq.A17. Otherwise, in all (h r , w ϕ ) parameter space reheating end is defined via explicit coupling. In such cases, T re and A re is defined as : when w ϕ < 5/9 T re = 6M 2 p (1 + w ϕ )m end ϕ H end 8πϵ(5 − 9w ϕ ) (h r ) 2

With thermal effect
Here we are mainly focused on the situation where T rad > m ϕ and the analysis for T rad < m ϕ will be the same as without thermal effect. In the limit, T rad > m ϕ , the thermal effect significantly impacts the reheating dynamics due to the explicit temperature dependence on the decay rate. The expression for decay rate in this limit can be written as (see, for instance, Eq.2) Inserting the above decay rate in Eq.6, one can find the evolution of the radiation energy density for fermionic coupling as ζ(w ϕ ) = 15(1+w ϕ ) 64π(7−15w ϕ ) . The aforementioned equation clearly suggests that the radiation energy density behaves differently depending on the inflaton equation state, whether it is greater than or less than 7/15. If the coupling strength h r > H c , fermionic coupling always defines reheating temperature and eventually equating ρ ϕ = ρ r f , we find reheating temperature as The above equation is valid when w ϕ < 7/15. However, for w ϕ > 7/15 we have . (B13) Moreover, in the limit of h r > H c , the maximum radiation energy density is also defined by the non-gravitational sector, which can be expressed as where, A max = 10 3+15w ϕ 2/7−15w ϕ . On the other hand, when h r < H c and w ϕ > 0.65, the gravitational sector drives the reheating dynamics and the reheating temperature to be the same as Eq.A17.
Appendix C: Analytical expressions of co-moving number density for FIMP which are produced from thermal bath For freeze in production of DM from the thermal bath, the DM can never reach thermal equilibrium i.e. n r x ≪ n r x,eq , hence the co-moving number density N r x = n r x A 3 follow the following simple equation, For this case, the DM production continue to happen untill the radiation temperature equals the DM mass T rad ≃ m x . Therefore, in order to solve the above equation, DM can be safely assumed to be relativistic, and in the limit m x < T rad , the equilibrium number density becomes, n r x,eq = j x T 3 rad π 2 (C2) The Hubble parameter can also be written as where H(A re ) = √ ϵT 2 re √ 3Mp hubble parameter at the end of reheating. So, combining above three equations, one can gate .
• For w ϕ < 5/9: The bath temperature can be written as Upon substituting the above equation in Eq.C4 along with Eq. (C3), the co-moving number density at any point A during reheating takes the following form When w ϕ > 3/7, the DM production happens instantaneously just at the end of inflation, as a result, the co-moving number density is constant.

With thermal effect :
When w ϕ > 7/15, the same situation will be happen as we discuss for without thermal effect for w ϕ > 5/9. Now for w ϕ < 7/15, the bath temperature is Upon substituting the above equation in Eq.C4 along with Eq. C3, the co-moving number density at any point A during reheating takes the following form When w ϕ > 9/25, the DM production happens instantaneously just at the end of inflation. As a result, the co-moving number density is constant.
Appendix D: Analytical expressions of minimum critical mass where both freeze-in and freeze-out mechanism coincides When the DM mass m x is smaller than T re , the DM abundance follows simple relation Ω x h 2 ∝ m x ⟨σv⟩ for the freeze-in scenario for a fixed T re . Therefore, the cross-section becomes inversely proportional to the DM mass. This suggests the existence of a critical ⟨σv⟩ and the associated mass for which the DM equilibrates with the thermal bath and freeze out happens. At that critical ⟨σv⟩ crit , DM number density must equate with the equilibrium number density at the end of reheating. Corresponding to these critical ⟨σv⟩ crit , there exists a minimum critical mass m x,min which satisfies present-day abundance. Below this mass, no mass will be available, which can give the correct abundance. Equating the solution of the DM number density n r x with the equilibrium number density n r x,eq at the end of reheating A re , one can find following expressions of ⟨σv⟩ crit without thermal effect for h r ϕf f with 3/7 < w ϕ < 5/9 3−7w ϕ Tre for h r ϕf f with w ϕ < 3/7 (D1) and with thermal effect for h r ϕf f with 9/25 < w ϕ < 7/15 2(9−25w ϕ ) 5Tre for h r ϕf f with w ϕ < 9/25 (D2) Above this critical ⟨σv⟩ crit , the DM can produce only from the freeze-out mechanism. Using the above critical ⟨σv⟩ crit in the DM abundance expressions, one can find the following expression of minimum critical mass The minimum mass is independent of background reheating dynamics, i.e., reheating temperature and inflation equation of states.  Appendix F: The re-annihilation of DM and its effects on DM parameter space Near the freeze-out temperature, the thermal production of DM receives Boltzmann suppression, and during this time the gravitational scattering may become the dominant channel for the DM production. For such case following conditions must satisfy Γ ϕϕ→SS/F F (t)ρ ϕ (t) m ϕ (t) ≥ ⟨σv⟩(n r x,eq (t)) 2 .
In our analysis, we ignored such an effect, and it is also important to note that such an effect is not important when the DMs freeze out during the RD era. However, when the DM freeze-out takes place during reheating, the DM re-annihilation turns out to be important within a very narrow DM mass range (10 11 − m end ϕ ) GeV, and its evolution till the Freeze-out point governed by, d(a 3 n r x ) dt = −a 3 ⟨σv⟩(n r x ) 2 + a 3 Γ ϕϕ→xx ρ ϕ (t) m ϕ (t) . (F2)