Two phase reheating: CMB constraints on inflaton and dark matter phenomenology

We propose a two-phase reheating scenario where the initial preheating dynamics is described by an effective dynamics followed by the standard perturbative reheating. Some of the important universal results of lattice simulation during preheating have been considered as crucial inputs in our two-phase dynamics. In this framework, detailed phenomenological constraints have been obtained on the inflaton couplings with reheating fields, and dark matter parameters in terms of CMB constrained inflationary scalar spectral index. It is observed that the conventional reheating scenario generically predicts the maximum reheating temperature $T_{re}^{max} \simeq 10^{15}$ GeV, corresponding to an almost instantaneous transition from the end of inflation to radiation domination. This fact will naturally lead to the problem of non-perturbative inflaton decay, which is in direct conflict with the perturbative reheating itself. Taking into account this by incorporating effective non-perturbative dynamics as the initial phase, our model of two-phase reheating scenarios also predicts model-independent maximum reheating temperature, which does not correspond to the instantaneous process. Further, $T_{re}^{max}$ is predicted to lie within $(10^{13}, 10^{10})$ GeV if CMB constraints on inflaton couplings with different reheating field are taken into account. We have further studied in detail the dark matter phenomenology in a model-independent manner and show how dark matter parameter space can be constrained through CMB parameters via the inflaton spectral index. Considering dark matter production during reheating via the Freeze-in mechanism, its parameter space has been observed to be highly constrained by our two-phase reheating than the constraints predicted by the conventional reheating scenarios, which are believed to theoretically incomplete.


I. INTRODUCTION
The inflationary universe [1] is currently the leading paradigm to explain the inhomogeneities in the Cosmic Microwave Background (CMB) [2], which plays the crucial role of seed perturbations for the large scale structure of the universe [3]. Within the present setup, the inflationary phase must be followed by a phase known as reheating 1 when the energy stored in the inflaton field is released to defrost the universe. Unlike inflation, the reheating phase is not constrained by direct observables; however, the modified expansion history of the universe due to the presence of the reheating phase prior to the hot big bang evolution influences the relation between physical scales of the CMB mode today and that at the time of their Hubble exit during inflation (See, Fig.2 in this context). This was the basic idea of reheating constraints to inflationary models from CMB [6]. Despite, the thermalization process erasing many micro-physical details of this phase, a better understanding of this phase is necessary and can shed light on how inflationary mechanism is connected to the rest of the universe dynamics [7][8][9], the production mechanism of baryonic asymmetry in the early universe called baryogenesis [10][11][12][13][14][15][16], the origin of dark matter [17], and the generation of primordial gravitational waves [18][19][20][21][22][23][24][25][26][27] and constraining the etc. There exist two approaches that can constrain the reheating phase through the inflationary models. We can either model the expansion during the reheating phase using an effective equation of state parameter [6,28], or we can solve the Boltzmann equation system supplemented with the background expansion [29]. Both the description have their limitations and are not theoretically complete.
However, the latter approach's advantage is that one can further generalize it by including matter components in addition to radiation, which could be physically motivated. For example, in the original work, the production of dark matter has been studied. The study revealed a very interesting link between the dark matter with CMB through inflation [29]. Later this formalism has been extended in considering various models of inflaton with general power-law types potential [30] and non-perturbative effect from numerical lattice simulation has also been considered [31]. In this paper, we take up this issue of the non-perturbative phase known as preheating and formulate an effective approach that will be shown to lead qualitatively different results than that of the usual reheating constraint analysis. Preheating is the phase when the occupation number of field quanta for both the inflaton and daughter field(s) grows exponentially due to parametric resonance [32].
The equation of state during steady state of the preheating phase is crucial for model independent reheating constraint analysis. Further, some scenarios can lead to non-standard case such as sudden blocking reheating process due to Higgs field [33] or breakdown of coherent oscillation without thermalization [34]. Considering those into our present scenario would interesting to consider. Although, the preheating phenomenon depends on the inflation model and its interaction with the daughter fields, certain universal behaviors have been observed to emerge irrespective of inflation models [31,35,36]. Namely, i) The preheating phase is episodic with at-least three distinct phases.
These results indicate that while non-perturbative processes dominate the initial stage, the inflaton decay in later stages should be described by perturbative channels. A systematic study of reheating constraints incorporating the non-perturbative phase for various models and interactions is missing in the literature.
In this work, we will extend the formalism developed in [29,31] to include the non-perturbative effects in reheating constrains analysis. As we just mentioned, reheating happens in multiple stages, with the initial non-perturbative stage followed by the perturbative one. Let us now briefly describe the main idea of the present work: endowed with the above few universal features, we will consider reheating as a two-stage process. We will model the initial non-perturbative stage (henceforth, phase-I) governed by an effective fluid with an effective equation of state (w ef f ). One of the boundary conditions of the effective non-perturbative dynamics will be set by the inflation model potential expressed in terms of the scalar spectral index (in this regard, the reader may find this in parallel with the conventional works in [6]). This phase is assumed to be continued until the inflaton energy decaying into 50% of its initial energy and the subsequent perturbative stage (henceforth, phase-II) follows. While evolving through phase-I and connecting to phase-II, we allow the system to satisfy an important consistency relation associated with the total energy conservation comprising inflaton and various daughter fields such as radiation, dark matter (which we considered separately). We will see that this consistency relation will restrict the possible values of w ef f during phase-I as opposed to the conventional analysis [6]. Furthermore, the perturbative decay in phase-II will be restricted by the CMB constrains [29].
We have structured our paper as follows: In section, II, we discuss the general analysis of singlestage perturbative reheating, and in section III, we will try to specify the possible limits on perturbative reheating considering some specified form of the interaction between inflation and radiation field. Finally, In section IV, we briefly describe our proposed two-phase reheating analysis and, in section V, illustrate the strategy of our numerical study. After that, in section VI, we will try to find out an analytical estimation of the maximum radiation temperature and reheating temperature.
Next, in section VII, we consider different inflationary models of and analyze in the context of the two-phase scenario and compare it with conventional reheating dynamics. In section VIII, we analyze the possible constraints on the coupling parameter correspond to different inflaton-radiation field interactions. Furthermore, in section IX, we include additional dark matter components and discuss the viable restrictions on the dark matter parameter space.

II. REHEATING CONSTRAINT ANALYSIS FOR PERTURBATIVELY DECAYING IN-FLATON
Before we directly jump into constructing the two phase reheating model, let us first elaborate on widely studied single phase perturbatie reheating with decaying inflaton following [29]. This not only explains the methodology of our analysis, but also helps us to identify the regime of its validity which will further motivate the reader, the need for considering physically more acceptable two phase reheating process mentioned in the introduction. While discussing this we will see one of the important results that is the existence of maximum reheating temperature. Subsequently the generalization to two phase reheating will show how the aforesaid maximum reheating temperature reduces depending upon the initial condition. Let us start with the following Einstein's equation for the cosmological scale factor and conservation of energy, ρ φ + 3ṅ re (ρ φ + p φ ) +ρ rad + 4ṅ re ρ rad = 0, with the following Freedman-Roberson-Walker (FRW) form of the metric ds 2 = −dt 2 + a 2 (dx 2 + dy 2 + dz 2 ) Where, "ρ"s are the energy densities of two different components. At any instant of time during reheating, we parametrize the duration of reheating by e-folding number n re (t) = ln(a/a i ), where "a" is the cosmological scale factor. The time derivative of n re is the Hubble expansion parameteṙ n re = H during reheating. During reheating we assume the effective equation of state of the inflaton w = p φ /ρ φ to be approximately constant. The fundamental difference between our present analysis followed from [29] and that in [6] is the consideration of Eq.1, where we consider the multiple dynamical components. Considering the evolution of ρ φ + ρ rad = ρ ef f together, the effective equation of state during reheating can be defined as Hence, w ef f will essentially interpolates between two values (w, 1/3) through non-trivial time dynamics for decaying inflaton ρ φ and the growing radiation field ρ rad . However, in [6], the authors have taken it to be constant during their analysis. Therefore, we not only employ realistic decay dynamics into the reheating constraint analysis but also provides a new framework to go beyond which is our main purpose of the present paper.
Keeping the above points in mind, let us express the total energy density as, which is followed from the conservation eq.1. The index "i" stands for the initial stage of reheating, which also marks the end of inflation. At the beginning of reheating we set ρ rad (t i ) = 0. For solving the above set of equations, the boundary condition is set by the inflaton energy density aṡ n re (t i ) = H(t i ) = ρ i φ /3M 2 p . The physical quantity of our interest is the ratio of the radiation energy density and the inflaton energy density. From eq.4, one gets Where, "f" corresponds to the final value of radiation density. We define total e-folding number during reheating as N re = n re (t f ).
The main goal of this whole program of reheating constraint analysis is to understand the relation among early universe inflaton dynamics, the intermediate reheating dynamics and late time dynamics. A particular cosmological scale k going out of the horizon during inflation will re-enter the horizon during late time cosmological evolution. This fact will provide an important relation among different phases just mentioned as follows The usual approach is to define the effective equation of state of the total energy density during reheating and study its evolution. However, we consider only the radiation part during reheating and try to understand the evolution of its temperature T rad as a function of scalar spectral index, and finally connect the temperature with CMB one on large scale [29]. The reheating temperature T re is identified with radiation temperature T rad at thermal equilibrium between the decaying inflaton and the radiation. From the entropy conservation of thermal radiation, the relation among T rad = T re at equilibrium, and (T 0 , T ν0 = (4/11) 1/3 T 0 ), temperature of the CMB photon and neutrino background at the present day respectively, can be written as Using eq. (6,7), one arrives at the following well known relation Where, g re ∼ 100 is the effective number of relativistic degrees of freedom during radiation phase.
In our subsequent study we identify the cosmological scale k as the pivot scale set by PLANCK, k/a 0 = 0.05M pc −1 and compare our result with the corresponding estimated scalar spectral index n s = 0.9682 ± 0.0062 [44].
A. Example I: Exactly solvable case As has already been discussed in [29], one of the important outcome of our formalism is the existence of maximum possible reheating temperature. In this and next section, we will elaborate on this considering simple ansatz of decaying inflaton. We first consider an analytically solvable case where the inflaton is decaying asρ Γ φ is a dimensionless constant, which parametrizes the decay of inflaton. This form of decay essentially modifies the Hubble friction term for the dynamics of inflaton during reheating. With aforementioned ansatz for the decaying inflaton, radiation density is analytically solved as The second term in the parenthesis is quantifying the fractional amount of inflaton energy left after the reheating process is over. Expressing ρ f in term of radiation temperature as ρ f rad = π 2 (g re /30)T 4 rad , the Eq.10 leads to the following maximum radiation temperature [47] for a given Γ φ as where, (x = 4/(3 + 3w +Γ φ ), P =Γ φ /(Γ φ + 3w − 1). This also can be clearly seen from the Fig.1 for each value ofΓ φ . From the perturbative point of view, the value of Γ φ should be ≤ 1. However, if we naively extrapolate the above result for largeΓ φ most important result turned out to be the existence of a maximum possible temperature, 2.9 × 10 15 GeV.
However, the numerical value of this maximum temperature turns out to be of the order of same as the limiting perturbative value forΓ φ = 1 as shown in the figure Fig.1. Therefore, above temperature can be naturally identified as maximum possible reheating temperature. This also corresponds to the maximum possible value of scalar spectra index n max s . Identifying associated temperature of the produced radiation in eq.10 with eq.8, we arrive at the following exact expression for (N re , T re ), In the fig.1, we have considered three possible values ofΓ φ for quadratic inflaton potential. The special value isΓ φ = 1, for which the equilibrium condition between the inflaton and the radiation can be achieved at the maximum temperature shown as a black dot. The maximum value of scalar spectral index turned out to be n max s 0.9654. This analysis motivates us to subsequently analyze more general case, and we will show that this conclusion still holds.

B. Example II: Standard pertrubative case
In this section we will consider the standard perturbatively decaying inflaton parameterizing by decay constant Γ φ as follows,ρ Where Γ φ is effective time independent inflaton decay constant. It is the phenomenological term which acts as a damping force during the oscillating inflaton. This term can be related to the total decay rate of inflaton to radiation. However, we believe our conclusion will remain same for time dependent Γ φ , which we will study later. Before doing any numerical analysis, let us examine the approximate solution which has already been discussed in the literature [47]. During the early stage of evolution, approximating ρ i φ e −Γ φ t ρ i φ , the radiation density can be calculated as Similar to the exactly solvable case in eq.11, the above equation also leads to a maximum radiation temperature [47] for a given Γ φ , Where, the relation T re = 0.45 (200/g re ) 1/4 Γ φ M p has been used. In the same way as our earlier exactly solvable case, maximum possible reheating temperature could be obtained, if one identifies a special point where two temperature meets, T max rad = T re . Our numerical analysis also shows the maximum reheating temperature at the aforementioned special point, Interestingly, the maximum reheating temperature T max re can also be computed for another exactly solvable case with w = 1/3. Corresponding result is as follows, . This expression is exactly the same as previously discussed. For this special value of w = 1/3, we also have exact expression for all the reheating parameters (T re , N re ) as follows, At this point let us again emphasize the fact that as long as we are in the perturbative regime, the relation among the scalar spectral index n s and the reheating temperature T re can be understood from our detail analysis above. However, existence of maximum reheating temperature will come if we extrapolate all our formulas for large Γ φ > 2ρ i /(3M 2 p ). For low scale inflation, Γ φ could always be in the perturbative regime. For large scale inflation, this could lead to non-perturbative regime, which will be discussed in our subsequent section. We will discuss about possible limits on the value of Γ φ below which our analysis will be valid. To this end, it is important to point out that in the effective reheating equation of state description [6], the maximum temperature can be explained in the limit of zero reheating e-folding number N re . Therefore, large Γ φ limit in our analysis can be thought of as equivalent to the zero N re limit of the previously studied reheating constraint analysis. However, it is important to remember that those two facts are certainly not identical.
Corresponding to our maximum temperature, we have a minimum reheating efolding number. Our prediction of maximum reheating temperature ∼ 10 15 GeV, and its model independence could be robust and they are intimately connected with the observed CMB scale.
Never the less the main point of our study is to understand the effect of decaying inflaton into the reheating constraint analysis. We think this is the appropriate procedure to understand the relation among (T re , n s ). Another advantage of our procedure is that we can easily generalize our analysis to include any other decay products during reheating such as dark matter which is observed to be dominant matter component of our universe [81]- [84], and that can shape the observed pattern in the CMB. Before, this, our main motivation would be to incorporate the non-perturbative aspects of reheating into our formalism.

III. REGIME OF VALIDITY OF PERTURBATIVE REHEATING
In this section we will try to mention the possible limits on the inflaton decay constant assuming some specific form of the interactions among the inflaton and the reheating field. As emphasized throughout the present work, we have assumed that the inflation decay to other components( for the present work the radiation component) is effectively described by a phenomenological decay term Γ φ . In fact, this was the first attempt to reheat the universe [53]. However, it was soon realized that once the particle production initiates, the inflaton decay is subject to various non-perturbative resonance production and feedback mechanisms. Those processes can change the reheating scenario dramatically. Though it has been argued in [7] that all such feedback mechanisms will have no effect on the CMB. Depending upon the coupling the parametric resonance can be very efficient which may complete the reheating era within a few efolding and in such cases the CMB will have a very little to tell about the reheating phase. Despite that the situation may not be such helpless as noted in [55] that the interactions amongst the produced particles can delay the parametric resonance extending the efolding number of reheating. This will eventually improve the situation of CMB constrain on reheating phase. It must also be noted that we can always choose the coupling constant small enough to evade the parametric resonance. Below, we will briefly mention the space of parameter region in which the perturbative treatment of reheating will be valid over the parametric resonance.
A. Inflaton decaying into scalar particle

Scalar φχ 2 interaction
First let us consider the case when inflaton decays into another scalar particle φ → χχ with the following interaction term L = −gφχ 2 . In this case the vacuum decay width for the decay process φ → χχ is given by [54] where, m φ , m χ are the mass of the inflaton and produced particle respectively, and g is the coupling constant. The mode function χ k of the decay product can be cast into the following Mathieu Where, z = (m φ t − 2z − π/2), A k = 4k 2 /m 2 φ , q = 4gΦ/m 2 φ . Φ is the initial amplitude of the inflaton during oscillations. The Mathieu equation is known to show resonance solutions of the form χ k ∝ exp(µ k z). The condition for the resonance to be efficient is formulated as This can be transformed into the following condition 2 on the dimensionless coupling constant g = g/m φ indicating the regime of perturbative validity [42] g ≤ V which can further expressed in terms of decay constant as Therefore, we see that if the decay width satisfies aforementioned condition, the perturbative reheating will be the only mechanism and our perturbative analysis will be at work. Given a model Γ cri φ (model) is the point which qualitatively separates the perturbative and non-perturbative effect of inflaton decay.

Scalar φχ 3 interaction
In this case, inflaton couples to another light scalar via interaction where y is the coupling constant. The vacuum decay rate of the inflaton field into three bodies φ → χχχ can be determined by Dalitz plot [40] as, At the tree level, the mode function χ k following the same Mathieu equation and the consideration from the previous scalar φχ 2 interaction can be correlated if one replacesgm φ Φ → h 2 Φ 2 . Therefore, the condition to treat the dynamics of reheating perturbatively is roughly To estimate the lower bound on the coupling for the resonance, we make a substitution Φ → φ end .
The above condition for the effectiveness of perturbative reheating can be written in terms of decay rate as, Thereafter, in our analysis, we want to examine whether this above condition consistent with our analysis or not for the different inflationary models.

B. Inflaton decaying into a pair of fermions
Let us now consider the case when the inflaton decays into a pair of massless fermions with the following Yukawa interaction where h is the dimensionless coupling constant. Now the vacuum decay rate is given by The condition for the validity of perturbative reheating in this case, as shown in [56], can be written as, Hence, in connection with decay rate, the equation (33) is rewritten as, In our proposed effective two phase dynamical scenario we will observe the existence of similar critical inflaton decay constant associated with the reheating e-folding number. We will see how aforementioned three different interacting model dependent critical decay constants restricts in initial parameter space of the reheating dynamics. In the following sections our attempt will be to build up a formalism which can effectively incorporate the non-perturbative dynamics at the initial stage of the reheating.
IV. REGIME OF EFFECTIVE NON PERTURBATIVE AND PERTURBATIVE REHEAT-

ING
The standard and well-studied mechanism to consider the non-perturbative effect during reheating is called preheating. This stage is essentially the combination of a highly non-linear process of parametric resonance and subsequent thermalization. The well-known fact that generically nonperturbative preheating mechanism does not completely decay inflaton into the radiation field.
Therefore, subsequent perturbative decay will be necessary to complete the reheating process. To the best of our knowledge, the reference [29] has considered this issue for the first time and studied perturbative reheating, followed by the preheating considering a specific model of chaotic type inflation. However, generically the preheating mechanism is model dependent. Hence, combining the end of preheating and subsequent model-independent perturbative reheating is somewhat irreconcilable. Therefore, our objective in the following sections would be to make these two phases reconcilable.
Instead of dwelling into explicit non-perturbative computation during the preheating stage, we will adopt an effective model-independent approach following the reference [6]. The basic idea is to assume the dynamics of preheating to be solely governed by an effective equation of state ω ef f supplemented with the total energy conservation law in terms of its constituents. As already emphasized in the introduction, the information about the actual non-perturbative dynamic will be encoded though considering its universal features into our effective dynamics.
As has been pointed out already, during the non-perturbative dynamics, inflaton decay is not complete, and typically it is around 50% of its total comoving energy, which is being transferred into the daughter fields. Furthermore, inflation models with quadratic potential near the minimum, the non-perturbative reheating does not lead to the equation of state, ω = 1 3 , which is expected at the end point of reheating [31,63]. Our essential idea would be to correctly utilize those results as the end point conditions of our proposed effective dynamics in place of preheating. After the end of this, the usual Boltzmann perturbative reheating process will follow. The second phase completes the reheating process by leading to the correct state equation with relativistic degrees of freedom as the dominant components collectively called radiation. The Fig.2 illustrates our methodology of calculation. Throughout this paper, we call this as two-phase reheating process.
phase-I:(Effective non-perturbative phase) During the early stage of reheating, the phase will be described by total energy density ρ T = ρ R + ρ φ and the constant effective equation state w ef f .
Hence the evolution will be described by, where ρ T e is the total energy density at the end of the inflation. a end is the scale factor at the end of the inflation. In this section we will build up our formalism considering two matter component with ρ φ , ρ R as the inflaton energy density and radiation energy density respectively at any instant of time. In the subsequent section we will add dark matter as a third component as an extension.
Nonetheless, from the total energy expression one can write down the following equation follows from Eq.35,ρ To reduce the number of unknown parameters, to this end we will also utilize total energy conservation relation considering individual equation of state of the inflaton (ω φ ) and the radiation field Given the aforementioned constraint relation one obtains the possible restricted value of the effective equation of state, ω ef f . To find those restrictions we combine above two eqs. (36) and (37), and obtain the following consistency relation, Right at the end of inflation or the beginning of the preheating phase, the energy density of the radiation part will naturally be closed to zero ρ R 0. As time evolves, ρ R increases due to decaying inflaton. This initial condition automatically restricts the possible values of w ef f to be very closed to that of the inflaton equation of state w φ . Of course more appropriate approach would be to assume the inflaton equation of state evolving from the value very closed to ω φ to the value required in the next phase. We will comment on this issue at the appropriate place during our discussion.
However, from the usual numerical lattice simulation, the preheating phase's significant duration is dominated by the inflaton. Hence, the equation of state will naturally be closed to that of the inflaton field. Thus, our effective dynamics approach towards preheating truly captures all these necessary properties of the non-perturbative dynamics. As w ef f turned out to be no longer a free parameter and constrained by the above consistency relation (38), our following analysis will be based on this important result. Throughout our study we consider models for which the inflaton equation of state during phase-I ω φ = 0, and consequently following Eq.38 we choose two values of effective equation of state ω ef f = (10 −6 , 10 −3 ). This choice will automatically fixes the initial radiation densities during phase-I as ρ R /(ρ φ +ρ R ) = (3×10 −6 , 3×10 −3 ). We will see the maximum reheating temperature crucially depends upon these initial conditions.
Phase-II: (Perturbative phase) Once the preheating dynamics ends the usual Boltzmann perturbative reheating follows. During this period various components of the total energy density satisfy the following standard Boltzmann equations [62], where the inflaton field φ decays into radiation with the decay rate Γ φ . ω 1 φ represents inflaton equation of state during perturbative reheating. Important to note that inflaton equation of state during phase-I, ω φ is taken to be different than that of the phase-II, ω φ . This is where we will again consider lattice simulation results as another important input. Now that we have identified the full reheating phase in terms of two distinct stages, we will numerically solve all those equations self consistently. With the appropriate dimensionless rescaled variables for the inflaton and radiation energy densities, governing equations for the effective dynamical preheating phase turn into the following form, and the associated governing equations for the subsequent perturbative phase will reduce into Phase-II The rescaled scale factor is defined as A = a a end . "Prime"( ) represents derivative with respect to A. The constant C 1 and redefined variables are, m φ is the mass of the inflaton. In next section we will describe the methodology for solving the above set of equations numerically.

V. PROCEDURE FOR NUMERICAL ANALYSIS AND BOUNDARY CONDITIONS
Let us describe the strategy of our numerical study. We first identify the inflation modeldependent input parameters as N k , H k , V end for a particular CMB scale k. For a given a canonical inflaton potential V (φ), the inflationary e-folding number N k and Hubble constant H k can be expressed as where, the field values at a particular scale k, (φ end , φ k ) are computed form the condition of end of inflation, and equating a particular value of scalar spectral index with n s (φ k ). Therefore, we will get explicit relations between (N k , n k s ) and (H k , n k s ). The well known inflationary input parameters can be found out from the following equations which are expressed in terms of slow-roll parameters The reheating parameters N re , T rad will implicitly depend upon the scalar spectral index n k s for a given scale. The above expression can be inverted to find φ k in terms of the scalar spectral index. After identifying all required parameters from inflation, we will set the initial conditions for subsequent reheating dynamics. Using all these relations among those parameters, one can establish the connection between CMB anisotropy and reheating through inflation.
Phase-I initial condition: The initial conditions for phase-I of the reheating dynamics (effective non-perturbative era) are set by the end of inflation at A = 1 and the equation 38. Those are as follows, Where V end (φ) which is defined at the end of inflation, fixed by φ end . The initial Hubble expansion rate is expressed as H I = ρ end φ /3M 2 p . Subsequent perturbative dynamics will now crucially depend on the end point of the first phase of reheating namely the phase-I. On this issue we rely on the actual non-perturbative lattice simulation results [31,35,63] considering specific model of reheating where the inflaton field is assumed to couple with the reheating field. This system has been studied quite extensively [66,67] in the literature by using the publicly available numerical package LATTICEEASY [68] and its parallelized version CLUSTEREASY [69]. The non-perturbative analysis for different inflationary models has been proved to yield some universal results which will be our important input for the numerical analysis. Extensive works on non-perturbative reheating analysis yields an important fact that only the 50% of the total comoving inflaton energy density is getting transferred into the daughter field. Additionally the inflaton equation of state tends to achieve a steady state value depending upon the power law form of the inflaton potential near its minimum. For example if one assumes the inflaton potential to be of power law form V ∼ φ n , for chaotic type model namely n = 2, non-perturbative phase ends with steady value of the equation of state ∼ 0.2.
However, for other value of n ≥ 4, the equation of state approaches ω = 1 3 at the end point of the non-perturbative reheating. These are the crucial quantitative results from non-perturbative preheating dynamics we will be utilizing in our analysis for the phase-II dynamcis.
Phase-II initial condition: After the phase-I dynamics, the pertubative dynamics will automatically follow. However, important point would be to identify the appropriate boundary conditions. The starting moment of phase-II will be set by the normalized scale factor A npre = a npre /a end which is the ratio between the scale factor at the end of the effective non-perturbative epoch namely phase-I a npre , and the end of the inflation. The initial conditions for the dimensionless comoving densities are It is important to realize that the initial condition is determined by 50% decay of the total comoving energy density ρ T . Further, we have numerically checked that our results do not seem to depend both qualitatively as well as quantitatively much on the amount of decay within 40% − 60% of the total energy at the end of phase-I. For the analysis, we further assume the inflaton equation of state ω 1 φ 0.2 irrespective of the models under consideration. This approximate value is again another important input from the Lattice simulation. For comparison, we also consider the cases where either phase-I or phase-II evolution completely governs the reheating dynamics.
Determining the reheating parameters: Once we numerically solve the reheating dynamics, we define one of the important parameters called reheating temperature T re , which is generically identified as the radiation temperature T rad when the condition H(t) = Γ φ is satisfied, where A re is the normalized scale factor at the end of the reheating. Accordingly, the reheating temperature in terms of radiation temperature (T rad ) is expressed as, Furthermore, the e-folding number during reheating N re consists of two contributions born out of two distinct phases as where N pre and N npre are the e-folding number during perturbative and effective nonperturbative region respectively. Combining equations (8) and (53), we obtain the most important modification of Eq.8 relating the reheating and inflationary parameters, Now connecting equations (51), (52) and (55), we can establish one to one correspondence between T re and Γ φ .
As described before, we will consider three possible cases and compare the results Kamionkawski et al 2014 [6] Case-III N npre = 0 ; N pre = 0 Phase-II, Discussed in the previous section (56) To this end let us specifically mention about the case-II, when perturbative dynamics ceases to exist. This particular procedure proposed in [6], has been studied quite extensively in the literature [32]. In this particular phase, dynamics is solely governed by the effective equation of state ω ef f .
The explicit decay of inflaton does not appear in the computation. However, information about the decay constant Γ φ is extracted from the equilibrium condition Γ φ = H, where the reheating of relativistic degrees of freedom. However, not to ignore an important difference between the phase-I described before and the approach devised in [6] or case-II for the present study is the additional conservation equation 38. This essentially differentiates the regime of applicability of these two approaches. Phase-I dynamics is assumed to be applicable in the early non-perturbative regime. Whereas, since condition 38 does not exist, the original Kamionkowski et al [6] approach is effectively applicable throughout the full period of reheating without any microscopic details.
Further, the value of ω ef f is no longer constrained to be very closed to the inflaton equation of state during phase-I. This relation essentially helps us to compare the results for various scenarios we consider. To avoid symbol confusion whenever we study case-II, we use the symbol ω K ef f instead ω ef f which we reserve for two-phase reheating dynamics.

VI. MAXIMUM RADIATION TEMPERATURE AND REHEATING TEMPERATURE: ANALYTIC STUDY
Before moving on to a particular model, let us analytically estimate the maximum reheating temperature and its dependence upon the initial condition following the same line as before. Considering the standard definition of the radiation temperature T rad = 30 π 2 g * ρ R

1/4
, and computing the radiation energy density during phase II following the Eqs. (39), (40), (36) and (37) the approximate radiation temperature assumes the following form (see appendix A for details calculation) where x, β, c and H in express as In the above expression ρ in φ , ρ in R represent inflaton and radiation energy density respectively at the end of phase-I or the beginning of phase-II The maximum radiation temperature defined at the point x max = a max /a npre , where dT rad dx = 0, which gives us the maximum radiation temperature for two phase reheating expressed in terms of dimensionless comoving densities, One particularly notices the correction term in the maximum radiation temperature due to initial comoving radiation density R(A npre ) at the beginning of phase-II. It boils down to well know expression T max re = D 1/4 in the R(A npre ) = 0 limit same as Eq.17. In this above expression, we ignored the contribution of dark matter. However, generically during the reheating period, dark matter is not the dominant component, therefore, the numerical value of the reheating temperature will not be affected. The analytic expression for the dimensionless comoving density during phase II related to the density at the end of inflation will be Where A npre is the normalized scale factor at the end of the effective dynamics (phase I) One of our important results from the above expression for the maximum reheating temperature is the highest radiation temperature, which is defined at T max rad = T max re corresponding to a given n max s .
As we change the value of ω ef f = (10 −3 → 10 −6 ), the maximum reheating temperature changes as T max re = (10 13 → 10 10 ) GeV. Once we set R(A npre ) = 0, the maximum reheating temperature becomes T max re ∼ 10 15 GeV as expected (see Eqs. 17,18). Proceeding further, we can also obtain the approximation expression for the reheating temperature itself. Utilizing the expression of Hubble constant at the equilibrium point (H re = Γ φ ), and subsequent entropy conservation, one arrives at the following expression Where, x re = a re /a npre can be recognize as Here The detailed derivation of all the aforementioned equations for the reheating temperature in the appendix B. Now we will consider a class of inflationary models of inflation and analyze our proposal of two-phase reheating scenario.

VII. INFLATION MODELS AND NUMERICAL RESUTLS
Based on our methodology discussed above, we will now consider a class of inflationary models for which the inflaton potentials assume quadratic form. We will also point out the regime of validity of the effective non-perturbative and perturbative era for the different inflationary models.
After the inflation, the inflaton field generically oscillates around the minimum of its potential V (φ). Reheating fields coupled with the oscillating inflaton is generically prone to non-perturbative particle production. Our objective is to replace this non-perturbative dynamics by an effective dynamical equation, which is solely governed by the effective equation of state, ω ef f supplemented with the additional constraint relation eq.38. We have already observed that during phase-I, ω ef f is closed to that of the inflaton equation of state, ω φ . Near the minimum of the potential if the form is taken to be power law as ∝ φ n , over multiple oscillations, the average inflaton equation of state is expressed as [101] For n = 2 model, ω φ assumes dust like equation of state (ω φ = 0). Throughout the subsequent study, we consider those inflationary models which have quadratic potential near their minimum.
Therefore, during phase-I of reheating, we set ω φ = 0. To this end, let us emphasize again that during phase-II, when the reheating dynamics enter into the perturbative phase, we assume corresponding to dotted green and dotted black curves. For all cases, ω φ = 0. One of the most important outcomes of our analysis is the emergence of a critical inflaton decay constant Γ φ = Γ cri φ denoted by red dots associated with each particular ω ef f . This indicates the fact that for Γ φ > Γ cri φ , the reheating period will be dominated by phase-I, effective non-perturbative dynamics, otherwise it is perturbative dominated. The critical value of inflaton decay constant increases with the decreasing ω ef f . This can be understood from several interconnecting physical effects. First of all most important Eq.38, which not only fixes the approximate value of ω ef f but also sets the initial condition for phase-I dynamics. Further, larger the value of ω ef f , higher will be initial radiation density R(A = 1) which automatically leads to smaller value of phase-I e-folding number N npre . Therefore, a particular Γ φ will naturally lead to larger N pre , as associated with each Γ φ there exits a reheating temperature which follows from T re ∝ e −Nre = e −(Npre+Nnpre) . On the other hand critical Γ cri φ is a point where N pre = N npre in N vs Γ φ space. From these two conditions one can argue that transition from perturbative to non-perturbative reheating phase would occur for larger critical value Γ cri φ for larger ω ef f value. Given a reheating model with specific inflaton-daughter field interaction, we have also discussed about the existence of critical inflation decay constant Γ cri φ (model) which were shown by vertical red lines in the plots. From the theoretical values of the critical inflaton decay constant, (Γ φ vs N k ) plots indicates that the value of ω ef f must lie within (10 −3 , 10 −6 ) irrespective of the inflationary models considered. At this point let us understand the physical meaning of nonvanishing initial radiation density ρ R (A = 1) (10 −3 , 10 −6 )ρ in φ considering ω ef f = (10 −3 , 10 −6 ). We replace the full non-perturbative dynamics by an effective dynamics, which naturally does not capture the complete picture. Typically non-perturbative phase contains three distinct phases: parametric resonance phase, thermalization phase and steady state phase. And this is the initial parametric resonance phase, where explosive particle production can naturally give raise to required initial radiation density ρ R (A = 1) 10 −6 ∼ 10 −3 in unit of total density ρ φ almost instantly.
To see whether our proposed phase I dynamics is justified or not, we compare our result with actual non-perturbative results. In order to do that, we use non-perturbative lattice simulation during the preheating, considering a specific inflaton-reheating field interaction 1 2 g 2 φ 2 χ 2 . In all the lattice simulation results, the initial radiation density typically assumes ρ R (A = 1) 10 −4 in units of initial inflaton energy density, which essentially lies within what we have considered. Once the preheating phase reaches the steady-state condition, we again solve perturbative dynamics, and found that the reheating ends at around the same value of A re (shown by the dashed black line) where our two-phase reheating ends for ω ef f = 10 −3 (solid pink line) and for ω ef f = 10 −6 (dashed pink line) accordingly. Therefore, our effective two-phase reheating approach seems to capture the essential properties of non-perturbative lattice results, except the non-perturbative e-folding number, which will be taken up in the future.
Nevertheless, for comparison, in the same plot, we also have drawn total reheating e-folding number for other two cases: dotted pink lines for case-II and solid blue lines for case-III mentioned before. It turns out that total number of reheating e-folding number for case-II, case-III, and the case-I, N re = (N npre + n pre ) are almost the same for all different values of the equation of state.
In an another class of plots in (n s vs T re ) space, we describe the variation of reheating tempera-ture T re with respect to the scalar spectral index n s . From these plots, we can read that two-phase reheating process (case-I) is crucially dependent upon the value of ω ef f . Furthermore, case-I results are qualitatively similar to that of the case-II ω K ef f = 0.212 (equation of state at the starting point of phase-II in two-phase analysis). On the other hand, perturbative reheating (case-III) results are qualitatively similar to that of the case-II for ω K ef f = 0. For usual perturbative reheating scenario (case-III) the semi-analytic approach discussed before reveals the existence of maximum possible reheating temperature ∼ 10 15 GeV. Our numerical computation also indicates the same thorough for ω ef f = 10 −3 , and N npre ∼ 12 for ω ef f = 10 −6 . Smaller the effective equation of state during phase-I, larger will be its duration N npre and consequently T max re will be reduced. As expected from our earlier analytical calculation the important results are the values of maximum reheating temperature T max re ∼ (10 13 , 10 10 ) GeV for ω ef f = (10 −3 , 10 −6 ) respectively. Physical origin of this two different limiting temperature is clear from the fact that increase of T re is directly connected with the increase of Γ φ . Hence with the increasing temperature reheating dynamics undergoes a transition from perturbative to non-perturbative regime at particular critical temperature T cri re associated with Γ cri φ , leading to a distinct value of N npre which is different for different ω ef f value. This leads to different T max re . Therefore, an important conclusion we can arrive at is that given the approximate estimates of model specific critical decay width Γ cri φ (model), the maximum reheating temperature T max re should be within (10 10 − 10 13 ) GeV, irrespective of the dynamics of the second phase-II and inflationary model under consideration. However, we must note that the associated maximum values of the n max s are model dependent, which will be discussed for each model. .
A. Chaotic inflation [48] Even though chaotic inflation is observationally disfavored, we consider this potential for its simple nature. For usual chaotic inflation the potential looks like, Where n = 2, 4, 6 . . . . If we consider only the absolute value of the field, n = 3, 5, . . . can also be included. m is parameter of mass dimension. For the purpose of our study, we only consider n = 2 mainly because ω φ = 0.
Initial conditions for phase-I: The initial densities to solve dynamical equation during phase-I can be calculated as where A δφ ∼ 10 −9 is the amplitude of the inflaton fluctuation which is measured from CMB observation. m φ is defined as second derivative of the inflaton potential. To establish the connection among inflationary and reheating parameters the inflationary e-folding number, N k and tensor to scalar ratio, r k are similarly calculated as, Initial conditions for phase-II: Additionally the initial conditions for the phase-II will be set at the normalized scale factor A npre where phase-I ends. The conditions are To establish the relation between reheating temperature (T re ) and inflationary index (n k s ), we follow the methodology explained in the previous section.
Observations: Important results for chaotic inflation are depicted in Fig.(4). As stated at length, This entails the fact that If Γ φ > Γ cri φ (T re > T cri re ), the reheating phase will be dominated by non-perturbative process.
For concreteness, let us bring specific reheating models into consideration. We have discussed three different interaction models with associated non-perturbative constraints equations (30), (26) and (34). Associated with those we have theoretical values of the critical inflaton decay constants Γ cri φ (model) = (0.003, 0.5, 11.8) GeV respectively. The first two values correspond to inflaton decaying into the scalar particle, and the third one corresponds to decaying into a pair of fermionic particles. Interestingly, comparing those numerical and theoretical values of Γ cri φ , one can observe that the initial effective equation of state ω ef f during phase-I must lie within (10 −3 , 10 −6 ). This essentially suggests that all the three models of inflaton interaction will lead to initial radiation density within the value (10 −3 , 10 −6 ) instantaneously, which we can immediately read off from the The variation of the reheating temperature as a function of the spectral index for the different reheating mechanism is shown in the fig.4. The behavior of reheating temperature with respect to n s appears to be model-independent.
B. Axion inflation [49,79] The potential for the axion/natural inflation is where, (Λ, f ) are the scale of inflation and axion decay constant of this present model. By tuning the value of the decay constant, this model marginally consistent with the recent observation [78].
To be consistent with CMB data, we consider two sample super-Planckian values of the axion decay constant, f = (10, 50)M p . The scale of this inflation, Λ fixes by the CMB normalization.
Initial conditions for phase-I: The initial conditions to solve the differential equations for effective non-perturbaive era are set at the end of inflation to be, . In addition the inflationary e-folding number, N k and tensor to scalar ratio, r k for natural inflation model are expressed in terms of scalar spectral index and model parameters as The initial condition for the phase-II dynamics will be the same as chaotic inflation has given in Eq.73.     (30), (26) and (34) accordingly. Let us point out again that the first two values correspond to inflaton decaying into the scalar particle, and the third one corresponds to decaying into a pair of fermionic particles.
Here again, from the left panel of Fig.5, one concludes that if the universe undergoes two-phase reheating, considering the specific interaction during reheating the initial ω ef f during phase-I must lie within (10 −3 , 10 −6 ).
The lower limit of n s has been set by the minimum possible reheating temperature due to of reheating parameters. Therefore, the more we decrease the error of the inflationary parameter more precisely, we will be able to fix the reheating parameters.  C. α−attractor model [52] This is a new class of models that unifies many of the existing inflationary models in a single framework and was first proposed in [52]. This is currently the most favored model from the observational point of view. A class of α− attractor potential, known as the E−model, is given as Where the mass scale Λ is fixed from the CMB power spectrum. An important feature of this class of potential is a large plateau region for the large field value. It also predicts a very low value of the scalar-to-tensor ratio for different n and α. However, it is worth noting that for n = 1, α = 1, this model reduces to the Higgs-Starobinsky model. So the form of the potential for the Higgs-Starobinsky model is as follows,   where the dimension full parameter β takes the following forms, Prefixes, S, H stand for Starobinsky and Higgs model, respectively. The aforementioned coupling parameters appear in the non-canonical Lagrangian are as follows, where, R J is the Ricci scalar in the Jordan frame. For the Higgs inflation model one assumes Initial conditions for phase-I: Initial coditions to solve the differential equations for the effective non-perturbative era in the context of present model can be expressed as, where Λ = M p 3π 2 rA s 2 2n(1 + 2n) + 4n 2 + 6α(1 + n)(1 − n s ) 4n(1 + n) The inflationary e-folding number, N k and tensor to scalar ratio, r k can be written interms of inflationary spectral index(n s ) as, D. Minimal plateau inflation model [57] The minimal plateau inflationary model is a non-polynomial modification of the power-law chaotic potential. The potential for this inflation is given by, the latest PLANCK data [2]. For numerical purpose, we consider φ * = (0.001, 0.1)M p with n = 2.
Initial conditions for phase-I: The initial conditions are set as, where We set Λ = 1, except for n = 4. Constraining the parameter Λ for n = 4 has been studied in the context of minimal Higgs inflation in [64]. The inflationary parameters N k and r k can be written as,   Similar to the other inflationary models, the initial conditions the second phase boundary condition is set at the normalized scale factor at A npre thought the equation Eq. 73. As we mentioned earlier, our main intention is to see the modification in reheating parameters (T re , N re ) in comparison with the usual analysis.
On the other hand our numerical analysis estimates the value of Γ cri φ = (2.3 × 10 3 , 4.8 × 10 −7 ) GeV for φ * = 0.01M p , and Γ cri φ = (394.7, 2.7 × 10 −8 ) GeV for φ * = 0.001M p . As discussed for other inflationary scenarios, for each model parameter value of φ * two bracketed values of Γ cri φ are calculated for ω ef f = (10 −3 , 10 −6 ) respectively. The reheating temperature linked with the decay width Γ cri φ , can be found to be T cri temperature, the bound on the inflationary parameters are given in the table III.

VIII. CONSTRAINING THE INFLATON COUPLING PARAMETERS
So far, we have discussed mainly understanding the reheating parameters and their constraints from reheating. In this section we qualitatively translate those results into constraints on coupling parameters (g = g/m φ , y, h) corresponding to specific inflaton-scalar interactionsgm φ φχ 2 , yφχ 3 , and inflaton-fermion interaction hφψψ respectively. So far, our analysis was independent of the specific inflaton interaction model. Therefore, the inflaton decay width was a free parameter with one-to-one correspondence with the reheating temperature. Constraining reheating models is very challenging from the perspective of its observational limitations. Therefore, indirect constraints on the inflaton coupling parameters through reheating dynamics would be significant from the model building point of view. Reheating temperature directly estimates the allowed ranges of dimensionless coupling parameter via the inflaton decay constant Γ φ . In this section for illustration only considers two observationally viable inflationary models: Higgs-Starobinsky and minimal plateau models with n = 2, which are consistent with the current observational bound on r < 0.064 [2].
Bounds on couplings: The constraints on different coupling constants are shown in figure 8.
Plots show how the dimensionless coupling parameterg, h, and y are intimately linked with CMB anisotropy via the inflationary observables such as (n s , r k ) for different types of reheating dynamics.
The mapping T re → Γ φ → (g, y, h) are directly followed from Eqs. (22,28,32,52). From these equations, we obtain the constraints on the coupling parameters with respect to the inflationary parameters. Any realistic scenario of reheating should include all possible inflaton coupling based on underlying symmetry. Therefore, the assumption of a specific inflaton coupling's contribution to be the dominant one throughout the entire reheating period may not be relevant. Hence, a more pragmatic approach would be to construct particle physics motivated models which we left for our future study. However, as a toy model analysis, the present study may guide us in for ω ef f = 10 −3 and the associated perturbative constraints for Yukawa interaction is h cri (model) 1.33 × 10 −5 . Therefore, we can infer from this observation that our two-phase reheating scenario essentially captures the necessary features of the non-perturbative phase.
So far, we have discussed the reheating dynamics considering inflaton and radiation as the two dynamical components. However, as we all know, dark matter is another important constituent of our present universe. One of this component's important properties is that it's coupling with the standard model fields must be very weak. Apart from this, not much is known about its other fundamental properties, such as charge, mass, and coupling. Experimental searches of this particle are going on across the globe without much success till now. The searches include both directly as well as indirectly observing the properties of this object and, finally, jointly constrain the parameter region. This paper will study the dark matter phenomenology based CMB parameter space following our previous work [30]. We essentially generalize our two-phase reheating formalism and include the dark matter as the third dynamical matter component.

IX. UNIFYING THE DARK SECTOR
In the previous section, we discussed the two-phase reheating process, where inflaton decays only into radiation. In the present discussion, we add additional dark matter components and discuss the impact on dark matter phenomenology. The assumption is that inflaton decays into radiation and then radiation to dark matter. The methodology of the analysis will be the same as before, except the new additional dynamical equations for dark matter.
phase-I:(Effective non-perturbative phase) dynamics is governed by where the new component ρ X is the energy density of the dark matter particle with mass M X and energy of the dark matter is expressed as E X = M 2 X + 9T 2 [70]. T is the temperature. The above equation can be written in differential form as, Besides the above equation, we consider additional conservation equation characterizing the dynamics of every individual energy components during this phase as, To solve the above equations (91) and (90), we need one more condition. We define the ratio of the dark matter and the radiation energy density as γ = ρ X ρ R . After combining the above two equations one finds, At the initial stage of the reheating, radiation energy density must be very small ρ R 0. Hence, as discussed extensively for the two component reheating, here also ω ef f must assume the value very closed to the inflaton equation of state ω φ , at least near the beginning. In terms of dimensionless variable this phase can be written as here the dimensionless dark matter density X = ρ X E X a 3 . phase-II (perturbative phase) The subsequent perturbative phase will now be governed by two more parameters related to the dark matter component. Apart from the inflaton equation of state ω 1 φ and the inflaton decay constant Γ φ , we have a thermal average of dark matter annihilation cross-section σv , and the dark matter mass M X . The corresponding dimensionless comoving energy densities' dynamics will be governed by the Boltzmann equation [29].
The equilibrium number density of the dark matter particle can be described in terms of the modified Bessel function of the second kind [70] n eq X = and the constants c 1 and c 2 are delineate as, We consider fermionic type dark matter particles with internal degrees of freedom g.
Initial conditions: The general form of the initial conditions during the first phase of reheating (phase-I) are, The initial values of the energy densities for the phase-II will be set at the normalized scale factor As described in detail in section V, radiation energy density is again assumed to be 50% of the total comoving energy density right after the completion of phase-I. Therefore, the dark matter number density will automatically be fixed for a given γ value. All the required equation of states for two different phases are assumed to take the same approximate values ω 1 φ 0.2 and ω φ 0. The methodology of solving the dynamics will be the same as before except some additional constraints in the dark sector after the end of reheating.
Boundary condition from observations : The condition for ending the reheating dynamics is set by the following equation, supplemented with the observational constraint namely the relation between reheating temperature follows from the above equation and present CMB temperature T 0 = 2.7K 2.35 × 10 −13 GeV though the relation Eq.55. Further additional observational constraint is the observed value of the dark matter abundance defined as Ω X [85,86] which is expressed in terms of radiation abundance Ω R (Ω R h 2 = 4.3×10 −5 ). T F is the temperature at very late time when both dark-matter and radiation energy components become stationary.
While solving the Boltzmann equations during perturbative reheating (phase-II), these condition will constrain the dark matter parameter σv (thermal average of the cross-section times velocity) for a fixed value of the dark matter mass, M X , and the inflaton decay constant in terms reheating temperature. The detailed analysis only on phase-II has already been done in [29] including the dark matter phenomenology. Nevertheless we only consider dark matter production via freeze-in mechanism. This mechanism indicates that the dark matter will never reach equilibrium with the thermal bath. This kind of dark matter is known as FIMP (feebly interacting dark matter) [87]- [90]. We can illustrate the production of dark matter via Freeze-in mechanism through the heavy mediator during reheating is sensitive to the early history of the universe before the UV dominated era [91, 99, 109-111, 113-116, 118].
Physical constraints: Further constraints on the dark matter parameter space will be inherited if one considers various theoretical limits on the scattering cross-section. Cross-section can not be arbitrarily large. Perturbative unitarity usually limits the cross-section σv in term of mass, σv max = 8π 73], which are shown by pink solid lines in Figs.(9, 10). On the other hand, we will also have another bound on the cross-sections coming from the fact that during reheating dark matter production peaks around the temperature of T * = M X 4 [70]. This provides a natural condition on the dark matter number density n X (T ) < n eq X (T * ) as for T < T * the dark matter production would be frozen, and it must be diluted subsequently due to the expansion of the universe. Aforementioned condition on the dark matter number density sets an upper bound on the cross-section σv ≈ σv T =T * [70] [100] σv * ≤ 7 × 10 −14 2 g g * (T * ) 10 We call it as reheating bound in the plot. This condition is depicted by black solid lines and black dotted lines for perturbative and two-phase reheating scenarios respectively in Figs. (9,10). We have shown both of these bounds in the subsequent plots for different inflation model once we fixed the dark matter mass. Interesting observation that can be made from this theoretical constraint is that given a dark matter mass, the perturbative unitarity bound and the dynamical condition Eq.104 modify the possible range of allowed n s values obtained from the previous analysis.
This can directly shed light on the inflationary model building. Conversely, one can state that for a given inflationary model, CMB can shed light on the possible nature of the dark matter candidate via reheating phase.
Nevertheless, unifying the dark sector into a single reheating framework is the primary motivation of this section. The basic philosophy is to look into further constraints on the dark-matter parameter space in a more realistic framework of two-phase reheating dynamics and compare with that of the usual perturbative reheating analysis Maity:2018dgy. An important outcome is the constraints dut to CMB temperature anisotropy. Particularly, constraints imparted on the dark matter and inflationary parameter space ( σv − n s ) by the CMB anisotropy could enable us to constrain the viable inflationary models though dark matter observable. Conversely, given a viable inflationary model, CMB can potentially shed light on the possible properties of dark matter.
keeping this in mind, we study dark matter phenomenology considering two observationally viable inflationary models: Higgs-Starobinsky and minimal plateau models, which are consistent with the current observational bound on r < 0.064 [2].

Higgs-Starobinsky model [71][72] and dark matter phenomenology
We have already discussed about the model in the previous section VII C, and the constraints on the reheating parameters (N re , T re ) in terms of spectral index (n s ). The inclusion of dark matter does not affect much on those parameters. Therefore, the main constraints will be on the thermally averaged cross-section times velocity ( σv ), and the dark matter mass M X . The first two plots of fig.9  From first two plots of fig.9, we read the variation of the cross-section for two different effective equations of state ω ef f . As we decrease the value of the ω ef f from 10 −3 → 10 −6 , the e-folding number N npre , which is nearly independent of inflationary parameter changes from 5.8 → 12.2.
Another interesting consequence of the Phase-I dynamics is the maximum possible value of dark matter mass M max X . To understand the underlying reason behind the origin of M max X , we have computed analytic expressions considering relativistic dark matter. The dark matter number density at the point of freeze-out n f X (see appendix C) is expressed as where expressions of various symbols are given the appendix.
x f = A f /A npre and A f is the normalized scale factor when both comoving dark matter and radiation component become constant.
By using the above expression, we can obtain dark matter abundance as The above expression indicates that the dark matter abundance increases with increasing dark matter mass. Moreover, at a particular value of the dark matter mass, the dark matter component's initial number density (n in X ) will also play in the final value of the observed dark matter abundance, Ω X h 2 = 0.12. It can be observed from the equation (106), if M X > M max X , then Ω X h 2 always ≥ 0.12. Therefore the maximum possible dark matter mass can be obtained from the above equation considering Ω X h 2 = 0.12 as , which is dependent on the initial dark matter number density for the phase-II evolution, For a given γ, the initial value of X is clearly set by the value of ω ef f . Therefore, for a given value of phase-I dynamics parameters ω ef f and γ, a particular value of the dark matter mass B. Minimal plateau inflation model [57] The details of this model are discussed in section VII D. As was already the case for Higgs's inflation, for the minimal inflation model also the reheating parameters such as (T re , N re ) will not be modified much because of dark matter dynamics. The reason being, the contribution of dark matter in the background evolution during the reheating phase is insignificant. Throughout this analysis, we consider φ * = 0.001M p , which satisfies the CMB observation. Another motivation is that as one increases φ * value, the models assume simple power law. Details constraints on the dark matter parameter space can be read off from the Fig.10.   In this paper, we propose an effective two-phase reheating scenario. After inflation, reheating has been studied extensively in the literature, either through a perturbative or non-perturbative approach. However, it is believed that both approaches independently should not capture the complete picture of the complicated dynamics. In this paper, we, for the first time, study this phase to the best of our knowledge, taking into account both the approaches together motivated by our previous work [31]. However, instead of considering explicit non-perturbative decay of the inflaton field through parametric resonance, we model the initial phase by effective dynamics governed by the standard conservation laws and parametrized by a constant effective equation of state (ω ef f ).
The combined form of conservation laws and the initial condition of the reheating dynamics put constraints on the effective equation of state during the effective non-perturbative process calling it as phase-I. However, during perturbative analysis due to explicit decay of the inflaton field into radiation, we obtain the non-trivial time-dependent effective equation of state. At this stage, let us remind the reader that in all the PLANCK analysis [58] on constraining the inflationary models w ef f is assumed to be a constant free parameter during reheating, which follows from the proposal described in [6]. What we argue is that those assumptions should not be correct. After inflation, every inflationary model has its own characteristic oscillatory period, which contributes to the equation of state during reheating. Therefore, considering w ef f as a free parameter loses some of the fundamental characteristic properties of the inflaton potential itself. Furthermore, if reheating occurs for a longer period of time, the time-dependent w ef f during the perturbative process should also be very important to get precise constraints on any inflationary model. This is where our analysis not only can play an important role in better understanding the inflationary models but also opens up the possibility of understanding the micro-physics of the reheating process through CMB physics. As we can clearly see how the CMB power spectrum constrains the value of inflation-radiation coupling parametrized by Γ φ through reheating temperature T re . The usual connection between Γ φ and T re will not be correct any more once we consider the decaying inflaton as it is a well-known fact that during the reheating process, even at the end of reheating time, Γ φ = H, inflaton does not decay into radiation completely. Therefore, one certainly needs to take into account this fact while calculating T re and its connection with the scalar power spectrum n s in the analysis. However, all the previous theoretical as well as in PLANCK analysis, complete decay of inflaton is assumed while relating the cosmological scales exiting and re-entering the horizon at two different time scales. Therefore, based on the two-phase reheating scenario, our prediction of reheating temperature corresponding to the inflationary power-spectrum is more accurate than the previous analysis.
At first, we analyzed the viable constraints on the decay width as well as reheating parameters (N re , T re ) considering the decay of the inflaton field in the perturbative Boltzmann framework. Perturbative dynamics have been shown to give rise to a maximum reheating temperature T max re 10 15 naturally, which essentially corresponds to almost instantaneous reheating. As long as the decay width is in the perturbative regime, the result from the only perturbative process is trustworthy. However, because of the straightforward relation between T re , and Γ φ , high reheating temperature limit can correspond to non-perturbative phenomena. This fact motives us to include non-perturbative aspects of reheating through effective dynamics. In our present scenario, the universe passes through two distinct phases during reheating. Combining the inflation and subsequent standard big-bang evolution with the intermediate two-phase reheating, our approach predicts the critical value of the inflaton decay constant Γ cri φ depending upon the phase-I equation of state ω ef f . The critical point naturally defined at N npre = N pre . Therefore, if Γ φ < Γ cri φ , the reheating phase will be dominated by perturbative one and vice versa. We also compare our numerical results of Γ cri φ with the critical decay width obtained from the theoretical consideration for different type of inflaton-reheating field interactions gφχ 2 , yφχ 3 , hφψψ. It turns out that all the theoretical values re ω ef f = 10 −3 3.5 × 10 10 7.2 × 10 10 4.0 × 10 10 9.0 × 10 10 2.2 × 10 10 8.8 × 10 9 T cri re ω ef f = 10 of Γ cri (model) correspond to an effective phase-I equation of state ω ef f within 10 −3 ∼ 10 −6 . Our actual lattice simulation results also appeared to be compatible with this conclusion (see Fig.3).
A summary table-VII for Γ cri φ is given for three observationally viable model. The inclusion of the initial non-perturbative phase naturally changes the maximum reheating temperature value because of its perturbative definition. T max re is no longer defined at the point of instantaneous reheating N re 0, rather is defined at N re ≈ N npre , which is equivalent in saying the phase-II e-folding number N pre 0. At the end of phase-I, approximately 50% of the total comoving energy density remains in the form of the inflaton, which naturally leads to different T max re defined in the perturbative phase-II dynamics. This phase further sets the final equation of the state of the system to 1/3. All these results have been shown to be crucially dependent upon the phase-I effective equation of state ω ef f . As one changes the value of ω ef f from 10 −3 → 10 −6 , the phase-I e-folding number N npre changes from 6 → 12. The maximum reheating temperature T max re accordingly changes from (10 13 → 10 10 ) GeV. Therefore, the conclusion that can be emphasized upon is that the value of reheating temperature may encode the information about the non-perturbative phase. Furthermore, all the inflationary models which are compatible with the observed CMB anisotropy, predict the same maximum reheating temperature (T max re ) for a given ω ef f . This is reminiscent of the maximum reheating temperature T max re obtained for purely perturbative reheating dynamics irrespective of the inflation model. Keeping this point in mind, we have performed a comparative analysis of different existing reheating formalisms such as conventional reheating dynamics (case II) and purely perturbative analysis (case III) with our proposed two-phase (case-I). For both cases II and III, the model-independent maximum value of the reheating temperature turns out to be T max re 10 15 GeV, which is reduced to 10 13 ∼ 10 10 GeV when considering two-phase reheating for ω ef f = 10 −3 ∼ 10 −6 . Further, two-phase reheating  M max X (maximum) 6.8 × 10 8 6.7 × 10 8 6.8 × 10 4 6.7 × 10 4 6.8 6.7 scenario constraints the inflation model within a very narrow range of allowed scalar spectral index compatible with CMB anisotropy.
Further generalization has been analyzed by including the dark matter component as one of the decay products of the inflaton. Depending upon the mass dark matter annihilation cross-section versus scalar spectral index parameter space has been shown to be reduced because of two-phase reheating as compared to that of standard reheating dynamics, which can be observed from Fig.9 and 10. Details of the allowed parameter space for various models can be obtained from the tables V and VI. Because of the non-trivial initial condition for two-phase dynamics, there exists a maximum possible mass M max X above which dark matter turned out to be overproduced no matter how small the annihilation cross-section is assumed. In the summary table VIII, we provide numerical values of maximum possible dark matter mass allowed for different viable models under consideration. summary As just stated, the value of M max X is directly connected to γ (γ = ρ X ρ R ), which is defined during phase I. Once we fixed the spectral index for a particular inflation model and these values of M max X is nearly independent of the choice of ω ef f . Nonetheless, one important point we should understand that the existing reheating scenarios, either perturbative or non-perturbative, are not the complete description of this phase. A unified description that connects both non-perturbative and perturbative dynamics is more appropriate.
In our present study, we, for the first time, try to construct such a unified description. As a first attempt towards this goal, we describe non-perturbative preheating dynamics by effective dynamics. Our present formalism is particularly suited for the class of inflation models with quadratic potential near its minimum. For inflaton potential with a power greater than two, lattice results generically predict the equation of state 1 3 after the end of non-perturbative dynamics [31]. So for those models, our two-phase reheating is no applicable. We will be considering this case in our future work. Instead of considering an effective non-perturbative approach, actual nonperturbative dynamics integrated with perturbative one would be more appropriate. Recently an interesting approach has been proposed to describe preheating phenomena in the Boltzmann framework [80]. In our present two-phase reheating dynamics, the aforementioned non-perturbative Boltzmann framework could be natural to integrate with the perturbative Boltzmann equations.
Another important fact we have not considered is the temperature dependency of the effective numbers of relativistic degrees of freedom (g * ). Constant effective degrees of freedom is reasonably good approximation for a wide range of temperature [102] [103] till the QCD hadronic transition happens at around 10 2 MeV scale, around which the vale of effective degrees of freedom changes as g * = 100 → 10 [104] [105]. So our eventual plane in the future to calculate dark matter and reheating parameter space accurately by acknowledging the precise evolution of those degrees of freedom in the thermal bath [106]- [108].
In order to solve analytically, we assume the inflaton energy density to follow the equation, .

(A4)
Here Γ φ is the time-independent inflaton decay constant. Notice that the effect of decay constant is being ignored assuming the fact that at the initial stage of perturbative reheating inflaton energy is the dominant one. ρ i φ and t i are initial density and initial time during the perturbative era respectively. Using the above equation the radiation energy can be solved as follows, d ρ R a 4 = Γ φ ρ φ (1 + ω 1 φ )a 4 + 2 E X σv n 2 X − n 2 X,eq a 4 dt φ da H + 2 E X σv n 2 X − n 2 X,eq a 3 da H . (A5) Using the following expression for the Hubble parameter, where ρ in R is the initial radiation density at the beginning of perturbative phase. For the reheating temperature computation we ignore the effect of dark matter whose contribution has been verified to be negligible in our full numerical computation. By solving Eq.A5 we obtain, In the above expression, we neglected higher-order terms of ρ in R /ρ in φ . Additionally in terms of radiation temperature T rad = 30 π 2 g * ρ R 1/4 the above equation transforms into following expression, Here x, β, c and H in defined as The maximum radiation temperature can be found by taking derivative of the above equation (A8) with respect to x and set it to zero In the limit of ρ φ /ρ R 1 (perturbative approximation), the values of x at the point of maximum radiation temperature appears as (A11) In our present analysis, the expression of x associated with maximum radiation temperature leads to the following relation where . Now after replacing the expression of x max into the above Eq.A8, the maximum radiation temperature turns out as . (A14) In the above expression we have neglected higher order terms of ρ in R /ρ in φ . Next, we will try to express all initial densities in terms of the inflaton energy density at the end of the inflation ρ end φ . The effective non-perturbative phase-I dynamics solves the radiation and inflaton energy density in terms of ρ end φ . Therefore, during phase I the dimensionless radiation energy density R I (A) can be correlate with inflaton energy density Φ I (A) (using (36), (37) and (41)) as The initial densities during phase II (perturbative era) in terms of dimensionless comoving energy densities are identified as where, A npre is the normalized scale factor at the end of the effective dynamics. A npre is defined when the dimensionless comoving radiation energy density becomes 50% of the total comoving energy density, R(Anpre) Φ(Anpre)+R(Anpre) 1 2 =⇒ Φ(A npre ) R(A npre ). Using Eq.A15 one can find A npre and corresponding e-folding number N npre as From our analytic expression above, we obtain N npre ∼ (5.8, 12.7) for two values of ω ef f = (10 −3 , 10 −6 ) accordingly. These values of the e-folding number during phase I almost exactly match with our numerical result.
The final expression for the maximum radiation temperature in terms of comoving energy densities is given by where

(B4)
Reheating temperature is defined when the inflaton field comes in thermal equilibrium with the radiation bath at the point, In the above equation, we ignore the contribution of inflaton energy density to be negligible. Using the expression for the radiation energy density we can obtain the decay width as follows, In the earlier expression, we can ignore the second term in the third bracket since x re 1 for most of the values of the spectral index. As a result, the Γ 2 φ can now be written as The average energy of the single component dark matter at the point of freeze-in can be expressed as