Dark Matter Freeze-In with a Heavy Mediator: Beyond the EFT Approach

We study dark matter freeze-in scenarios where the mass of the mediator particle that couples dark matter to the Standard Model is larger than the reheat temperature, TRH, in the early Universe. In such setups, the standard approach is to work with an effective field theory (EFT) where the mediator is integrated out. We examine the validity of this approach in various generic s- and t-channel mediator frameworks. We find that the EFT approach breaks down when the mediator mass is between one to two orders of magnitude larger than TRH due to various effects such as s-channel resonance, a small thermally-suppressed abundance of the mediator, or decays of Standard Model particles through loops induced by the mediator. This highlights the necessity of including these contributions in such dark matter freeze-in studies. We also discuss the collider phenomenology of the heavy mediators, which is qualitatively different from standard freeze-in scenarios. We highlight that, due to the low TRH, the Standard Model-dark matter coupling in these scenarios can be relatively larger than in standard freeze-in scenarios, improving the testability prospects of these setups.


Introduction and Motivation
The freeze-in mechanism [1,2] has been extensively studied in the literature as a viable production mechanism for the observed relic abundance of dark matter (DM). In contrast to thermal freeze-out scenarios, the freeze-in mechanism is characterized by couplings between the dark matter particle and the Standard Model (SM) thermal bath that are so feeble that the two populations never thermalize. The dark matter abundance is instead built up gradually over the cosmological history through these feeble interactions. Such feeble couplings suppress most experimental probes of dark matter such as indirect or direct detection. The most promising detection avenues consist of producing dark matter parent particles at high energy colliders and observing their displaced decays, see, e.g. [3][4][5][6][7][8].
Dark matter freeze-in scenarios can broadly be classified into two categories. If dark matter is produced via renormalizable interactions or decays of heavier particles, production dominantly occurs at temperatures close to the mass of the decaying or annihilating particles (in such cases, obtaining the correct relic density requires extremely small couplings); this scenario is dubbed infrared (IR) freeze-in. In contrast, if higher dimensional operators are involved, dark matter production primarily occurs at the highest temperatures, leading to scenarios referred to as ultraviolet (UV) freeze-in [9], where the relic abundance is sensitive to the highest temperature reached by the thermal bath, the reheat temperature, T RH , at which the radiation dominated evolution of the Universe commences after inflation. Although this dependence on UV physics is unattractive, such scenarios are inevitable in several beyond the Standard Model (BSM) theories. If these mediator particles are heavier than T RH 1 , integrating them out gives rise to higher dimensional operators in an effective field theory (EFT) relevant for early Universe cosmology, leading to UV rather than IR freeze-in of dark matter [10]. In such cases, the feeble interactions required to produce the correct dark matter relic density are a natural consequence of integrating out the heavy mediators of mass M , which leads to the effective couplings getting suppressed by powers of T RH /M . In this paper, we explore various cosmological and collider aspects of UV freeze-in scenarios in the presence of heavy mediator particles with M med > T RH . Our focus is to study the limitations of the EFT obtained by integrating out the heavy mediators, and how this affects the relic density calculation. We find that considerations of the existence of the heavy mediators can give rise to important effects not present in the EFT treatment, such as enhanced dark matter production from resonant effects, modification of the DM momentum distribution, and loop induced decays of SM particles into DM. In this paper, we will explore such scenarios in frameworks with s-as well as t-channel mediators. For s-channel mediator scenarios, we will consider scalar and vector mediators that mix with the SM Higgs and Z boson respectively, as well as a heavy Z that couples directly to both SM and dark matter. For the t-channel mediator scenario, we will consider a new heavy scalar that couples to SM fermions and the DM candidate. Several papers in the literature have studied DM freeze-in with s-channel scalar [11][12][13][14][15][16][17][18], vector [19][20][21], as well as spin-2 [22] mediators. Likewise, tchannel mediators have been discussed in [5,23]. For more general studies of freeze-in via portal mediators, also see [24][25][26]. Most of these studies consider reheat temperatures above the mediator mass, so that the mediator is part of the thermal bath. Our paper, which studies the opposite regime, is therefore complementary to these studies. It is interesting to note that our scenario interpolates between the standard UV freeze-in scenario (T RH M med ), and the standard IR freeze-in scenario (T RH > M med ).
In this paper, we also study the collider phenomenology of such mediators. In general, collider-accessible particles that can produce the correct dark matter abundance via IR freezein have very long lifetimes due to the feeble couplings involved, resulting in decays outside collider detectors. Therefore, IR freeze-in scenarios generically require non-standard cosmological histories in order to obtain mediators with modified lifetimes that can be directly tested at colliders, as well as larger production cross sections [3][4][5][6][7][8]. This does not apply to the frameworks of interest to us, which feature much larger couplings of the heavy mediator with SM and dark matter particles. Hence we find that collider signatures of such setup are qualitatively different from those expected from standard freeze-in scenarios.
The paper is organized as follows: In Section 2, we provide an overview of the simplified models and cosmological history we base our study on. Cosmological aspects of dark matter freeze-in for the cases of s-and t-channel mediators are discussed in Sections 3 and 4, respectively. Section 5 is devoted to a discussion of various phenomenological aspects of such heavy mediator frameworks. We end with a summary of our main findings in Section 6. A brief discussion of the impact of the epoch before radiation domination on dark matter production is presented in Appendix A. Finally, in Appendix B, we present some specific model realizations of our s-and t-channel mediated freeze-in scenarios.

Simplified Models
We perform our analyses in the framework of simplified models. We consider a Dirac fermion dark matter particle, X, with mass m X , that is stable and singlet under the SM gauge group. X interacts with the SM fermions through a heavy mediator. We choose SM fermions rather than gauge bosons since they are lighter and therefore more abundant in the thermal bath at low temperatures. Integrating out the mediator therefore gives rise to an EFT with fourfermion interactions between a pair of DM particles and a pair of SM fermions, ff ↔XX, which will produce DM via UV freeze-in. We consider this setup under two broad categories: an s-channel (scalar or vector) mediator (Sec. 2.1) and a t-channel mediator (Sec. 2.2).

Scalar mediator
We consider a Higgs-portal model, consisting of a SM-singlet real scalar mediator,Ŝ, with the following interactions: with v = 246 GeV, the SM-like Higgs h and the mass eigenstate S can be written as where m S and m h are the S and h masses, respectively. Note that the mixing angle θ h can be chosen relatively independently of the value of m S . This mixing gives rise to ff ↔XX interactions mediated by both S and h in the s-channel. The couplings involved in these interactions are given by where y f = m f /v is the SM fermion Yukawa. The DM mass can arise from the Yukawa interaction with S, from additional interactions with a broader range of dark sector particles, or be vector-like. We will therefore treat it as a free parameter in our numerical investigations.

Vector mediator Kinetic Mixing
A second possibility for the s-channel mediator is a dark U (1) gauge boson (the dark photon [27]),Ẑ , that mixes with the SM hypercharge gauge boson,B. The corresponding Lagrangian is given by where D µ = ∂ µ − ig D q DẐ µ , with g D and q D the dark U (1) gauge coupling and charge of the X particle, respectively. TheẐ mass, mẐ , can come either from the interaction with a dark Higgs or from the Stueckelberg mechanism [28]. The kinetic mixing parameter, , induces a mixing of theẐ with the SM hypercharge gauge boson. In particular, if we define Z 0 = √ 1 − 2Ẑ , the dark Z and the SM Z mass eigenstates are given by where θ W is the Weinberg angle,Ẑ is the would-be SM Z boson, and the last approximated expression only holds for 1 and mẐ m Z . In the same limit, the mass of the physical Z is mẐ (corrections arise at the 2 order). The mixing in Eq. (2.5) is responsible for the Z coupling with SM fermions, and for the SM Z boson with DM: where the couplings are given by T 3 and Y are, respectively, the third component of the isospin and the hypercharge of the SM fermion, where we use the convention Q = T 3 + Y , with Q the fermion electric charge. The corresponding couplings of the Z boson (g L , g R , g X ) are obtained via the exchange − sin α → cos α and cos α → sin α. These couplings give rise to the ff ↔ XX interactions mediated by both the Z and the SM Z boson.
In both of the simplified models discussed above, the s-channel mediator obtains couplings to SM fermions via its mixing with a SM particle (h or Z). As a consequence, both the mediator and the particle it mixes with can mediate interactions between SM fermions and dark matter. Here we instead consider a setup where the mediator couples directly to SM fermions without necessarily mixing with any SM particle. Such a scenario is realized, for instance, in models obtained by gauging one of the anomaly-free global symmetries of the SM, such as L µ − L τ or B − L. If DM is charged under this symmetry, it can interact with SM fermions via the new Z without the SM Z boson featuring in the process.
The minimal L µ − L τ DM model is defined by the following Lagrangian: 8) where q and q χ are free parameters quantifying the charge of the SM leptons and DM under the new U (1) µ−τ gauge symmetry, and g the gauge coupling strength of U (1) µ−τ . In this model, a kinetic mixing between the new Z and the SM Z boson is induced via loops of SM taus and muons, which are charged under both U (1) µ−τ and the hypercharge U (1) Y . The kinetic mixing coefficient is given by where g 1 is the hypercharge gauge coupling. We will see later that this tiny mixing is inconsequential for dark matter production and phenomenology.

T-Channel Mediator
Another possibility is the existence of a t-channel mediator. In this paper, we focus on a scalar mediator, S T , that interacts as where f is a SM fermion. We assume that the scalar potential is such that S T does not get a VEV, and therefore does not mix with the SM Higgs boson. The above coupling gives rise to the ff ↔XX interactions via S T in the t-channel. Such a setup can be realized, e.g. in supersymmetric frameworks, where S T is identified as a sfermion and X as the bino or axino (see e.g. [23]). Note that S T cannot be a SM singlet but carries the same charges asf . For example, if f is a SU (2) L doublet, S T must likewise be a doublet with multiple degrees of freedom. While this would give rise to a greater variety of collider signatures, in this paper we focus on the simpler case where f is a right-handed fermion, so that only one, SU (2) L singlet, mediator is involved. In particular, we will discuss two cases, (1) f = e R , and (2) f = t R , for which the cosmology is qualitatively different. In both cases, S T couples to the SM γ, Z bosons; in the latter case, S T also couples to gluons.

General Features of the Cosmological History
In the remainder of the paper, we will explore in detail various aspects of dark matter freezein in each of the scenarios outlined above. In this section, we provide a brief overview before delving into the details in the following sections.
We assume that the early Universe after the end of inflation is radiation dominated, and denote the reheat temperature at the beginning of this era as T RH . We assume that the epoch before radiation domination, when the Universe is dominated by the energy density of the decaying inflaton field, contributes negligibly to the dark matter abundance (the conditions for the validity of this assumption are discussed in Appendix A). Therefore, in all scenarios we are interested in, all of the dark matter is produced by freeze-in processes at temperatures below T RH .
The process common to all (s-or t-channel) scenarios is dark matter production from the annihilation of SM fermions, ff → XX. The dark matter yield from this process can be calculated as [19] , M Pl is the Planck mass, dΩ is the integral over the solid angle, s is the Mandelstam variable, K 1 (z) is the 1st-order modified Bessel function of the second kind, g ρ * and g S * are the effective numbers of degrees of freedom of the thermal bath for the energy and entropy densities, respectively, and M is the spinaveraged matrix element for the process. This equation is applicable to the case of s-or t-channel mediators, as well as to the higher dimensional, effective four-fermion interactions, with appropriate specifications of the matrix element M.
In addition to this fermion annihilation process, scenario-specific annihilation and decay processes not captured in the EFT, such as annihilation processes that produce one or two mediator particles in the final state, can also contribute to the dark matter abundance.
Here we comment on the production of dark matter through the decay of a mediator, concluding that this process should not be added to the aforementioned fermion annihilation process. In the s-channel case, when the mediator mixes with a SM particle (such as the Higgs or Z bosons), decays of the SM particle population in the bath can in principle also contribute to the dark matter abundance. The contribution from the mediator decay can be estimated as Here, g SM , m SM are the number of degrees of freedom and mass of the SM mediator, respectively. The second approximation holds for m SM > T RH , with Γ [x, y] the incomplete Gamma function. If T RH m SM , this abundance is Boltzmann suppressed, as encapsulated in K 1 (x). We can have m DM > T RH as long as kinematically accessible in the decay. Likewise, the decay of the heavy (non-SM) mediator particle also contributes, with the contribution given by Eq. (2.12) with appropriate replacements. For all mediator and SM decay contributions, the width of the decaying particle is greater than the Hubble rate for all T RH we consider in this paper. Thus these particles should be thought of as resonances rather than long-lived particles during the cosmological epoch of interest for dark matter production, which nevertheless maintain an "equilibrium" thermal abundance due to rapid inverse decays of SM states (for s-channel mediators), justifying the use of Eq. (2.12) for such scenarios.
Note, however, that the s-channel 2 → 2 process ff → XX contains a resonance regime, where the center of mass energy of the incoming particles matches the mass of the mediator, resulting in an enhancement of the cross section. This resonance regime corresponds to the intermediate particle being produced on shell, then decaying into dark matter particles. This is precisely the decay contribution given by Eq. (2.12). Therefore, including both contributions amounts to double counting. In such scenarios, it is therefore sufficient to only consider the annihilation contribution in Eq. (2.11), which includes both on-and off-shell contributions from the mediator (see e.g. [29,30] for related discussions). Nevertheless, in scenarios where the resonant regime dominates the production process, the decay contribution calculated from Eq. (2.12) will match the yield from Eq. (2.11). In this case, the decay contribution, which is simpler to calculate, can be used to estimate the dark matter abundance.

Dark Matter Freeze-In: S-Channel Mediator
In this section, we study freeze-in scenarios mediated by s-channel mediators, exploring the interplay between the various production channels in the dark scalar (Sec. 3.1) and dark vector (Sec. 3.2) cases.

Contributions to Dark Matter Freeze-In
The dark matter abundance receives contributions from several ff → XX processes. These freeze-in processes have been considered in a EFT framework after integrating out the heavy mediator [9]. While hh, SS → XX annihilations also contribute, these contributions are negligible for low reheat temperatures T RH m h as the DM yield is proportional to the square of the abundance of the Higgs or the scalar S, which are Boltzmann suppressed.
As we will show, among the SM fermion annihilations, the SM Higgs exchange processes dominate since we consider the regime m h m S . The identity of the SM fermion that dominates the process depends on m X and T RH . The fermions that couple appreciably to the Higgs and have (relatively) unsuppressed thermal abundance at T RH m h are the bottom and charm quarks, and τ leptons.
The relevant part of the Lagrangian, before integrating out the SM Higgs boson and scalar S, is given by where the two couplings can be expressed in terms of the Lagrangian parameters and Higgs mixing angle in Eq. (2.3) as y hXX = y S sin θ h , y hf f = −y f cos θ h , and analogously y SXX = −y S cos θ h , y Sf f = −y f sin θ h . In the EFT framework, we can compute the DM yield from freeze-in via a four-fermion dimension-6 operator 1 where N c is the number of colors of the fermion f , and we take g * S = g ρ * = 100 as the effective number of relativistic degrees of freedom around the GeV scale. Note that this formula assumes that the mass of the annihilating fermion is negligible compared to T RH .
We now turn to a treatment of this process that takes into account the physical nature of the mediator rather than treating the interaction as an EFT. The full matrix elements in the annihilation process are given by where  Dark matter differential yield from the annihilation processbb → XX as a function of √ s, as obtained using the approximate EFT formula (dashed curves) and the full calculation (solid curves) for two temperatures T = 5 GeV and T = 20 GeV. We show the SM Higgs (h) and heavy scalar (S) mediated contributions separately to highlight their relative sizes (the interference term is not shown); the total contribution (including interference) is shown with a black curve. For this plot, we use m S = 500 GeV, m X = 2 GeV, y hXX = 10 −6 , and sin θ h = 0.1.
Here Γ h (= 4.1 MeV) and Γ S are the widths of the SM Higgs and heavy scalar, respectively.
S is the width of the corresponding SM-like Higgs at the same mass, and Γ XX S the width into DM particles. Substituting these matrix elements into Eq. (2.11) enables us to calculate the DM abundance from these annihilation processes. Since the combination of couplings entering the Higgs and S contributions are the same (|y hXX y hf f | = |y SXX y Sf f |), and the mixing angle θ h is treated as a parameter independent of the mass m S , we expect the Higgs exchange to generically dominate over the S exchange or interference term. This is confirmed by Fig. 1.
In Fig. 1, the solid curves show the differential DM yield contributions from the two mediators, h and S, from thebb → XX process as a function of the center of mass energy at two different temperatures, T = 5 GeV (in blue and green) and T = 20 GeV (in red and orange). We choose a benchmark scenario with m S = 500 GeV, m X = 2 GeV, y hXX = 10 −6 , and sin θ h = 0.1. The Higgs exchange contribution (in blue and red) dominates at most energies, but at higher temperatures the S resonance can produce the dominant effect when √ s ∼ m S , as can be seen from the orange curve for T = 20 GeV. For comparison, we also show (as dashed curves) the contributions obtained by dropping the s dependence in the denominators of the matrix elements in Eq. (3.4), which leads to the EFT approximation in Eq.
S (Decay) Figure 2. Relative contributions of various annihilation channels, ff → XX, as a function of the reheat temperature for m X = 1 GeV, m S = 500 GeV, sin θ h = 0.1. For each channel, we plot the effective Yukawa coupling, y hXX , needed to obtain the measured relic abundance. Lower curves correspond to the more dominant channels. The dotted curves are obtained using the approximate EFT formulae. For illustrative purposes, we also include the contributions from h and S decays (black solid and dashed curves, respectively).
regions, which is a crucial aspect of considering the physical nature of the heavy mediator rather than treating the annihilation in an EFT framework. Such resonant contributions grow with T , as a larger part of the thermal distributions of the SM fermions can reach the resonant regime.

Cosmological History
We now study the interplay between various contributions in producing the measured dark matter relic abundance. In Fig. 2, we show the relative sizes of the several fermion annihilation channels for a representative choice of heavy scalar mass (m S = 500 GeV), heavy scalar-Higgs mixing (sin θ h = 0.1), and dark matter mass (m X = 1 GeV). The y-axis shows the size of the Higgs-DM coupling, y hXX , needed for a particular channel to fully provide the observed relic density of dark matter as a function of the reheat temperature. The lower a curve, the more efficient the channel is in producing dark matter. Therefore, the lowest curve represents the most dominant contribution. The solid and dashed curves show the full calculation of the annihilation contributions mediated by the Higgs and S, respectively. We have checked that interference effects are unimportant. For comparison, the dotted curves show the contributions as computed from the EFT. For illustrative purposes, we also include the contributions from h and S decays (black solid and dashed curves, respectively).
As we can see, the dominant contribution comes from bb → XX through the exchange of a Higgs (solid green curve). For very low reheat temperatures, T RH ∼ < 5 GeV, the EFT calculations match the full calculation (see solid vs. dashed green curves), hence the EFT language appropriately captures the dark matter production. In this regime, the correct dark matter abundance is obtained for y hXX ∼ 10 −7 − 10 −5 . On the other hand, for T RH 5 GeV, the curves depart from the EFT results as the resonant behavior due to the presence of the physical s-channel mediator (SM Higgs) becomes relevant (see also Fig. 1). In this region, the value of y hXX required to produce the observed DM relic abundance can be more than two orders of magnitude smaller than that predicted from the EFT treatment. In this regime, T RH is sufficiently high that production is dominated by the contribution from the resonance region √ s m h , which is effectively captured by calculating the contribution from h decay (solid black curve in the figure, as obtained from Eq. (2.12) with Γ SM→XX → Γ h→XX and m SM → m h ). Here, we find that couplings as small as y hXX ∼ 10 −11 (for T RH ∼ > 20 GeV) can produce the correct dark matter relic abundance. Such numbers are comparable to the feeble couplings associated with traditional IR dominated freeze-in scenarios. These patterns also hold for the contributions from the heavy Higgs exchange 3 . While the total contribution from S-mediated processes is always subdominant to that from the SM Higgs at these low T RH , it is interesting to note that tt → S * → XX is more efficient than tt → h * → XX (see the solid vs. dashed orange curves) as the latter does not get any resonant enhancement.
In Fig. 3, we show contours of the Higgs-DM coupling y hXX that produce the measured DM relic abundance as a function of T RH and m X . The sharp features at m X ∼ m h /2 are due to the fact that much larger couplings are needed at higher DM masses as the Higgs resonant enhancement is no longer possible. Relatively large Higgs-DM couplings are needed (y hXX ∼ O(0.01) and above 4 ) for m X O(10 GeV) and T RH ∼ O(GeV). We also show the boundary between regions where dark matter production is dominated by non-resonant annihilation (where the EFT approach provides a good approximation of the yield), and regions where production is dominated by resonant annihilation (where either the h or S decay approximation from Eq. (2.12) appropriately captures the dark matter yield). In the figure, we also show the bound from LHC searches for the Higgs decaying invisibly (red curve) and the direct detection bounds from XENON1T [31] (orange curve), see Secs. 5.1.1, 5.2 for more details. These searches already probe part of the parameter space of the model.

Vector Mediator
We now discuss dark matter freeze-in in the simplified model with a vector mediator, Z . We first discuss the case of a Z that mixes kinetically with the SM hypercharge, followed by the case of a Z arising from the gauged L µ − L τ symmetry. The results in this subsection will  be qualitatively similar to those in the previous (scalar) subsection, but with some crucial differences. In particular, in the kinetically mixed scenario, due to the dependence of the mixing angle on the ratio of the Z, Z masses (see Eq. (2.5)), we will find that the heavier mediator as well as the interference term play a more important role. In the L µ − L τ scenario, only the Z contributes to the dark matter abundance.

Contributions to Dark Matter Freeze-In
In the EFT framework, the relevant interactions are derived from the four-fermion dimension-6 operators 1 , with P L,R the left-handed and right-handed projection operators, respectively, and the coefficients Λ L,R given by with the relevant couplings defined in and below Eq. (2.7). The dark matter yield from these operators can be estimated as This estimate is similar to that from the scalar case (Eq. (3.2)), except for different prefactors due to a different Lorentz structure of the operators. The full matrix element for this annihilation process can be written as where sin θ W for 1 and m Z m Z , all of these squared matrix elements scale as ∼ 2 m 4 Z /m 4 Z and therefore are of comparable importance. Note that this is in contrast to the scalar mediator case, where the mixing angle sin θ h can be set independent of the heavy mediator mass, and we chose to fix it to a constant value (see the discussion below Eq. (3.4)).

Cosmological History
In Fig. 4, we plot the size of the coupling combination q D g D needed to produce the observed dark matter relic abundance as a function of T RH for various channels for m Z = 500 GeV, m X = 1 GeV. As in Fig. 2, the lowest curve represents the dominant production channel. The T RH /GeV ν Figure 4. Relative contributions of various channels for a Z mediator that kinetically mixes with the SM hypercharge, for m Z = 500 GeV, m X = 1 GeV. We include both the Z and Z mediated diagrams. We only show a subset of the EFT (non-resonant) contributions, for the bottom quark (dashed green) and the electron (dashed orange). For illustrative purposes, we also include the contribution from the Z decay. The decay contribution from a thermal Z population does not feature on this plot as it is independent of at leading order.
solid curves represent the result of the full calculation. For illustrative purposes, we show the EFT calculation for bb and eē annihilations only (dashed curves). At low T RH 7 GeV, nonresonant annihilation is dominant, with the largest contributions coming from the electron and the up quark 5 . At higher T RH 7 GeV, resonant effects start to become important, leading to departures from the EFT approximation. At these temperatures, the behavior is instead matched by the decay of a thermal Z population (blue curve obtained from Eq. (2.12)). Overall, couplings ∼ O(10 −7 ) are needed to produce the correct relic abundance from nonresonant annihilation at low T RH 7 GeV, whereas ∼ O(10 −9 ) couplings are sufficient for resonant annihilation at T RH 10 GeV.
In Fig. 5, we show contours of the values of g D q D needed to obtain the measured relic abundance as a function of T RH and the heavy mediator mass m Z for m X = 1 GeV. A large range of couplings, O(10 −11 − 10 −5 ) can give the correct dark matter abundance. We also show the boundary between regions where non-resonant annihilation dominates, so that the EFT approach gives a good approximation of the yield, and where resonant effects become dominant, and the full calculation must be performed. This boundary occurs at T RH ∼ 5 GeV and is essentially insensitive to the exact value of m Z . We also show bounds from LHC searches for a heavy Z (see caption of the figure and Sec. 5.1 for details) for two different  . Value of g D q D needed to obtain the measured relic abundance as a function of m Z and T RH for m X = 1 GeV. The white curve separates regions of parameter space where different contributions dominate dark matter production, as specified by the labels. We also show bounds from LHC searches for a heavy Z decaying into a lepton pair: the grey curve represents the CMS bound [32]; the black curves represent the stronger bound between the ATLAS [33] and CMS [34] searches. For the LHC bounds, we fix g D q D = 3 × 10 −6 (solid curves) and 3 × 10 −5 (dashed curve). sets of parameters, g D q D = 3 × 10 −6 (solid curves) and 3 × 10 −5 (dashed curve). Thus LHC constraints can be quite strong in the region of parameter space of interest if g D q D is small.

Modified Vector Mediator
In this section, we study dark matter freeze-in in the L µ − L τ model. The main difference between this framework and the kinetically mixed Z model is that the heavy Z mediator here couples directly to both dark matter and SM particles (muon and tau leptons and neutrinos) without requiring mixing with the SM Z boson. We have checked that processes mediated by the SM Z boson, which acquires a coupling to DM via loop-suppressed mixing effects (see Eq. (2.9)), are suppressed and negligible. Therefore, Z mediated processes dominate the dark matter production and phenomenology. 10 -4 ν μ/τ Figure 6. Value of the gauge coupling g needed to achieve the correct relic abundance of m X = 1 GeV dark matter from individual channels in the L µ − L τ model for q X = q l = 1 and m Z = 500 GeV. Solid curves represent the full calculation, the red dashed curve is the EFT (non-resonant) result for the muons and taus, and the blue curve denotes the contribution from Z decays. The Z decay contribution (green curve) is suppressed due to kinetic mixing only being induced at the loop-level.
The various contributions to the DM abundance are shown in Fig. 6, where we set q X = q l = 1 for simplicity. The solid curves denote the full calculation, while the dashed curve denotes the EFT (non-resonant) treatment (shown only for the muons and taus). The solid blue and green curves represent the contributions from the decays of thermal populations of Z and Z bosons, respectively. We see that contributions from Z decays are always subdominant due to the loop suppression of the kinetic mixing that gives rise to such decays. Otherwise, in line with previous results, non-resonant contributions dominate at low T RH 15 GeV, where the correct abundance is obtained for g ∼ 10 −4 . For higher T RH , the resonant behaviour is important, and the result is instead captured closely by considering decays of a thermally suppressed abundance of Z bosons, which gives the correct relic abundance for much smaller couplings. The couplings involved in producing the correct relic abundance in this model are larger than those involved in the kinetic mixing case due to the smaller couplings of the mediator with the SM particles as well as the absence of the Z-mediated interactions.

Salient Features
We now discuss some salient features common to all s-channel mediator frameworks.
As we saw in the previous subsections, the main feature is the s-channel resonance, which can dominate the dark matter production process and invalidate the EFT approach. The importance of this effect depends on the couplings involved as well as on the nature of the mediator (whether it couples to both dark matter and the SM directly or via mixing with some SM particle). In Fig. 7, we show this feature for the various mediators we have discussed. We plot the ratio Y f ull /Y non−res , where Y f ull is the dark matter relic abundance from the full calculation (including resonance effects), whereas Y non−res is the abundance obtained from the EFT approximations (obtained by dropping the s dependence in the denominators of the matrix elements). We use the values of the couplings that give the correct relic density for m X = 1 GeV with the full calculation. This ratio is plotted as a function of T RH /M med , where M med is the mass of the mediator that provides the largest contribution. From left to right, the curves are for the L µ − L τ model with M med = m Z = 0.1, 1, 10 TeV (red, blue, orange curves respectively), the scalar model (green curve, with M med = m h ) and Z kinetic mixing model (purple curve, with M med = m Z ).
At low values of T RH /M med , Y f ull /Y non−res ≈ 1 for all curves, showing that the EFT gives the correct dark matter relic abundance for sufficiently low T RH in all cases. As T RH /M med increases, an increasingly larger fraction of the thermal distribution of SM fermions can access the s-channel resonance regime, resulting in deviations from the EFT calculations. Recall that the matrix elements scale as |M| 2 ∼ 1/M 4 med in the non-resonant EFT limit but as |M| 2 ∼ 1/M 2 med Γ 2 med in the resonant regime. Therefore, the deviation from the EFT result is controlled by the magnitude of Γ med relative to M med : smaller widths, i.e. smaller values of Γ med /M med result in earlier deviations from the Y f ull /Y non−res ≈ 1 limit. This is indeed visible in the plot: the Z boson in the L µ − L τ model has tiny couplings to the dark and SM particles, therefore a very narrow width, and thus begins to deviate from the EFT calculation already at T RH /M Z ≈ 0.025. The scalar curve deviates next at T RH /M h ≈ 0.035, since the SM Higgs also has a relatively narrow width. Finally, the Z/Z curve deviates from the EFT result at T RH /M Z ≈ 0.05, since the Z width is larger. This implies that for dark matter freeze-in in s-channel models, the EFT approach already breaks down when the mediator mass is one or two orders of magnitude above the reheat temperature. The relative steepness of the curves is also determined by the size of the width relative to the mass of the mediator: the smaller the width, the faster the rate of departure from the EFT result. Finally, for the L µ −L τ model, we have shown three curves, for m Z = 0.1, 1, 10 TeV (red, blue, orange curves respectively). Consistent with the discussion above, a lighter Z departs from the EFT result at lower T RH /M Z since it requires smaller couplings to obtain the correct relic density, hence Γ med /M med is smaller.
The prominence of the s-channel resonance with heavy mediators significantly changes not only the relic abundance of dark matter but also its momentum distribution. This is illustrated in Fig. 8, where we plot the DM momentum distribution from SM Higgs decays (representative of resonant annihilation) 6 for various values of T RH = 5, 20, 40 GeV. In general, dark matter momentum distributions peak at p/T ∼ 1 since dark matter particles are produced, on average, with energy corresponding to the temperature of the thermal bath. However, when the process is dominated by resonant annihilation via the s-channel Higgs mediator, dark matter particles get produced dominantly at p ∼ m h /2. Since the temperature of the thermal bath is much lower, the DM momentum distribution peaks in this case at p/T ∼ m h /(2T RH ). This ratio gets further suppressed by an O(1) factor as entropy released from subsequent decoupling and decays of the SM particles heats the thermal bath but not the DM population. The various curves in the figure agree with these considerations. For comparison, we also show the momentum distribution expected in standard IR freeze-in scenarios (purple curve, adapted from Ref. [17]), where T RH M med . In this scenario, most of the dark matter production occurs at T M med , and consequently the IR freeze-in distribution is colder than the other distributions with smaller T RH , as well as broader with a larger lower momentum counterpart. Such modified momenta distributions provide another observable difference between standard IR freeze-in scenarios and freeze-in with a heavy mediator and a low reheat temperature.

Dark Matter Freeze-In: T-Channel Mediator
We now focus on dark matter freeze-in with a t-channel mediator. This scenario contains several qualitative differences from the s-channel mediator framework discussed in the previous section. In particular, the t-channel framework does not feature the resonant behavior of the annihilation processes that resulted in large deviations from the EFT approach in the schannel scenario. Instead, as we will see, deviations from the EFT arise due to on-shell production of the t-channel mediator through annihilation processes, as well as new loopinduced decays of SM particles. As mentioned earlier, we study the t-channel scenario in two cases: (1) mediator coupling to the right-handed electron, e R ; (2) mediator coupling to the right-handed top quark, t R .

Contributions to Dark Matter Freeze-In
The various processes contributing to dark matter production in the t-channel framework are shown schematically in Fig. 9. In the low energy EFT, DM freeze-in arises via four-fermion interactions obtained by integrating out the heavy mediator S T (diagram (a) ), analogous to the s-channel case. However, since DM is now no longer Z 2 -symmetric, and since S T carries SM charges, additional contributions exist and must be taken into account. Since the S T − X system carries an effective Z 2 -symmetry, there are three classes of annihilation processes: XX, X S † T , or S T S † T (top row in Fig. 9). Since S T →f X is the only decay channel available for S T at tree-level, each S T particle produced in the early Universe will result in a DM particle, hence the latter two diagrams contribute to secondary DM production. In addition to these annihilation processes, the mediator S T also gives rise to novel loop-level processes that produce DM via decays of SM particles. This is shown in the second row of Fig. 9. We now discuss these various contributions in more detail.

(a) ff → XX annihilation
As in the s-channel case, we can estimate the DM yield from freeze-in via a four-fermion dimension-6 operator 1 Λ 2 (fX)(Xf ) obtained from integrating out the t-channel mediator from assuming m f T RH . For the full calculation, we use the full amplitude for the process, where t is the standard Mandelstam variable. Note that there is no resonance effect for a t-channel mediator, as was the case for the s-channel mediator, as the propagator cannot go on-shell. Thus, Eq. (4.1) is expected to remain a good approximation of the annihilation contribution.
(b) (Z/γ/g/h) f → S † T X coannihilation Fermion coannihliation with an electrically neutral SM boson (Z/γ/g/h)f → S † T X proceeds via the two diagrams shown in Fig. 9 (b). In contrast to the ff → XX process, the cross section for this class of diagrams scales as y 2 T instead of y 4 T . However, since these processes involve the mediator being produced on-shell and we have assumed T RH < m S T , the probability for such interactions to occur is Boltzmann suppressed. Note that all initial state particles, including Z, h, t, can be approximated to be massless since their masses are negligible compared to the kinematic threshold required for this process.
These interactions, if sufficiently rapid, will produce a thermal abundance of dark matter. Hence we require that these interactions remain slower than Hubble at all times, n σv < H, where the sum is over all processes that contribute to the production of the XS † T final state. For the case where S T couples to t R , the dominant process is coannihilation with a gluon, g t R → S † T X (the h t L → S † T X contribution is of the same order of magnitude due to the large top Yukawa), which enforces the approximate condition where x = m S T /T RH , and we have ignored O(1) factors in the estimate. For the case where S T couples to e R , coannihilations with Z/γ are the most important: summing these contributions, the approximate condition to avoid thermalization is the same as in Eq. (4.3) with α s → α.
We focus on regions of parameter space where the non-thermalization conditions are satisfied, so that dark matter is produced from freeze-in. In this case, freeze-in formulae analogous to Eq. (2.11) (see e.g. discussions in [9]) can be used to calculate the freeze-in abundance. For instance, for the process h t L → S † T X, the first diagram in Fig. 9 For this process, g 2 T = 3(m t /v) 2 7 . For the other processes, both diagrams in Fig. 9 (b) contribute, and the DM yield does not admit a simple closed form as above, but they are straightforward to evaluate numerically.
The final diagram in the top row of Fig. 9 represents annihilations of SM particles that pair-produce the mediator. There are several classes of diagrams contributing with the most important being annihilations of gauge bosons. Such dimension-4 operators involve no small couplings and can therefore be efficient enough to thermalize the S T population if T RH is sufficiently high. This occurs if the reheat temperature is greater than the S T freeze-out temperature T f.o. , which can be very roughly estimated as T f.o. ≈ m S T /20 for weak scale interactions. In this regime, the dark matter abundance is a result of thermal freeze-out of the mediator, whose decays then populate dark matter, rather than freeze-in processes.
For simplicity, we will focus on the T RH < T f.o. ≈ m S T /20 regime, and assume that S T is never in equilibrium with the SM bath 8 but instead can be produced via freeze-in. Note that the DM abundance from this contribution is independent of y T . However, given that the S T abundance from freeze-in has to be smaller than its thermal freeze-out abundance, 7 In the case of the Higgs we do not include the contribution from the second diagram in Fig. 9 (b), as it depends on the coupling hST S † T , which is an independent parameter of the Lagrangian. We neglect this coupling, for simplicity. However, in general, this diagram will contribute to the dark matter abundance, especially at relatively sizable values of TRH . However, even in this regime, we do not expect a sizable change in the value of yT needed to achieve the correct relic abundance compared to what is shown in Fig. 10. 8 If ST interacts with gluons, as in the case where it couples to tR, thermal equilibrium is maintained to far lower temperatures. However, in this case the ST freeze-out abundance is also significantly smaller, and the subsequent contribution to dark matter abundance becomes negligible.
such freeze-in contributions are very small and generally do not contribute significantly to the dark matter relic abundance.

(d,e) Z, h decay
The SM Higgs and Z bosons can decay into a pair of dark matter particles through oneloop diagrams with f, S T in the loop, as shown in the bottom panel of Fig. 9. We compute the decay widths with the help of Package-X [35,36].
In the model in which S T couples to e R , the Z → XX decay width in the m e , m Z m S T limit can be written as (4.5) where g T = e tan θ W is the Z coupling to right-handed electrons, and we have dropped terms suppressed by powers of m 2 e /m 2 S T . In the model in which S T couples to t R , the corresponding width receives additional comparable contributions that consist of powers of m 2 t /m 2 S T instead of m 2 Z /m 2 S T and depend on the Z coupling to the SM left-handed top quark. We do not present the full expression for this decay width since it is quite lengthy 9 , but we use it in our numerical calculations.
We evaluate the Higgs diagram in Fig. 9 (d) to be for both models (f = e R , t R , with N c = 1, 3, respectively). The dark matter abundance from such decays can be estimated by substituting the above decay widths into Eq. (2.12). Note that while the Z decay width can be enhanced by large logarithms, the Higgs decay width is suppressed by powers of both dark matter and loop fermion masses, as the process requires a helicity flip since the mediator couples only to righthanded fermions. 10 Due to these additional suppression factors, as well as a smaller number density of Higgs bosons compared to Z bosons in the early Universe at T RH < m h , m Z , the contribution from Higgs decay is several orders of magnitude smaller than the contribution from Z decay in all cases we consider (see Fig. 10). 9 For Z decays involving top quark loops, the part of the amplitude leading to Eq. (4.5) The additional contributions are proportional to , where g L T is the SM Z coupling to tL. These terms lead to contributions comparable to the ones in Eq. (4.5). 10 For simplicity, for the Higgs decay, we ignore the possible coupling between the Higgs and ST of the form λvhST S † T , which would introduce a new diagram analogous to Fig. 9 (e). This contribution would scale as ∼ (λv 2 ) 2 instead of (m 2 f ) 2 in Eq. (4.6), and can be important. However, since the Higgs contribution to the relic abundance turns out to be subdominant in all cases (see Fig. 10), this additional contribution would not change our conclusions.  Figure 10. The size of the coupling y T needed to produce the correct relic density from various processes as a function of the reheat temperature, T RH , for S T coupling to the right-handed electron e R (top panel) and to the right-handed top-quark t R (bottom panel), for m X = 1 GeV and m S T = 1200 GeV. The dashed curve in the electron plot is the EFT approximation. In the shaded regions, coannihilation processes are rapid enough to produce a thermal DM population.

Cosmological History
We now turn to a discussion of the interplay between the above contributions in setting the correct dark matter relic density from freeze-in. Fig. 10 shows the relative contributions of the various channels as a function of T RH for m X = 1 GeV and m S T = 1200 GeV for scenarios where the mediator couples to e R (top panel) or t R (bottom panel). Here, we choose m S T = 1200 GeV for the mediator mass in view of strong constraints from LHC searches (see Sec. 5.1.3).
For the electron case (top panel of Fig. 10), the EFT approximation from Eq. (4.1) (dashed blue curve) closely matches the full calculation for e + e − → XX (solid blue curve) throughout, as no strong resonance features are present for a t-channel mediator. The contribution from loop-mediated Z decays (purple curve) is seen to be subdominant to this annihilation contribution throughout. Higgs decay contributions (orange curve) are even weaker and do not feature in this plot. At low T RH (< 60 GeV), fermion annihilation dominates dark matter production, giving the correct relic density for y T ∼ O(10 −4 ). At higher T RH , coannihilation with a SM boson, in this case γ and Z, become increasingly important, as an increasingly greater fraction of the thermal population gains enough energy to produce S T via (γ/Z)e → S † T X. These coannihilation processes overtake fermion annihilation as the dominant dark matter production channel around T RH ∼ 60 GeV. In the shaded region, the coannihilation processes thermalize the DM with the SM bath as determined by Eq. (4.3) with α s → α. As anticipated, the curves for the correct relic density from freeze-in lie away from this region.
Some of these features change when we consider the scenario where the mediator couples to t R (bottom panel of Fig. 10). Here, the fermion annihilation curve (solid blue) is significantly weaker (i.e. requires significantly larger couplings) than for the electron case, as the number density of top quarks in the thermal bath is severely Boltzmann suppressed at such low temperatures. 11 For this reason, Z loop decays, although suppressed by several factors, provide the dominant contribution to dark matter freeze-in for T RH 50 GeV, producing the correct relic density for y T ∼ 10 −4 − 10 −3 . The contribution from Higgs decays, enhanced by the large top Yukawa coupling, is also visible, but remains subdominant to the Z decay contribution. As in the e R case, coannihilation processes with SM bosons become dominant at larger T RH . In this case, coannihilation with a gluon, gt → S † T X (brown curve in the figure), dominates for T RH 50 GeV, and significantly smaller couplings ∼ 10 −5 can produce the correct relic abundance.
In Fig. 11, we show contour plots of the value of the coupling y T needed to obtain the correct relic abundance as a function of the mediator mass m S T and the reheat temperature T RH . As explained earlier, we terminate the x-axis at T RH /m S T = 1/20, beyond which the mediator is likely in thermal equilibrium with the SM bath. We also delineate regions of parameter space where different processes dominate dark matter production. For the case of S T coupling to e R (top panel), the annihilation process e + e − → XX dominates in large regions of parameter space, giving the desired relic abundance with y T = O(10 −4 ). This changes at T RH /M S T ≈ 0.05, beyond which the coannihilation process (γ/Z)e → S † T X dominates. This boundary is delineated by a red curve in the plot. Thus, the EFT approach where S T is integrated out breaks down at T RH /M S T ≈ 0.05.
The bottom panel of Fig. 11 shows the analogous contour plot for S T coupling to t R . As explained before, in most of the low T RH parameter space, Z decays provide the dominant contribution as tt annihilation is Boltzmann suppressed for T RH < m t . Since the Z-decay width into dark matter is parametrically suppressed by the heavy mediator mass as well as loop factors (see Eq.  Figure 11. Contours of the coupling y T needed to produce the correct relic density as a function of the reheat temperature, T RH , and mediator mass, m S T , for m X = 1 GeV, for scenarios where the mediator couples to the right-handed electron (top panel) or the right-handed top quark (bottom panel). Red curves separate regions where different processes dominate dark matter production, as denoted by the labels (see text for more details). The dashed grey curve in the bottom panel denotes the approximate lower bound on m S T from the LHC. coupling is > 1. As T RH increases, top-gluon coannihilation, gt → S † T X, grows to dominate, and smaller couplings y T ∼ 10 −4 are sufficient. In the plot, this occurs to the right of the red curve at T RH /m S T ≈ 0.04. For m S T 1200 GeV, we also see the emergence of a third region, where tt → XX dominates: in this region, T RH is sufficiently lower than m S T that the coannihilation processes involving on-shell production of S T are suppressed, but high enough that the Boltzmann suppression of the thermal abundance of top quarks is no longer too severe, so that dark matter production from tt annihilations are dominant. Hence all three processes -fermion annihilation, fermion-boson coannihilation, and boson decays -can be the leading dark matter production mode in this t-channel scenario.
It is instructive to compare the nature of EFT breakdown in the t-channel scenario with those from the s-channel framework (Fig. 7). In both cases, the EFT breaks down due to the emergence of processes where the mediator is produced on-shell. In the s-channel scenarios, the mediator is produced on resonance via the inverse decay process ff → h, S, Z, Z , whereas in the t-channel case it is a product of coannihilation between a SM fermion and a boson.
We end with a brief comment regarding the dark matter momentum distribution in the t-channel mediator scenario. The various contributing processes have distinct energy scales associated with them, and will therefore produce dark matter with different momenta. For the electron case, fermion annihilations produce DM with p ∼ T , whereas coannihilations produce DM with p ∼ 0 directly, as well as with p ∼ m S T /2 from the subsequent decays of S T . Likewise, the top quark case features annihilations (p ∼ m t ), coannihilations (p ∼ 0, m S T /2), as well as Z decays (p ∼ m Z /2). Thus, the dark matter momentum distribution carries imprints of the dominant production process. We leave a detailed investigation of such features for future study.

Phenomenology
Having discussed the early Universe history, we now turn to a discussion of the phenomenological aspects of the various frameworks discussed in this paper. It is well known that indirect detection of dark matter annihilations is extremely unlikely in freeze-in models due to the small effective couplings involved, which remains true in the setups we studied in this paper. However, while the effective couplings are small, the real couplings involved in SM-dark matter interactions, relevant for direct experimental probes such as direct detection and colliders, can be relatively sizable for T RH M med , improving detection prospects on these fronts. We discuss collider prospects below, and follow with a short discussion of direct detection.

Probing mediators at collider experiments
Collider phenomenology of standard freeze-in setups often involve displaced decays of the mediator particles due to the feeble couplings involved, see, e.g. [3][4][5][6][7]. For the setups we have considered, which can involve larger couplings, the signatures can be qualitatively different.

S-channel scalar mediator
The DM-Higgs coupling, y hXX , needed for DM freeze-in will induce an exotic decay channel of the SM Higgs into dark matter with the decay rate The combination of the most recent ATLAS searches for invisible Higgs decays, performed with 5, 20, 139 fb −1 of 7, 8, and 13 TeV data, sets a bound of BR(h → invisible) ∼ 0.11 at 95% C.L. [37]. Similarly, the latest CMS combination of Higgs invisible decay searches performed with 5, 20, 36 fb −1 of 7, 8, and 13 TeV data provides the bound BR(h → invisible) ∼ 0.19 at 95% C.L. [38]. The ATLAS bound translates into a bound on the Higgs coupling to DM y hXX 0.01 for m X m h . This bound is shown in red in Fig. 3, and is seen to probe our freeze-in scenario at very low reheat temperatures T RH ∼ GeV and for DM masses 10 GeV. Future projections show that a bound on the Higgs invisible branching ratio at the level of ∼ 2% can be achievable at the HL-LHC [39,40], which translates into y hXX 4 × 10 −3 , slightly extending the coverage in parameter space.
The phenomenology of the heavy scalar, S, depends not only on the coupling y hXX but also on its mixing with the SM Higgs boson, sin θ h . It can be singly or pair produced at the LHC via its Higgs portal coupling, with the production cross section given by the corresponding cross section for a SM Higgs with the same mass, suppressed by sin 2 θ h . Thus, non-negligible production at colliders requires this mixing to be sizable. For example, for the parameters used in Sec. 3.1, m S = 500 GeV and sin θ h = 0.1, the heavy scalar production cross section calculated at NNLO+NNLL is ∼ 45 fb [41]. This value of sin θ h will lead to prompt decays of S into SM fermions and gauge bosons with width Γ SM S ∼ 0.6 GeV. Its decay width into dark matter particles depends on y hXX and sin θ h . For the values needed to obtain the correct relic abundance via freeze-in, this width is generally negligible, except at very low values of T RH ∼ GeV and m X 10 GeV (in and around the white region in Fig. 3, which corresponds to y hXX > 1). LHC searches for new scalar resonances do not yet constrain this heavy scalar with m S = 500 GeV, but it can be probed at the HL-LHC via its decays to ZZ [39,42]. More exotic decays of S are possible if the dark sector contains additional structure. This is a model dependent question, and a specific example that realizes such possibilities is discussed in Appendix B.1.

Kinetic Mixing
If the Z − Z kinetic mixing is significant, the Z can be produced copiously at the LHC, and its decays can provide observable signatures. Depending on the parameters, the Z can decay dominantly into dark matter or to SM states.
The strongest LHC limits on Z resonances are derived from CMS and ATLAS searches for narrow dilepton resonances. The most important searches are: CMS search [32] for dimuon resonances in 110 GeV < m Z < 200 GeV; CMS search [34] and ATLAS search [33] for dilepton resonances at higher masses up to 6 TeV. Ref. [43] shows that in a kinetically mixed Z model where the Z decays 100% into SM states, the bounds are at the level of ∼ 10 −2 across the range of masses we consider in our paper. In Fig. 5, we set bounds on the Z parameter space for two different values of g D q D = 3 × 10 −6 , 3 × 10 −5 . For small values of g D q D , these bounds are seen to be significant for low T RH models.
Invisible decays of Z into dark matter, Z → XX, could be probed, e.g. by monojet searches [44,45]. We have checked that the monojet cross sections predicted from the parameter space in Fig. 5 are several orders of magnitude smaller than what is currently probed by LHC searches.
Indirect constraints from electroweak precision measurements [46,47] or measurements of Z invisible decay width (Γ inv Z = 499.0 ± 1.5 MeV [48]) are also very weak in the m Z 100 GeV mass range that we focus on in this paper.
The Z gauge boson arising from gauging L µ − L τ is only mildly constrained by collider data if m Z > m Z . The Z can be produced at the LHC through its coupling to muons, taus, and neutrinos via the processes pp → Z µ + µ − , pp → Z νν, and pp → Z µν µ , where the muons can be replaced with taus. However, so far, LHC searches have only been performed in the mass range (5 − 70) GeV, where the Z is produced from Z decay (Z → Z µ + µ − , Z → µ + µ − ) [49]. Additional bounds can be obtained recasting ATLAS and CMS multilepton analyses. This has been done in [50], showing that the Z masses up to 550 GeV are probed for g = O(1).
Additional constraints arise from high intensity experiments. The most stringent constraints come from the measurement of the neutrino trident process ν µ N → ν µ µ + µ − N [51,52] by the CCFR experiment [53], but this only constrains light Z masses. Finally, the L µ −L τ Z can address the (g − 2) µ anomaly. However, for couplings g O(1), this requires m Z 200 GeV [51].
All these bounds can in principle be affected by the mixing of the Z with the SM hypercharge gauge boson. In the L µ − L τ model, this mixing is generated at one loop (see Eq. (2.9)). For the values of g needed to produce the measured relic abundance, this mixing is very small and does not appreciably affect the collider bounds on Z .

T-channel Mediator
The t-channel mediator S T is a scalar with the same quantum numbers as the antifermion f (f = t R or f = e R ) that it couples to, whereas the S T − X system shares an effective Z 2 -symmetry, as can be inferred from the Lagrangian in Eq. (2.10). Consequently, the S T phenomenology is very similar to that of a right-handed stop or slepton in the MSSM, with X being the Bino LSP (lightest supersymmetry particle) and the Z 2 -symmetry being the R-symmetry. The LHC searches for stops and sleptons pair production, followed by prompt decays into the corresponding fermion and the LSP (missing energy). The range of values of the coupling y T needed to obtain the correct dark matter abundance results in S T decays that are prompt for the purpose of LHC searches. Therefore, the most stringent bounds on S T arise from the ATLAS and CMS searches for promptly decaying stops [54][55][56][57] and sleptons [58,59]. The LHC bounds on such particles, with an essentially massless LSP (recall that in Sec. 4, we focused on benchmarks with m X = 1 GeV), are approximately m S T 1200 GeV ( 400 GeV) for S T coupled to the top quark (electron), based on ∼ 140/fb LHC Run II data.

Direct Detection
Direct detection signals are generically suppressed in dark matter freeze-in models due to the feeble couplings involved. However, as we discuss here, our low reheat temperature freeze-in scenarios have better direct detection prospects due to generically larger SM-dark matter interactions.
For s-channel Higgs mediated models, the direct detection spin-independent DM-nucleon cross section can be approximated by where we have neglected the contribution of the heavy scalar, and µ Xn is the dark matternucleon reduced mass. The XENON1T result [31] constraints some parts of the parameter space of our freeze-in framework for dark matter masses above ∼ 10 GeV (see the orange curve in Fig. 3). For Dirac dark matter in the kinetically mixed scenario, the spin-independent scattering cross section receives contributions from both the Z and the Z and is given by For light dark matter masses as we consider in the plots of Sec. 3.2, the most relevant bounds are from the CRESST-III experiment [60], which constrain σ SI ∼ 10 −37 cm 2 . Future experiments such as SuperCDMS and NEWS-G will improve on this by several orders of magnitude [61]. However, even these more stringent projections are unable to probe the parameter space relevant for the freeze-in scenario discussed in this paper. For larger values of the dark matter mass and relatively large values of g D q D , the cross sections in Eq. (5.3) could be tested with XENON1T data [31] or with future LZ [62] or DARWIN [63] data. For the t-channel mediator framework, different processes can play the leading role for direct detection. For the model with S T coupled to right-handed electrons, dark matter can scatter with electrons at tree-level (s-channel analog of Fig. 9 (a)) with an approximate cross section For y T ∼ 10 −5 as roughly needed for the correct relic abundance, this cross section is far too small to be probed experimentally. Scattering with nuclei are mediated by one-loop penguin diagrams (analogous to Fig. 9 (d,e)) with the Z/h mediating the scattering with nuclei. For the model with S T coupled to right-handed top quarks, analogous diagrams mediated by gluons are also relevant. Ref. [64] finds that the direct detection cross section mediated by gauge bosons features a further suppression by |q| 2 /m 2 f , where q is the momentum transferred in the scattering process. The Higgs penguin is also negligible since the hXX effective vertex is much smaller than the ZXX effective vertex due to additional mass suppressions (see the discussion below Eq. (4.6)). Therefore, it will be quite challenging to probe the t-channel mediator framework considered in this paper at direct detection experiments.

Summary
In this paper, we have studied scenarios of dark matter freeze-in where the mediator particle that gives rise to interactions between the dark matter (X) and the Standard Model is heavier than the reheat temperature, T RH , in the early Universe. In such setup, the standard approach is to integrate out the mediator and focus on an effective field theory (EFT) with higher dimensional interactions between the SM and DM. We examined the validity of this approach in the regime M med T RH .
We studied three classes of s-channel mediator frameworks: (i) a heavy scalar that mixes with the SM Higgs boson, (ii) a heavy Z that mixes kinetically with the SM hypercharge, and (iii) a Z gauge boson from a gauged L µ − L τ symmetry that couples directly to both SM particles and dark matter. In all cases, the EFT approach (integrating out the mediatorsincluding the SM Higgs / Z boson -and focusing only on the resulting ff → XX processes) captures the correct dark matter freeze-in abundance at very low T RH , where dark matter is dominantly produced through non-resonant annihilation of SM fermions: T RH 5 GeV for the scalar case, T RH 7 GeV for the kinetically mixed Z case, and T RH 0.025 m Z for the L µ − L τ Z model. However, at higher T RH , the resonant enhancement of the schannel annihilation cross section, not captured within the EFT, becomes important. The predicted dark matter freeze-in abundance from fermion annihilation in the EFT can deviate from the full result by several orders of magnitude (see Fig. 7). In this regime, the dark matter abundance is instead appropriately captured by considering decays of exponentially suppressed thermal abundances of the mediators in the bath.
Similarly, we studied t-channel mediator scenarios with couplings to top quarks or electrons (both right-handed). For the e R case, we found that the EFT calculation reproduces the correct dark matter abundance for T RH 0.05 M med . At higher T RH , dark matter is dominantly produced through coannihilation processes (Z/γ) e → S † T X. Similarly, for the t R case, gt → S † T X dominates at high temperatures, T RH 0.04 M med . In contrast, for low T RH scenarios, given the suppressed abundance of top quarks in the thermal bath, loop decays of the SM Z boson induced by the t-channel mediator are the dominant source of dark matter abundance.
We thus find that in both s-and t-channel scenarios, novel channels that are not captured by the EFT dominate dark matter production even when T RH is more than an order of magnitude below the mass of the mediator. It is therefore important to include such contributions when studying dark matter freeze-in production with heavy mediators. Furthermore, these new production channels also change the momentum distribution of dark matter, peaking it towards a warmer distribution than standard freeze-in scenarios.
Finally, we discussed the collider phenomenology of the heavy mediators and the prospects of testing these scenarios at direct detection experiments. We find that, in contrast to mediators in standard freeze-in scenarios, the mediators in our setup can have large couplings with both DM and SM particles, leading to qualitatively different collider phenomenology compared to standard freeze-in setups. Parts of the parameter space of our low T RH scenarios are already probed by LHC searches for Higgs invisible decays, searches for prompt dilepton resonances and for SUSY stops and sleptons, as well as by dark matter direct detection experiments.
temperature signifying the temperature of the Universe at the onset of radiation domination can be written as The maximum temperature the radiation bath reaches during this phase can be written as where ρ φ0 is the initial energy density of the inflaton. Thus, temperatures prior to the radiation dominated era can be higher than T RH if ρ φ0 > Γ 2 φ M 2 P l , i.e., if the initial inflaton energy density is large and transferred to the radiation bath at a very slow rate. Parameterizing we expect T M AX > T RH if α φ < M φ /M P l . In this case, the higher temperatures during this era can result in greater dark matter production from decays of the heavy mediator(s), whose abundances are no longer Boltzmann suppressed, as well as annihilation processes, which can proceed faster. We numerically solve the above differential equations and calculate the production of dark matter from various processes, and find, indeed, that the dark matter abundance from this era dominates over the abundance from the subsequent radiation dominated era for α φ M φ /M P l for which T M AX > T RH , as discussed above. Thus, for this epoch before radiation domination to contribute negligibly to the dark matter abundance in the Universe, we must assume α φ M φ /M P l .

B Specific Models
In this Appendix, we discuss some specific realizations of the simplified models discussed in this paper, and explore how model-specific details can affect phenomenological aspects.

B.1 S-channel Mediator: Sterile Neutrino Dark Matter
Here we discuss a specific model for the s-channel simplified model presented in Sec. 3.1. We consider an extension of the canonical (type I) seesaw mechanism, where singlet right-handed neutrinos, N j , act as portals to a hidden sector (see e.g. [13,15,17,[67][68][69][70] for details): Here L i (i = 1, 2, 3) and h are the SM lepton doublet and Higgs fields, respectively, y ij , y kj are dimensionless Yukawa couplings, and L k , h are hidden sector fermion and scalar states, charged under a hidden sector U (1) symmetry. Integrating out the N j , which we assume to be much heavier than other scales in the model (these will henceforth be ignored, and the notation N k will refer to the light sterile states L ) and dropping indices for simplicity, we have: As before, we also assume a quartic term λ h 2 h 2 that leads to mixing between the two scalars. Here, the scalar h serves as the heavy mediator, S, while one of the sterile neutrinos, e.g. N 1 , is dark matter. If the hidden sector scalar obtains a VEV, v , the dark matter mass is given by y 2 v 2 /M and is naturally much smaller than the mediator mass m S ∼ v , if M v . If the U (1) is gauged, this model also includes the vector (Z ) mediator with mass m Z ∼ g v , which is again significantly heavier than dark matter.
The effective Higgs-sterile neutrino couplings are approximately Here, h, h denote the SM-like and heavy Higgs mass eigenstates, and ν, N denote the SM and hidden sector sterile neutrinos. Note that all of the couplings are suppressed by the heavy scale, M , and are therefore expected to be small. Considering a specific model also gives rise to a novel phenomenology that is not captured by the simplified model discussion. While freeze-in production of dark matter N 1 proceeds via both annihilation of SM fermions and decays of h, h , in this specific model there are also additional production channels, due to the presence of the additional sterile neutrinos N i . These can lead to both annihilation N i N i → N 1 N 1 and decay N i → N 1 N 1 ν contributions via the Higgs and neutrino portals. While the decays of h, h to N 1 are invisible, the presence of the additional sterile neutrino states also gives rise to new collider signatures: h, h → N i N i , which can then further decay into SM fermions, e.g. as N i → νe + e − with displaced vertices (see e.g. [71]).

B.2 T-channel Mediator: Axino Dark Matter
Here we discuss a specific model for the t-channel simplified model presented in Sec. 4. The model contains axino dark matter and is based on Ref. [23] (the interested reader is referred to this paper for more details). Models that solve the strong CP problem using the Peccei Quinn (PQ) symmetry contain a new particle, the axion. Supersymmetric extensions also contain its superpartner, the axino. The couplings of the axino to the MSSM particles (as is the case for the axion) are suppressed by the PQ scale, f a . Therefore, the axino only has feeble couplings to the remainder of the MSSM field content and, as we discuss below, can be a dark matter candidate with an abundance set by the freeze-in mechanism.
Since the axino is part of a chiral multiplet, it does not acquire a tree-level Majorana soft mass term from supersymmetry breaking, but can obtain a small effective mass due to mixing with neutral heavy states by virtue of the presence of Dirac mass terms, or from loop effects.
The axino can therefore naturally be much lighter than the MSSM particles. We henceforth treat the axino mass as a free parameter. We also confine ourselves to scenarios where T RH is significantly below the scales of both PQ and supersymmetry breaking, so that the axino is the only relevant SUSY particle in the early Universe.
Following Ref. [23], we focus on the KSVZ-type axion models [72,73], where SM particles are not charged under the U (1) P Q symmetry. In this setup, the axino couples to MSSM states via loops of heavy PQ states (with masses at the PQ scale, f a ). The Lagrangian contains a dimension-5 axino-gluon-gluino vertex, which is given by Likewise, an axino-quark-squark vertex arises at two loops and is given by [23] g ef f ≈ α 2 Depending on the mass hierarchy between the squarks and gluinos, either of these two interactions can dominate dark matter production in the early Universe through the t-channel processes gg → aa, qq → aa. If we assume that gluinos are much heavier than squarks and can be neglected, then the only relevant interaction for dark matter production is the axinoquark-squark interaction. Thus, this well motivated supersymmetric framework maps onto our t-channel simplified model Lagrangian in Eq. (2.10), with the squarks acting as the heavy mediator, S T , and g ef f ( 1, as mg f a ) as the coupling y T .