Freezing In, Heating Up, and Freezing Out: Predictive Nonthermal Dark Matter and Low-Mass Direct Detection

Freeze-in dark matter (DM) mediated by a light ($\ll$ keV) weakly-coupled dark-photon is an important benchmark for the emerging low-mass direct detection program. Since this is one of the only predictive, detectable freeze-in models, we investigate how robustly such testability extends to other scenarios. For concreteness, we perform a detailed study of models in which DM couples to a light scalar mediator and acquires a freeze-in abundance through Higgs-mediator mixing. Unlike dark-photons, whose thermal properties weaken stellar cooling bounds, the scalar coupling to Standard Model (SM) particles is subject to strong astrophysical constraints, which severely limit the fraction of DM that can be produced via freeze-in. While it seems naively possible to compensate for this reduction by increasing the mediator-DM coupling, sufficiently large values eventually thermalize the dark sector with itself and yield efficient DM annihilation to mediators, which depletes the freeze-in population; only a small window of DM candidate masses near the $\sim$ GeV scale can accommodate the total observed abundance. Since many qualitatively similar issues arise for other light mediators, we find it generically difficult to realize a viable freeze-in scenario in which production arises only from renormalizable interactions with SM particles. We also comment on several model variations that may evade these conclusions.


Introduction
Despite overwhelming evidence demonstrating the existence of dark matter (DM) on galactic and cosmological scales [1], identifying its possible non-gravitational interactions remains one of the greatest challenges in contemporary physics. Historically, much of the DM search program has been guided by the weakly interacting massive particle (WIMP) hypothesis in which ∼ TeV-scale DM carries electroweak charge and acquires its abundance by freezing out of equilibrium from the Standard Model (SM). However, decades of null searches [2] have motivated a new, broader paradigm in which DM is the lightest stable SM singlet in a "hidden sector" with its own dynamics and cosmological history (see [3,4] for a review). In these scenarios the viable DM mass range can span many orders of magnitude, so there is strong motivation for new detection techniques.
Recently, there has been great progress in proposing new direct-detection targets for the ∼ keV-GeV DM mass range. Some representative ideas involve scattering off atomic electrons [5][6][7], superconductors [8,9], semiconductors [5,[10][11][12], Helium evaporation [13], Fermi-degenerate materials [14], Dirac metals [15], Graphene targets [16], magnetic bubble chambers [17], color centers [18], and chemical bonds [19] (for a review, see [4]). Unlike traditional direct detection experiments, designed for WIMP scattering off heavy nuclei with ∼ keV recoil energy thresholds, many of these techniques exploit observable transitions between small internal energy levels (e.g. atomic ionization) to probe much lighter DM, which typically deposits ∼ meV − few eV per scatter. Furthermore, these targets are particularly sensitive to interactions mediated by light particles φ whose recoil spectra are sharply peaked towards low recoil energies, so the scattering cross section can be greatly enhanced and may be observable even for extremely small couplings. Here g χ/f are the φ couplings to dark/visible particles where f is a SM fermion and, for illustration, we have taken m χ m e and normalized to the characteristic momentum transfer q ∼ αm e for DM scattering off atomic electrons [5]. However, if g f is sufficiently large to thermalize φ with the SM in the early universe, cosmological bounds from BBN and ∆N eff generically require m φ MeV [20], which is outside the range in which Eq. (1.1) holds. Thus, in order to exploit this huge enhancement, we demand g f 1 to avoid thermalizing with the photons, which excludes thermal DM production with a sub-MeV mediator, 1 but still allows for non-equilibrium alternatives.
Aside from thermal freeze-out, the only other predictive, UV insensitive production mechanism is "freeze-in" [22,23], in which DM is not initially populated during reheating, but produced later through sub-Hubble interactions with SM particles. Since DM never achieves a thermal number density, its Boltzmann equation is integrable and the final abundance is Hs σv SM→χχ (traditional freeze in), (1.2) where T is the SM temperature, H is the Hubble expansion rate, s is the entropy density, ρ cr is the critical density, n SM is the number density for some SM species, and a 0 subscript represents a present day value. As with freeze-out, the abundance is fixed by the DM-SM interaction and is insensitive 2 to the details of earlier, unknown cosmological epochs (e.g. inflation, reheating). Although the couplings required for freeze-in production are generically very small for light mediators, it may nonetheless be possible to exploit the low thresholds of these experiments to bring the cross section in Eq. (1.1) into the observable range. It has been shown that many new direct-detection targets are potentially sensitive to the freezein parameter space for DM with interactions mediated by an ultra-light ( keV) dark-photon [12,24]. However, this particular model has many special features, which do not generalize to other mediators (or even to heavier dark-photons). Most notably, the in-medium thermal suppression of kinetic mixing at finite temperature counterintuitively makes the nearly-massless dark-photon safer from astrophysical and cosmological bounds because all on-shell production vanishes in the massless mediator limit [25,26]. Furthermore, in this limit, the cosmological dark-photon production is also suppressed, so it can be neglected in the Boltzmann system and Eq. (1.2) can be straightforwardly applied. However, in models for which mediator production is not parametrically suppressed in the early universe, Eq. (1.2) is no longer valid and it is necessary to solve the full Boltzmann system that tracks both χ and φ densities. Similarly, in such models there are nontrivial astrophysical bounds on the mediator-SM coupling g f [25,26]. Thus, it is important to identify other viable light-mediator scenarios that enjoy the enhancement in Eq. (1.1) and understand the generic issues they present.
Pioneering earlier work [24] studied DM production through renormalizable portals and its subsequent self-thermalization via dark force interactions. However, this study primarily 3 considered a massless dark-photon mediator, with the above-mentioned special properties. Furthermore, [27] and [28] recently studied the models and constraints for light mediators accessible to new direct detec-1 A loophole around this argument involves sub-MeV DM with a comparably light mediator, both of which thermalize with neutrinos after they decouple from the photons [21]. However, such models must still avoid thermalization with charged particles.
2 However, if χ is produced through nonrenormalizable higher-dimension operators whose suppression scale exceeds the highest temperature ever realized in the early universe Tmax, then the abundance in Eq. (1.2) will depend on UV physics through the ratio Tmax/Λ. However, this is never the case for the light mediators considered in this paper. 3 This paper also includes a discussion of heavy scalar DM S with a Z 2 stabilizing symmetry. In this variation, DM is produced directly through the Higgs portal interaction S 2 H † H and there is no light mediator. tion techniques; however, these analyses leave open the question of whether there exist predictive non-thermal production mechanisms beyond the ultra-light dark-photon mediated scenario.
In this paper we extend this discussion by studying a representative benchmark of relevance for low-mass direct detection: the predictive freeze-in production of dark matter through a light, feebly-coupled scalar mediator φ with Higgs portal mixing. We find that even when g f 1 and all SM-DM interaction rates are slower than Hubble, there is generically a sizable φ population, which can thermalize with χ if g χ is sufficiently large. For fixed, viable choices of g f , there are three qualitatively distinct regimes: • Small Coupling: If g χ is sufficiently small, the hidden sector never thermalizes with itself, so hidden-annihilation is always sub-Hubble Γ(χχ → φφ) H, and production approximates traditional freeze-in. In this regime, the χ abundance increases with g χ and Eq. (1.2) is approximately valid.
• Large Coupling: If g χ is sufficiently large, hidden-annihilation becomes efficient Γ(χχ → φφ) H and the hidden sector thermalizes with itself (but still not with the SM). In this regime DM production initially resembles freeze-in, but eventually the χ-φ system reaches thermal equilibrium (at a lower temperature than the SM) and χ undergoes an additional phase of freeze-out. In this regime, increasing g χ merely converts more of the χ into dark radiation.
• Intermediate Coupling: At intermediate values of g χ for which the maximum hidden-annihilation rate is briefly comparable to Hubble Γ(χχ → φφ) ∼ H, production interpolates between the above regimes and yields a maximum DM abundance for each choice of mass.
Intriguingly, we find that for viable mediator-SM couplings g f 1, safe from astrophysical and cosmological bounds, it is generically difficult to account for the total freeze-in abundance with only renormalizable couplings to SM particles. Nonetheless, for choices of g f that saturate existing limits, the maximum abundance that can be produced is significantly greater than the naive expectation from Eq. (1.2) once mediator-DM interactions are included in the Boltzmann system. This paper is organized as follows: section 2 introduces our benchmark model; section 3 describes the early universe production mechanism; section 4 discusses model variations that evade some of our general conclusions; and section 5 we summarize our findings and suggest future directions of inquiry.

Benchmark Model
Our benchmark scenario involves a scalar singlet mediator φ and mixes with the SM through the Higgs portal operators where A and κ are the renormalizable portal couplings. After electroweak symmetry breaking (EWSB) the neutral component of the doublet h mixes with φ and, in the small mixing limit, the system is diagonalized with the shift h → h + sin θ φ, so φ acquires mass-proportional couplings to SM fermions. Including a φ Yukawa interaction with a fermionic DM candidate χ, the relevant Lagrangian becomes 2) h f f Figure 1. Leading Feynman diagrams giving rise to non-thermal χ freeze-in production from SM annihilation during the early universe (left) and reannihilation into dark sector mediators φ (middle, right) once the hidden sector thermalizes with itself. Although the SM initiated process must be slower than hubble to avoid thermalizing the dark sector, the χ−φ interactions can achieve chemical equilibrium at a separate temperature set by the energy density transmitted to the dark sector via the left diagram and those depicted in Fig. 2.
where f is a SM fermion of mass m f and v 246 GeV is the SM Higgs vacuum expectation value. Although this is only one of several possible scalar mediated scenarios, as we will see, it captures much of the essential physics and many of the issues encountered here apply to a much broader class of variations on this simple setup (see Sec. 4 for a discussion and [29,30] for a complementary studies of thermal freeze-out with the same field content).
In this model, the cross section for nonrelativistic χe scattering is where q is the three-momentum transfer, µ χe is the reduced mass, and we take the light mediator limit to highlight the parametric enhancement from Eq. (1.1) in terms of q ∼ αm e , the characteristic momentum transfer for χ scattering off atomic electrons [5]. Our main results (described in Sec. 3 and displayed in Fig. 5, right panel) are presented in terms of the convention in Eq. (2.3) and are largely independent of the mediator mass so long as m φ αm e , m χ . Although our emphasis in this paper is on the freeze-in production for this scenario, it has also been shown that models with light scalar mediators m φ ∼ 10 MeV can yield sufficiently strong DM self interactions σ χχ /m χ ∼ cm 2 g −1 to reduce tension between DM-only N-body simulations and the observed small scale properties of haloes [31]. It may be interesting to explore whether DM with a scalar mediator in the mass range considered here m φ αm e could preserve some of these features (as in [29] for heavier scalars), but this question is beyond the scope of the present work.

Boltzmann System
Our initial condition assumes a hot, radiation dominated universe in the broken electroweak phase. 4 We further assume that n χ = n φ = 0 after reheating so that sub-Hubble interactions with SM fields are solely responsible for populating the hidden sector. 5 The Boltzmann equation for χ production is Figure 2. Representative Feynman diagrams giving contributing to φ freeze in production in the early universe. These processes determine ρ φ which (in addition to ρχ) determines the hidden sector's temperature TD. The effective vertex in the right diagram arises from integrating out the top quark for energies in between the weak scale and the QCD confinement scale (see Appendix B).
where H ≡ 1.66 √ g * T 2 /m Pl is the Hubble rate during radiation domination, m Pl = 1.22 × 10 −19 GeV is the Planck mass, and g * is the number of relativistic degrees of freedom. Here an eq superscript denotes an equilibrium quantity, and we do not assume equilibrium between visible and hidden sectors (T = T D ); all relevant SM species are in equilibrium with each other, so we omit their superscript. As discussed above, existing bounds on the mixing angle require that the dark sector never thermalize with the SM in the early universe, so n χ n eq χ (T ), but we allow for the possibility that χχ ←→ φφ annihilation (shown in the middle and right diagrams of Fig. 1) equilibrates the hidden sector so that χ eventually achieves a thermal distribution n χ → n eq χ (T D ), where T D T . Thus, we can simplify the Boltzmann equation where the first two collision terms are independent of n χ and serve merely as sources for χ production.
In the limit where the hidden-sector annihilation rate is also sub-Hubble throughout the production process Γ(χχ → φφ) H, Eq. (3.2) recovers the usual freeze-in form in Eq. (1.2).

Dark Temperature Evolution
In this scenario, the hidden sector temperature T D is defined in terms of the evolving dark sector energy density, which is also produced through SM interactions. The energy transfer Boltzmann equation for direct χ production is and the contributions from mediator production (see Fig. 2) are and we define the thermally averaged quantity where ∆E is the energy transferred to the hidden sector and v mø is the Møller velocity (see Appendix A). Note that the definition of the energy transfer differs depending on whether we are considering χ 10 -10 Equilibrium Yields, g χ = 10 -3 , sin θ = 10 -10 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -5 Maximum Abundance, eV ≲ m ϕ ≪ αm e , sin θ = 10 -10 Figure 3. Left: Comoving equilibrium yields Y eq χ = n eq χ (TD)/s(T ) plotted against the dimensionless time variable z = mχ/T for a various mχ. Note that the dark temperature is determined primarily by the φ energy density, which is independent of mχ, so each mass point has comparable TD during production; thus, sufficiently large masses are highly Boltzmann suppressed. Right: Maximum DM abundance for a range of mχ with fixed sin θ = 10 −10 , which is still viable, but near the Red Giant exclusion bound (see left panel of Fig. 5). Below mχ 200 MeV, the blue curve represents the coupling for which the maximum abundance (red curve) can be achieved; larger values merely annihilate away more of the DM as the number density tracks the equilibrium distributions (left) out to larger values of z (see Fig 4,left). In the green band between the dashed vertical lines, max Ωχh 2 > 0.12 so it is possible to achieve the full abundance for two choices of couplings, as shown for a representative mass point in Fig. 4.
pair production in Eq. (3.3) for which ∆E tracks both final state particles, or single φ production in Eq. (3.4) where it tracks only the energy transferred to φ. The combined dark sector energy density ρ D = ρ φ + ρ χ defines a dark temperature where ξ D counts the relativistic degrees of freedom in the hidden sector. Note that the dark temperature initially satisfies T D = 0 and depends nontrivially on the visible temperature through Eqs. (3.3) and (3.4), which source its growth in Eq. (3.6). The dark temperature also determines the equilibrium χ number density in the hidden sector which enters into the RHS of Eq. (3.2). In our Boltzmann system we have omitted all processes involving lighter SM fermions, which only contribute negligibly to χ and φ production in the early universe; the dark sector is mainly populated while heavy electroweak states are in equilibrium with the SM thermal bath with corrections of order ∼ (m f /m t ) 2 for all lighter SM fermions. For this reason, we can safely neglect longitudinal plasmon mixing effects which can resonantly enhance the φ production at lower temperatures for which m φ is near the plasma frequency [26]; this enhancement does not compensate for the parametric suppression in the coupling relative to continuum production off heavy electroweak states. However, such an effect may be important in model variations in which φ production arises mainly from its coupling to electrons (see Sec. 4 for a discussion). Although freeze-in is insensitive to the precise value of the reheat temperature, we assume it to be above m t for all SM electroweak states to be relativistic when dark sector production begins; if the reheat temperature is lower, DM production will proceed through lighter fermions and the resulting abundance will be significantly suppressed, but the qualitative features of our discussion are unaffected.

Numerical Results
For the relevant parameter space considered in this paper m φ αm e , the hidden sector temperature is primarily set by ρ φ whose growth rate is proportional to sin 2 θ in Eq. (3.4), whereas the collision terms for ρ χ in Eq. (3.3) all depend on g 2 χ sin 2 θ and are further suppressed. Furthermore, ρ φ continues to accumulate through gg → gφ processes long after direct χ production processes (e.g. tt → χχ) have frozen out. Thus, the cosmological evolution of T D is insensitive to the DM mass, so sufficiently large values for which m χ T D encounter hidden-sector Boltzmann suppression throughout their production, φ-thermalization, and re-annihilation phases.
The left panel of 3 presents the equilibrium χ yield Y eq χ ≡ n eq χ (T D )/s(T ) normalized to the total entropy, where n eq χ (T D ) is computed using Eq. (3.7) and depends on the visible temperature implicitly through Eq. (3.6). Here we see that, as m χ increases, the heights of these curves begin to fall, which ultimately limits the total abundance that can be generated for a given mass point; increasing g χ increases the abundance in the freeze-in regime, but once the dark sector equilibrates with itself, the yield tracks these equilibrium distributions and eventually diminishes as more of the DM annihilates away into φ radiation. Thus, for each choice of mass there is a maximum abundance as represented by the red curve in Fig. 3 (right panel), which corresponds to the value of g χ from the blue curve.  ] Higgs-Portal Mediator, sin θ = 10 -10 , eV ≲ m ϕ ≪ α m e Figure 5. Left: Available parameter space for a light, Higgs-mixed scalar particle plotted alongside constraints from 5th forces [10,[35][36][37][38][39], stellar cooling [26], and cosmology [40]. The vertical line represents the approximate boundary where the χe scattering rate undergoes a parametric shift for scattering off atomic electrons, which is taken to be the reference value for presenting different techniques in [4]. For other materials with different characteristic momenta, this shift takes place for different values of m φ . Right: Constraints and projections for the χe cross section for sub-GeV DM. Some representative projections for the SENSEI 2e − , Silicon 1e − , and superconductor targets are taken from [4], and the Dirac metals projection is from [15]; all projections assume ΩχΩDM. For points outside the green vertical band, the blue curve depicts (max Ωχ/ΩDM)× σe where sub-Hubble production as described in Sec. 3 generates the maximum fractional abundance (max Ωχ < ΩDM) for each mχ. Inside the vertical green band, the blue curves represent σe evaluated at gχ for which the full DM abundance can be accommodated; note there are two such values as depicted in Fig. 4 and on the right panel of Fig. 3. Although the Higgs portal mixing sin θ = 10 −10 may seem ad-hoc, this is chosen as a representative viable value near the stellar cooling exclusion region depicted on the left panel of this figure.

Cosmology
Note that inside the green band, the blue curve splits into two segments corresponding to the distinct values that yield the observed DM abundance Ω χ h 2 0.12, and not the maximum abundance (as it does outside this band).
This counterintuitive behavior can be understood from the plots in Fig. 4 which shows the production history for a representative mass point m χ = GeV and sin θ = 10 −10 as the coupling g χ is varied. On the left panel, the red curves are solutions to the Boltzmann equation in Eq. (3.2) for two choices of coupling that yield the observed DM abundance; the blue curves at the bottom of the plot also show the "naive" freeze-in result for the same choices of g χ ; as in Eq. (1.2), the blue curves are computed by merely integrating the SM collision terms in Eq. (3.2) without including the additional χχ ←→ φφ dynamics. The right panel of Fig. 4 presents the DM abundance as a function of g χ as computed with the full Boltzmann system (red curve) and with the naive freeze-in integration (blue curve).
On the right panel of Fig. 5 we present our main result: the predictive parameter space for χ production in the early universe plotted in terms of the fiducial χe scattering cross section where A(q 0 ) is the matrix element for a free-electron target evaluated at reference momentum q 0 = αm e following the convention in [4], which also applies to the other projections shown. For the blue curve in this plot, this cross section is given by Eq. (2.3) in the m φ αm e limit and evaluated at g χ values, which either yield the full DM abundance or the largest possible subdominant fraction where this is not possible as shown in Fig. 3 (right); in latter case the cross section is rescaled by the fractional abundance. The vertical green band in plot highlights the mass window for which the total abundance can be achieved and, as in Fig. 3 (right), the blue curve splits into two segments in the shaded green band to reflect the two values of g χ that can accommodate the total abundance in; the top segment corresponds to χ that yield freeze-out in the hidden sector, whereas the bottom segment presents the freeze-in like value. Note that for most of the parameter points shown here max Ω χ 10 −2 Ω DM (see Fig. 3, right), so DM self interaction bounds do not constrain the viable parameter space; for a discussion of this issue see [27]. Relatedly, there is no perturbative unitary violation in χχ scattering for the g χ along this curve since the values that yield the maximum possible abundance for viable choices of sin θ are all well below unity (see Fig. 4, right). For completeness, we also show the viable parameter space for a Higgs portal mediator on the left panel of Fig. 5 to motivate our benchmarks.

Model Variations
Thus far, we have restricted our attention to the scenario in which both χ and φ are produced through h-φ mixing during the early universe. We have found that the bounds on sin θ 5 × 10 −10 (see left panel of Fig 5), severely limit the total energy density transferred to the hidden sector, thereby making it impossible to achieve the total DM abundance for m χ few 100 MeV for all values of g χ . Here we consider some possible solutions to this problem that still feature sub-Hubble DM production and light scalar mediators.

Different SM Current
Since much of the discussion above depends on the Higgs-portal coupling relations, it may be possible to evade some of the above complications with a different pattern of couplings. For the scenario studied in this paper, the main impediment to generating the observed DM abundance over the eV m φ αm e parameter space of interest is the limit sin θ 5 × 10 −10 , which is driven primarily by resonantly enhanced φ production through its coupling to electrons in Red Giants (RG) [26] -see Fig. 5 (left). For a minimal Higgs-mixing scenario, the same mixing angle θ rescales the φ coupling to electroweak states and thereby also suppresses the total energy energy/number density transferred to the hidden sector via electroweak production in the early universe. However, the early universe production is dominated by the heaviest SM states (t, h, W/Z) at T ∼ v, so if the φ coupling to these particles can be enhanced relative to the electron coupling in a minimal Higgs-mixing scenario it is possible to increase the overall abundance at late times.

Two Higgs Doublets
One simple extension, which can accommodate such a variation involves mixing φ with the the neutral CP-even states of a two-Higgs doublet model (2HDM) to break the usual SM relations between quark and lepton couplings. For instance, in a Type II 2HDM, φ could mix with the predominantly uptype doublet, thereby increasing the relative ratio of t/e couplings and enhancing dark sector particle production in electroweak processes. However, it is not clear whether a concrete realization of such a scenario can simultaneously accommodate the observed DM abundance while satisfying the existing limits on 2HDMs, which typically force the 125 GeV CP even scalar close to alignment limit, which approximately reproduces the usual SM coupling patters [41,42], so a careful study is necessary to see if this is possible.

Vectorlike Fermions
Alternatively, the scalar mediator can couple to heavy vectorlike quarks, which are integrated out to induce effective interactions with SM particles (e.g. φ G µν G µν , φ F µν F µν , or φ qq if vectorlike states mix with SM quarks). If the vectorlike quarks are sufficiently decoupled, all dark sector production proceeds through effective SM interactions and therefore this class of models makes a firm prediction for the DM-SM scattering rate. Furthermore, in this variation there is no tree-level electron coupling, so Red Giant bounds can be significantly weakened [26] and more DM can be produced via freeze-in, but this scenario is harder to test with new detection techniques, most of which probe DM-electron scattering in various materials.
If, instead, the mediator couples to heavy vectorlike leptons, the IR theory will contain the φ F µν F µν operator (and also φ¯ if the new heavier states mix with right handed leptons). However, this variation does not ameliorate the Red Giant bounds, which now constrain all early-universe production. Unlike in the Higgs-mixing scenario, where production off heavier electroweak states in Eqs. (3.3) and ( 3.4) is parametrically enhanced relative to the electron scattering cross section, here the production and detection depend on the same couplings, so dark sector production is greatly suppressed. Although the dominant φ production in this leptophilic scenario occur at lower temperatures T ∼ m φ where there is a resonant enhancement from φ-longitudinal-plasmon mixing, this effect has been found to be modest (order-few) relative to non-resonant continuum production in astrophysical contexts [26], so it is unlikely to compensate for the large reduction in the overall rate. Nonetheless, a full calculation of resonant φ production at these lower temperatures is worth further investigation. 6

Asymmetric Freeze-In
For the the scenario considered in earlier sections, the abundance is ultimately limited by the additional χχ → φφ annihilation phase that becomes efficient at depleting the χ population for sufficiently large g χ . However, if the DM carries a conserved global quantum number and has an asymmetric population at late times, then increasing g χ will not lead to further depletion once all antiparticles are annihilated away. In order to maintain the qualitative features of the scenario considered here (light mediator and sub-Hubble production), the freeze-in interaction that produces the DM must also satisfy the Sakharov conditions [44]. Although models of asymmetric freeze-in have been proposed in the literature [45], the asymmetry they produce is suppressed by subtle cancellations [46] and it is not clear whether these limitations can be overcome, but this problem deserves further study (see [47] for a discussion).

Additional Production Source
One simple way to enhance DM production is to introduce a new source of χ or φ production beyond just the Higgs portal coupling to SM particles in Eq. (2.2). For instance φ could also mix with a new, heavy singlet scalar S which thermalizes with the SM in the early universe. Since the S-φ mixing angle is a new free parameter, it is possible to increase the χ and φ abundances via S initiated scattering and decay processes, which still realize freeze-in production if this mixing is suitably small. However, this modification is no longer predictive and manifestly depends on presently unknown UV physics, thereby detracting somewhat from the appeal of the freeze-in mechanism.
In addition to changing the SM-φ flavor structure, the SM extensions discussed in Sec. 4.1 can provide additional production sources if the new heavier particles are populated on shell after reheating. New interactions from these states can enhance φ and χ production in the early universe, but the magnitude of this enhancement now depends on the details of the extended sector(s), which adds UV sensitivity to the production mechanism.

Concluding Remarks
The recent burst of creativity in devising new techniques sub-GeV direct detection may soon make it possible to probe feebly coupled "freeze-in" dark matter coupled to a very light keV mediator. However, the only predictive, UV complete example in which this has been demonstrated (DM coupled to an ultra-light dark-photon) has subtly special properties that may not generalize to other mediators. To understand how robustly such scenarios can be probed, we have studied the freeze-in production of dark matter through a light scalar mediator with Higgs portal mixing.
In our analysis we find that the recently updated stellar cooling bounds on light scalars have severe, generic implications for any such freeze-in scenario. Most notably, for any sufficiently light ( keV) scalar mediator, compensating for the suppression from this bound requires a large, orderone coupling to dark matter, which quickly thermalizes the dark sector with itself. Even though all interactions between dark and visible sectors remain slower than Hubble expansion, the dark sector's separate thermal bath enables efficient dark matter annihilation into mediators, thereby depleting its abundance once the mediator-dark-matter coupling increases beyond the threshold for hidden sector thermalization. Thus, for each choice of dark matter mass and Higgs mixing-angle, there is a maximum cosmological abundance that can be accommodated if the dark sector is only produced through its interactions with the Standard Model. We find that for Higgs portal mixing-angles near the limit imposed by astrophysical bounds, it is impossible to generate the observed dark matter abundance except for a narrow window in the ∼ 100 MeV -few GeV range. Furthermore, even inside this range, the cross sections for direct detection off electron targets are several orders of magnitude below future sensitivity projections for proposed experimental techniques.
We have also considered some possible model variations, which may alter these conclusions and allow for successful freeze-in with a light ( keV) scalar mediator. In order to produce more DM during the early universe, while satisfying astrophysical bounds (which mostly constrain the mediatorelectron coupling), it may be possible to generalize the minimal Higgs-portal scenario sector to mix the mediator with a larger electroweak sector (possibly a 2HDM) and thereby enhance the top-mediator coupling. This modification increases early universe freeze-in production while heavy electroweak states are still in equilibrium, but maintains an appreciable coupling to electrons for direct detection at late times. Similarly, it may be possible to enhance DM production by preferentially coupling the light mediator to a heavy vectorlike fourth generation of quarks, which are integrated out to generate higher dimension operators with SM quarks and gluons. Such a may alleviate astrophysical bounds on the mediator coupling and thereby enable viable freeze in production. Other possible variations include either additional (beyond Standard Model) sources of mediator or DM production and/or additional CP violation to realize a light-mediator variation on asymmetric freeze-in, but these possibilities are beyond the scope of the present work and deserve further study.

Appendix A: Thermal Averaging
In this appendix we apply the methods used in [48] to compute thermal averages for annihilation cross sections and energy transfer rates in the early universe.

Number Density Collision Terms
The full Boltzmann equation for DM depletion via χ(p 1 )χ(p 2 ) → A(p 3 )B(p 4 ) into SM final states A and B is H ≡ȧ/a = 1.66 √ g * T 2 /m P l is the Hubble expansion rate during radiation domination and the collision term for the number density is where |M| 2 is the spin averaged squared matrix element and we have omitted bose(pauli) enhancement(blocking) factors. Using detailed balance f eq χ (p 1 )f eq χ (p 2 ) = f eq A (p 3 )f eq B (p 4 ) and following the derivation in [48], the Botlzmann equation becomes where the thermally averaged annihilation cross section is defined to be where K 1,2 are Bessel functions of the first and second kinds, respectively. In the traditional "Freeze In" regime, the initial condition is n χ (0) = 0 and the χ production rate in SM annihilation processes is slower than Hubble expansion, so the n 2 χ term can be neglected and the Boltzmann equation becomes integrable and we can define comoving yields Y i ≡ n i /s and change the time variable to z ≡ m χ /T to obtain which can be integrated to obtain the asymptotic abundance at late times where s 0 2969cm −3 is the present day CMB temperature, ρ cr = 8.1h 2 × 10 −47 GeV 4 is the critical density.

Energy Density Collision Terms
Here we generalize the argument in [48] to calculate the thermally averaged energy transferred to particle 1 (which we will later identify with χ or φ) in a 1 + 2 ←→ 3 + 4 process governed by the Boltzmann equation Ignoring quantum statistical factors and assuming the other particles are in equilibrium with the radiation bath, the collision term is where we make no assumption about particle 1 being in thermal equilibrium with the others. Using detailed balance f eq 3 (p 3 )f eq 4 (p 4 ) = f eq 1 (p 1 )f eq 2 (p 2 ) and the definition of the cross section where F = (p 1 · p 2 ) 2 − m 2 1 m 2 2 is the usual flux factor, so we get Since the production/absorption rate is always slower than Hubble expansion for our purposes, the f 1 term can be neglected, so we get C ρ = d 3 p 1 2E 1 (2π) 3 d 3 p 2 2E 2 (2π) 3 4F σ 12→34 E 1 f eq 1 (p 1 )f eq 2 (p 2 ) = P 12→34 n eq 2 n eq 1 , (A.12) where the thermally averaged power transfer rate is and where the Møller velocity is v mø ≡ F/E 1 E 2 . the thermally averaged power is and denominator factor can be evaluated analytically (A.14) To evaluate the numerator of Eq. (A.13), we follow the procedure in [48] we perform all trivial angular integrations and perform the coordinate transformation (E 1 , E 2 cos θ) → (E + , E − , s) where θ is the CM angle between the incoming three-vectors, E ± = E + ± E − , and s = (p 1 + p 2 ) 2 to get where we have used ε = 1 − s 0 /s E 2 + − s, so the power transfer becomes 16) which is valid even in the m i → 0 limit provided that the Bessel functions are properly expanded to cancel the mass dependence in the prefactor. Using detailed balance, we can relate the resulting collision term to the reverse process P 12→34 n eq 1 n eq 2 = P 34→12 n eq 3 n eq 4 , (A. 17) which describes the SM+SM → SM φ processes in Eq. (3.4) where we identify φ with particle 1 of this argument. For processes in which 2χ particles are produced via SM annihilation, a similar argument yields the expression which tracks the energy transfer to both hidden sector particles produced in this process. whose form is appropriate for the collision terms in Eq. (3.3).