Dark matter from mediator decay in early matter domination

We study dark matter production from mediator decays in scenarios with an epoch of early matter domination. Particles that mediate interactions between dark matter and the standard model particles are kinematically accessible to the thermal bath as long as their mass is below the reheating temperature of the Universe after inflation. Decay of on-shell mediators can then lead to copious production of dark matter during early matter domination or a preceding radiation-dominated phase. In particular, for mediators that are charged under the standard model, it can exceed the standard freeze-in channel due to inverse annihilations at much lower temperatures (often by many orders of magnitude). The requirement to obtain the correct relic abundance severely constrains the parameter space for dark matter masses above a few TeV.


Introduction
While there are many lines of evidence for the existence of dark matter (DM) in the Universe [1], the identity of DM is an outstanding problem at the interface of particle physics and cosmology.Explaining the observed DM abundance is another challenge that in addition depends on the details of the thermal history of the early Universe.Thermal freeze-out in a radiation-dominated (RD) Universe is a simple and attractive mechanism that can yield the correct relic abundance if the (thermally averaged) DM annihilation rate takes the specific value ⟨σ ann v⟩ = 3 × 10 −26 cm 3 s −1 .However, in nonstandard cosmological histories, the correct DM abundance can be obtained for much larger or smaller values of ⟨σ ann v⟩ [2,3].
It is known that well-motivated classes of models arising from string theory generically lead to nonstandard histories that involve one or more epochs of early matter domination (EMD) [4].In general, an EMD phase arises when a matter-like component dominates the energy density of the Universe and eventually decays to establish a RD Universe prior to big bang nucleosynthesis (BBN).The matter equation of state can be due to coherent oscillations of a scalar field, like a string modulus, that is displaced from the minimum of its potential during inflation [5,6,7,8].It can also arise from long-lived nonrelativistic quanta produced in the postinflationary Universe [9,10,11,12].In fact, it has been recently shown that a visible-sector long-lived particle (LLP) with a weak-scale mass may achieve this [13,14].
The indirect DM searches have set stringent constraints on the DM annihilation rate, most notably the recent Fermi-LAT results from dwarf spheroidal galaxies [15] and Milky Way satellites [16].An analysis of these results has pushed the upper limit on ⟨σ ann v⟩ below the nominal value of 3 × 10 −26 cm 3 s −1 for DM masses below 20 GeV in a model-independent way [17] (larger masses can be excluded for specific annihilation channels).The increasingly tighter experimental bounds therefore warrant studying cases with a small DM annihilation rate.If ⟨σ ann v⟩ is very small, inverse annihilation will never be in thermal equilibrium and DM production occurs in the freeze-in regime.
Freeze-in production during EMD has been studied in the literature [18,19,20,21,22] (for a detailed review, see [23]).These studies have mostly focused on the latter part of an EMD epoch where evolution of radiation is nonadiabatic.It has been noted that earlier stages of the thermal history can significantly contribute to freeze-in production and thereby dominate the DM relic abundance (for example, see [24]).This is an important consideration because an epoch of EMD is typically preceded by other phases in the postinflationary history (notably the RD phase established at the end of inflationary reheating), especially for the last epoch of EMD that can end just before BBN.
A period of EMD starting with a non-negligible initial amount of radiation goes through an adiabatic phase during which T ∝ H 2/3 , as opposed to T ∝ H 1/4 in the nonadiabatic phase.In fact, as we will see, the Universe is generally in this adiabatic phase for most of the EMD period, in terms of the temperature evolution.The faster redshift of temperature in this phase then implies that particles that mediate DM interactions with the standard model (SM) particles can be kinematically accessible to the thermal bath during EMD (or in a preceding phase) even if they are very heavy.The mediator particles can indeed reach equilibrium if their couplings to the SM particles are not too small, and hence their decay can be an important source of DM production.
In this work, we present a detailed study of the contribution from on-shell mediator decays to the DM relic abundance in the freeze-in regime.By paying special attention to the adiabatic phase of EMD, we will show that these decays can easily dominate over the standard freeze-in production from inverse annihilations at much lower temperatures by many orders of magnitude.In fact, mediator decays alone can easily lead to DM overproduction in large parts of the parameter space.As we will see, the parameter space is tightly constrained for DM masses above a few TeV.The resulting constraints are milder and the allowed parameter space will extend to smaller DM masses if DM and the mediator belong to a hidden sector with its own gauge symmetry.
The rest of this paper is organized as follows.In Section 2, we introduce the general picture for mediator coupling to DM and SM particles.In Section 3, we discuss some important details of EMD epochs.In Section 4, we present our main results and identify the allowed region of the parameter space in order for mediator decays not to overproduce DM.We also comment on how our results are modified in the case of hidden sector DM models.We conclude the paper in Section 5. Some details of our calculations are discussed in the Appendix.

The set up
The main particles that are of interest in this work are the DM candidate χ and the mediator of its interactions with the SM particles X.Our focus is on DM production from X decay, and hence we consider dimension-four operators that are linear in X and involve χ and the SM particles (collectively denoted by ψ): Here, h i are dimensionless couplings.We assume that X is charged under the SM gauge group, while χ may or may not be a SM singlet.All operators in Eq. ( 1) are invariant under the SM gauge symmetry and other symmetries that a given model may possess (for example, the symmetry that is responsible for stability of DM).For clarity, we consider some specific forms of operators in Eq. ( 1) that arise in well-motivated particle physics models: 1 -X a scalar and χ a fermion.In this case, one can have the following types of interactions (using the four-component notation for fermions): where ψ are the SM fermions.Supersymmetric extensions of the SM with R-parity conservation give rise to the first term in Eq. ( 2), while R-parity implies a 2 = b 2 = a 3 = b 3 = 0.In these models, the lightest supersymmetric particle (LSP) is stable and hence a DM candidate (for a review, see [25]).Then X can be any scalar superpartner that is coupled to the neutralino DM.We note that both X and χ are charged under the SM in this case.
Another example is the model with GeV-scale DM proposed in [26], which gives rise to the first and third terms in (2) with a 1 = b 1 = a 3 = b 3 = 1/2.In this case X is a color-triplet and iso-singlet scalar that is coupled to the up-type and down-type RH quarks (as well as DM), and χ is a SM singlet fermion.
2 -X a gauge boson and χ a fermion.In this case, one can have interactions of the following type: h 1 X µ χ(a 1 +b 1 γ 5 )γ µ ψ+h.c., h 2 X µ χ(a 2 +b 2 γ 5 )γ µ χ+h.c., h 3 X µ ψ(a 3 +b 3 γ 5 )γ µ ψ+h.c., (3) with ψ denoting the SM fermions.This happens, for example, when the SM gauge symmetry is extended and X is a gauge boson of the new symmetries.An example is the U (1) B−L extension of the SM [27] (also see [28]) that leads to the second and third terms in (3) with b 2 = b 3 = 0.In this case, the Z ′ B−L is coupled to all SM fermions as well as the RH neutrinos, the lightest of which can be the DM candidate [29] 1 .
Eqs. (2,3) are meant as illustrative examples, rather than a full list, of dimension-four operators in Eq. (1).The important conclusion is that, when m χ ≪ m X , the first two operators in Eq. (1) generally result in the following partial decay width: where , and C 1 is a factor whose exact value depends on the nature of X and χ and the details of their couplings.In the following, we will take C 1 ∼ 1.Given that X is in equilibrium with the thermal bath at T ≥ m X , because of its SM charges, its decay can lead to copious production of DM even if h is very small.
The interactions in Eq. ( 1) also lead to annihilation/inverse annihilation χχ ↔ ψ ψ and scattering χψ ↔ χψ processes mediated by X.The former contributes to χ production from the thermal bath.However, X decay typically dominates χ production as long as the Universe starts at a temperature T ≫ m X so that on-shell X particles exist in sufficient abundance.This is because Γ ann ≪ Γ X→χ , where Γ ann = ⟨σ ann v⟩n χ .
We note that if Γ ann ≪ H(T ∼ m χ ), the comoving number density of DM particles produced from X decay will remain frozen.The annihilation rate at energies E ∼ m χ depends on the dimension of the effective operator χχψψ and its Lorentz structure.For example, assuming S-wave dominance, at E ≪ m X we have: where C 2 is a multiplicity factor.One final comment is in order.In general, nonrenormalizable operators mediated by X can significantly contribute to production of DM in the freeze-in regime.Such contributions, however, are typically subdominant to X decay due to the higher order of corresponding operators (see Appendix B).Our main goal here is to demonstrate that mediator decays set tight constraints on the parameter space in order to not overproduce DM.Including additional contributions can further restrict the allowed regions of the parameter space.

Early matter domination: generalities
In this work, we consider EMD that starts at some point after the completion of inflationary reheating.The matter component ϕ carries an energy density ρ ϕ , and its decay at the rate Γ ϕ feeds the radiation energy density ρ r .The evolution of ρ r and ρ ϕ is governed by the following system of Boltzmann equations: where the Hubble rate H is set by the total energy density ρ tot = ρ ϕ + ρ r .We take the moment when ρ ϕ = ρ tot /2 to be the onset of the EMD epoch.We denote the Hubble rate and the temperature at this time by H O and T O respectively.Given that ρ r = ρ tot /2, we have: where M P is the reduced Planck mass and g * counts the relativistic degrees of freedom at the time indicated by the subscript.Note that in the most extreme case we have T O ≃ T reh (where T reh is the inflationary reheat temperature).The EMD epoch that starts at H ≃ H O consists of two phases: Adiabatic phase -In this phase, the initial radiation dominates over that produced from ϕ decay.As a result, T ∝ a −1 , with a being the scale factor, similar to RD even though the Universe is in an EMD epoch where a ∝ t 2/3 .During this stage, the relation between H and T follows: This phase eventually ends when the temperature drops to the following value: Nonadiabatic phase -There is a transition to the nonadiabatic phase at T tr where radiation produced by ϕ decay becomes dominant over the preexisting amount.During this stage, we have: which implies that T ∝ a −3/8 .This stage extends all the way from T ≃ T tr to the end of the EMD epoch when H ∼ Γ ϕ and T ≃ T R .
This discussion holds even if some non-radiation component with an equation-of-state parameter 0 < w ≤ 1 makes the dominant contribution to ρ tot /2 at the onset of EMD, provided that it is not converted into radiation.Note that a component with w > 0 is redshifted faster than matter, and hence will become increasingly subdominant during EMD without any need to decay to radiation.If it feeds radiation, on the other hand, entropy will increase at a faster rate during EMD thereby leading to a shorter adiabatic phase and a higher T tr than what is mentioned above.
To illustrate the relative duration of the two stages of EMD, we consider two generic possibilities regarding the origin of such an epoch: coherent oscillations of a scalar field that is displaced from the minimum of its potential, and nonrelativistic particles produced from the thermal bath.
(1) Coherent oscillations -Consider a scalar field ϕ with mass m ϕ (a notable example is string moduli [4]).If m ϕ ≪ H inf , where H inf is the Hubble rate during inflation, then ϕ can be significantly displaced from the minimum of its low-energy potential during inflation.It typically starts to oscillate about this minimum when H ≃ m ϕ and its coherent oscillations behave like non-relativistic particles of mass m ϕ .If the Universe is RD at this time, and the initial fractional energy density of ϕ is α, coherent oscillations start to dominate when the Hubble expansion rate is: This signals the start of an EMD epoch when the temperature reaches the following value: On the other hand, if the Universe is still dominated by inflaton oscillations when H ∼ m ϕ , we will instead have: where H reh < m ϕ is the Hubble rate when a RD Universe is established after inflationary reheating.Therefore, given that α ≲ 12 , the temperature at the onset of the EMD epoch follows: (2) Nonrelativistic quanta -Consider a (bosonic or fermionic) particle ϕ from decay or scattering processes in the thermal bath.The produced quanta become nonrelativistic at T ∼ m ϕ .If the fractional energy density of ϕ at this time is α, its quanta dominate the Universe when the Hubble expansion rate is: We note that α ≲ g −1 * m ϕ , with the maximum occurring when ϕ reaches thermal equilibrium.Thus, in this case, we have: To get a feeling of the relative duration of the adiabatic and nonadiabatic phases, let us consider the m ϕ ≃ 10 2 − 10 6 GeV mass range.The upper and lower limits correspond to the cases when ϕ is a string modulus [6,7] and a weak-scale visible sector particle [13,14] respectively.Eqs.(15,17) then give us: where we have used g * ∼ O(100) as in the case of the SM.We therefore see that T O can be within a vast range spanning many orders of magnitude.The important point is that T tr is in general not extremely larger than T R .It follows from Eq. ( 10) that T tr /T R scales as (T O /T R ) 1/5 .For example, for T R ∼ 10 MeV, we find: We see that while T O is typically much higher in case 1, see Eq. ( 18), in both cases the majority of the EMD epoch is spent in the adiabatic phase T tr ≲ T ≲ T O .

Results
We consider the evolution of the χ number density produced from the decay of X in the bath, as described by the following Boltzmann equation [18,21,33,34]: where K i are the modified Bessel functions of the second kind.The superscript 'eq' denotes the corresponding equilibrium number density.Note that X reaches equilibrium with the thermal bath at T ≫ m X due to its SM charges.n X will follow its equilibrium value at T < m X if annihilation to the SM particles is efficient or in the presence of an efficient decay mode via the hXψψ term in Eq. ( 1).
In general, the χ number density will become frozen at some time in the cosmological history.Depending on the values of the parameters, this can occur in any of the phases described in the previous section.Furthermore, the freezing process may proceed as a freezein, where χ never achieves equilibrium with the bath, or a freeze-out, where the production rate is sufficient to reach equilibrium, with χ subsequently decoupling.
We can analytically estimate the abundance of χ produced from X decay, expressed as the number density normalized by the entropy density n χ /s, in the following way (see Appendix A for details).We first define the temperature of the background bath at the time when the χ number density freezes as T f .Beyond this point the number density is only redshifted with expansion as n χ ∝ a −3 .We express the number density n χ at this time as: where the function F (γ χ ) accounts for the possibility of not reaching equilibrium, with γ χ ≡ Γ X→χ /H(T = m X ).We will typically consider T f ≫ m χ such that the equilibrium number density above is given by the relativistic expression proportional to T 3 f .In general, the χ number density produced from X decay can become frozen at different stages of the cosmological history described in the previous section: (1) T R < m X ≲ T tr -In this case, n χ freezes during the nonadiabatic phase of EMD, and the final DM abundance is given by: where g χ is the number of degrees of freedom in χ.
(2) T tr < m X ≲ T O -In this case, n χ freezes during the adiabatic phase of EMD and the DM relic abundance follows: (3) T O < m X ≲ T reh -In this case, n χ freezes during the RD phase preceding the EMD epoch.Given that the frozen n χ and s are both redshifted ∝ a −3 during this RD phase and the adiabatic phase of EMD alike, we have: If γ χ ≪ 1, production of χ will not be strong enough to reach equilibrium and the number density of χ saturates to n χ ≈ 7γ χ n eq χ .The temperature at the time of χ freeze-in in this case falls roughly within the range m X /10 ≲ T f ≲ m X /5.On the other hand, if equilibrium is reached, we have γ χ ≳ 1 and F (γ χ ) = 1, with the temperature at χ freeze-out not being limited to the previous range.
In the case that χ reaches equilibrium in the nonadiabatic phase (or maintains it through the transition from the adiabatic phase), we can estimate the decoupling temperature using the RHS of the Boltzmann equation Eq. ( 20): where we have made use of n eq χ (T ) = (g χ /2π 2 )m 2 χ T K 2 (m χ /T ), and similarly for n eq X (T ) 3 .Using this to find T f , which we identify as the decoupling temperature in this case, we can obtain the abundance at the end of EMD (n χ /s) R for cases where χ remains relativistic and in equilibrium in the nonadiabatic phase at temperatures below m X .We note that if equilibrium is reached during the adiabatic phase, but equilibrium is not maintained through the transition (i.e., decoupling occurs before T tr ), then there is no dependence on the decoupling temperature, as seen in Eq. ( 23) above.
With T f determined, we can now compare the frozen χ abundance to the observed DM relic abundance in order to constrain such scenarios: Here, we only require that DM is not overproduced from X decays.Underproduction may be compensated for by other sources that create χ (for example, direct decay of ϕ or inverse annihilations).We numerically solve the Boltzmann equations for the background energy densities, given by Eq. ( 7) together with the equation for the DM number density produced from X decay, Eq. ( 20).We would like to emphasize that while Γ X→χ is essentially an input in (20), it can be calculated within the particle physics models mentioned in Section 2 according to Eq. ( 4).In particular, the model in [14] provides an explicit example where decay of a scalar X is the main source of GeV-scale fermionic DM and also produces long-lived particles that are responsible for a period of EMD [13].
In Figs. 1, 2, and 3 we show contours in the T O -m X plane that achieve the observed DM relic abundance for the values of the parameters shown in each figure.The region above and to the right of a given contour corresponds to underproduction of DM.In Fig. 1 we vary T R for fixed h and m χ .In Fig. 2 we vary the DM mass m χ while keeping h and T R fixed.In Fig. 3 we vary h for fixed values of T R and m χ .
The vertical dashed line denotes the current LHC limit on new particles charged under the SM, which we loosely take to be m X ≳ O(TeV) 4 .One could also constrain T O by considering the absolute upper bound T O < T reh where T reh ≲ 10 −1 (H inf M P ) 1/25 .The current Planck bounds [37] on the tensor-to-scalar ratio correspond to H inf ≲ 10 13 GeV, which results in a very weak limit T reh ≲ 10 15 GeV.Other considerations can significantly tighten this bound.For example, the recently proposed trans-Planckian censorship conjecture of the swampland program sets an upper limit of (H inf M P ) 1/2 ≲ 10 9 GeV in order for modes with sub-Planckian wavelength to not exit the horizon during inflation.This would result in T O < 10 9 GeV thereby removing the very top of the T O − m X plane.Model-dependent bounds on H inf lead to similar restrictions.
The diagonal dashed line shows where T O = m X .Extending the contours to the region below this line assumes that the Universe is in a RD phase from T = T O all the way up to T ≃ m X .This will not be the case if m X exceeds the reheating temperature after inflation T reh , or if there is another epoch of EMD at H > H O .Thus, the dependence on the details of the thermal history above T O requires the m X > T O half-plane to be treated with some care.
A general feature seen in all figures is that each contour starts as a vertical line at sufficiently small m X .Along the vertical segments the bulk of DM production takes place in the nonadiabatic phase corresponding to m X ≪ T tr where, see Eq. ( 22), the DM relic abundance is independent from T O .The vertical segment starts to turn at m X ∼ T tr .We note that in Fig. 2 the turning points are stacked on top of each other around values of m X within the same order of magnitude (unlike Fig. 1).This is because T tr is much more sensitive to T R than T O as can be seen in Eq. (10).Since T R is fixed in Fig. 2, the change in T tr for different contours is rather minimal, while this is not the case in Fig. 1.The more significant dependence of T tr on T R also explains the apparently missing vertical segment in the contours for T R = 10 MeV and 100 MeV in Fig. 1.In this case, T tr is smaller than the range of m X shown in that figure .Another visible feature in all figures is the presence of a plateau beyond m X ∼ T tr where the contours continue as a horizontal line for some range of m X .Along this segment, DM production from X decay occurs in the adiabatic phase.According to Eq. ( 23), the relic abundance has no dependence on m X as long as F (γ χ ) = 1 (happening for γ χ ≳ 1).Note that in the adiabatic phase, see Eq. ( 9), γ χ ∝ (m X T O ) −1/2 .This implies that for sufficiently small values of T O , we have γ χ ≳ 1 all the way up to m X = T O .Assuming a RD phase prior to EMD, γ χ ∝ m −1 X for m X > T O .It will therefore drop below one at some value of m X , beyond which DM production is not in equilibrium and the contour has a negative slope.This behavior is seen for contours with smaller values of T R (m χ ) in Fig. 1 (Fig. 2).Increasing T O results in a smaller value for γ χ during the EMD epoch, in which case DM production can drop out of equilibrium during the adiabatic phase of EMD.The contour Figure 1: The T O -m X plane for varied EMD end temperature T R with other parameters held fixed.Solid contours show parameter combinations that achieve the observed DM relic abundance.The region above (below) a given contour corresponds to underproduction (overproduction) and is therefore allowed (constrained).The diagonal dashed line shows T O = m X for reference, while the vertical dashed line sits at m X = 1 TeV.then starts to fall at m X < T O and keeps doing so for m X > T O (but with a different slope due to the different dependence of γ χ on m X ).This is the behavior seen for contours with larger values of T R (m χ ) in Fig. 1 (Fig. 2).
An additional feature seen in Fig. 3 is the asymptotic behavior of the contours as h increases.In fact, the flat horizontal segment of contours with sufficiently large h starts at the same T O and m X .This is because once γ χ ∼ O(1) at T tr , χ reaches thermal equilibrium during the adiabatic phase of EMD.Further increasing h only extends the plateau toward larger values of m X when γ χ falls below one again.
Given that the region below and to the left of each contour is ruled out as it gives rise to DM overproduction, we can draw some important conclusions about the allowed parts of the T O − m X plane: • The constraint m X ≳ O(TeV) does not severely restrict regions with DM underproduction unless T R ≪ 1 GeV and/or h ≪ 10 −5 .It mainly affects regions where DM is overproduced, which are already ruled out.
• The parameter space significantly opens up in the lower half-plane T O < m X .However, in this part the main contribution to the DM relic abundance comes from the pre-EMD phase that requires detailed knowledge of the postinflationary thermal history.
• Cases with m χ > O(TeV) are essentially ruled out unless T O and m X are (well) above 10 10  GeV (see Fig. 2).This is only marginally possible in models with large T O (like modulusdriven EMD), and not viable in models with intermediate values of T O , see Eq. (18).
We would like to note that throughout the allowed parameter space shown in Figs.1,2. the expressions given in Eqs.(5,6) result in ⟨σ ann v⟩ ≲ 10 −43 cm 3 s −1 .The most extreme case in Fig. 3, corresponding to h = 10 −3 and m X = 1 TeV, results in ⟨σ ann v⟩ ≃ 10 −37 cm 3 s −1 .This is still small enough to be in the freeze-in regime [20].It also results in an annihilation rate for χ that satisfies Γ ann ≪ H at temperatures T ≲ m χ , which ensures that the comoving number density of χ remains frozen upon production from X decay.In Appendix B, we discuss DM production from inverse annihilations and find that in large regions of the parameter space it can be completely neglected.In cases with m χ < T R and m X ≫ T tr , this contribution becomes significant and may even give rise to the correct abundance while X decays lead to DM underproduction.

Scenarios with hidden sector DM
In this section, we consider the case where X and χ are part of a hidden sector with its own gauge symmetry and temperature that can be distinct from those of the visible sector.The Boltzmann equations in this case are slightly modified to include the radiation energy density of the hidden sector: where superscripts denote the relevant sector and Br ϕ→VS is the branching fraction for ϕ decays to the visible sector.In order to guarantee a RD Universe that is compatible with observations, we need to have ρ VS r ≫ ρ HS r at the end of EMD.This implies that ϕ must decay predominantly to the visible sector, and hence Br ϕ→VS ≃ 1.
There are two major differences between this case and that considered in the previous sections: (1) m X does not need to be restricted to being larger than O(TeV).This is because X is a hidden-sector particle, and hence it can easily evade the LHC bounds.All that is required for our analysis to be valid in this case is that m X ≫ m χ .
(2) Given that ϕ decay must mainly reheat the visible sector, the amount of radiation injected to the hidden sector does not necessarily dominate over the initial radiation therein.As a result, the nonadiabatic phase of EMD may be shortened or even non-existent as far as the hidden sector is concerned.
In order to solve the Boltzmann equations in Eq. ( 27), we need to set the initial values of ρ HS r and ρ VS r as well as the value of Br ϕ→VS .Here, we consider a case where radiation is almost entirely in the hidden sector prior to EMD.This, for example, is the case in the model of high scale SUSY discussed in [38,39].We also use the current observational limits from PLANCK 2018 on dark radiation in the form of ∆N eff using the full TT,TE,EE+lowE+lensing+BAO data [40] to set a lower bound on Br ϕ→VS .These conditions yield a maximally constraining case for the configuration of the HS in that it is initially dominant and reheated to the maximum amount allowed by observations.An initially subdominant HS, as well as a shorter or even absent nonadiabatic phase, results in weaker constrains in the T O -m X plane.
In Fig. 4 we show contours for different values of m χ that yield the correct DM relic abundance for the benchmark with T R = 1 GeV and h = 10 −5 .We note that the vertical segments of the contours have shifted to the left, compared with Fig. 1, or completely disappeared.This is a consequence of having Br ϕ→HS ≡ 1 − Br ϕ→VS ≪ 1.Combined with the removal of the m X ≳ O(TeV) constraint, this results in a significant opening up of the allowed parameter space for small values of m X in comparison with Fig. 1.

Conclusion
We presented a detailed study of mediator production and decay and their contribution to the DM relic abundance in scenarios with an epoch of EMD.Special attention was paid to mediators that are charged under the SM gauge group and decay during early stages of EMD (or before its onset).
We showed that decay of on-shell mediators can totally dominate over the standard freezein contribution from inverse annihilations that happen at much lower temperatures.The requirement of not overproducing DM then leads to stringent constraints on the parameter space, demonstrated in Figs.1-3, for DM masses above a few TeV.Additional information on the scale of inflation H inf , or an upper bound on it, will further restrict the allowed parameter space in favor of small DM masses.
The situation will be more relaxed if the DM candidate and mediators belong to a hidden sector with its own gauge symmetries.The absence of collider bounds on the mediator mass combined with a subdominant reheating of the hidden sector (as required by observational bounds) opens up new regions of the parameter space in this case as seen in Fig. 4.
While mentioning some specific examples of cosmological histories with EMD and DM candidates, we remained largely agnostic about the particle physics origin of DM and the underlying model governing the cosmological history in this work.Our results can therefore be used by both phenomenologists and early Universe model builders to constrain DM models for a given thermal history and vice versa.They can also be used to pinpoint the parameters

A DM abundance from X decay
To estimate the final abundance of χ particles produced by decay of X during EMD, we consider the redshifted number density of χ from the time that its comoving value becomes frozen, marked by T = T f , to the end of EMD when T = T R , after which point entropy injection ceases: with assuming T ≫ m χ .Taking T R to be when ρ ϕ = ρ r at the end of EMD, to better match numerical results, we have (1) T tr < T f < T O -In this case the comoving χ number density freezes during the adiabatic phase and the Hubble rate at that time is therefore given by: where we have used a factor of 2 for ρ ϕ (T O ) = ρ r (T O ), and g * T 3 ∝ a −3 .After normalizing by the entropy density at T R , this gives: (2) T R < T f < T tr -The number density freezes during the nonadiabatic phase and the corresponding Hubble rate is given by: where we have taken Γ ϕ = √ 3H R to match our numerical results.After normalization, this gives:

B Production from inverse annihilations
Freeze-in production of DM from inverse annihilations can happen in three different regimes: (1) m χ > T tr .In this regime, inverse annihilations during the adiabatic phase of the EMD epoch (or whatever phase precedes it) will be the dominant source.In fact, as shown in [24], the main contribution comes from the highest temperature at which the expressions

Figure 2 :
Figure 2: The T O -m X plane for varied DM mass m χ with other parameters held fixed.As in Fig. 1, solid contours correspond to the observed DM relic abundance, with the region above (below) a given contour being allowed (constrained).

Figure 3 :
Figure 3: The T O -m X plane for varied coupling h with other parameters held fixed.As in Fig. 1, solid contours correspond to the observed DM relic abundance, with the region above (below) a given contour being allowed (constrained).

Figure 4 :
Figure 4: The T O -m X plane for varied DM mass m χ , with other parameters held fixed, for the case that X and χ are in a hidden sector.The hidden sector is taken to be dominant in the RD period before EMD, and is partially reheated during EMD to the amount allowed by PLANCK 2018 constraints on dark radiation from ∆N eff .