Prethermalization Production of Dark Matter

At the end of inflation, the inflaton field decays into an initially nonthermal population of relativistic particles which eventually thermalize. We consider the production of dark matter from this relativistic plasma, focusing on the prethermal phase. We find that for a production cross section $\sigma(E)\sim E^n$ with $n>2$, the present dark matter abundance is produced during the prethermal phase of its progenitors. For $n\le 2$, entropy production during reheating makes the nonthermal contribution to the present dark matter abundance subdominant compared to that produced thermally. As specific examples, we verify that the nonthermal contribution is irrelevant for gravitino production in low scale supersymmetric models ($n=0$) and is dominant for gravitino production in high scale supersymmetry models ($n=6$).


I. INTRODUCTION
The non-thermal universe at the end of inflation is often assumed to be hidden from observations due to subsequent thermalization and radiation domination of the universe prior to Big Bang Nucleosynthesis. However gravitational effects and stable non-gravitational relics, if produced in this period, are sensitive to/provide potential probes of this era. Examples of gravitational effects include production of gravitational waves, changes in the expansion history of the universe, and production of primordial black holes. Non-gravitational relics include the matter-antimatter asymmetry and the dark matter (DM) abundance [1]. In this paper, we provide conditions under which the present DM abundance can be sensitive to the non-thermal phase immediately following inflation, and carry out a detailed calculation of this abundance. More specifically, we study concrete DM models where non-thermal conditions change the predictions from the thermal scenario by many orders of magnitude.
We consider the case where DM is primarily produced via freeze-in from a relativistic plasma, which in turn is produced via a small coupling to the inflaton. 1 In this paper we restrict ourselves to perturbative processes during reheating after inflation. Under these conditions we arrive at a rather simple conclusion: If the total cross-section for DM production from the relativistic plasma has a sufficiently steep dependence on energy: σ(E) ∝ E n with n > 2, and the thermalization of the inflaton decay products is mediated by the Standard Model (SM) gauge interactions, then it is the production of DM from the non-thermal phase of the relativistic plasma, at the earliest stages after the end of inflation, that predominantly determines the DM abundance. * marcos.garcia@rice.edu, mustafa.a.amin@rice.edu 1 We do not consider the case of DM produced non-thermally through direct decays; this scenario is studied in [2,3].
For gauge-mediated thermalization, our results supersede other ultraviolet (UV) sensitive DM production scenarios considered in the literature. In [4], one of us found that the dependence of the DM abundance on the highest temperature during reheating T max is obtained for n ≥ 6, assuming a relativistic plasma in thermal equilibrium. 2 . In the present paper, by including non-thermal effects, we find UV-sensitivity at the milder power n > 2, which makes our result relevant for a broader range of models. [5][6][7][8][9][10] The rest of the paper is organized as follows. We briefly review the DM production mechanisms in the literature in section I A. Therein we also introduce relevant notation and definitions, and we schematically summarize our results. We carry out the detailed calculation of nonthermal DM production in section II. We then specialize to the application of the formalism to gravitino production in section III, in the case of weak scale supersymmetry breaking (n = 0), as well as that of very high energy supersymmetry breaking models, where the population of almost the whole spectrum of superpartners is kinematically forbidden during reheating (n = 6). Section IV is reserved for our summary and future applications.

A. Dark Matter Production Scenarios (Review)
Particle DM production scenarios can be roughly divided in two classes. 3 In thermal freeze-out models, the DM particles (hereafter denoted by χ) interact sufficiently strongly at high energies with the primordial relativistic plasma (denoted by γ) keeping γ and χ in thermal equilibrium. Once χ becomes non-relativistic and the expansion rate overtakes the annihilation rate, the DM yield 45 freezes-out of equilibrium. Y χ is therefore insensitive to the early thermal history of the Universe up to the decoupling temperature, T dec 20 m DM for processes mediated by the weak interaction [13]. That is, freeze-out is infrared (IR) dominated.
Let us assume that DM production takes place via 2 → 2 scattering processes in the presence of a thermalized background, with the integrated effective cross section where √ s = E is the center of mass energy of the scattering. Fig. 1 schematically shows freeze-in scenarios for different values of n, which we discuss further below. For non-negative n, this type of energy dependence will generically arise if the DM sector, decoupled at low energies, communicates to the visible sector via high dimensional operators. The mass scale M in (2) can be thought to be parametrically related to the mass of a heavy mediator in the UV theory.
In the absence of entropy production (for example, during the radiation dominated era that follows after the end of inflationary reheating) freeze-in is UV-dominated for n > −1 [3,4,15]. This result then holds under the assumption that the inflaton φ instantaneously decays into a thermal bath of initial temperature T reh = 40 at the decay time t reh = Γ −1 φ , where g reh denotes the number of effective degrees of freedom in the thermal bath, Γ φ is the (perturbative) inflaton decay rate, and M P = (8πG) −1/2 is the reduced Planck mass.
However, the reheating process is not instantaneous, but involves a continuous transferal of inflaton energy 4 We denote by γ a single degree of freedom of the relativistic plasma; for a species in thermal equilibrium, nγ → n T γ = ζ(3) π 2 T 3 . 5 In the absence of entropy conservation, this ratio is better suited for calculations than the DM-to-entropy ratio [3].
density into its relativistic decay products, until the former is depleted around the time when the universe becomes dominated by radiation. Entropy is not conserved during reheating, and DM is continuously produced as the thermal plasma is populated. A more detailed calculation which takes into account the finite duration of the reheating process shows that under the assumption of an instantaneously thermalized background at the maximum temperature during reheating T max ∼ (m φ /Γ φ ) 1/4 T reh [3,[57][58][59], the freeze-in DM yield takes the form for n ≥ 0, n = 6. The special n = 6 case yields a logarithmic dependence: Y T χ ∝ T 7 reh ln(T max /T reh ) [4]. Remarkably, this result implies that for n ≥ 6 the temperature dependence of the DM yield is dominated by T max , that is, DM production is UV-dominated during reheating. Only such a steep temperature dependence of the scattering cross section can be sufficiently competitive with respect to the dilution rate due to entropy production.
It is worth noting that in UV-dominated models, the decay of the inflaton determines the DM yield. In the present epoch, T T reh , the DM is non-relativistic, and its relic abundance will be given by Accounting only for Standard Model degrees of freedom, g reh = 427/4 and g(T 1 MeV) = 3.91, and a DM mass m χ = 100 GeV, the observed DM abundance Ω χ h 2 0.12 [60] requires the yield at the end of reheating to be Y χ 10 −9 1. The smallness of Y χ will constrain us to focus on Planck-suppressed decay rates or smaller for φ,

B. Non-thermal Effects in Dark Matter Production
The high-energy sensitivity discussed above invites us to consider non-thermal effects. Recall that when an inflaton decays, it produces high energy quanta with energy ∼ m φ > T max . Until elastic and inelastic processes within the plasma become sufficiently efficient to maintain kinetic and chemical equilibrium, these notyet-thermalized decay products can then copiously produce DM particles if the production cross section is UVsensitive. Therefore, if thermalization is significantly delayed compared to the time-scale attaining T max , the bulk of the DM abundance may be produced before the onset of thermal equilibrium in the universe. As we will show later, in scenarios for which thermal equilibrium is reached via gauge interactions, this delay is sufficient to lead to a non-thermally produced DM abundance for n > 2.  The evolution of the inflaton field φ (black curve) and its relativistic decay products γ (gray shading). The time t end denotes the end of inflation, t th is the time of full thermalization of γ, and t reh denotes the end of reheating when the inflaton decays almost entirely to γ. After t reh , the universe becomes radiation dominated. Dark matter χ can be produced from γ for t > t end in the non-thermal and thermal phases. Lower panel: Different paths to the observed dark matter abundance. For a production cross section σ(E) ∝ E n , the production of χ from thermal γ dominates for 2 ≥ n (gray), while the production of χ from non-thermal γ dominates for n > 2 (black). The n > 2 case is sensitive to the earliest stages at the end of inflation.
A schematic summary of our results is shown in Fig. 1. The upper panel shows the typical evolution of the inflaton field and its decay products (shown as a shaded background), starting from the slow roll of φ in the inflationary era, continuing through the reheating epoch, dominated by the pressureless oscillation of φ, and finishing in the radiation dominated era. Note the finite non-thermal epoch, spanning times between the end of inflation at t end and the time of formation of a fully thermalized background at t th . The lower panel shows the evolution of the DM yield in freeze-in models as a function of n in (2) For n < −1 the cross section favors lower energies, and the bulk of the DM is produced around the time it becomes non-relativistic. For 2 ≥ n > −1, the production is enhanced at higher energies, but this enhancement cannot compete against the dilution due to the decay of the inflaton, and as a result the population created around t reh dominates. For n > 2, the DM population produced by the non-thermal particles survives after the end of reheating and constitutes the majority of the DM particles.
We next consider the pre-thermalization epoch in detail. The treatment is mostly analytic. Numerical calculations for specific models are done in section III.

II. PRE-THERMALIZATION PARTICLE PRODUCTION
The transfer of the energy density stored in the inflaton field to its decay products during perturbative reheating can be approximated by the following set of equations, where we denote by ρ φ and ρ γ the energy densities of the inflaton condensate and that of its relativistic decay products, respectively. Since we are interested in a freezein scenario, we implicitly assume that the contribution of DM particles to ρ γ is negligible. We also assume that no direct decay channel of the inflaton to DM is available.
A. Before Thermalization of γ In the very earliest stages of reheating, Γ φ (t−t end ) 1 where the subscript "end" denotes the end of inflation, the thermalization of the inflaton decay products will be at best incomplete. In the limit in which the created particles are not interacting, we can then approximate their number density at a given time by counting the number of inflaton quanta that have decayed up to that instant. Straightforward integration of (6) leads to where a denotes the scale factor, m φ is the inflaton mass, and is a natural dimensionless time variable, v reh = 1. Let us parametrize the inflaton decay rate as Γ φ = m φ |y| 2 /8π, and denote byg the effective number of relativistic degrees of freedom contributing to the number density of the thermal plasma. It is then easy to verify that the non-interacting decay products are in a regime of under-occupation with respect to their thermal counterparts,g n T γ /n γ > 1, for |y| O(10 −5 ), until number-increasing processes become efficient [3,61].
This "weak reheating" regime will be generically realized for Planck-suppressed inflaton decays, and, as we will discuss in the following sections, it corresponds to the regime allowed by gravitino production.
In order to compute the DM production rate we need not only (10), but the phase space distribution function f γ (p, t) of the decay products. It is illustrative to derive it directly from the Boltzmann equation. In a comoving frame, the inflaton particles sit in a condensate at rest, with a distribution given by f φ = (2π) 3 n φ δ 3 (p). The resulting Boltzmann transport equation for the process φ → γ + γ, assuming a two body decay for simplicity, is then given by At very early times the exponential in (9) may be approximated as a constant. The transport equation can then be integrated to yield As expected, the energies of these "hard" decay products are centered at much higher values, p ∼ m φ , than their would-be thermal counterparts [62,63]. Let us consider now the DM creation process γ + γ → χ + χ. Recalling that the progenitors are under-occupied with respect to their thermal distributions, the Boltzmann equation for the DM density n χ can be simply written aṡ where σ(s) ≡ σ(s) γ+γ→χ+χ , s the Mandelstam variable, and where we have accounted for the number of degrees of freedom involved in the process. Upon substitution of (13), and following the same steps as in the usual thermal calculation (see for example [64,65]), we can rewrite the previous expression aṡ where the effective non-thermal cross section σv NT is defined by the expression above, in analogy with the thermal case.

B. Thermalization of γ
As reheating progresses, number-conserving and violating interactions eventually bring the plasma into kinetic and chemical equilibrium. For processes mediated by a gauge interaction with coupling α, inelastic 2 → 3 splittings at small angles become effective at early times and rapidly populate the relativistic plasma with "soft" (p m φ ) particles. In turn, these (over-populated) soft products eventually produce a fully thermalized background, that is capable of cooling down the more energetic decay products. A systematic tracking of the phasespace distributions of the hard and soft sectors reveals that the complete depletion of the hard non-thermal sector is achieved at [61,66] v th α −16 which occurs well before the end of reheating given that Γ φ M 3 P α 8 /m 2 φ , this being the case for Plancksuppressed decays and SM gauge interactions (α ∼ 10 −2 ). It is worth noting that the temperature associated with the time scale (17), namely, is smaller by a factor of ∼ α 4/5 (|y|M P /m φ ) 3/10 than the maximum temperature T max |y| −1/2 T reh assuming instantaneous thermalization.
Let us now go back to (15) and note that all the time dependence is contained in the non-thermal radiation number density n γ . For simplicity let us assume that the transition from the non-thermal distribution f γ to its thermal form occurs suddenly at the time (17). We can then integrate with respect to time to obtain the DM number density at v th , Well before the end of reheating, when the equationof-state parameter w = p/ρ 0, the scale factor evolves as [3] where Here the O(1) factor in the second equality is approximately equal to 2.8 for Starobinsky inflation, and to 1.3 for a quadratic potential. Since v max ∼ A v th v reh , the integral in (19) can be evaluated to give with v th given by (17). After the onset of thermal equilibrium, eq. (15) becomeṡ where σv T denotes the thermally averaged production cross section, and for which (22) should be taken as initial condition. Alternatively, given the linearity of the previous equation with respect to the DM number density, we can treat the thermal and non-thermal abundances as independent populations. Thus, for v > v th , the non-thermal occupation number simply redshifts. At the end of reheating, when entropy production has ceased, we can estimate the non-thermal yield as For t > t reh , this quantity is conserved (up to changes in the number of relativistic degrees of freedom) assuming that the universe has a standard thermal history after reheating. 6

C. Comparison of Thermal vs. Non-thermal Yield
With the previous results at hand, we can estimate the ratio of the non-thermal to the thermal abundance if we, for simplicity, assume that for a given model the total scattering cross section can be written as where M is a mass scale relevant to the model, λ is the dimensionless coupling and c n = 32ζ(3) 2 /2 4+n πΓ(2 + n 2 )Γ(3 + n 2 ). This convoluted numerical factor is chosen so that, assuming Maxwell-Boltzmann-distributed progenitors [64] σv T g 2 In turn, where Analogous expressions for f (n) may be written if we assume the correct Fermi-Dirac or Bose-Einstein thermal distributions for γ. At n = 0 the difference between these expressions and (29) is 31%, and decreases exponentially for larger n, it being 3% for n = 6. This correction is negligible compared to the dependence on n and we will disregard it for the present discussion, but not in the detailed calculations of the following section.
Substitution of (28) into (24) gives For comparison, note that the purely thermal yield at the end of reheating is given by [4] Y T χ 96ζ(3)g 2 χ g 2 γ λM P T 7 reh √ 40g Contours of the non-thermal to thermal yield ratio at the end of reheating, Rχ = Y NT χ /Y T χ , in the (n, |y|)-plane. The parameter n determines the energy dependence of the scattering cross section, σ(E) ∝ E n , and y is the effective dimensionless coupling that parametrizes the decay rate of the inflaton, Γ φ = |y| 2 m φ /8π. Thermalization is assumed to proceed via the strong gauge interaction, including y-dependent radiative corrections. In the shaded region non-equilibrium particle production is non-thermally dominated. The width of the level curves corresponds to the difference in gauge coupling running and in degrees of freedom between the Standard Model and its minimal supersymmetric extension. The blue six-pointed star signals the upper bound for light gravitino production, while the red five-pointed star corresponds to the observed DM abundance for the EeV gravitino.
We can write the ratio of the non-thermal yield to the thermal one, R χ ≡ Y NT χ /Y T χ ; for n < 6 it has the form R χ 0.2 (6 − n)f (n) |y| By itself this expression is rather obscure, which is why we show selected level curves for R χ in the (n, |y|)-plane in Fig. 2. For definiteness, we assume thermalization through the strong interaction, with a y-dependent running coupling, and we consider both a pure Standard Model spectrum, or an MSSM-like spectrum; the difference is accounted by the width of the contours. In both scenarios, we arrive to the conclusion that the DM abundance produced before thermalization is the dominant constituent for n > 2. Note that the expression (25) is an over-simplification, as it assumes that the same production channels are available before and after the plasma reaches thermal equilibrium.

III. APPLICATION TO TWO MODELS
We now consider two applications of the formalism developed in the previous section. We will first revisit the non-thermal production of light gravitinos, and we will recover the results of [3] in a simpler way. We then study non-thermal production of an EeV gravitino [5][6][7][8]. We will not only focus on the final abundances, but we will also track numerically the evolution of the thermal and the non-thermal abundances during reheating, with the help of eqs. (15) and (17).
In order to avoid discontinuities in the numerical solution, we include the effect of the depletion of the nonthermal distribution and the population of its thermal counterpart for v < v th . Following [61,66], the depletion of the non-thermal γ population can be accounted by substituting where denotes the splitting momentum below which a particle can deposit an order one fraction of its energy and participate in the thermal plasma. In turn, the soft component of the plasma can be considered to be partially thermalized for v v th , and an effective temperature can be associated to it. This effective temperature is lower than the naive estimate assuming instantaneous thermalization (v th → v max ), and can be used to determine the corresponding thermal number density n T γ . We do not present here the explicit form of this effective temperature, 7 since the two corrections mentioned above are negligible, at least compared to the correction to the time dependence of the scale factor (20) arising from the continuous transition from matter to radiation domination around v reh .

A. Light Gravitino
Gravitino production in weak scale supersymmetry breaking models has been extensively studied in the literature, the spin-3/2 particle being the poster child for the cosmological moduli problem [3,. The production cross section is simply proportional to M −2 P , and it is dependent on the instantaneous temperature only through radiative corrections to the gauge and Yukawa coupling. In the presence of a thermalized plasma, the average cross section for the production of gravitinos can be written as [50,51,53] σv = σv top + σv gauge , with σv top = 1.29 where A t is the top-quark supersymmetry-breaking trilinear coupling, and where the mg i are the gaugino masses and the constants c i , k i depend on the gauge group (see [3] for details). The first term in the gaugino mass-dependent factors (1 + m 2 gi /3m 2 3/2 ) corresponds to the production of the transversally polarized gravitino, while the second term is associated with the production of the longitudinal (Goldstino) component. We therefore have n 0. For thermally produced gravitinos the estimate derived from the assumption of instantaneous reheating is a good approximation, as the population produced at T > T reh is rapidly diluted by the entropy produced in subsequent inflaton decays [3,4,[54][55][56]. The yield at low temperatures T 1 MeV can be computed to give Overproduction of DM is then averted for |y| 10 −5 [3]. Let us now for definiteness consider a specific scenario in which the inflaton decays predominantly into gauge bosons, φ → g + g, as is the case in no-scale supergravity models of inflation with a non-trivial gauge kinetic function [70][71][72]. The pre-thermalization gravitino abundance can then be computed as that given by the channel g + g →g +G, with amplitude [50] where t, s are Mandelstam variables. The resulting nonthermal production cross section is then given by Since the weak-scale gravitino corresponds to n = 0 in (25), we expect that the early produced non-thermal abundance will be overwhelmingly drowned by the thermal population produced during the final stages of reheating; this is a known fact [3]. In Fig. 3 we show Ratio of the instantaneous yield to its final value after reheating assuming instantaneous thermalization, , for the weak-scale gravitino. Here v = Γ φ (t − t end ). Shown are the instantaneous reheating approximation (green, dotted), the instantaneous thermalization approximation (orange, dashed), the pure non-thermal population (blue, dot-dashed), and the combined thermal + non-thermal result (black, solid). The peak in the black curve is caused by the definition of Y 3/2 in terms of a single γ degree of freedom, instead of all degrees of freedom present in the plasma.
the results for the numerical integration of the transport equation for n 3/2 assuming (i) instantaneous reheating and thermalization at v = 1 (dotted, green), (ii) instantaneous thermalization at T max with no non-thermal production (dashed, orange), (iii) pure non-thermal production stopping at v th 10 −6 (blue, dash-dotted), and (iv) combined thermal and non-thermal production (black, continuous). In the later case the yield is defined as Y NT 3/2 = n 3/2 /n NT γ for v < v th . For simplicity we suppress the top quark contribution and the longitudinal component. The model parameters are chosen to coincide with the upper bound for DM production, m φ 3×10 13 GeV, g reh = 915/4 and |y| 2.7 × 10 −5 . For these, we obtain a final non-thermal population which is a factor of ∼ 3 × 10 −7 less abundant than the thermally produced one, in agreement with the simple estimate (32), which predicts R 3/2 ∼ 10 −6 . Pre-thermalization effects are negligible for the light gravitino, and its production is IR-dominated with respect to the reheating process, although it is UV-dominated with respect to the subsequent radiation-dominated epoch.

B. Heavy Gravitino
Let us consider now gravitino production in high scale supersymmetry models, where the population of the spectrum of superpartners is kinematically forbidden during reheating, since their masses are assumed to be larger than m φ . Only the longitudinal component of the gravitino, whose mass is significantly above the weak scale, m 3/2 ∼ 1 EeV, may be produced during reheating [5][6][7][8]. 8 In this scenario, the dominant production channel is mediated by effective SM + SM →G +G vertices that are suppressed by the supersymmetry breaking scale F = √ 3m 3/2 M P . The total averaged production rate can be written as [5] σv and it includes production from scalar, fermion and vector initial states. For these models n = 6. Under the assumption of an instantaneously thermalized plasma, the final yield acquires the T max -dependent correction of the form [4] where Y 3/2, instant ∝ T 7 reh M P /F 4 . We now evaluate the non-thermal yield. Similarly to the previous case, let us assume for simplicity an inflaton that decays into gauge bosons. The amplitude for the process g + g →G +G has then the form [5] with t, s Mandelstam variables. Including the contributions of all 12 SM gauge bosons, the total non-thermal cross section can be computed from (15) to be If we assume that the late-time DM abundance is primarily populated by these non-thermally produced relics, we can substitute the previous expression in (24), with g reh = 427/4, to obtain after reheating, and from (5), Ratio of the instantaneous yield to its final value after reheating assuming instantaneous thermalization, , for the EeV gravitino. Here v = Γ φ (t − t end ). Shown are the instantaneous reheating approximation (green, dotted), the instantaneous thermalization approximation (orange, dashed), the pure non-thermal population (blue, dot-dashed), and the combined thermal + non-thermal result (black, solid). The dot-dashed blue curve and the continuous black curve overlap for all v. Despite the huge enhancement for the non-thermal yield relative to the thermal yield, its maximum value is Y NT 3/2 4 × 10 −6 1, consistent with the non-equilibrium assumption.
at T 1 MeV. Substitution of (17) assuming again thermalization through the strong interaction finally gives Therefore, the observed DM relic density is obtained for |y| 2.9 × 10 −7 , or T reh 2 × 10 8 GeV. In comparison, the assumption of instantaneous reheating and thermalization yields [6] Ω inst 3/2 h 2 0.11 0.1 EeV which underestimates the final abundance by a factor of ∼ 10 −14 given the same y. Fig. 4 shows the results for the numerical integration of the transport equation for n 3/2 under the same assumptions discussed in the previous section. The model parameters are chosen to yield the observed DM abundance; they coincide with their analytical estimates in (48) within ∼ 5%. Here the gauge coupling running accounts only for SM degrees of freedom, with v th 10 −7 . It is clear from the figure that, by the end of reheating, the overwhelming majority of the DM yield is produced before thermalization, with the instantaneous thermalization assumption underestimating it by a factor of 4 × 10 12 , and the instantaneous reheating approximation by ∼ 10 14 . It is also worth noting that, despite the enormous enhancement of n 3/2 relative to the thermal result, one can check that Y 3/2 1 at all times, consistent with the non-equilibrium approximation.

IV. CONCLUSION
In this paper, we considered the production of DM from a relativistic plasma produced perturbatively at the end of inflation. We provided a simple criteria to determine when DM production from the highly energetic, non-thermal population of the relativistic plasma provides the dominant contribution to the present DM abundance. For the production cross section σ(E) ∼ E n , if n > 2, the non-thermal contribution dominates, whereas for n ≤ 2 the production of entropy dilutes away the nonthermal contribution. In the former case, the time-scale for the onset of thermalization is crucial to determine the particle abundances at late times, due to the higher energies of the inflaton decay products compared to their thermally distributed counterparts. In the latter case, the usual calculations assuming a thermal plasma suffice.
We tracked the (non-thermal) phase space distributions of the relativistic plasma (produced by the perturbative decay of the inflaton) as well as the DM distribution function (sourced by the relativistic plasma) using coupled Boltzmann equations with certain simplifying assumptions. Using this analysis, we provided analytic estimates of the non-thermal yield of DM. When compared with the thermal yield, the non-thermal yield can be larger by many orders of magnitude (see Fig. 2). Under the assumption of gauge-mediated thermalization and Planck-suppressed inflaton decay, eq. (24) provides a simple analytical estimate for the non-thermal DM yield at the end of reheating.
Our determination of the non-thermal production rate benefits from several simplifying approximations. We made use of the non-interacting decay distribution function to integrate the transport equation to determine n χ . However, the relativistic plasma slowly but steadily evolves into a thermalized fluid and its constituents interact sufficiently often to create DM. For UV-dominated models, it may then be necessary to consider a more careful account of the production mechanism, such as the backreaction of DM production on the thermalization rate.
In terms of concrete models, we considered thermal and non-thermal creation of gravitinos in supersymmetry breaking models with a superparticle spectrum lying at the weak scale, or above the inflationary scale. In the former case, the early primordial abundance is diluted by entropy production, and the instantaneous reheating and thermalization assumption is a remarkably good estimate (as expected, since n = 0 here). In contrast, for the EeV gravitino, the steep dependence of the scattering cross section on the momenta of the progenitors (n = 6) leads to a relic density that is very sensitive to the pre-thermalization epoch; the observed abundance is obtained for a reheating temperature two orders of magnitude smaller than that assuming instantaneous decay at t reh . The derivation of the DM production rate presented here relies on four ingredients: (i) a slowly decaying, oscillating scalar field that dominates the energy budget of the universe, (ii) whose coupling to the dark sector is negligible, (iii) and which injects entropy into the universe in the form of relativistic degrees of freedom, (iv) which in turn do not thermalize instantaneously. Although these ingredients can be clearly present in the reheating epoch that follows the end of inflation, they will also generically be present during a possible moduli-dominated era. Their effect must be accounted in addition to other "nonthermal" phenomena [35,68,69,[73][74][75][76][77][78][79][80][81][82][83][84]. Our results also open the door for consideration of non-perturbative, preheating effects in determining the DM abundance in UVsensitive scenarios, which will be pursued elsewhere. It must be noted that while non-thermal effects can impact the DM abundance, this abundance alone cannot be a model independent probe of the non-thermal processes. Additional observables, or a definite model are required to make this inference possible.
Finally, while we have focussed on DM production, a similar analysis can be applied to the production of feebly interacting relativistic relics (dark radiation) through their contribution to the effective number of light relativistic species N eff and its effect on the cosmic microwave background [85][86][87][88][89][90][91]. We leave this study for a future publication.