Measuring Inflaton Couplings via Primordial Gravitational Waves

We investigate the reach of future gravitational wave (GW) detectors in probing inflaton couplings with visible sector particles that can either be bosonic or fermionic in nature. Assuming reheating takes place through perturbative quantum production from vacuum in presence of classical inflaton background field, we find that the spectral energy density of the primordial GW generated during inflation becomes sensitive to inflaton-matter coupling. We conclude, obeying bounds from Big Bang Nucleosysthesis and Cosmic Microwave Background, that, e.g., inflaton-scalar couplings of the order of $\sim\mathcal{O}(10^{-20})$ GeV fall within the sensitivity range of several proposed GW detector facilities. However, this prediction is sensitive to the size of the inflationary scale, nature of the inflaton-matter interaction and shape of the potential during reheating. Having found the time-dependent effective inflaton decay width, we also discuss its implications for dark matter (DM) production from the thermal plasma via UV freeze-in during reheating. It is shown, that one can reproduce the observed DM abundance for its mass up to several PeVs, depending on the dimension of the operator connecting DM with the thermal bath and the associated scale of the UV physics. Thus we promote primordial GW to observables sensitive to feebly coupled inflaton, which is very challenging if not impossible to test in conventional particle physics laboratories or astrophysical measurements.


Introduction
Inflation stands as one of the most fundamental pillars of contemporary cosmology, explaining several puzzles of the early Universe, for example, the horizon or the flatness problem [1,2].In its simplest form, inflation can be described by a single, slowly-rolling scalar field with an approximately flat potential that dominates the energy density of the primordial Universe.The flatness of the potential is generally expressed in terms of the socalled slow roll parameters, which are connected to cosmological observables measured by the spectrum of the cosmic microwave background (CMB).Cosmic inflation dilutes any preexisting matter and radiation and thus requires a reheating mechanism to eventually result in the Universe dominated by radiation.Perturbative reheating can be realized through interactions between the classical, coherently-oscillating inflaton and thermal bath.In this framework, the Standard Model (SM) particles are produced in a quantum process from the vacuum in the homogeneous inflaton background.Hereafter, this process is dubbed as the inflaton decay.The common assumption of perturbative reheating, a constant decay width of the inflaton, cannot be justified in a generic scenario, e.g., if inflaton φ oscillates in a potential of the form V (φ) ∝ φ 2n , with n = 1 [3,4] or a daughter field acquires a vast mass due to its coupling with φ, which generates kinematic suppression [5,6].
The primordial gravitational wave (GW) is one of the most crucial predictions of inflationary paradigm.During inflation, quantum fluctuations inevitably give rise to a scale-invariant spectrum of tensor metric perturbations at super-Hubble scales.In a standard post-inflationary scenario, tensor modes become sub-horizon during the radiationdominated (RD) stage.It is, however, well established that the presence of non-standard cosmologies, e.g., an early stiff era before radiation domination, breaks such a scale invariance.In this case, the GW spectrum becomes significantly blue-tilted in the frequency range corresponding to the modes crossing the horizon during the stiff period [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] 1 .In addition, the amplitude of the tensor power spectrum becomes considerably enhanced, as compared to the standard scenario, where the Universe immediately transitioned from the inflationary into radiation-dominated phase.For instance, the onset of the RD epoch can be delayed if the reheating process is not instantaneous.This happens in a class of reheating models with a time-dependent inflaton decay rate.
The fact that Dark Matter (DM) constitutes about 24% of the matter-energy budget of the Universe, has been unequivocally established from several cosmological and astrophysical observations [24][25][26].As an alternative to the weakly interacting massive particle (WIMP) paradigm [27,28], feebly interacting massive particles (FIMPs) have gained quite an attention [29][30][31].Freeze-in, as opposed to freeze-out, requires very suppressed interaction rates between the dark and visible sectors, which can be achieved either via small couplings, known as IR freeze-in [30], or via non-renormalizable operators, suppressed by a high mass scale, called UV freeze-in [30,32].The latter scenario is particularly interesting, as the DM yield is sensitive to the highest temperature reached by the SM plasma, controlled by the dynamics of the inflaton decay.The highest temperature of the bath, on the other hand, depends on the nature of the inflaton coupling or equivalently, the decay rate.
Motivated by the above arguments, we explore a generic reheating scenario with the time-dependent inflaton decay width whose dependence on the scale factor is parameterized by a power-law dependence [4,5,33,34].The non-standard evolution of the inflaton decay rate results in the non-trivial behavior of the inflaton and SM energy densities in the post-inflationary phase.We consider several types of inflaton interactions with bosonic and fermionic fields, and for each case, we obtain the size of the inflaton-matter coupling that could be probed by the proposed GW detectors, satisfying bounds from the Big Bang Nucleosynthesis (BBN) and CMB on GW energy density.We finally discuss DM production via UV freeze-in under the influence of time-dependent inflaton decay width and show that it is possible to have DM over a wide mass range ∼ O(MeV) − O(PeV), satisfying the PLANCK observed relic abundance.Interestingly, since both DM dynamics and GW energy density are controlled by the inflaton decay width, the proposed GW detectors can probe inflaton couplings that can give rise to the correct relic abundance for the DM, making a direct correspondence between DM and GW detection prospects.
The paper is organized as follows.In Sec. 2, we discuss the thermal history of the Universe in the post-inflationary regime, elaborating on the underlying inflaton-matter interaction.The derivation of the spectral energy density of the primordial gravitational wave from inflation is concerned in Sec. 3. Subsequent results of the GW analysis are furnished in Sec. 4. The impact of time-dependent inflaton decay on dark matter abundance is discussed in Sec. 5. Finally, in Sec.6, we summarize our findings.Detailed calculations of inflaton decay are presented in Appendix.A, while Appendix.B discusses the bound on the scale of inflation from inflationary observables.

The Framework
Let us consider the following action of a system composed of the inflaton φ and the SM sector: where L φ ,SM ,int are respectively the Lagrangian densities corresponding to the inflaton, the SM field, and relevant interactions, which we will elaborate in a moment.We assume that all matter fields are minimally coupled to gravity so that their kinetic terms are canonically normalized.We denote the Ricci scalar via R, and M P ≡ 1/ √ 8πG = 2.435×10 18 GeV is the reduced Planck mass.The background metric is assumed to be in the FLRW form.In the above expression, g denotes the determinant of the metric tensor, whose line element is with a(t) being the scale factor.The inflaton Lagrangian reads with V (φ) being the inflaton potential.In this work, we focus on the α-attractor T-model of inflation [35], where V (φ) is given by where Λ is the scale of inflation, and M = √ 6 α M P .Hereafter the parameter α will be fixed at 1/6.In the following discussion of the inflationary period V (φ) = Λ 4 will be adopted, while during reheating we will use V (φ) = Λ 4 (|φ|/M ) 2n .The Lagrangian L int , describing the interaction between the inflaton and matter shall be specified shortly.

Post-inflationary evolution of the Universe
We first discuss the evolution of energy density of inflaton and radiation during the era of reheating, when the massive scalar field φ (inflaton) governs the evolution of the SM radiation energy density.The dynamics of such a system is controlled by the following set of coupled Boltzmann equations (BEQ) for time-averaged (over one period of inflaton oscillations [6]) energy densities of inflaton (ρ φ ) and radiation (ρ SM ) as where ±(1 + w)Γ φ ρ φ accounts for the energy gain (loss) of the SM bath (the inflaton field).Note that above, we have assumed that φ mainly transfers its energy into the SM sector.Here the dots denote derivatives with respect to the cosmic time t, and w is the time-averaged equation-of-state (EoS) parameter for the inflaton where ... denotes the time average, with p being the inflaton pressure.For a given inflaton potential, the EoS parameter is calculable, e.g., for monomial potentials, it can be related to the index n through the following formula: w = (n − 1)/(n + 1).At the onset of the reheating phase, the total energy density of the Universe is mainly in the form of the inflaton, while its end is defined as a moment when the inflaton and radiation energy densities become equal.In this work, we limit ourselves to the class of reheating models characterized by a stiff EoS parameter, i.e., w > 1/3.The evolution with w > 1/3 is referred to as the stiff epoch.In Fig. 1, we have schematically shown different periods of the evolution of the Universe and their corresponding equation of states.We adopt the following [5] parametrization of the inflaton width Γ φ : where a e is the initial value of the scale factor, indicating the end of inflation, Γ e φ denotes the inflaton width at a = a e , while β is assumed to be a constant parameter.The above parametrization of the inflaton width can be systematically derived for a homogeneous inflaton condensate with a power-law evolution of the inflaton energy density, assuming massless final states, as elaborated in Appendix.A adopting the α-attractor T-model of inflation.More generically, such a parametrization is valid for osciallation of the inflaton in any monomial potential during reheating.In addition, as shown in [3,6], the generic power-law form of the inflaton width is also well justified if the time dependence of Γ φ arises due to kinematical suppression effects for non-zero final state masses.Further, a general time (temperature) dependent dissipation rate can have several physical origins as discussed in Sec.3 of Ref. [33].It is worth to realize and emphasize that corrections that originate from kinematical mass effects do not change the form of a-dependence of the inflaton width, it turns out they merely modify the power β in (2.7).In the following part, we will show how β can be connected to the underlying inflationary potential.The Hubble rate, H, evolves according to the Friedmann equation The initial conditions for the inflaton and SM radiation are given by ρ φ (a e ) ≡ ρ e , ρ SM (a e ) = 0 . (2.9) At early times, the term associated with the expansion, i.e., 3 (1 + w) Hρ φ , typically dominates over the interaction term (1 + w) Γ φ ρ φ .The right-hand side of the first Boltzmann equation (2.5) can then be neglected, and straightforward integration gives a solution for the inflaton energy density during the oscillatory phase , a e < a ≤ a rh .
(2.12) which can further be written as where we have used the value of the radiation energy density at the end of reheating, defined by the inflaton-radiation equality There is a comment here regarding continuity of ρ RH SM (a) as a function of n at β = (n + 4)/(n + 1).Within any given inflaton model, β is a fixed function of n, e.g., in the case of the inflaton Yukawa interactions β = β ψ = 3(n − 1)/(n + 1).Then, ρ RH SM (a) should be a continuous function of n, including the values of n that follows from the condition β = (n+4)/(n+1).Indeed ρ RH SM (a) given by Eq. (2.12), does satisfy this condition.To prove the continuity, it is essential to keep both terms in the bracket of Eq. (2.12).However, to perform further analytical calculations, we will be dropping one of those terms at a time.Therefore cases β (n + 4)/(n + 1) and β (n + 4)/(n + 1) will appear, and it is important to note that such approximations are not applicable in the close vicinity of β = (n + 4)/(n + 1).
As the inflaton energy density dominates the early stages of reheating, the background dynamics is determined by the value of w, whereas the behavior of the SM sector depends on both w and β.Note that, for β < (n + 4)/(n + 1), the first term in the square bracket in (2.12) dominates, whereas if β > (n + 4)/(n + 1) the energy density of the radiation bath decreases as a −4 .Thus, depending on the value of β + 3 n/(n + 1), one obtains a very different evolution of the radiation sector, which implies a non-trivial scaling of a thermal bath temperature T with the scale factor during reheating.In particular, the temperature of the SM sector is measured by the radiation energy density where g (T ) counts the effective number of relativistic degrees of freedom at temperature T .Utilizing (2.12), one finds that during the non-standard phase of reheating, T behaves for a e < a ≤ a rh as (2.16) The scale factor at the end of reheating can be determined from Eqs. (2.10) and (2.12) as follows where W[z] denotes the Lambert W-function.Note that in order to obtain the above result, we had to drop either the first or the second term in the bracket in Eq. (2.12).During non-instantaneous decay of the inflaton, the bath temperature can rise to several orders of magnitude above the reheating temperature [36].This happens at a = a max , defined via so that the maximum temperature of the radiation bath is From the discussion so far, it is understandable that once one specifies the inflaton potential V (φ) and the form of the inflaton-SM interactions (the "model"), w, H e , β, Γ e φ can be determined.
Before proceeding further, let us briefly discuss the details of inflaton dynamics during the oscillatory phase.In the FLRW background, the classical equation of motion (EoM) for the inflaton is given by where we have assumed that φ is spatially homogeneous2 .During the period of reheating, the oscillating inflaton field with a time-dependent amplitude can be parametrized as [3,[37][38][39][40][41].
φ(t) = ϕ(t) .P(t) . (2.21) Here, P(t) is a quasi-periodic, fast-oscillating function, and ϕ(t) is a slowly-varying envelope.It is instructive to introduced the effective mass of the inflaton The slow-roll evolution of the inflaton field during cosmic inflation can be parameterized by the so-called potential slow-roll parameters that are related to the cosmological observables, namely the spectral index n s and the tensor-to-scalar ratio r via At the very end of inflation, ä = 0, and the first potential slow-roll parameter is V 1.This condition determines the field at a = a e , i.e., Moreover, at a = a e , the potential energy of φ matches with its kinetic energy, so that where we have assumed that P(a e ) = 1.From the above expression, we see that ρ e has a very mild dependence on the index n for a fixed α and Λ.For example, with α = 1/6 and Λ = 3 × 10 −3 M P , we find ρ e 2 × 10 63 GeV 4 for n ∈ [2 , 10].

Models of inflaton-matter interactions
Let us now specify interactions between the inflaton and other matter particles.Here, we assume that perturbative expansion is justified while calculating the quantum production of radiation (reheating) from the vacuum in the classical inflaton background.Once the inflaton couples to matter fields, its oscillations are severely damped by the decay.Without going into the details of the UV completion, we consider four possible interaction vertices 3 and parametrize them as where S , V , ψ, a are generic scalar field, vector boson, Dirac fermion, and particles with derivative interactions (e.g., axion-like particles), respectively.We refer to scenarios with an inflaton-matter coupling specified as to models and consider them individually.It is important to stress that in the adopted parametrization, the coupling g Sφ and g V φ have dimensions of mass, g ψφ is dimensionless, while h aφ has inverse mass dimension.Let us also notice that once the model of reheating is specified, Γ e φ can be calculated, see Eq. (A.21), in terms of the model coupling g iφ (i ∈ SS , ψψ , VV , aa) and parameters of the inflaton potential, which also determines the value of the β parameter.Moreover, the time-averaged EoS parameter w is predicted solely by the shape of the inflaton potential Eq. (2.4).Within a given model of inflaton-matter interaction, we adopt the following set of independent parameters4 {Γ e φ→f , Λ, n} , where f = SS, ψψ, VV , aa.

Spectrum of Primordial Gravitational Waves
Here we briefly summarize formalism to calculate the spectrum of a stochastic GW background of primordial origin.The complete derivations can be found in, for example, Refs.[15,[42][43][44], hence without going into the details, here we highlight the salient points that help in obtaining the final expression for primordial GW spectrum.We shall mainly stick to the convention followed in Ref. [44] while elaborating on the definition of different relevant quantities.GWs are described as the tensor metric perturbations in a spatially-flat FLRW Universe Tensor fluctuations are assumed to be small perturbations, i.e., |h ij | 1, satisfying the transverse-traceless conditions,

The GWs equation of motion follows the
Einstein equations linearized to first order in h ij over the FRLW background with Π TT ij being the transverse and traceless part of the anisotropic stress tensor.It is instructive to decompose the perturbation h ij over two polarization states λ as where λ ij (k) are spin-2 polarization tensors satisfying orthonormality and completeness relations with P ij = δ ij − ki kj being the projection operator and k being the unit vector along the direction of k.
The mode function h λ (t, k) ≡ h λ k satisfies the following equation: where Π λ k is the Fourier components of Π TT ij .It is convenient rewrite the above equation in conformal coordinates where prime denotes a derivative with respect to the conformal time coordinate, defined by dt = adτ .Eq. (3.6) admits approximate analytical solutions in two extreme regimes as follows • Super-Hubble scale: For modes far outside the Hubble horizon, i.e., k aH, one can write Eq. (3.6) as where we have ignored the source term.This equation has a solution of the form (3.8) 5 In presence of a viscous background the mode equation receives correction due to (bulk and shear) viscosity, that causes damping of the GW amplitude [44,45].However, such a damping is inefficient in an expanding background as long as the interaction rate of particles in the viscous medium is much faster compared to the cosmic expansion rate [44], which is satisfied in the present analysis via instantaneous thermalization.
where C 1,2 are constants of integration.The second term in the above expression decays with time 6 .Therefore, ignoring the decaying term, one concludes that h λ k stays constant for super-Hubble modes [43,44].At some point after inflation, these modes become sub-Hubble, i.e., k > a H and re-enter the horizon.
• Sub-Hubble scale: After the end of inflation, modes eventually re-enter the horizon (k > a H) and start to oscillate.In this case, Eq. (3.6) can be solved by assuming a solution of the form where A and B are real functions.Substituting this in Eq. (3.5), and comparing the real and imaginary components, we obtain where again, we have dropped the source term.Now, considering the oscillation to be very rapid compared to the time variation of the amplitude, and the modes to be well inside the horizon k a /a, we obtain from the first equation which on substitution in the second equation brings us to The behaviour of the mode function h λ k is thus described by a WKB solution of the from where C is some arbitrary constant.
The energy density carried by the GWs is given [15,43,44,46] by where ... denotes the spatial average.As it could be seen from Eq. (3.5), once a particular mode re-enters (in other words crosses) the horizon, i.e., k > aH, the corresponding mode function obeys , therefore the energy density could be written as It is useful to define the spectral GW density per logarithmic wavenumber interval, normalized to the critical density Then, starting with Eq. (3.14), one can obtain [43,44] 7 where the primordial power spectrum has been defined as with the primordial mode function h λ k ,prim defined as the mode function h λ k (τ ) at some moment shortly after the end of inflation, when all the modes have already exited the horizon.
The primordial tensor power spectrum P T,prim is determined by the Hubble parameter at the time when the corresponding mode crosses the horizon during inflation (k = aH) [15,43,44] where we have used the mode solutions by matching the sub-Hubble modes (during inflation) with the super-Hubble modes (at the end of inflation) at k = a H [15,42].The transfer function, T (τ , k), adopted in Eq.(3.17), connects primordial mode functions with mode functions at some later time h λ k (τ ) as [42,44] (3.20) As shown earlier, the modes h λ k remain constant on super-horizon scales, while they decrease as a −1 once they re-enter the Hubble horizon.In other words, in the sub-horizon regime, i.e., k > aH, h λ k = a hc /a h hc,λ k ; therefore, the transfer function could be written [18,47,48] as where the prefactor 1/2 appears due to oscillation-averaging the tensor mode functions [44,49,50] and a hc is scale factor at the horizon crossing.As evident from Eq. (3.21), the transfer function characterizes the expansion history between the moment of horizon crossing τ = τ hc of a given mode k and some later moment τ > τ hc .From Eq. (3.17), one can see that the spectral GW energy density at the present time is given by where a 0 = 1 and H 0 are respectively the scale factor and the Hubble parameter measured today.
One anticipates that Ω (0) GW (k) has a power-law dependence on k.Using the fact that k = a hc H(a hc ), and the equation-of-state parameter after the horizon crossing is w, one obtains Hence, from Eq. (3.22) we find the following scaling applies (3.24) Let us now find the coefficient of proportionality in Eq. (3.24).This coefficient might, in particular, depend on parameters such as, e.g., the inflaton width Γ e φ at the end of inflation, which, in turn, is determined by the inflaton-matter couplings.The sensitivity of the spectral GW energy density upon the inflaton couplings is one of the central tasks of this study.We will determine the present value of the spectral GW energy density for the modes re-entering the horizon at different epochs; after the end of reheating, during RD, and prior to the onset of RD, that is, during inflaton domination (reheating).At first, let us assume that the horizon re-entry happens during the RD epoch after the end of reheating, i.e., a hc RD ∈ (a rh , a eq ).In this case, we reproduce the standard scale-invariant result for the spectral GW energy density where the function F(g ) ≡ 3 tracks the number of degrees of freedom from the end of reheating till present, and Ω (0) γ = ρ γ,0 /ρ c = 2.47 × 10 −5 h −2 is the fraction of the energy density of photons at the present epoch.Here we have assumed entropy conservation from the moment of horizon crossing (RD) till today, implying T ∝ a −1 g −1/3 s .On the other hand, if the horizon crossing happens during reheating, before the radiation era, i.e., a hc RH ∈ (a e , a rh ), one obtains where we have used the following relations: Now, for the re-entry of the modes happening during the period of reheating, redshifting the GW frequency f to present, we obtain [44] where ρ rh in the total energy density of the Universe at the end of reheating.
There exists an upper bound on frequencies, dictated by the modes that re-entered the horizon after the end of inflation with which shows the blue-tilted behaviour of the spectrum for w > 1/3, which is equivalent to n > 2. One can then, exploiting Eq. (2.12) and Eq.(2.17), rewrite ρ rh in Eq. (3.30) in terms of Γ e φ and remaining parameters (Λ, β, n), so that the spectral GW energy density as a function of the frequency measured today reads with The above expression explicitly shows how the proportionality coefficient of Eq. (3.24) depends on the inflaton decay width Γ e φ at the end of inflation.Since the only frequency dependence in Eq. (3.32) enters through f 2 n−4 2 n−1 , one can write the full spectrum today in an approximate piecewise form as where f c can be analytically computed as (3.35) Here we would like to emphasize that Eq. (3.32) is a general expression for GW energy density today, agnostic to the underlying inflaton coupling, and only relies on the fact that the inflaton decay width depends on the scale factor the way it has been assumed in Eq. (2.7).In Sec. 4 we present plots of Ω (0) GW (f ) calculated according to Eq. (3.34) as a function of f for different inflaton models and various choices of parameters {g iφ , n , Λ}.It is worth noting the presence of the constant contribution for f < ∼ f c .

Results and Discussion
So far, we have considered primordial gravitational waves (PGW) and the connection between the spectral energy density of PGW and the inflaton decay width at the end of inflation, Γ e φ , i.e., we have found a correspondence between the particle physics model and GW spectrum.Hereafter we assume that all the inflaton decay products belong exclusively to the thermal bath, and they thermalize instantaneously after being produced 8 .In particular, we emphasize that the inflaton does not decay into dark matter (DM) pairs.Below, first, we review existing experimental constraints on the parameter space.Then, we present and discuss results for GWs spectral energy density within each inflaton model defined in Sec.2.2.

Constraints on the model parameters
There exist experimental constraints on the parameter space {g iφ , Λ , n} or, equivalently, upon {Γ e φ→f , Λ, n}.They will be reviewed shortly below.

∆N eff
As it has been shown in Sec. 3, the energy density of GW background scales as ρ GW ∝ a −4 , i.e., the same way as that of radiation energy density.This implies that the GWs act as an additional source of radiation.Any observable capable of probing the background evolution of the Universe (and hence its energy content) has, therefore, the potential ability to constrain the GW energy density.In fact, it is possible to put an upper limit on the GWs energy density at the time of BBN and CMB decoupling.As it is well known, an upper bound on any extra radiation component, in addition to those of the SM, can be expressed in terms of the ∆N eff .The number of effective relativistic degrees of freedom N eff is defined through the expression for the radiation energy density in the late Universe as where ρ γ , ρ ν , and ρ GW correspond to the photon, SM neutrino, and GW energy densities, respectively, with T ν /T γ = (4/11) 1/3 .Note that, the relevant temperature here is the photon temperature after e + e − annihilation.Within the SM, the prediction taking into account the non-instantaneous neutrino decoupling is N SM eff = 3.044 [51][52][53][54][55][56][57][58][59], whereas the presence of GW implies Since ρ GW ∝ a −4 and ρ γ ∝ T 4 , the above relation can be utilized to put a constraint on the GW energy density redshifted to today via where, as before, Ω γ h 2 2.47 × 10 −5 .The bound in Eq. (4.3) is applicable for an integrated energy density as follows where f max is given by Eq. (3.29).Note that the upper limit of the integration is f max for modes that match the size of the horizon at the end of inflation, while the lower limit can be f BBN , which corresponds to the modes that match the horizon size at the time of BBN.One can then perform the integration analytically using Eq.(3.32), and for all [11,15,17,60], for n > 2. From Eq. (3.30), we note that, for n > 2, Ω (0),RH GW is inversely proportional to a positive power of the radiation energy density at the end of reheating.This implies, an increase in ρ rh should relax the ∆N eff bound on Ω (0) GW .In the following subsections we will see, for instance, a larger inflaton-matter coupling remains safe from ∆N eff constraint as it leads to more radiation production, diluting the GW energy density.
Hereafter we use orange color to denote regions of the parameter space excluded by ∆N eff measurements.Within the framework of ΛCDM, the Planck legacy data produces N eff = 2.99 ± 0.34 at 95% CL [61].In our plots, this is shown by the solid orange horizontal line.Once the baryon acoustic oscillation (BAO) data are included, the measurement becomes more stringent: N eff = 2.99±0.17.As computed in Ref. [62], a combined BBN+CMB analysis shows N eff = 2.880 ± 0.144.We denote this by the dashed horizontal line.Upcoming CMB experiments like CMB-S4 [63] and CMB-HD [64] will be sensitive to a precision of ∆N eff 0.06 and ∆N eff 0.027, respectively.These are denoted by dot-dashed and dotted lines in subsequent figures.The next generation of satellite missions, such as COrE [65] and Euclid [66], shall improve the limit even further up to ∆N eff 0.013, as shown by the large dashed line.The orange-shaded region is thus disallowed, depending on the sensitivity of a particular experiment.We collectively named them "∆N eff bounds" in all the subsequent figures.

BBN: T rh
The reheating temperature can be obtained by evaluating Eq.(2.16) at a = a rh and utilizing Eq.(2.17) where 2) .Note that when the reheating scenario is specified, β becomes a given function of n.For instance, for inflaton-fermion Yukawa coupling, β = β ψ = 3(n − 1)/(n + 1).However, regardless the explicit form of β(n), in the middle case T rh does not depend on β.When a particular model is chosen, corresponding n for the middle case is fixed, e.g., for Yukawa interaction, this corresponds to n = 7/2, while the most upper and most lower cases correspond to n 7/2 and n 7/2, respectively.To avoid conflict with BBN predictions, we require T rh 4 MeV [67][68][69][70][71][72][73].As one can understand from Eq. (4.5), this requirement puts a constraint on the inflaton couplings (through Γ e φ ) and potential parameters Λ and n.Here we would like to mention that for w 0.65 (equivalently n 5), purely gravitational reheating (PGR) becomes important [40,[74][75][76] and may dominate over perturbative reheating as shown in Ref. [74,76].Namely, for w 0.65, PGR alone cannot successfully reheat the Universe because i) the energy density of the SM radiation redshifts faster than the inflaton energy ( w ≤ 0), and ii) the reheating temperature is below the BBN bound ( w 0.65).Since our goal is to predict constraints on perturbative inflaton-matter coupling, hereafter, we consider n ∈ [2 , 5], where the gravitational reheating could be neglected.For completeness, we would like to mention that there exist several other kinds of reheating mechanisms.For example, in minimal preheating scenario, discussed in [77], the oscillating inflaton condensate directly transfers its energy to the SM Higgs field, even in the absence of a direct inflaton-Higgs coupling.Such a process, which relies typically on the tachyonic instablities [78], turns out to be very efficient compared to the perturbative reheating scenario.Also, preheating arising from tri-linear vertex can feature tachyonic resonance [79], that can lead to nonperturbative copious particle production immediately after inflation.As it has recently been shown in [80], the presence of inflaton self-interactions can lead to inflaton fragmentation that can severely influence the reheating dynamics.A complete analysis including all such collective effects requires a dedicated lattice simulation [80][81][82], which we plan to address in a future work.

Inflation: Λ
CMB measurement of the inflationary observables r (the tensor-to-scalar ratio) and n s (the spectral index) puts an upper limit on the Hubble scale at the end of inflation, that can be recast into a limit on the scale of inflation Λ so that Λ 10 16 GeV, see Appendix.B for details.In addition, adopting the end-of-inflation condition V (φ e ) 1, one can determine the inflaton field strength at the end of inflation φ e .

Scenario I: φ → SS
Now, we are ready to discuss various reheating models.We begin with the tri-linear scalar interaction involving the inflaton and a pair of real scalars S. Adopting Eq. (A.16), the following result for φ → SS decay width has been obtained with β S = 3(1−n) n+1 .Also note that β S = (n + 4)/(n + 1) for n = 1/2, so it happens to be outside of our region of variation for 2 < n < 5. Unlike the model-independent analysis CMB disallowed Γ ϕ e /GeV=1.7×10 -2   Γ ϕ e =1.7×10 -35   Γ ϕ e =1.7×10 -45   Γ ϕ e =1.7×10 -51 Top Left: Reheating temperature T rh as a function of n.The red part of the curve is discarded from present bound on ∆N eff due to PLANCK [61].The grey-shaded region is in conflict with the BBN limit on T rh .Top Right: Ω GW as a function of frequency, where different curves correspond to different choices of Λ for a fixed g Sφ and n.The orange-shaded regions are discarded from ∆N eff due to overproduction of GWs.Bottom Left: Same as top right, but for fixed Λ and n, with different choices of g Sφ specified below the curves.Bottom Right: Same as bottom left but for a fixed Λ and g Sφ , different curves correspond to different n values.In the bottom and upper-right panels the parameters Λ, n and g Sφ are chosen so that the BBN limit T rh > 4 MeV is satisfied [cf.Fig. 5].Moreover, in the bottom panels, we show sensitivities of future GWs experiments.before, here, for each n, the exponent β S of the decay width is fixed by construction.Now, assuming reheating takes place through the decay of inflaton into scalars, we can compute the temperature T rh at the end of reheating from Eq. (4.5).This is, of course, a function of the decay width itself, which in turn makes it dependent on the inflatonscalar coupling g Sφ .Note that the same coupling strength also controls GW spectral energy density [cf.Eq. (3.32)].The question that one can therefore ask is, for what size of g Sφ is it possible to have GW energy density that may have some detection prospects satisfying ∆N eff constraint, together with successful reheating, keeping the BBN predictions unharmed.
In the top left panel of Fig. 2, we choose different coupling strengths and compute corresponding T rh as a function of n for fixed inflationary scale Λ = 10 15 GeV.From this figure, one can see that for inflaton-scalar coupling g Sφ 10 −12 GeV, the reheating temperature exceeds the BBN bound T rh 4 MeV for all considered values of n.Moreover, the strongest the g Sφ coupling, the highest T rh .This can be understood from Eq. (4.5) noting that T rh ∝ [(g Sφ ) n /Λ] 1/(2 n−1) , hence a large coupling leads to higher reheating temperature for n > 2. Following Eq. (3.32), the GW spectral energy density goes as . Therefore, for n > 7/8 and given g Sφ , a large scale of inflation Λ implies the overproduction of GWs.This is evident from the top right panel of Fig. 2, where we see that a large Λ can easily be in conflict with ∆N eff constraint (albeit there exists an upper bound from CMB).Also, in the top left panel, the red curves correspond to the spectral energy density of GW that is ruled out from the present bound on ∆N eff due to PLANCK [61].As expected, this bound is stronger for smaller g Sφ , following the argument adopted above.The same feature is also visible from the bottom left panel.Finally, in the bottom right panel, we show the shape of Ω (0) GW (f ) for different choices of n, where we see the spectral energy density corresponding to n = 3 , 4 lies within reach of the proposed GWs detectors while remaining beyond the reach of the future ∆N eff constraints.We show existing and expected sensitivity reaches 9 from LIGO [84][85][86][87][88][89], LISA [90,91], CE [92,93], ET [94,95], BBO [96,97], DECIGO [98][99][100][101][102], µ-ARES [103], THEIA [104], AEDGE [105,106] and AION [106][107][108][109] in the bottom panel, and in the subsequent figures in this section.

Scenario II: φ → ψψ
For the fermionic final state, the inflaton decay width at the end of inflation is given by following Eq.(A.17), with β ψ = 3 (n − 1)/(n + 1).Note that, in this case, β ψ ≷ (n + 4)/(n + 1) for n ≷ 7/2.We see the effect of this transition in the top left panel of Fig. 3, where the evolution of the reheating temperature has different behaviour around n = 7/2.This can be understood from Eq. (4.5), where, for a fixed Λ and g ψφ , one finds that This shows, T rh decreases with n for 2 < n 7/2, while for 7/2 n < 5, an increase is observed.As before, we see from the bottom left panel, for fixed n = 3, small couplings  are disallowed from the present ∆N eff constraint on Ω (0) GW .This is understandable from Eq. (3.32), as for n = 3, the GW spectral energy density behaves as Ω (0) GW ∝ Λ 6 /g ψφ 3/5 , implying, a too small coupling may result in Ω (0) GW overproduction for a given Λ.On the other hand, it also shows that large Λ may overproduce GW energy density for a fixed n and coupling strength, as seen in the top right panel.In the bottom right panel, we show different values of n that are within reach of several future experiments for g ψφ = 10 −5 and Λ = 10 16 GeV.

Scenario III: φ → aa
The derivative interaction between inflaton and a pair of scalar bosons give rise to a decay width of with β a = 9(n−1) 1+n .In this case β a = (n + 4)/(n + 1) for n = 13/8, so it is beyond our range of variation for n.Note that the h aφ coupling has an inverse mass dimension.We find that Γ ϕ e /GeV=1.5×10 10 / G e V -1 for a fixed Λ = 10 16 GeV, couplings weaker than h aφ 10 −20 GeV −1 are entirely forbidden from present ∆N eff constraint on Ω (0) GW .This is shown via the red curve in the top left panel of Fig. 4. One can understand this behaviour from the fact that for n > 2, smaller couplings lead to GW overproduction, as Ω (0) GW ∝ h show the dependence of Ω (0) GW (f ) on the scale of inflation Λ for a given h aφ in the top right panel.It is interesting to note here that the cut-off frequency f max is independent of the choice of the scale of inflation Λ.As mentioned, for n > 2, one cannot go to a very small coupling since it results in GW overproduction and contradicts ∆N eff bounds.This can be realized from the bottom left panel.However, for h aφ = 10 −17 GeV −1 the GW amplitude falls within reach of the future detector sensitivities, still satisfying ∆N eff constraint, as shown in the bottom right panel.

Combined limits
In Fig. 5, we present the parameter space that remains allowed after imposing the limits discussed in the above subsections.The red shaded regions are disallowed from BBN bound on reheating temperature, that requires T rh > 4 MeV [cf.subsection.4.1.2],for different choices of the couplings.The "CMB bound", shown via the gray shaded region parallel to the horizontal axes, implies bound on the scale of inflation Λ 10 16 GeV, from CMB observables [cf.subsection.4.1.3].The cyan shaded regions are forbidden from ∆N eff bound from PLANCK on overproduction of GW energy density [cf.subsection.4.1.1].Note that, in case of derivative interaction, the cyan region corresponds only to h aφ = 10 −20 GeV −1 , as for other values of h aφ this bound is absent.Before moving on we would like to clarify that limits on tri-linear inflaton couplings from CMB observables have been derived in, for example, Ref. [110,111], considering α-attractor potential for the inflaton and taking non-perturbative effects into account.These analyses assume n = 1 reheating scenario.We, on the other hand, are typically interested in n > 2 in order to include bounds from ∆N eff , while considering only perturbative reheating.

UV freeze-in with time-dependent inflaton decay
In this section, we discuss Dark Matter (DM) production via UV freeze-in in a modelagnostic way, i.e., adopting the time-dependent inflaton decay rate as shown in Eq. (2.7).The evolution of the DM number density n DM is governed by the Boltzmann equation (BEQ), which can be written in a generalized form as where γ(T ) corresponds to the DM production rate out of SM particles as a function of the bath temperature T .Note that the equation (5.1) is coupled through H to the Friedmann equation, so, in this case, we are solving the full set of equations together with Eqs.(2.5).
When the SM entropy is conserved, it is instructive to rewrite Eq. (5.1) in terms of the DM yield defined as a ratio of DM number to entropy density of the Universe, Y ≡ n DM /s, where s(T ) ≡ 2π 2 45 g s (T ) T 3 , with g s (T ) being the number of relativistic degrees of freedom contributing to the SM entropy.Then Eq. (5.1) reads10 Figure 5. φ → SS (top left), φ → ψψ (top right) and φ → aa (bottom) scenario.In all cases the red shaded region is disallowed from BBN bound on reheating temperature T rh < 4 MeV, the cyan dashed region is forbidden from PLANCK observed ∆N eff bound on GW overproduction at f = f max , and the "CMB bound" discards scale of inflation Λ > 1.34 × 10 16 GeV from constraint on tensor to scalar ratio.
For the UV freeze-in, we assume that DM communicates with the visible sector through non-renormalizable operators suppressed by an appropriate power of the scale of new physics (NP) Λ UV and that the maximal temperature of the bath is well below Λ UV , so T max Λ UV .If, in addition, we assume that the temperature of the thermal bath at the end of reheating is large enough to neglect radiation and DM masses (m SM , m DM T rh ), then the production rate of DM from the SM bath can be parametrized as [32,[113][114][115] where κ = 2 (d − 2) with d being mass dimension of the relevant effective operator, (d ≥ 5).Note that this scale of new physics is uncorrelated with the scale of inflation Λ and can be larger, smaller, or equal to Λ. Finally, to match the whole observed abundance, the present DM yield Y 0 has to be fixed so that GeV, where m DM is the DM mass in GeV, ρ c 1.1 × 10 −5 h 2 GeV/cm 3 is the present critical energy density, s 0 2.9 × 10 3 cm −3 is the entropy density at present, and Ω DM h 2 0.12 [61].

Case-I: T rh m DM
Utilizing the solutions for T (a) and H(a), obtained in the previous section, one finds that the DM comoving number N ≡ n DM a 3 during the non-standard reheating evolves as which, after on integration between a e and a rh with the standard freeze-in initial condition N (a e ) = 0, leads to DM number at the end of reheating where and is the exponential integral function defined as whereas Γ[a, z] is the incomplete gamma function.Note that the above result (5.5) is valid for light DM, i.e., with mass lower than the reheating temperature.Moreover, the UV freeze-in dark species are mainly produced just after the cosmic inflation, when the thermal bath temperature is the highest.Hence, the contribution from the late times, i.e., a > a rh , is negligible.Although the SM entropy density is not conserved when the inflaton is decaying, one can calculate the DM yield at a rh as Y (a rh ) = n DM (a rh )/s(T rh ) = N (a rh )/ a 3 rh s(T rh ) .After the end of reheating, the SM entropy is conserved and therefore Y (a rh ) remains constant.Thus we obtain where a e /a rh can be read from Eq. (2.17), whereas T rh is given by Eq. (4.5).To obtain the right abundance for DM, m DM Y (a rh ) needs to match the present DM abundance, which leads to which shows that for a given decay channel, and n, Λ the ratio Λ κ−4 UV /m DM is fixed.

UV
(5.12) While the DM yield is conserved after reheating and not during reheating due to the entropy injection, the final yield at the end of reheating reads where a DM /a rh is given by Eq. (5.11).As before requiring the right abundance for the DM, one can show, for a given channel, m DM /Λ κ−4 UV is a constant for a given inflaton decay channel for a fixed n , Λ [cf.Eq. (5.10)].Note that, in determining the UV freeze-in yield, we introduce three more free parameters here, namely, the DM mass m DM , the scale of effective interaction Λ UV , and the mass dimension of the DM-SM operator, κ.

Dark matter production
In Fig. 6, we show the contours that reproduce the observed DM relic abundance.We consider different reheating scenarios, e.g., φ → SS , ψψ , aa and project the parameter space in Λ UV − m DM plane by fixing the scale of inflation Λ = 10 16 GeV, n = 5 and κ = 6 (or equivalently, d = 5).We choose the strength of the inflaton-matter coupling in such a way that they give rise to observable GW spectrum as discussed in Sec. 4. From the top left panel, we first note that in all cases, the parameter space shows the existence of a cut-off for certain m DM , beyond which DM can no longer be produced.This cut-off depends on the two cases discussed before.In case where m DM T rh , DM species heavier than T rh cannot be created from the thermal bath due to the Boltzmann suppression.On the other hand, for dark particles with mass T rh m DM T max , the production ceases beyond m DM = T max .In both cases, UV freeze-in production becomes inefficient when the SM temperature drops below dark matter mass.The slope of each contour is exactly given by 1/(κ − 4), which can easily be obtained from Eq. (5.10) once all parameters are fixed except for m DM and Λ UV .Next, we note that in all cases, a larger inflaton coupling g iφ requires larger Λ UV to produce the correct DM abundance for a given DM mass (keeping all other parameters fixed).This can be understood from Eq. (5.9), where we find, for n = 5 , κ = 6 , Λ = 10 16 GeV, the DM yield approximately reads Y (a rh ) ∝ g

5/2
Sφ /Λ 2 UV in case of inflaton decaying into scalars, while for fermionic decay Y (a rh ) ∝ g 3/2 ψφ /Λ 2 UV and for the derivative interaction Y (a rh ) ∝ h 3/2 aφ /Λ 2 UV (similar dependence can also be found from Eq. (5.12), for a fixed DM mass).Therefore, irrespective of the inflaton decay products, a larger coupling leads to a larger abundance that can be compensated by a higher Λ UV for a given DM mass.Moreover, for our choice of parameters, we find that the decay width into pairs of scalars is small compared to the fermionic or derivative scenario.Since the DM yield in all cases is proportional to (Γ e φ ) x /Λ κ−4 UV (x is some function of n, the exact form of which depends on the decay mode), hence we see the scale Λ UV is more suppressed for the scalar case compared to other two cases.Finally, in the bottom right panel, we show the effect of DM-SM operators of different dimensions on the DM parameter space.As an illustrative example, we choose the fermionic reheating scenario and fix all the relevant parameters.Now, we see, a larger d (≡ κ) relaxes the bound on Λ UV for a given DM mass since Y (a rh ) ∝ m DM Λ 8−2 d UV .Hence, operators of larger dimensions lead to more suppression, allowing lower Λ UV .Since the same inflaton-matter coupling controls DM abundance and GW spectrum, it is, therefore, possible to make a one-to-one correspondence between DM abundance and observable PGW spectrum through such interactions.This is what is portrayed in Fig. 7.In the left panel, we show the parameter space that produces right relic abundance considering all possible inflaton-matter interactions and choosing {g Sφ , g ψφ , h aφ } = {10 −15 GeV , 10 −7 , 10 −18 GeV −1 }, for a fixed Λ UV , n, and considering DM-SM operator of mass dimension d = 6.For the same choice of parameters, we show the GW spectrum in the right panel, where each curve corresponds to a range of m DM that gives rise to the right abundance.This shows, for the fermionic reheating scenario, that it is possible to probe the relic density allowed parameter space with GW detectors like BBO or DECIGO, for m DM {10 −5 − 10 11 } GeV.This is also true for the bosonic (scalar) reheating scenario, where the DM mass lies in the range m DM {10 −5 − 10 0 } GeV.Although derivative interactions can provide the right relic abundance, the prediction remains beyond the reach of future GW detectors, as shown by the black dashed curve.However, for this set of parameters, φ → SS coupling could potentially be ruled out from future constraints on ∆N eff from CMB.It is understandable that by choosing the inflaton-matter coupling and the scale of inflation appropriately, it is always possible to have a detectable GWs spectrum, satisfying BBN constraint, while accommodating the proper abundance for the DM requires tuning Λ UV and m DM .Thus, a future GW detection can be regarded as a potential freeze-in DM signal, originating from an underlying new physics leading to a DM-SM operator of a given dimension.

Conclusions
This work investigates the possibility of probing inflaton couplings to the visible sector by primordial gravitational waves (GWs) of inflationary origin.We consider the inflaton field oscillating with a time-dependent amplitude assuming the α-attractor T-model potential.Such a scenario inevitably implies a time-dependent decay width of the inflaton for potential different from quadratic (for small field strength).The inflaton is assumed to couple to a pair of bosonic or fermionic matter within a perturbative regime that can give rise to successful reheating.We then compute the primordial power spectrum and the spectral energy density of primordial GW.The energy density of the Universe turns out to be a function of the inflaton decay width at the end of inflation, which, in turn, is related to the inflaton couplings to matter.Satisfying bounds from the requirement of successful BBN and CMB, we find the inflaton coupling as small as g Sφ ∼ 10 −15 GeV becomes sensitive to future gravitational wave detectors in case the inflaton decays into a pair of scalar bosons via an operator of dimension three.For inflaton, Yukawa interactions with fermions, the corresponding coupling that falls within the detector sensitivity turns out to be g ψφ ∼ 10 −5 .For derivative interaction involving a mass dimension five operator, coupling strength h aφ ∼ 10 −17 GeV −1 becomes sensitive to the experiments.It is important to note that our result strongly depends on other free parameters of the model, for example, the shape of the potential or scale of inflation which is currently bounded from PLANCK data.However we envisage that once we got more information on inflation from next generation CMB experiments, like measuring the scale of inflation, uncertainties of our prediction would be reduced.
We finally discuss the production of dark matter (DM) during reheating from the radiation bath, assuming the DM communicates with the SM via non-renormalizable operators suppressed by large-scale Λ UV .In this scenario, the UV freeze-in from the thermal bath is the primary mechanism of DM production.The freeze-in scenario requires tiny couplings between the bath and DM particles.Therefore it is very challenging if not impossible to probe such very feebly coupled dark sector.Here we promote GW as a novel observable to probe such dark sectors and identify the parameter space that satisfies the observed relic density.Because of the time-dependent inflaton decay width, the standard DM yield also gets modified in this framework.It turns out to be possible to tune the scale of new physics Λ UV and the DM mass to satisfy the PLANCK observed relic DM density.The interesting point here is that depending on the nature of the inflaton decay channel, one can produce DM over a large range of mass.As the same inflaton decay width influences DM production and controls the GW spectrum, hence a detectable GW may be interpreted as an indication of the DM detection and provide information about the scale of new physics.
Finally we comment on the fact that in our analysis we have only considered one GW detector running at a time.However once more than one GW detector starts operating then sensitivities to GW spectral shapes are expected to improve so that much larger regions of the parameter space could be probed.Such an analysis is beyond the scope of the current manuscript but we envision precision measurements that GW astronomy promises due to the planned global network of GW detectors, which can fulfill the dream of testing high-scale and weakly-coupled fundamental BSM scenarios a reality in very near future.the form (A.12) The total probability can then be found by summing over each outgoing momenta.In the continuum limit, this reduces to multiplying P (p 1 , p 1 )/T by a factor V d 3 p 1 /(2π) 3 V d 3 p 2 /(2π) 3 .The energy gain from created particles in volume V and time dt is

B Constraints on reheating from inflation
Here we derive constraints of the α-attractor T-model from the combined WMAP, Planck, and BICEP/Keck data [118][119][120].Let us start with revisiting some of the necessary inflationary parameters.The so-called potential slow-roll parameters are defined as Note that at the end of inflation ä = 0, which roughly corresponds to V (φ e ) 1.This condition allows us to find the inflaton field value at the end of inflation Moreover, the inflaton potential at a = a can be expressed in terms of the scalar powerspectrum and tensor-to-scalar ratio as Since, at a = a the field value φ > M P , the inflaton potential can be approximated by a constant value V (φ ) ≈ Λ 4 , which in turn sets an upper bound on the inflationary scale, Λ ≤ 1.34 × 10 16 GeV, which has a feeble dependence on n, α.

Figure 1 .
Figure1.Evolution of cosmological comoving length scales with the scale factor at different epochs of the evolution of the universe.Here "RD" stands for radiation domination, "MD" implies matter domination, and "DE" indicates dark energy domination.In each case corresponding equation of state is also mentioned.

5. 2 2 2 P κ 4 H e m 3 DM Λ κ− 4
Case-II: T rh m DM T max Production of heavy DM particles with mass exceeding the reheating temperature, is Boltzmann suppressed in the region where m DM > T rh .The present DM number density in that case could be estimated by integrating Eq. (5.4) from a e to a DM ≡ a(T = m DM ), where a DM = a rh Y (a DM ) = N (a DM )/ a 3 DM s(m DM ) then reads g ρ (n+1) Γ e φ H e M

Λ 7 Figure 6 .
Figure 6.Dark Matter parameter space producing right relic density in Λ UV − m DM plane, considering reheating via φ → SS (top left), φ → ψψ (top right) and φ → aa (bottom left), for a fixed Λ = 10 15 GeV, n = 10 and DM-SM operator of dimension d = 5.In the bottom right panel we compare parameter space for DM -SM operators of dimension d = {5 , 6 , 7}, considering fermionic reheating as a benchmark scenario.