Nonthermal Production of Dark Matter Prior to Early Matter Domination

Freeze-out or freeze-in during a period of early matter domination can yield the correct dark matter abundance for 〈σannv〉f < 3 × 10−26 cm s−1. However, in a generic non-standard thermal history, such a period is typically preceded by other phases. Here, we study nonthermal production of dark matter in a simple post-inflationary history where a radiation-dominated phase after reheating is followed by an epoch of early matter domination. Focusing on the freeze-in regime, we show that dark matter production prior to early matter domination can dominate the relic abundance in large parts of the parameter space, including weak scale dark matter masses, and the allowed regions are highly dependent on the entire post-inflationary history. Moreover, for a very broad range of 〈σannv〉f spanning over several decades, dark matter particles can start in chemical equilibrium early on and decouple during early matter domination, thereby rendering the relic abundance essentially independent of 〈σannv〉f . We briefly discuss connections to different observables as a possible means to test the elusive freeze-in scenario in this case.


I. INTRODUCTION
Despite various lines of evidence for the existence of dark matter (DM) [1], its identity remains a major problem at the interface of cosmology and particle physics. Weakly interacting massive particles (WIMPs) have long been the focus of direct, indirect, and collider searches for DM. Thermal freeze-out in a radiation-dominated (RD) universe can yield the correct DM abundance if the annihilation rate takes the nominal value σ ann v f = 3 × 10 −26 cm 3 s −1 , which is the ballpark value for weak scale masses and interactions (hence called the "WIMP miracle"). However, this scenario has come under increasing scrutiny by recent experiments, most notably the Fermi-LAT results from observations of dwarf spheroidal galaxies [2] and newly discovered Milky Way satellites [3]. A recent analysis [4] has ruled out thermal DM with a mass below 20 GeV in a model-independent way (unless there is P -wave annihilation or co-annihilation), while masses up to 100 GeV can be excluded for specific annihilation channels.
The situation can change if the universe is not in a RD phase at the time of freeze-out [5]. An important example is an epoch of early matter domination (EMD), which is a generic feature of early universe models arising from string theory constructions (for a review, see [6]). In this context, an EMD era is driven by long-lived scalar fields (for example, string moduli) that are displaced from the minimum of their potential during inflation, come to dominate the energy density of the post-inflationary universe, and eventually decay to form a RD universe prior to big bang nucleosynthesis (BBN). Both thermal freeze-out and thermal freeze-in during EMD can produce the correct DM relic abundance for small annihilation rates σ ann v f < 3 × 10 −26 cm 3 s −1 [7,8].
In general, however, an EMD epoch is only one of the stages in the post-inflationary universe. Unless it is driven by oscillations of the inflaton itself, it is typically preceded by a RD phase or a period with a more general equation of state. Thermal freeze-out or freeze-in at these early stages can also contribute to the DM relic abundance (for example, see [9][10][11]) and thereby affect the allowed parameter space. Since freeze-out occurs at temperatures below the DM mass m χ , pre-EMD production will only be relevant for (very) large DM masses in this case. On the other hand, in the case of freeze-in, the main contribution to the DM abundance can arise when DM particles are (ultra) relativistic. Thus, pre-EMD production can have a significant affect on the freeze-in side of the parameter space for values of m χ at the weak scale or below. In such cases, an exact calculation of the DM relic abundance requires knowledge of the earlier stages of the non-standard thermal history.
In this work, we perform a detailed study of freeze-in within a non-standard thermal history that involves a RD phase after inflationary reheating followed by a period of EMD. Such an EMD phase is characterized by two distinct periods: a transition during which the initial radiation energy density redshifts away, which we call the memory phase, and the usual EMD behavior once entropy production becomes significant. We calculate the contributions to the DM relic abundance from production during the prior RD phase as well as the memory phase for the case that σ ann v f is constant over the temperature range of interest. This early contribution to DM production depends on the temperature of the universe at the onset of the EMD epoch T 0 and at the completion of inflationary reheating T reh , in addition to that at the end of the EMD epoch T R . We show that pre-EMD production can totally dominate the DM relic abundance in large parts of the m χ − σ ann v f plane, and the allowed parameter space is highly sensitive to T 0 and T R . A particularly notable observation is that the relic abundance is virtually independent of σ ann v f for a very broad range of σ ann v f , spanning over many decades, where DM particles start in chemical equilibrium in the pre-EMD epoch and decouple later on. In this case, measurement of m χ at collider experiments, in combination with other cosmological implications of an EMD epoch, may be used as a potential probe of the elusive freeze-in scenario.
The rest of this paper is organized as follows. In Section II, we discuss a simple non-standard thermal history and calculate the pre-EMD contributions to DM production via freeze-in. In Section III, we present our main results including the allowed regions in the m χ − σ ann v f plane and their sensitivity to the history prior to EMD. We elaborate on the correlation between m χ and T 0 /T R in large parts of the parameter space in Section IV, and discuss connections to observables as well as implications for a modulus-driven EMD. We conclude the paper with some discussions in Section V. Some of the details of our calculations are included in the Appendix.

II. DARK MATTER ABUNDANCE IN THE FREEZE-IN SCENARIO
As mentioned, a period of EMD naturally arises in a well-motivated class of early universe models. However, it is typically only one of the stages in the post-inflationary history. We consider a simple scenario that starts with a RD phase at the end of inflationary reheating (for reviews, see [12,13]), followed by an EMD epoch driven by oscillations of a long-lived scalar field, or non-relativistic quanta produced in the post-inflationary universe [14][15][16][17]. A standard RD universe is established at the end of EMD and before BBN.
Here, we are mainly interested in the evolution of temperature and the freeze-in production of DM without delving into the details of inflationary reheating, the specific particle physics origin of EMD, or the explicit models for DM freeze-in (as done, for example, in [18]). This can be done by introducing three parameters: the largest temperature in the RD phase after inflationary reheating T reh , the temperature at the onset of EMD T 0 , and the highest temperature in the subsequent RD phase T R . The corresponding Hubble expansion rates are denoted by H reh , H 0 , and H R respectively.
In calculating the DM relic abundance, we consider the case where σ ann v f is constant over the temperature range of interest (namely, T R T T reh ). This can be the case, for example, if DM interaction with the standard model (SM) particles is mediated by particles whose mass is above T R . Also, without loss of generality, we assume that the DM particle (denoted by χ) represents one degree of freedom. In Section V, we will briefly comment on cases with T dependence of σ ann v f and more general thermal histories involving multiple epochs of EMD separated by RD phases. The details of our calculations can be found in the Appendix.

A. The Standard Lore
Let us briefly recap the well-known case of freeze-in during the late phase of EMD, where decay of the matter component determines the temperature evolution. Deep inside the EMD era, H R H H 0 , there is a subdominant radiation component due to decay of the species that drive EMD. Assuming that the decay products thermalize promptly, a thermal bath forms with the instantaneous temperature [7,8]: where g * is the number of relativistic degrees of freedom at temperature T , and M P is the reduced Planck mass 1 .
In the freeze-in scenario, σ ann v f is so small that DM particles do not reach thermal equilibrium during EMD, and are produced from annihilations of the SM particles. The main contribution to the DM abundance occurs when T ∼ m χ /4 [7,8]. The reason being that particles produced at higher temperatures are quickly diluted by the Hubble expansion, while production at lower temperatures is Boltzmann suppressed. The relic abundance from freeze-in during EMD is given by [7,8]: where g * ,mχ/4 is the number of relativistic degrees of freedom at T = m χ /4. For a given DM mass, the maximum value of σ ann v f in the freeze-in regime can be approximately found by setting the DM number density n χ equal to its thermal equilibrium value n χ,eq at T ∼ m χ /4. For larger values of σ ann v f , DM particles reach equilibrium with the thermal bath and hence production transitions to the freeze-out regime. For T R ∼ (1-10) GeV, the maximum value of σ ann v f in the freeze-in regime lies within the 10 −33 − 10 −32 cm 3 s −1 range [8].

B. Production Prior to Early Matter Domination
In general, an EMD period can be more complicated than that in the previous Subsection. Particularly, if the abundance of radiation during EMD is ever significantly larger than the contribution from decay, then the relation between T and H will deviate, for a time, from that in Eq. (1). The temperature evolution throughout the full EMD era may thus have multiple phases. A large abundance of radiation can arise from the presence of additional decaying components (see, for example, [19]), or, more simply, from a RD phase that precedes EMD. In the post-inflationary history we are considering, a period of RD is present before EMD resulting in a large abundance of radiation at the onset of EMD 2 . This radiation then redshifts away until the decay contribution becomes dominant, recovering Eq. (1). The temperature evolution during the full EMD period therefore has two phases: an initial phase during which the memory of the prior radiation is being erased, followed by the usual decay-driven phase.
We now consider the evolution of temperature in the two periods prior to the entropy-producing phase of EMD.
RD phase prior to EMD: The universe is in a RD period for H 0 H H reh , during which radiation simply redshifts and temperature is inversely proportional to the scale factor: T ∝ a −1 . This results in the standard relation between T and H: Memory phase of EMD: The radiation and matter components have comparable energy densities at H = H 0 , which signals the beginning of the EMD era 3 . Although H ∝ a −3/2 for H R H H 0 , the existing radiation dominates over the contribution from the decaying matter component(s) driving EMD for some time, and continues to redshift. As a result, T ∝ a −1 for H tran H H 0 , where: and the relation between T and H is 4 : T ≈ 90 Once H H tran , the memory of the initial radiation is completely erased and the T dependence on H transitions to that of Eq. (1) until the end of EMD.
We show the evolution of T in terms of H in Fig. 1, derived from numerical calculations, through the entire nonstandard thermal history considered here for T reh = 10 12 GeV, T 0 = 10 10 GeV, and T R = 10 GeV. The evolution during the different stages is in very good agreement with the relations given in Eqs. (1,3,5).
As shown in the Appendix, the relation between T and H in these two phases implies that the main contribution to DM production occurs at the highest temperature in each phase. Hence, with T reh T 0 , the DM relic abundance is set in the pre-EMD epoch. This pre-EMD component of the relic abundance has two distinct regimes based on the value of σ ann v f . If it is very small, DM particles will not be able to establish chemical equilibrium throughout the post-inflationary history. Production in this decoupling regime will dominantly occur at T T reh . If σ ann v f is large enough, then DM particles will reach chemical equilibrium in the pre-EMD phase but can decouple at H H tran , hence the early-equilibrium regime. The condition on σ ann v f to be in either of these regimes is given in Eqs (16,20) of the Appendix: The pre-EMD contribution to the DM abundance follows from Eqs. (19,22): Let us underline some important differences between the pre-EMD and EMD components of the relic abundance given in Eqs. (7) and (2) respectively: • The pre-EMD component depends on T 0 and T reh in addition to T R . Freeze-in during pre-EMD is most efficient at the highest temperature in that era, which explains the appearance of T reh in the first expression in Eq. (7). Since DM particles start in chemical equilibrium in the decoupling regime, the second expression in (7) is independent from T reh . The factor T R /T 0 appears in both cases indicating dilution by entropy generation during the EMD epoch.
• The pre-EMD component is proportional to m χ . This can be understood by noting that the bulk of DM production in the pre-EMD phase occurs when T m χ . Therefore, the resulting DM number density, see Eqs. (19,22), is independent of m χ .
• The pre-EMD component has no explicit dependence on σ ann v f in the early-equilibrium regime, the second expression in Eq. (7), because DM particles start in chemical equilibrium. The only role of σ ann v f in this case is to determine the decoupling temperature T dec , which results in an implicit dependence through g * ,dec .
• The pre-EMD component has a much milder dependence on T R and m χ than the EMD component. As a result, moderate changes in these parameters can render the latter totally negligible, significantly affecting the allowed parameter space. This will become clear when we present our results in the next Section.

III. RESULTS
In this Section, we present our results. We numerically solve the set of Boltzmann equations governing the evolution of radiation, r, the decaying component driving EMD, φ, and DM particles, χ: where Γ φ is the decay rate of φ, E χ ≈ m 2 χ + 9T 2 1/2 is the average energy per DM particle 5 , and n χ,eq denotes the thermal equilibrium value of the DM number density. The energy density in φ is tracked until it is sufficiently small to be unimportant for the subsequent evolution, and is then dropped to facilitate faster numerical calculation. We have taken the detailed temperature dependence of the g * factor into account down to well below T R . In order to calculate the DM relic abundance, we normalize the DM number density with the entropy density long after the end of the EMD epoch.
We investigate the parameter space in the m χ − σ ann v f plane that yields the correct DM abundance for various values of T R , T 0 , and T reh . In each case, we use the contribution to the relic abundance from the entropy-producing phase of EMD as a baseline for comparison.

FIG. 2:
Values of m χ and σ ann v f that yield the DM relic abundance of Ω χ h 2 ≈ 0.12, obtained by numerical solution of the Boltzmann equations (8). The solid curve corresponds to production during the entire thermal history shown in Fig. 1, while the dashed curve depicts the freeze-in side of the corresponding baseline curve from production in the entropy-generating phase of EMD alone. In region 1, the DM relic abundance is dominated by the late EMD contribution, while the early-equilibrium and decoupling regimes of the pre-EMD contribution dominate in regions 2 and 3, respectively. The transition to freeze-out of the baseline curve is seen to begin at the top-right corner.
In Fig. 2, we show the curve that represents points for which Ω χ h 2 tot = Ω χ h 2 EMD + Ω χ h 2 pre-EMD reproduces the observed DM abundance for T reh = 10 12 GeV, T 0 = 10 10 GeV, and T R = 10 GeV. The curve consists of three distinct regions, 1, 2, and 3, as follows: (1) Region 1 starts at small DM masses and initially follows the baseline curve, but moves above it as m χ and σ ann v f increase. This behavior can be understood from the analytical approximations in Eqs. (2,7). The EMD component dominates at small masses due to its scaling ∝ m −5 χ . As m χ and σ ann v f both increase along the baseline curve, the pre-EMD component (7) becomes more relevant. Obtaining the correct DM abundance then requires a larger m χ than the EMD component alone, and hence the Ω χ h 2 tot curve goes above the baseline. (2) Region 2 starts at the turning point (the abrupt departure from the baseline curve) and extends to very small values of σ ann v f . It is essentially horizontal for the following reason. In this region, the pre-EMD component dominates and σ ann v f is large enough that we are in the early-equilibrium regime. The relic abundance is then given by the second expression in Eq. (7), which is almost independent from σ ann v f . However, region 2 is not exactly horizontal as Ω χ h 2 EMD brings in a very mild m χ dependence that is too small to be noticeable in the figure. Figs. 3 and 4 depict the sensitivity of the allowed parameter space on the post-inflationary history. In Fig. 3, we show variation of the curve from Fig. 2 for different values of T 0 for fixed T R and T reh (left panel), and for different values of T R when T 0 and T reh are kept constant (right panel). Fig. 4, shows the change in the curve when T reh is varied for fixed T R and T 0 . We observe the following main features in the figures: • In Fig. 3, region 2 moves up with increasing T 0 and down with increasing T R , while there is no change when T reh varies in Fig. 4. Since pre-EMD production in the early-equilibrium regime dominates in this region, we expect m χ ∝ T 0 /T R , independent of T reh , which agrees with what we observe in the figures.
• The turning point is highly dependent on T R and T 0 , as seen in Fig. 3, but does not change with T reh in Fig. 4. Its position can be estimated by setting the sum of the EMD component (2) and the second expression in Eq. (7) equal to the observed DM abundance and finding the local maximum of σ ann v f in terms of m χ . As it turns out, σ ann v f ∝ T 5 0 T −12 R at the turning point. This is in agreement with the considerable horizontal movement (specially when T R changes) in Fig. 3. As expected, the vertical shift follows that of region 2 mentioned above.
• The left end of region 2, where it meets region 3, moves horizontally with changing T reh in Fig. 4, and vertically when T 0 and T R are varied in Fig. 3. This point divides the decoupling and early-equilibrium regimes of the pre-EMD contribution, where we have σ ann v f ∝ T −1 reh , which explains the horizontal movement. As expected, the dependence on T 0 and T R follows that of region 2.
A very important point is that pre-EMD production opens up vast regions of the parameter space at very small σ ann v f that are not allowed when the EMD component alone is considered. In fact, the freeze-in side of the baseline curve does not extend below a certain value of σ ann v f . This is because lowering σ ann v f results in a smaller m χ in order to obtain the correct relic abundance. However, once T f ∼ m χ /4 drops below T R , the contribution to the relic abundance from RD after EMD becomes important. This sets a lower bound on σ ann v f , for a given T R , beyond which the baseline curve does not extend. Nevertheless, the pre-EMD contribution can still dominate for such small m χ and σ ann v f , especially for combinations of the parameters that lower region 2, such as T reh = 10 12 GeV, T 0 = 10 8 GeV, and T R = 10 GeV.
For the values of T 0 and T R in Figs. 2-4, the pre-EMD component leads to a separation of the freeze-in and freezeout parts of the allowed parameter space by a horizontal gap. Though the freeze-out part is not affected as much as the freeze-in part, and generally lies to the right of the figures, we include a short segment at the top right corner of Fig. 2 for reference. Decreasing T R and/or increasing T 0 results in a growing overlap between region 1 and the baseline curve, and a significant movement of the turning point to the right, making region 2 larger. However, the turning point cannot go beyond the peak of the baseline curve because larger values of σ ann v f actually give rise to freeze-out. Similarly, the freeze-out side moves toward the left and up with decreasing (increasing) T R (T 0 ).

FIG. 5:
Variation of T R , with constant T reh and T 0 , for values that display the transition between the horizontal gap and vertical gap behavior described in the text. For the larger two values of T R , the solid curves are separated into a freeze-in part on the left and a freeze-out part on the right (which we do not show apart from the top right corner). For the two smaller values, the curves are split into a lower part that follows the baseline curve and an upper part that connects the freeze-in and freeze-out sides.
At some point, the freeze-in and freeze-out parts join and region 1 coincides with the freeze-in side of the baseline curve, while region 2 splits away. By further decreasing T R and/or increasing T 0 , the allowed parameter space is again divided into two disjointed parts that are now separated by a vertical gap. We clearly see this in Fig. 5 where the smaller two values of T R have fully split from the baseline curve. For these two, region 1 follows the full baseline curve, including the freeze-out side which is not shown. The upper segments include regions 2 and 3, discussed above, that make a smooth transition to the freeze-out regime rising up on the right side. In summary, at higher T R and/or lower T 0 , the effect of pre-EMD production of DM is to disconnect the allowed parameter space into the freeze-in and freeze-out parts, while for lower T R and/or higher T 0 , it gives rise to a new curve that smoothly interpolates between the freeze-in and freeze-out regimes but is situated at larger values of m χ compared to the baseline curve.

IV. EARLY-EQUILIBRIUM REGIME: A CLOSER LOOK
A remarkable feature observed in Figs. 2-5 is that the DM relic abundance is essentially independent from σ ann v f in large parts of the parameter space. This is due to the fact that DM particles start in chemical equilibrium during the pre-EMD phase for a broad range of σ ann v f . This range, corresponding to the early-equilibrium regime, spans over many orders of magnitude extending from a minimum value determined by T reh , see Eq. (20), to a maximuum value that depnds on T 0 and T R and can be as large as that at the peak of the baseline curve. This feature can be considered as the freeze-in analogue to the WIMP miracle in a complementary way: while the relic abundance mainly depends on σ ann v f for the latter, it is mostly dependent on m χ in this case.
The DM abundance in this case is given by the second expression in Eq. (7). This is at most equal to the observed relic abundance as, in general, there exist other sources that contribute to the total abundance. Most notably, DM particles may be directly produced in the decay of the field(s) driving the EMD epoch [21][22][23][24]. As a result, we find the following inequality in order to not overproduce DM: In Fig. 6, we show the contours corresponding to the maximum allowed value of m χ in the T 0 − T R plane. The thick line segments are from full numeircal calculations, while the thin lines represent the RH side of (9). In the latter, we have taken g * ,dec = 106.75, the value for the SM, which is a very good approximation for the range of T R and T 0 shown in the figure. The numerical and analytical results agree very well for the range of parameters chosen.

A. Connection to Observables
If the two sides of (9) can be connected to experiments, then the inequality can be considered as a consistency relation of the pre-EMD early-equilibrium regime. The LH side is the DM mass, which can in principle be measured at the LHC (or future colliders) through the standard missing energy signal for the range of m χ shown in Figs. 2-5. One may also impose an upper bound on the RH side, which is a direct measure of the duration of the EMD era, in the context of inflationary cosmology.
In the non-standard history we have considered here, the number of e-foldings of inflation between the time when cosmologically relevant perturbations, corresponding to the pivot scale k * = 0.05 Mpc −1 , left the horizon and the end of inflation can be written as [25,26]: where: Here, H inf is the Hubble rate duting inflation, r is the tensor-to-scalar ratio, and w reh represents the equation of state during reheating. Using the relation between H and T at the onset and the end of the EMD epoch, we can write: Theoretical arguments and numerical simulations suggest that generally 0 ≤ w reh ≤ 1/3 [27,28], which implies that ∆N reh ≥ 0. Combined with the experimental bound r < 1 [29], and the typical range of values for g * ,0 and g * ,R , we find: In general, N k * is related to the scalar spectral index n s that is constrained by the cosmic microwave background (CMB) experiments. In particular, there are simple relations between N k * and n s [30] in two important universality classes of single field models of inflation that include a large number of models compatible with the latest Planck results [29]. One can therefore use the experimental bounds on n s to impose a lower limit on N k * and, through (13), an upper limit on T 0 /T R [31].
In addition, one may use the spectrum of primordial gravity waves to further constrain the RH side of (9). The initial spectrum of gravitational waves produced during inflation depends on r and the tensor spectral index n T . However, their subsequent evolution depends on the post-inflationary thermal history and is in particular affected by an epoch of EMD [32,33] (also, see [34][35][36]). Tensor modes that enter the horizon during EMD, and the modes that are already at subhorizon scales, experience a suppression, compared to a standard thermal history, due to entropy generation in this epoch. As a result, the shape of the tensor spectrum is sensitive to the beginning and the end of the EMD phase, equivalently T 0 and T R . This is complementary to the information from the scalar spectral index, mentioned above, which only depends on the duration of the EMD period encoded in T 0 /T R . Therefore, if r is not much smaller than the current experimental bound r max 0.064 [29], a future detection by (or limits from) the gravitational wave detectors could further constrain the allowed regions in the T 0 − T R plane.
We note that because of the very small value of σ ann v f in the early-equilibrium regime, see Figs. 2-5, there is no realistic prospect for a detectable signal from indirect detection searches. For the same reason, one can expect that this regime will also escape direct detection. However, with the help of (9), a combination of collider experiments and cosmological observations could be used as an indirect test of this otherwise elusive scenario.

B. An Example
We now discuss a specific particle physics scenario of EMD to see how the inequality in Eq. (9) is translated into constraints on the underlying model parameters. As mentioned before, moduli that arise in string theory are natural candidates for driving a period of EMD [6]. Consider a modulus field φ with mass m φ . It has gravitationally suppressed coupling to other fields resulting in a decay width: where typically c ∼ O(0.1). The modulus φ gets displaced from the minimum of its potential during inflation, and starts oscillating about it when H m φ . Generic arguments based on effective field theory estimates [37][38][39][40], or explicit calculations [41], suggest that the initial amplitude of oscillations is φ 0 O(0.1M P ). This implies that φ oscillations, which behave like matter, dominate the universe after their start, hence H 0 m φ . Oscillations eventually decay when H Γ φ and establish a RD universe with temperature T R 6 . Since the universe is approximately RD at the onset and the end of the EMD phase, and after using Eq. (14) with c ∼ 0.1, the inequality in (9) can now be written as: m χ 120 g * ,dec g * ,R g * ,0 where the RH side is basically controlled by m φ . In Fig. 7, we show the allowed region of the m χ − m φ plane accordingly. The curve depicts the upper bound of (15), obtained from full numerical calculations. We have taken g * ,0 = 106.75 corresponding to the high values of T 0 . Including new degrees of freedom beyond the SM (for example g * ,0 = 228.75, as in its minimal supersymmetric extension) will shift the curve slightly. Note that the slope of the curve does change, as expected, due to the implicit dependence of g * ,R on m φ when T R drops below the electroweak scale.
FIG. 7: Upper bound on the DM mass m χ for a given modulus mass m φ in the case of modulus-driven EMD. The curve is obtained numerically from the early-equilibrium regime.

V. DISCUSSION AND CONCLUSION
We now turn to discussing how our results may be extended to more general situations that involve a non-standard thermal history.
Throughout the paper, we have considered a non-standard thermal history where the EMD epoch is preceded by a RD phase established at the end of inflationary reheating. However, for a very slowly decaying inflaton, the field(s) driving EMD may have comparable energy density to the inflaton before reheating completes. An example is a modulus-driven EMD scenario with H reh m φ 7 . In this case, the universe does not enter a truly RD phase at the end of inflationary reheating but rather a phase where the energy density of the radiation and matter components are roughly equal. One may approximate this case by taking T reh = T 0 in the non-standard thermal history considered here. We have checked that full numerical calculations are in very good agreement with this approximation. Such a scenario is also of the "two-field" type, discussed in [19], where the two matter-components have equal initial energy densities. We note that the freeze-in and freeze-out sides of the curves in Figs. 2-5 are thus expected to merge at very high DM masses, m χ > T 0 , mimicking the shape of the baseline curve, though offset to smaller σ ann v f .
A general post-inflationary thermal history may include multiple epochs of EMD with respective parameters T 0,i and T R,i , 1 ≤ i ≤ n, separated by intermediate RD phases. Because of the temperature dependence discussed in the Appendix, the EMD component of the DM relic abundance is mainly due to production in the last (i.e., n-th) bout of EMD. This implies that T R must be replaced by T R,n in Eq. (2). The pre-EMD component, in the decoupling regime, is dominated by production at temperature T reh , while in the early-equilibrium regime, it is set at temperature T dec when chemical decoupling of DM particles occurs. Hence, the main modification to the expressions in (7) is replacing T R /T 0 with the product of T R,i /T 0,i to take all (relevant) EMD periods into account.
As mentioned at the beginning of Section II, we have assumed σ ann v f to be constant within the temperature range of interest. In general, however, σ ann v f can have a temperature dependence. Then, any allowed σ ann v f in Figs. 2-5 must be considered as the value of σ ann v f (T ) at the relevant temperature. I.E., T ∼ m χ /4 for region 1, T T dec for region 2, and T T reh for region 3 in those figures. In cases where σ ann v f (T ) ∝ T n with n > 0 (such as models studied in [45,46]), the pre-EMD production of DM particles will be enhanced compared to our calculations. The pre-EMD component of the relic abundance can be larger in this case thereby affecting the allowed parameter space even more prominently than shown in Figs. 2-5. In conclusion, early stages in non-standard thermal histories that include a period of EMD can significantly affect non-thermal production of DM via freeze-in for weak scale DM masses. We have demonstrated this in a postinflationary history involving an early RD phase followed by an EMD epoch that reheats the universe to a temperature T R before BBN. In the case where σ ann v f is constant over the temperature range of interest, the pre-EMD component of the DM relic abundance depends on the temperature at the onset of the EMD phase T 0 and the reheating temperature after inflation T reh , in addition to T R . This opens up vast regions of the m χ − σ ann v f plane where pre-EMD production can totally dominate the relic abundance as shown in Figs. 2-5. Moreover, DM particles reach chemical equilibrium in the pre-EMD era for a very broad range of σ ann v f spanning over many decades. The relic abundance in this case is virtually independent of σ ann v f , and avoiding DM overproduction yields an inequality between m χ and T 0 /T R . This brings an interesting possibility of combining collider searches (to measure m χ ) with CMB and gravitational wave detector experiments (to constrain T 0 /T R ) to test an elusive scenario that escapes indirect and direct detection. We leave a detailed study of this for a future investigation.