Momentum distribution of dark matter produced in inflaton decay: effect of inflaton mediated scatterings

Post-inflationary reheating is a widely discussed mechanism for non-thermal production of dark matter (DM). In this scenario the momentum distribution of the produced DM particles is usually taken to be the one obtained at reheating, red-shifted at later times due to the expansion of the Universe. However, since in such a scenario both the DM and the standard model (SM) fields couple to the inflaton, the DM particles necessarily undergo self-scatterings, as well as elastic and inelastic scattering reactions with the SM bath, all of which proceed through $s-$channel or $t-$channel inflaton exchange. We compute the momentum distribution of the DM particles including the effect of these scatterings, and find that the distributions can be significantly altered, even though DM remains non-thermal throughout the cosmological evolution. We observe that if the inflaton dominantly couples to the SM Higgs boson through a renormalizable interaction, then reheating temperatures and inflaton masses at the TeV scale lead to a large effect from the scattering processes, with the DM-inflaton coupling constrained by the DM density. The scattering effects are found to be sensitive to the duration of the reheating process -- larger the duration, more momentum modes are filled at reheating, leading to an enhanced scattering probability. We also obtain the free-streaming length of such DM using the resulting non-thermal momentum distribution, which can be used to estimate the implications of the Lyman-$\alpha$ constraints on the DM mass. It is observed that in the scenarios considered, including the scattering effects can reduce the DM average velocity at matter-radiation equality, and its free-streaming length, by upto a factor of $40$, thereby making the constraints on light DM produced in inflaton decay significantly weaker.


Introduction
Inflation is a leading candidate for a theory of the pre-radiation dominated Universe, in which the horizon and flatness problems in big-bang cosmology can be addressed, and which also makes predictions about properties of fluctuations in the cosmic microwave background and large-scale structure [1][2][3][4][5][6][7][8].
At the end of the slow-roll phase of inflation, the inflaton field φ undergoes a damped oscillation around the minimum of the potential φ 0 . When φ is close to φ 0 , the potential can be approximated by a quadratic form: This is equivalent to the field theory of spin-0 particles of mass m φ , and negligible velocity. Furthermore, since inflation needs to end eventually, giving rise to a radiation dominated Universe, the inflaton field must also couple and decay to other fields including standard model (SM) matter and radiation. This epoch of inflaton oscillation and decay is known as the period of reheating [1,6,[9][10][11][12][13][14][15][16]. If the dark matter (DM) field couples to the inflaton as well, DM may be produced non-thermally during the reheating epoch. Such a possibility has been extensively considered in the literature, see for example, Refs [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. Light, cold DM may also be produced during reheating, and with suitable choices of the reheating temperature T R and the inflaton mass m φ , the DM particles from inflaton decay can become cold enough due to the red-shift of the momenta [30].
Since the DM field in this scenario couples to the inflaton, inflaton mediated DM self-scatterings necessarily take place. In addition, since the SM fields also couple to the inflaton, DM-SM scattering processes become relevant as well. These scattering reactions, if they have appreciable rates, can modify the momentum distribution of the DM particles. Furthermore, once the inflaton mass and its couplings to the SM and DM fields are fixed to reproduce a given reheating temperature and a desired DM density, the DM-DM and DM-SM scattering rates are completely determined. The impact of the scattering processes on the DM momentum distribution should therefore be studied in detail. This is the primary objective of the present study. For light DM, the momentum distribution is directly relevant in determining the constraints from structure formation. We have, to this end, made a simple estimate of the free-streaming length of the DM particles, which can be used to bound the DM mass using constraints from the Lyman-α forest data [34][35][36]. This work, therefore, should be relevant for more accurate determinations of the structure formation constraints on light non-thermal DM produced in inflaton decay.
The rest of the paper is organized as follows. In Sec. 2 we discuss the production of DM and SM radiation during reheating, and obtain the DM number density and momentum distribution at the end of the reheating epoch, including the effects of quantum statistics and back-reaction. We then describe the framework adopted to obtain our primary new results on the effect of inflaton-mediated scatterings on the DM momentum distribution in Sec. 3. The computations of the DM average velocity and its co-moving free-streaming length are discussed in Sec. 4. We present our results on the modifications to the DM momentum distribution due to scatterings, and the impact of the reheating temperature, inflaton mass, the DM-inflaton coupling and the duration of the reheating process on it in Sec. 5. We also discuss the impact of the scatterings on the resulting DM average velocity at matter-radiation equality and the DM free-streaming length, along with the implications of the Lyman-α constraints in this section. Sec. 6 summarizes our results.

Dark matter production during reheating
We consider a minimal scenario in which the scalar inflaton field φ couples through renormalizable interactions to both the SM and the DM sectors. We take the DM field ψ to be an SM singlet fermion in the following discussion. The results can be straightforwardly extended to the case of a singlet scalar DM as well, with the essential modification being in the quantum statistics involved. The relevant interaction Lagrangian is given by Here, H is the SM Higgs doublet, the only SM field to which the inflaton has a renormalizable coupling 1 . Although here we have considered the reheating of the SM sector 1 The scenario considered here is a minimal one with the inflaton coupling to only the Higgs field through a renormalizable interaction. In general, the inflaton can couple to all the SM fields, and the broad qualitative conclusions in the subsequent sections of this study should continue to hold in such scenarios as well. Similar considerations apply to the R 2 (Starobinsky) inflation scenario, in which the reheating in the SM sector proceeds through the particle production by the curvature scalar, which is universally coupled to all types of elementary particles.
through the trilinear scalar interaction, the quartic interaction term φ 2 H † H can also be relevant [37].

DM number density and distribution in the instantaneous decay approximation
To understand the parametric dependence of the DM number density and its momentum distribution on the reheating process variables, we can start with the approximation that the inflaton field decays instantaneously, giving rise to a radiation dominated Universe. In that case, the DM number density at the reheating time t R , n ψ (t R ), can be computed as follows. Since the DM is produced from the decay of the inflaton quanta, φ → ψψ, we have where, we have added up the contributions to the DM particle and anti-particle number densities in n ψ (t R ), and BR(φ → ψ) denotes the branching ratio of φ decays to a DM pair.
We also have the inflaton number density is the inflaton energy density at t R . Let us also denote by Γ φ , Γ R and Γ ψ the total decay width of the inflaton, and its partial decay widths to the SM radiation and DM, respectively. In terms of the couplings in Eq. 2.1, we have Γ R = µ 2 H /(32πm φ ) and Γ ψ = λ 2 m φ /(8π). Assuming that the inflaton dominantly decays to the SM sector (which is necessary in this framework to have a radiation-dominated Universe at the big-bang nucleosynthesis (BBN) epoch), we can approximately write, is the radiation energy density at t R . Furthermore, the reheating time t R can be defined by the following relation Further assuming that the SM radiation bath thermalizes instantaneously with a temperature T R at time t R , we have where, g * (T R ) is the number of relativistic degrees of freedom at the temperature T R . Thus the DM number density at t = t R is proportional to the combination T 2 R /m φ , in addition to Γ ψ . If no DM particle is produced or annihilated after t = t R , n ψ (t R ) can be used to obtain the DM density Ω ψ = ρ ψ (t 0 )/ρ c at the present epoch t = t 0 , where ρ c is the critical density: .
Here, s 0 is the entropy density at the present epoch. Thus, keeping the partial decay width Γ ψ fixed, Ω ψ is decreased for higher values of T R m φ . Having obtained the DM number density n ψ (t R ), we can also determine its momentum distribution f ψ ( k, t R ) using the same instantaneous decay approximation, where the distribution function is normalized as g ψ being the internal degrees of freedom of ψ. At t = t R , all the DM particles are produced with momentum k = p ψ m φ /2, for m φ >> m ψ . Hence, where, A is a proportionality constant and by isotropy of the Universe f ψ is a function of | k| ≡ k only. The constant A can be determined using Eq. 2.7 with the n ψ (t R ) obtained in Eq. 2.4 leading to In the absence of any collisions, for t > t R , the momenta simply redshift as k(t) ∝ 1/a(t), where a(t) is the scale factor in the Friedmann-Robertson-Walker space-time metric.

Beyond the instantaneous decay approximation
To obtain a more accurate expression for the number density of the DM n ψ (t), we need to consider the evolution of the energy densities of the inflaton field (ρ φ ) and the SM radiation bath (ρ R ), and n ψ , which are governed by the following coupled set of approximate equations, along with the Friedmann equation for the Hubble expansion rate H In writing these equations, we have made the approximation that the back-reactions and the effect of quantum statistics are not significant. However, they will be important for the light DM to be considered here, and will be included in the subsequent discussion of the momentum distribution and number density of DM. We have also considered the DM to be a massive matter particle, whose number in a co-moving volume (n ψ a 3 ) is conserved in the absence of production and annihilation processes. The factor of 2 in the RHS of Eq. 2.12 accounts for the fact that 2 DM particles are produced in the decay of each inflaton.
Eq. 2.10 for ρ φ can be integrated to obtain where, the time t I is the beginning of the period of inflaton oscillations and decay, when the energy density in the inflaton field is ρ φ (t I ). We can also integrate Eq. 2.12 using the solution in Eq. 2.14 as an input, and obtain We can further rewrite Eq. 2.15 with the help of Eq. 2.14 to express n ψ (t R ) as follows: Comparing Eq. 2.16 with Eq. 2.2, we see that compared to the instantaneous decay approximation, there is an enhancement in n ψ (t R ), due to a larger inflaton number density n φ (t), for t I ≤ t ≤ t R . With Γ φ t R ∼ 1, and Γ φ t I → 0, this amounts to an enhancement by a factor of around (e − 1) 1.72. Finally, Eq. 2.11 gives the radiation energy density as In writing Eqs. 2.15 and 2.17, we have used the boundary conditions that n ψ (t I ) = 0 = ρ R (t I ). By solving the coupled system of equations 2.10-2.13 numerically, we can determine ρ R (t). Simple approximate solutions may also be obtained by assuming that the energy density of the Universe is dominated by a stable matter field φ during the time t I ≤ t < t R , where t R is the time when radiation domination starts. Using the fact that in a stable matter dominated Universe H = 2/3t, we obtain at t = t I and hence, where we have used a(t) ∝ t 2/3 in a stable matter dominated era. Using Eq. 2.19 in Eq. 2.11 we can solve for ρ R (t) in terms of an infinite series as [38] ρ Here again, we have imposed the boundary condition that ρ R (t I ) = 0.

Including the effects of back-reaction and quantum statistics
The computation so far has not included the effect of back-reaction and quantum statistics. As the DM number density builds up, both the effects may become relevant. In order to include these effects let us consider the Boltzmann kinetic equation for f ψ ( k, t), with f ψ ( k, t)d 3 k being the number density of DM particles in the momentum range k to k + d k: where, the collision term for the process φ( p) ↔ ψ( k) + ψ( k ) is given by In writing this equation we have assumed that the distribution functions for the DM particle and anti-particle are the same.
Since the inflaton field during reheating may be taken to be dominated by its zero mode, the decaying inflaton particles are nearly at rest. Thus, we can reduce Eq. 2.22 to the following form (2.23) In deriving Eq. 2.23, we have assumed that the phase-space density of the inflaton quanta is large, and hence 1 + f φ ( p, t) f φ ( p, t). Here, as before, p ψ is the momentum of the DM particle at production, and is given by For f ψ ( k, t) << 1, the effects of back-reaction and Fermi statistics are not important, and integrating Eq. 2.21 in that limit reproduces Eq. 2.15, as expected. The DM momentum k at time t is related to its production time t p , where t I ≤ t p ≤ t R , since at t p its momentum is always p ψ , which is then redshifted to k at time t. Hence, the k dependence of f ψ ( k, t) may be replaced by its t p dependence as follows: Here, we have used the relation which is obtained for a(t) ∝ t 2/3 valid in a stable matter dominated Universe. In writing the collision term in Eq. 2.25 from Eq. 2.23, we have also imposed the condition of isotropy of the Universe, which implies Integrating Eq. 2.25 and imposing the boundary condition that f ψ ( k, t I ) = 0, where t I is the time reheating begins, we obtain , Here, the Θ−function ensures the condition that k(t) < p ψ . It is useful to observe the difference between the distribution functions obtained (a) without back-reaction and Fermi statistics, (b) including back-reaction but with classical statistics, and (c) with both back-reaction and Fermi statistics. While scenario-(c) is already discussed above, the distribution in scenario-(a) is found to be and in scenario-(b) the distribution is obtained as follows The functionf (t p ) in Eq. 2.28 can be written in terms of the momentum where, we have set Γ φ t R ∼ 1. The distribution function in Eq. 2.27, along with Eq. 2.31 are used as the initial condition at t = t R for the evolution of f ψ ( k, t) due to the effect of collisions. We compare these three distribution functions and the resulting DM number densities in Fig. 1, for fixed values of the relevant parameters shown in the figure. From the left panel of Fig. 1, we see that scenario-(a) (red solid curve), for which a closed form analytic expression for n ψ (t R ) can be obtained, overpredicts the DM momentum distribution in all but the highest momentum values. In scenario-(b)(orange dashed curve), the back-reaction effect suppresses the distribution throughout − lower the momentum at reheating, higher the suppression. This is because the particles produced at the earliest times have the lowest momentum at reheating due to a large redshift effect. On the otherhand, it is these particles which have the highest probability of recombining to produce an inflaton through back-reactions, due to the longer timespan available. Including the effect of Fermi statistics in scenario-(c) (blue dotted curve) suppresses the distribution further due to Pauli blocking by a smaller ratio. For comparison, we also show the delta-function-like distribution in (d) the instantaneous decay approximation (green solid line). For our subsequent numerical studies, we have included the effects of both back-reaction and Fermi statistics.
Back-reaction,Fermi λ Ω ψ h 2 Back-reaction, Fermi In the right panel of Fig. 1 we have shown the DM relic abundance Ω ψ h 2 as a function of the DM-inflaton coupling in all the four cases (a)-(d). Following the distribution functions, it is straightforward to understand the Ω ψ h 2 curves. It is interesting to observe that Ω ψ h 2 reaches a constant value above a certain coupling of around λ 10 −8 , when back-reactions and Fermi statistics effects are included. This is essentially because for couplings larger than this value, the initial DM number density builds up quickly, such that the forward and backward reaction rates become the same, thereby freezing the change in the co-moving DM number density n ψ a 3 . We also observe that the instantaneous decay approximation differs from the scenario-(1) prediction by a constant factor (of about 1.72) as seen in the comparison between Eq. 2.16 with Eq. 2.2. We find that for T R = 1 TeV and m φ = 5 TeV, a 10 keV DM particle achieves the relic density of Ω ψ h 2 = 0.12 for a DM-inflaton coupling of around λ 10 −9 .

Dark matter momentum distribution: effect of inflaton mediated scatterings
The collision term f coll ψ ( k, t) in the kinetic equation 2.21 for f ψ ( k, t) receives further contributions from inflaton-mediated scattering processes, which are necessarily present in the reheating scenario studied here. In particular, the interaction Lagrangian in Eq. 2.1 implies the following reactions: 1. s−channel and t−channel inflaton mediated scatterings: 3. s−channel inflaton mediated scattering: Here, we have assumed the reheating temperature T R > m h , so that the Higgs boson h is present in the thermal bath. For T R < m h , the same interaction Lagrangian Eq. 2.1 would imply DM scatterings with W −bosons, or b−quarks or SM leptons, depending upon the relevant temperature.
In presence of the above 2 → 2 scatterings, the DM momentum distribution evolves according to the following Boltzmann equation where we have written the Liouville operator in the LHS in polar co-ordinates for the momentum p 1 , by imposing the condition of isotropy to set the derivatives with respect to the polar and azimuthal angles to zero. The collision term C(E 1 , t) encapsulates the effects of the 2 → 2 processes given above and can be written as follows: which can be conveniently broken up into two terms,

(3.4)
In contrast to the collision term for the φ ↔ ψψ process, in which the effect of Pauli blocking was included while determining the DM momentum distribution at reheating in Sec. 2, for the 2 → 2 processes considered above, we have made the approximation of omitting the effects of quantum statistics. This is a small effect here, and its inclusion makes the numerical evaluation of the collision integral more computationally intensive.
(3.5) With this, we can reduce the backward collision term to the following form (3.6) where, the angular integral is given by a cos 2 α + b cos α + c Θ(a cos 2 α + b cos α + c) . (3.7) In Eqn. 3.7 the parameters a, b and c are defined as follows: with E i = p 2 i + m 2 i . The Heaviside theta functions in Eqs. 3.6 and 3.7 put the relevant kinematic constraints on the phase-space variables.
(3.11) Using these parametrizations, we obtain the following form of the forward collision term (3.12) with the angular integral given by with E i = p 2 i + m 2 i . The matrix elements squared appearing in Eqs. 3.7 and 3.13 are obtained to be the following 1. ψψ → ψψ: 2. ψh → ψh or ψh → ψh: 3. ψψ → hh: where s = (p 1 + p 2 ) 2 = (p 3 + p 4 ) 2 and t = (p 3 − p 1 ) 2 represent the standard Mandelstam variables. We numerically solve the Boltzmann Eq. 3.1 with the collision term given in Eq. 3.2. The angular integrals in Eqs. 3.7 and 3.13 are evaluated using Monte-Carlo integration implemented in CUBA [43] to calculate the functions F (p 1 , p 3 , p 4 ) and F (p 1 , p 2 , p 3 ), respectively. These functions are then fed into Eqs. 3.6 and 3.12 to obtain the full collision integral, which are subsequently used to solve Eq. 3.1.
We have followed the evolution of the DM momentum distribution from reheating (t R ) to matter-radiation equality (t EQ ). The boundary condition is set by the distribution function at t = t R which is shown in Eq. 2.27. For numerical convenience we have traded the time variable in terms of the temperature of the SM plasma using the following relation: The Boltzmann equation is then solved using the backward difference formula in 200 bins spanning the momentum range 10 −10 GeV to 10 10 GeV. To solve the Boltzmann equation we have used a suitable differential equation solver from the SciPy library of Python [41,42]. We found that for all the cases considered in this work, the shape of the distribution freezes around T ∼ 100 GeV, and then subsequently redshifts till matter-radiation equality, implying that the scattering processes are no longer active below this bath temperature.

Dark matter average velocity and free-streaming length
Having obtained the momentum distribution of the DM particles, including the effects of inflaton-mediated scattering processes, we can easily estimate the average DM velocity at different epochs, and from there the DM free-streaming length λ FSH . We recall that until the DM perturbations become Jeans unstable and begin to grow at t = t EQ , t EQ being the time of matter-radiation equality, collisionless particles can stream out of overdense regions into underdense regions. This process can smooth out inhomogeneities, and would be constrained by the DM power spectrum deduced from cosmological observations [34]. This effect can also be approximately understood using the DM co-moving free-streaming length. The co-moving free-streaming length is defined as the co-moving distance travelled by the DM particles between the time of decoupling of scattering reactions t dec and the time of matter-radiation equality t EQ [34,35]: where, the average velocity at the time t, v(t) is given by Eq. 4.1 can be re-written by exchanging time t with the scale factor a(t) as If the DM decoupling takes place in a radiation-dominated Universe, then, we have the following relations between the scale factor and temperature: a dec = T 0 T dec g * ,s (T 0 ) g * ,s (T ) and a EQ = T 0 T EQ , as the relativistic degrees of freedom g * ,s at equality and the present epoch are the same. Here T 0 denotes the present temperature, with a(T 0 ) = 1. Including the contributions from matter and radiation energy densities, the Hubble parameter may be expressed as where, H 0 is the Hubble expansion parameter at the present epoch, and Ω R = ρ R (T 0 )/ρ c , ρ c being the critical density. Using Eq. 4.4, the free-streaming length can be expressed as

3/2 EQ
√ a + a EQ da. where a EQ = Ω R /Ω M = 2.91 × 10 −4 . It is convenient to re-write the momentum integral for v(a) in terms of the scaled

Results
The 2 → 2 scattering cross-sections discussed in the previous section are functions of the DM-inflaton coupling λ, the SM-inflaton coupling µ H , the inflaton mass m φ and the DM mass m ψ . Since the inflaton decay to the SM sector dominates its total decay width, we can trade the coupling µ H with the reheating temperature T R . Thus, for a fixed m ψ and λ, the scattering rates are controlled by T R and m φ . T R also controls the branching fraction of inflaton decay to DM, and hence its number density at reheating. In addition to these effects, T R and m φ affect the impact of scatterings on the DM momentum distribution in a number of ways, as discussed subsequently. Since the collision integrals and the solution of the Boltzmann equations for the DM momentum distribution are numerically highly demanding even for a single set of parameter values, a full scan of the {T R , m φ , m ψ , λ} parameter space is beyond the scope of this work. We shall therefore choose a particular set of these parameters to illustrate the effect of collisions. The effect of the scatterings is also controlled by how broad the DM momentum distribution is at reheating. This in turn is dependent on the duration of the reheating process. We quantify the duration of the reheating process by the scale factor ratio a(t I )/a(t R ), where t I denotes the beginning of the reheating phase, and t R denotes the beginning of radiation domination with Γ φ t R ∼ 1. Smaller values of a(t I )/a(t R ) would correspond to a longer duration of the reheating process. In Fig. 2, we show the momentum distribution of DM including the effect of the three 2 → 2 inflaton-mediated scattering processes discussed in Sec. 3 (red solid histogram). We have shown the distribution at the epoch with SM bath temperature T = 1 MeV. The corresponding distribution at the matter-radiation equality, T EQ ∼ 1 eV remains the same in shape, with each momentum getting red-shifted by the same factor. For comparison, we also show the distribution at reheating with T R = 1 TeV (green dot-dashed curve) obtained earlier, and the corresponding red-shifted distribution without the effect of collisions at T = 1 MeV (green dotted curve). As we can see by comparing the distributions at T = 1 MeV with and without collisions, the two-body scatterings lead to a significant modification over the whole momentum range. However, the effective change is more pronounced in the lower and medium momentum range. The modifications are driven by two factors − one being the energy exchange through the DM self-scatterings and the t−channel scattering of DM particles with the SM bath, and the other being the non-thermal freeze-in production of DM through the inflaton-mediated s-channel h(p 3 )h(p 4 ) → ψ(p 1 )ψ(p 2 ) process. For the DM coupling and inflaton mass chosen, the collisions do not thermalize the DM with the SM sector. In order to demonstrate this, we also show in Fig. 2 the thermal Fermi-Dirac distribution at T = 1 MeV for a particle of mass 10 keV with zero chemical potential (black dot-dashed curve). Clearly, in addition to the small difference in the normalization, the shape of the distributions, especially at lower momenta, differ significantly. All the results are shown with the inflaton mass m φ = 5 TeV, reheating temperature T R = 1 TeV, DM mass m ψ = 10 keV, DM-inflaton coupling λ = 6.6 × 10 −9 and the scale factor ratio chosen as a(t I )/a(t R ) = 5 × 10 −5 .  Figure 3: The average DM velocity v as a function of the scale factor ratio a/a EQ , with the effect of the inflaton-mediated 2 → 2 scatterings included (solid curve) and without them (dot-dashed curve). All the relevant parameters are kept fixed at the values shown in the figure, which are the same as in Fig. 2. The value of the DM free-streaming length λ FSH in both cases are also shown, which, along with v is significantly modified by the effect of collisions.
The impact of the scattering processes can be quantified using the observationally relevant quantities such as the average DM velocity discussed in Sec. 4. To understand this effect, we show in Fig. 3, v as a function of the scale factor ratio a/a EQ , with the effect of the inflaton-mediated 2 → 2 scatterings included (solid curve) and without them (dotdashed curve). All the relevant parameters are kept fixed at the values shown in Fig. 3, which are the same as in Fig. 2. There are two different regimes in both the curves. The first is the relativistic regime, in which v 1. For higher values of the scale factor we have the non-relativistic regime in which we observe a constant slope for the curve without collisions, as in this case, v ∝ | p| ∝ 1/a. The onset of the non-relativistic regime is found to be much earlier for the case with scatterings included. This is because, the scattering processes lead to the enhanced occupation of lower momentum modes, which in turn reduces the average velocity. As we can see from this figure, the difference between the two cases is highly significant, with the v at matter-radiation equality, where a/a EQ = 1, differing by a factor of 40. This effect can also be quantified by the free-streaming length of DM defined in Sec. 4, which differs in the two cases by a factor of around 37. Thus for the DM mass shown, m ψ = 10 keV, while the value of λ FSH without including the collisions is found to be 0.0734 Mpc, inclusion of the collision effects reduces it to 0.0020 Mpc. Therefore within this estimate, while the number without the effect of collision would imply strong constraints on such a DM from the Lyman−α forest data [36], including the collision effects shows that for fixed values of the other parameters, this DM mass is highly consistent with the structure formation requirements. This crucial effect of the scattering processes is the primary result of this paper. We observe that for fixed values of the couplings, the scattering effects are sensitive to the duration of the reheating process. This is because of two reasons. Firstly, as seen earlier, for example in Eq. 2.16, the duration of the reheating process affects the number density of DM produced at reheating, which in turn enters the reaction rates. Secondly, longer the duration, the momentum distribution at reheating is broader and more stretched towards lower momentum values. A larger range of DM momenta increases the probability of the scatterings further. Even if the number density remains very similar, the broadening of the distribution itself can make a large difference in the effect of the collisions. In order to demonstrate this effect, we show in Fig. 4 the DM momentum distributions with a(t I )/a(t R ) = 5 × 10 −3 . For this figure, all other parameters have been kept fixed at the same values as in Fig. 2, for which a(t I )/a(t R ) was taken to be 5 × 10 −5 . For the two a(t I )/a(t R ) values shown here, the DM number density at T R are nearly the same. However, by comparing Fig. 4 with Fig. 2, we see that larger the duration of reheating, the initial distribution at T R is broader, which is easily understood by taking into account the larger redshift of the DM momenta the earlier it is produced. By further comparison of Fig. 4 with Fig. 2, we find that the effective impact of the collisions on the distribution is larger for a smaller a(t I )/a(t R ).

Without Collision
With Collision a(tI )/a(tR) v(aEQ) λFSH (Mpc)   For a more detailed numerical comparison, we show the average velocity at matter radiation equality ( v(a EQ ) ), the co-moving DM free-streaming length (λ FSH ) and the scaled DM number density at the present epoch n ψ (T 0 )/T 3 0 , for different values of a(t I )/a(t R ) in Table 1. As we see from this table, these quantities without including the collision effects already converge for a(t I )/a(t R ) = 5 × 10 −3 . This implies that the average properties of the distributions without collisions converge by this value of a(t I )/a(t R ). However, since the broadness of the distribution affects the collision probabilities significantly, the values of v(a EQ ) , λ FSH and n ψ (T 0 )/T 3 0 are modified by large factors once collision effects are included. v(a EQ ) and λ FSH reduce more with longer duration of the reheating phase, since the momenta are distributed more towards smaller values, due to enhanced collision probabilities. On the otherhand, for this same reason, the freeze-in contribution increases, thus enhancing n ψ (T 0 )/T 3 0 . For our choice of model parameters, with the reheating temperature T R = 1 TeV and the inflaton mass m φ = 5 TeV, the scaled DM number density increases by upto a factor of around 60, depending upon the value of the scale factor ratio a(t I )/a(t R ). For 2T R < m φ as chosen by us, no resonant enhancement of the hh →ψψ process is possible. For the opposite hierarchy of 2T R > m φ , this effect can enhance the post-reheating production of DM pair from the primeval plasma substantially, and may thus be disfavoured by relic density considerations.
Given its important role, how small can r = a(t I )/a(t R ) be? From Eq. 2.14, we find that where we have used the conditions that the reheating time t R ∼ 1/Γ φ and ρ φ (T R ) = ρ R (T R ). We have further assumed that between t I ≤ t ≤ t R , the Universe is approximately dominated by a stable matter field, and therefore, the scale factor a(t) ∝ t 2/3 . Since the constraints from the cosmic microwave background data on inflationary scenarios restrict the energy density in the inflaton field at the end of the slow-roll phase to be around ρ φ (t I ) 10 16 GeV 4 , for T R = 1 TeV, we have r 10 −17 . If we increase the reheating temperature, the constraint on r becomes stronger. Even though we have chosen r values much larger than the lower bound allowed, as we observe in the above discussion, the physical quantities of interest in the DM sector, such as the average velocity at reheating, its free-streaming length and its density all reach their asymptotic values before collisions by r ∼ 5×10 −3 for the scenarios studied here, although they differ considerably on inclusion of the collisions.
As we can see from Fig. 1, the DM-inflaton coupling chosen for all the above results, λ = 6.6 × 10 −9 , overpredicts the DM relic abundance Ω ψ h 2 by a factor of 50. Thus, to be consistent with the measurements of the DM density, such scenarios would require a 50 times lighter DM particle of mass m ψ = 0.2 keV. While such a light DM would be ruled out by Lyman-α constraints on its free-streaming length if we did not include the effect of collisions, on including the scattering effects it may still be marginally allowed, with an λ FSH 0.1Mpc. The same 10 keV mass DM may also be consistent with the DM density if a small amount of entropy is produced after the decoupling of the collision processes but before the onset of the BBN epoch. We have checked that for the choices of T R and m φ made, if the DM-inflaton coupling is reduced by a factor of around 7 to match the DM abundance with m ψ = 10 keV, the effect of scatterings become less significant, due to the combined effects of the decrease in both the DM initial number density and the scattering cross-sections. This prompts us to vary T R and m φ over a broad range to determine in which region the scatterings become crucial − a numerically very costly task we intend to perform in a future study. Nevertheless, given the analysis presented here, it is clear that the collision effects will be important in a considerable range of the four unknown parameters, being consistent with the DM density requirements.
Can we understand qualitatively the scale of T R and m φ that led to a large effect from scatterings? Since we have considered the effects of on-shell Higgs boson scatterings with the DM particles, the reheating temperature T R should be larger than the Higgs mass m h . For the given interaction Lagrangian, a lower T R would also lead to scatterings with W − bosons, b−quarks or SM leptons, but with a further suppressed rate. Since the scattering reactions must take place over a length of time for redistributing the momenta considerably, we conclude that T R should be somewhat higher than m h , thus requiring T R to be around the TeV scale or higher. In most common inflationary scenarios with a reheating phase, we have T R < m φ . However, with the momentum scale of the SM particles being of the order of T R , m φ cannot be too high compared to T R for the inflaton-mediated scattering cross-sections to be relevant. Thus with T R being around a TeV, we can at most have m φ to be of the order of a few TeV.

Summary
One of the most well-studied non-thermal mechanism invoked to produce DM in the early Universe is through the decay of the inflaton in the post-inflationary reheating phase. In such scenarios, the resulting DM momentum distribution is taken to be the one produced at reheating, which is then red-shifted at later epochs due to the expansion of the Universe. However, in such scenarios since both the SM and the DM fields couple to the inflaton for successful reheating into these sectors, inflaton mediated 2 → 2 scatterings involving the DM and SM particles are necessarily present. These reactions include the self-scattering of a pair of DM particles, and the s− and t−channel inflaton mediated DM-SM scatterings. We address the question of the impact of these scattering processes on the DM momentum distribution in this study. The same parameters that enter the reheating dynamics, and hence the resulting reheating temperature of the SM bath and the density of DM produced, also completely determine the scattering rates. Since the DM momentum distribution directly enters the considerations of structure formation of the Universe, studying the effects of the scatterings is highly relevant observationally as well. This is especially true for light DM, for which constraints coming from, for example, the Lyman-α forest data on its average velocity, and hence the free-streaming properties are stringent.
We find that for a range of the inflaton mass m φ and the reheating temperature T R , including the scattering effects are crucial for determining an accurate DM momentum distribution. For the scenario considered here, in which the inflaton dominantly couples to the SM Higgs boson through renormalizable interactions, this range of T R and m φ is largely controlled by the scale of the Higgs mass m h . In particular, we argue, and find that a reheating temperature T R at the TeV scale and an inflaton mass scale of a few TeV lead to a large effect from the inflaton mediated 2 → 2 scatterings involving the DM and SM particles. We also show, by studying a number of example cases, that for fixed DM mass, DM-inflaton coupling, T R and m φ , the impact of the scattering processes also depend significantly upon the duration of the reheating process.
We show that the inclusion of the effect of scatterings can modify the average DM velocity at matter-radiation equality and the DM free-streaming length by upto a factor of 40 for the scenarios studied. This makes including the scattering effects crucial while considering constraints from the Lyman-α forest data. Since the scattering processes tend to redistribute the momenta in such a way that more lower momentum modes are filled, the DM average velocity goes down, and hence the free-streaming length. Therefore, inclusion of the scattering effects makes a much larger range of DM mass allowed from structure formation considerations.
A lot remains to be explored in this direction. First of all, within the scenario studied, a full exploration of the parameter space, being a numerically very challenging task, is left for a future work. As we emphasized above, the scale T R and m φ found to lead to the maximum impact of scatterings is essentially set by the scale of Higgs mass m h . However, this is specific to the interaction Lagrangian considered, which is motivated by the fact that the inflaton has a renormalizable coupling with the Higgs field in the SM. However, in different scenarios for inflation and reheating, the inflaton might have a dominant coupling to other SM fields, thus changing the relevant DM-SM scattering processes. In such a scenario, a different scale of T R and m φ will lead to the maximum impact of the scattering processes, which should be studied in detail. As shown by the results obtained in this study, precision determination of the momentum distribution of DM produced in inflaton decay including the effect of inflaton mediated scatterings is highly relevant in determining constraints from data on cosmological large-scale structure.