Freeze-in, glaciation, and UV sensitivity from light mediators

Dark matter (DM) freeze-in through a light mediator is an appealing model with excellent detection prospects at current and future experiments. Light mediator freeze-in is UV-insensitive insofar as most DM is produced at late times, and thus the DM abundance does not depend on the unknown early evolution of our universe. However the final DM yield retains a dependence on the initial DM population, which is usually assumed to be exactly zero. We point out that in models with light mediators, the final DM yield will also depend on the initial conditions assumed for the light mediator population. We describe a class of scenarios we call"glaciation"where DM freezing in from the SM encounters a pre-existing thermal bath of mediators, and study the dependence of the final DM yield on the initial temperature of this dark radiation bath. To compute DM scattering rates in this cosmology, we derive for the first time an exact integral expression for the Boltzmann collision term describing interactions between two species at different temperatures. We quantify the dependence of the DM yield on the initial dark temperature and find that it can be sizeable in regions near the traditional (zero initial abundance) freeze-in curve. We generalize the freeze-in curve to a glaciation band, which can extend as much as an order of magnitude below the traditional freeze-in direct detection target, and point out that the DM phase space distribution as well as the yield can be strongly dependent on initial conditions.


Introduction
The hypothesis of thermal contact between dark matter (DM) and the Standard Model (SM) is a powerful organizing principle for predictive and testable models of DM. The most common such paradigm is thermal freeze-out [1], where DM is in thermal equilibrium with the SM for some period before the expansion rate of the universe exceeds the DM annihilation rate. In this case the relic abundance of DM is a remnant of the original thermal population, and is "UV-insensitive" in the sense that it is only physics at late times that sets the DM abundance. An alternate paradigm is freeze-in [2][3][4][5][6], where the relic abundance of DM is gradually built up through rare processes that produce DM from the SM thermal plasma. In this scenario the DM never attains thermal equilibrium with the SM and thus predictions for its abundance necessarily retain some dependence on initial conditions; however, when the mediating interactions are renormalizable, the DM production rate peaks at late times, resulting in a weaker but still valuable form of UV-insensitivity. The freeze-in paradigm is especially appealing as a target for direct detection [7,8], because if the mediator of the interaction is light (for example, a kinetically-mixed dark photon), even weak DM-SM interactions are enhanced at low velocities, leading to promising sensitivity at current and future terrestrial experiments  (see Ref. [58] for a review of recent progress).
In models of "traditional" IR-dominated freeze-in where the DM does not interact appreciably after it is produced, the residual UV sensitivity amounts to a constant offset in the DM yield for a given parameter point, as we briefly discuss. However, producing DM via a light mediator necessarily implies that DM interactions with mediator particles can give rise to DM-number-changing processes at cosmologically-interesting rates. A cosmological population of dark mediators can therefore substantially affect the final DM number density that results from a particular coupling to the SM. Here we quantitatively assess the UV sensitivity that arises in models with different initial conditions for the light mediator, and demonstrate that different initial conditions for the dark sector can give rise to very different cosmological histories for the same couplings. We consider the simple and generic scenario when this population is thermal, i.e., in kinetic equilibrium at a temperatureT , which in general will differ from the SM temperature T .
There is a substantial body of literature studying the interplay between freeze-in and freeze-out processes in determining the final relic abundance of DM in hidden sectors with light mediators [59][60][61][62][63][64][65][66][67][68]. These studies consider the case where the energy density in the dark radiation bath is built up entirely from the energy injected from the SM. The novel point we focus on here is the qualitatively new sensitivity of the final DM relic abundance to the initial dark sector population (DM and light mediators), which we parameterize through an initial temperature ratio ξ i =T i /T i ; previous work corresponds to setting ξ i = 0. We call freeze-in into a pre-existing thermal bath "glaciation". Taking a Dirac fermion χ interacting with a light kinetically-mixed dark photon Z D (m Z D m χ ) as our benchmark model for the dark sector, we establish the regions of parameter space where both traditional freeze-in and glaciation are self-consistent descriptions of the theory. We demonstrate that for larger values of the model couplings, the energy injection from the SM overwhelms the initial conditions and predictions are UV-insensitive, while for values near the traditional freeze-in curve, the realized DM abundance can depend sensitively on the initial temperature ratio. Our results substantially clarify the theoretical status of the freeze-in curve as a target for direct detection experiments, and motivate an expanded "glaciation band" which can extend up to an order of magnitude below the freeze-in cross section. This paper is organized as follows. In Sec. 2, we review the traditional freeze-in paradigm with a light kinetically-mixed mediator and zero initial abundance, and show that the assumption of no self-interactions is valid up to a maximum value of α D . In Sec. 3, we introduce the thermalized dark sector population and set up and solve the Boltzmann equations relevant for the more general glaciation scenario. As a consequence of our analysis, we derive for the first time an exact expression for the collision term describing interactions between two populations at different temperatures. These Boltzmann equations assume that the injected DM achieves rapid kinetic equilibrium with the SM, and in this section we delineate the parameter space where this assumption is valid. We present our results in Sec. 4, including the implications for direct detection experiments searching for DM-electron scattering. We conclude in Sec. 5. Details of our solutions to the Boltzmann equations are given in Appendices A and B.
2 Freeze-in with a light dark photon mediator DM lighter than 10 GeV is strongly constrained by energy injection constraints from the cosmic microwave background (CMB) [69,70]. The freeze-in mechanism is a generic way to avoid excess late-time DM annihilation, since there is never enough DM for the annihilation process to be active, and thus there is no need to appeal to a velocity-suppressed annihilation cross section which implies constraints on the spin and parity of the DM or mediator.

Benchmark dark photon model
A standard benchmark model which realizes the "traditional" freeze-in scenario contains Dirac fermion DM χ that interacts with a dark photon, Z D , with dark gauge coupling g D . The dark photon communicates with the SM through kinetic mixing with SM hypercharge [71,72], (2.1) We take the dark photon to have a small but non-zero mass m Z D , which for simplicity we consider to arise from a Stückelberg mechanism [73,74]. We will typically be interested in χ masses below 1 GeV.
In the regime where DM never attains thermal equilibrium with the SM, the portal coupling is very small, and the couplings of Z D , Z to SM fermions f and DM are to an excellent approximation given by: while the coupling of the Z boson to SM fermions is to leading order unaltered. This DM model can thus be described at the Lagrangian level by four parameters, which we will take to be α D , , m χ and m Z D where α D = g 2 D /(4π). However, when m Z D m χ , the regime of greatest interest for direct detection, the dark cosmological history as well as the resulting direct detection signals are largely insensitive to the specific value of the dark photon mass.
In the limit m Z D m Z of interest, the couplings of the Z D reduce to the simpler expressions g Z D f ≈ − eQ f , g Z D χ ≈ g D and g Zχ ≈ g D tan θ W .

Traditional freeze-in: review
With the mass hierarchy m Z D m χ , and the absence of any additional dark sector species, freeze-in is UV-insensitive in the following sense. DM is produced from annihilation of SM particles in the thermal plasma, SM + SM → χ +χ. The DM abundance grows monotonically with time, reaching a maximum once the temperature drops below either m χ or m e , whichever is larger: in the former case, DM production becomes Boltzmann-suppressed at T = m χ , and in the latter case, the abundance of SM particles coupling to the dark photon becomes Boltzmann-suppressed after positron annihilation and plasmon decays become more important [75][76][77] for the production of DM. The lightness of the dark photon is crucial here, allowing the s-channel annihilation to be dominated by the lightest mass scale (or lowest temperature) in the problem, rather than (say) by the mass of a new heavy mediator. Since DM production originates from the thermal SM plasma and most of the DM is produced at late times, this mechanism is insensitive to the unknown early history of our universe. The parameters required to achieve the observed relic abundance are [8] Including DM production through the plasmon channel decreases the couplings required to achieve the freeze-in relic density by up to an order of magnitude for m χ < m e [75,77], but the above estimates are sufficient since we primarily focus on the regime m χ > m e in this paper. As mentioned in the Introduction, the hidden UV sensitivity in this nominally UVinsensitive scenario is the choice of initial DM abundance, which is customarily taken to be zero. We refer to this scenario as "traditional" freeze-in. We now show that the only effect of a nonzero initial χ abundance is a simple offset in the late-time relic abundance, rendering this residual UV-sensitivity rather trivial. Assuming for simplicity that e + e − → χχ is the only process which populates the dark sector aside from any primordial abundance, the Boltzmann equation relating the DM abundance n χ to the electron abundance n e iṡ n χ + 3Hn χ = 2 σv n 2 e , where H ≡ȧ/a is the Hubble parameter and σv is the thermally-averaged annihilation cross section. Changing variables to the comoving yield Y χ = n χ /s (where s is the entropy density) and to the dimensionless time variable x = m χ /T , we have Consider first the regime where x 1, and assume that m χ > m e . In that case, for annihilation through a light Z D , σv ∼ πα D 2 α/T 2 ∝ x 2 since there are no Boltzmann suppressions or kinematic endpoints. Electrons are always relativistic, so n e ∝ T 3 ∝ 1/x 3 , and similarly, s ∝ T 3 ∝ 1/x 3 . By assumption, freeze-in is taking place during radiation domination, so H ∝ 1/x 2 . Collecting the x dependence, we find which has the trivial solution Y χ = Y 0 + const. × x. So the effect of a primordial abundance Y 0 is simply to offset the linear growth of Y χ , which will push the slope of Y χ (x) to smaller values (in other words, smaller couplings 2 α D ) to achieve the same DM abundance when freeze-in turns off around x ∼ 1.

Self-consistency of traditional freeze-in
The above analysis has made an implicit assumption that dark particles, once produced from the SM, subsequently free-stream without further interaction. In some portions of our fourdimensional parameter space, this assumption does indeed hold. In other parts of parameter space, however, interactions of the injected dark matter particles, both with each other and with the light dark mediator particle, are important and can lead to sizable impacts on the DM phase-space distribution or even the final relic abundance. We estimate here the regime of validity of this traditional freeze-in treatment by requiring that a frozen-in particle does not undergo further scattering after production. Writing the elastic scattering rate between frozen-in DM particles as n χ σv el , we have from Eq. (2.6) that n χ ≈ n 2 f σv /H where f is an SM fermion. Therefore, by simply imposing that the elastic scattering rate be smaller than the Hubble rate, we have the condition Interestingly, the rate is given by the number density of the particles annihilating into DM and by an effective cross section which is the geometric mean of the scattering and the annihilation cross sections. We estimate the thermally-averaged scattering cross section as σv el ≈ πα 2 D /T 2 , and thus n f σv ≈ 2αα D 2 /(3π 3 )T 2 for T m χ . The bound for the combination of couplings α D and for negligible elastic scattering would be where g * ρ is the effective number of relativistic degrees of freedom related to the energy density.
To get the correct freeze-in abundance (Eq. (2.4)), we need 2 α D ≈ 3.5 × 10 −24 . Therefore, in this case our estimation for the maximum value of α D self-consistent with the traditional freeze-in mechanism is For other values of and α D , both self-scattering and self-annihilations are important, and the initial condition dependence becomes more involved. We turn to this region of parameter space in the following section.

Freeze-in into a pre-existing thermal bath
A more interesting type of UV sensitivity, with rich accompanying dynamics, arises when there is a pre-existing population of a dark sector containing χ, rather than simply a non-interacting primordial DM abundance. As DM is injected into this dark thermal bath, it will exchange kinetic energy with bath particles. Further, annihilations within the dark sector may begin to deplete the DM abundance, in sharp contrast to the monotonic increase in the traditional freeze-in scenario described above. We parameterize the initial conditions on this dark thermal bath through an initial temperature ratio ξ i ≡T i /T i . Our regime of interest is ξ i < 1, and therefore the Hubble rate is always dominated by the SM energy density, H(T,T ) ≈ H(T ).

Boltzmann equations
For a kinetically mixed Z D , the dominant source of energy injection into the hidden sector is through DM pair production. Since this injection can easily occur after DM has already departed from full chemical equilibrium, it is important to track how much of this energy is converted into the shared dark sector temperatureT and how much remains sequestered as rest mass. In other words, the energy density of the hidden sector, ρ HS = ρ Z D + 2ρ χ , as well as the number density of DM, n DM = 2n χ , are determined by the DM chemical potential µ as well as the hidden sector temperatureT . 1 The corresponding Boltzmann equations can be written aṡ where the sums run over SM fermions f and P HS = P Z D + 2P χ is the pressure of the hidden sector. The collision terms appearing in Eq. (3.1), C ρ ff →χχ (T ) = n 2 f (T ) σvE and C ρ Z→χχ (T ) = ΓE Z n Z (T ), govern the injection of energy into the HS from DM pair production, where E = E 1 + E 2 (E = E Z ) is the total energy of the annihilating fermions (decaying Z boson). In the regime of primary interest to us, the first term, describing production from SM fermion annihilations, dominates over the second term, which indicates the contribution from Z decays.
The collision terms appearing in Eq. (3.2) include the effect of DM annihilations within the hidden sector as well as the injection of DM from the SM. The specific expressions for the various thermally-averaged quantities appearing in the collision terms are given in Appendix A.5. The analogous Boltzmann equation for the SM temperature (including the effect of reverse annihilations), along with the Friedmann equation giving the dependence of H on T andT , provide a closed system of equations. We solve this set of equations numerically using the dimensionless time variable x = m χ /T with initial condition x i = 10 −2 ξ i . This initial condition defines the initial temperature ratio ξ i at the SM temperature T i = 10 2 m χ /ξ i , which ensuresT i = 10 2 m χ for all values of ξ i , and thus makes sure we set initial conditions early enough to capture the correct DM evolution for all cases.
These Boltzmann equations have made one major assumption: that the DM number density (and thus energy density and pressure) can be described entirely in terms of µ and T , or in other words, that DM can always be taken to be in kinetic equilibrium with the mediator bath. This is the opposite limit from traditional freeze-in, where after production the DM phase-space distribution evolves only through redshifting. The description in terms of µ andT is valid when DM produced via freeze-in rapidly reaches kinetic equilibrium with the dark radiation bath, which holds over the parameter space of primary interest to us; we demonstrate the self-consistency of this assumption in Sec. 3.2.
We can gain some intuition about this system of equations by first considering the situations where the dark sector is in internal chemical equilibrium, in which case Eq. (3.1) for ρ HS is the only necessary equation to solve. In this case the dark sector temperatureT evolves non-adiabatically with scale factor once the rate of energy injection from the SM becomes comparable to the rate of energy dilution owing to the expansion of the universe [59,61,63,64,78]. Once the energy injection from the SM shuts off, the energy density in both sectors resumes adiabatic evolution. Thus during the time that the dark sector is in chemical equilibrium, the temperature evolution during the non-adiabatic period, which we will refer to as the "leak-in" phase for clarity, follows a cosmological attractor solutionT LI (a) [64]: givenT T and a collision term C E (T ) ∝ 2 α D describing the rate of energy transfer into the HS,T LI (a) is entirely fixed in terms of the SM temperature, withT LI (a) ∝ ( 2 α D ) 1/4 . When C E ∝ T 5 , as is generic in the absence of mass thresholds, the resulting leak-in solution givesT LI (a) ∝ a −3/4 . Hidden sectors withT (a i ) >T LI (a i ) evolve adiabatically untilT (a) =T LI (a) and subsequently follow the leak-in solution, while hidden sectors withT (a i ) <T LI (a i ) see their temperature rapidly rise up to the attractor solution. The approximate scaling of the attractor solution, normalized to the SM temperature and written in terms of temperature instead of scale factor for future convenience, is where α = e 2 /(4π) is the QED coupling. The existence of this IR-dominated attractor solution helps mitigate the sensitivity of the DM relic abundance to the initial value of ξ i , since sectors with ξ i < ξ LI (T i ) will trend toward to the attractor temperature ratio ξ LI (T ). However as chemical equilibrium is lost within the dark sector, it is necessary to keep more careful track of how much the energy injected from the SM is distributed. The system will leave the attractor solution once any of the following conditions are met: (i) the energy injection from the SM shuts off; (ii) the HS departs from chemical equilibrium; (iii) the energy density in the HS is dominated by matter, rather than radiation. To understand in detail which of these conditions is most relevant for any given parameter point, we need to numerically solve the full system described by Eqs. (3.1) and (3.2), to which we turn in the next section.
Finally, sufficiently large portal couplings will thermalize the dark sector with the SM.
In other words, at a sufficiently large value of , the dark sector reachesT = T for a given α D . The attractor solution gives a quick way to estimate when thermalization occurs. On the attractor, the dark temperature is given byT 4 constant. Thus settingT = T = CM Pl 2 α D lets us estimate when the two sectors thermalize.
We are interested in the temperature range T > m χ , and therefore the value of at which the two sectors thermalize is

Kinetic equilibration
The Boltzmann equations given in Eqs. (3.1)-(3.2) are a good description of the system as long as the DM produced from out-of-equilibrium interactions with the SM rapidly reach kinetic equilibrium with the dark radiation bath. A χχ pair injected into a dark thermal bath of temperatureT can interact with both the DM and dark photons within the bath. Kinetic equilibrium can be obtained through scattering of injected DM with bath DM particles via a t-channel Z D , as well as the Compton scattering of injected DM from a Z D in the bath.
The injected DM can also approach chemical equilibrium through annihilating with bath particles, via the t-channel process χχ → Z D Z D . In the regime of interest T T , the Hubble rate is determined by the SM temperature, meaning that H ∝ T 2 /M pl . To attain kinetic equilibrium, the momentum loss rate of the injected DM due to scattering with some particle in the pre-existing dark thermal bath (χ,χ or Z D ) needs to be greater that the Hubble rate, i.e., where we have defined the fractional momentum loss rate with respect to the momentum p(T ) of an injected DM particle in a Lorentz-invariant way. To compute this rate we derive new exacts result for collision terms describing the scattering of particles at two different temperatures, given in Apps. A.1-A.3. First consider the case whenT evolves adiabatically, and therefore ξ is constant (up to mass thresholds). Fig. 1 shows the minimum values of α D for which the assumption of rapid kinetic equilibrium is satisfied for a given fixed ξ (solid lines). Notice as the hidden temperature gets closer to the SM temperature, smaller values for α D are needed to obtain rapid kinetic equilibrium in the dark sector. This can be understood from the fact that as the hidden temperature increases, the number density of bath particles increases as well, giving higher interaction rates. On the other hand, when the hidden temperature is significantly less that the SM temperature, there will be fewer interactions and a bigger interaction coupling is needed for the injected χ to efficiently lose its momentum. Finally, the gray dashed line shows conservative constraints on α D coming from measurements of halo ellipticities [81,82] or relaxation of the halo profiles of galaxy groups and clusters [83] (see Ref. [80] for a review of self-interacting DM constraints). This fixed-ξ estimate can be overly conservative, however, depending on the value of , as it neglects the effect of energy injection from the SM on the dark temperature. In the right panel of Fig. 1 we show the evolution ofT for a range of initial ξ i and compare to the attractor solution corresponding to a particular α D , pair (also used in Fig. 2 below). The larger two initial temperature ratios (red and blue lines) begin above the attractor solution (dotted black) and redshift adiabatically down, while the initially underabundant purple curve rises up rapidly to the attractor. Meanwhile the green curve redshifts down until it meets the attractor, after which it follows the attractor solution. The effect of the QCD phase transition is visible at x ∼ 0.07, where the approximate attractor solution does not account for this effect. Neglecting the SM energy injection is thus an excellent approximation for the red and blue lines but underestimates the HS temperature and therefore the scattering rate for the green and especially the purple lines, for which ξ < ξ LI (a; , α D ) for some a. The impact of this non-adiabatic evolution on the requirement of kinetic equilibration is illustrated with the dashed purple line in Fig. 1, which shows the minimum values of α D that give rapid kinetic equilibration, given the attractor solution corresponding to = 2×10 −8 .
Further details about the calculation of kinetic equilibration are given in Appendix A.4. For the values of (m χ , , α D ) of interest in this work, rapid kinetic equilibration is a good approximation in a substantial portion of parameter space, and in particular the portion of parameter space that displays interesting dependence on initial conditions.

Results
The DM number density is obtained after solving the system of equations (3.1) and (3.2). To develop some intuition for the strength of the couplings needed to obtain the correct DM relic abundance, we first explore the parameter space as a function of the initial temperature ratio ξ i . We show the results in Fig. 2 for m χ = 10 MeV and different initial temperature ratios. The left panel in Fig. 2 shows contours of Ω χ (normalized to the observed DM relic density) in the α D -plane. This plot illustrates two distinct regimes at small coupling (bottom left corner): I. For small ξ i (short-dashed curves), at small couplings there is not enough DM in the hidden sector to achieve the required relic abundance through hidden-sector freeze-out alone, and instead the relic abundance is obtained through freeze-in, which implies a minimum for a given α D .
II. For large ξ i (solid curves), obtaining the observed relic abundance is possible for arbitrarily small values of , since the DM can freeze out entirely within the hidden sector, decoupled from the SM.
We have checked that the approximation of rapid kinetic equilibrium, the conditions for which can be seen from Fig. 1, holds for all of the parameter points shown in colored points (curves) in the left (right) panel of this figure, except the pink point (curve); the brown point (curve) lies at the boundary of the rapidly equilibrated region of parameter space along the freeze-in line. For sufficiently large couplings, contours for different values of ξ i converge on the attractor solution described in Sec. 3. At these larger couplings, there are two qualitatively different scenarios for achieving the correct relic abundance, regardless of the initial temperature ratio. For above the gray dashed line, the hidden sector thermalizes with the SM and freeze-out obtains in the traditional way. For between the dotted blue and dashed gray lines, DM can obtain the correct relic abundance through leak-in (i.e., one phase of freeze-out during a period of non-adiabatic temperature evolution) and/or reannihilation (i.e., two distinct phases of freezeout). The right panel of Fig. 2 shows the evolution of the DM yield for the colored points marked in the left panel, showing the transition from freeze-out to leak-in/reannihilation 2 to freeze-in for ξ i = 10 −3 . The existence of Regime II demonstrates the UV sensitivity of freeze-in with a light mediator: these secluded freeze-out solutions are available only for some initial values of ξ i , and the specific value of α D that yields the correct relic abundance through secluded freezeout depends on the specific value of ξ i . At sufficiently small ξ i , however, secluded freeze-out does not occur, and the relic abundance is instead dominated by freeze-in processes. We can understand the division between Regimes I and II straightforwardly by looking at the initial DM abundance as a function of ξ i . First, let us define the comoving DM number density as which is equivalent to the standard, and more familiar, form Ω DM h 2 = 0.12. In the case of interest where the DM chemical potential is zero and its temperatureT is different from the SM temperature T , we have where ξ =T /T ,x = m χ /T , g * s counts the effective relativistic degrees of freedom contributing to the entropy density, g χ = 4 for a Dirac fermion, and K 2 is a modified Bessel function. Therefore, the initial DM yield (x 1) can be expressed as Notice that the initial yield is independent of the DM mass, which is just the statement that DM is relativistic forT > m χ . As a result, the initial yield is entirely fixed by the initial hidden-to-SM temperature ratio, ξ i . On the other hand, the late-time DM yield, Eq. (4.1), only depends on the DM mass. This leads to two possibilities: Initially Underproduced g * s = 100 g * s = 10 Figure 3. Initial DM density as a function of the initial hidden-to-SM temperature ratio ξ i and DM mass m χ . The boundary between initial overproduction and initial underproduction (which depends on g * s ) defines the parameter space for which the freeze-out or freeze-in mechanisms are viable.
• If Y i χ > Y DM , there is too much DM initially and DM needs to annihilate to reproduce the correct relic density, which is accomplished by freeze-out. 3 • If Y i χ < Y DM , there is too little DM initially and the DM abundance needs to build up over time, which is accomplished by freeze-in.
We show in Fig. 3 the values of ξ i and m χ for which the freeze-out or freeze-in mechanism is needed. We use g χ = 4 and we show two representative values of g * s : 100, when T i 200 MeV (i.e. above the QCD phase transition) and 10 when T i 20 MeV (below the QCD phase transition). Finally, we define the critical temperature ratio ξ * i , such that Y i χ = Y DM , meaning the initial yield precisely coincides with the observed relic abundance.
These two possibilities (underproduced vs. overproduced) map onto Regimes I and II discussed above. However, due to the attractor solution, the "true" initial temperature ratio at early times will not be ξ i but rather ξ LI (T i ), so long as ξ i < ξ LI (T i ), where the temperature evolution of the attractor solution is given in Eq. (3.3). Thus we need to check whether the boundary between the initially under-vs. overproduced regimes is robust against the attractor solution for values of α D and along the freeze-in curve. An initial temperature ratio ξ i < ξ * i will remain in the underproduced region so long as ξ LI (T i ) < ξ * i where T i is the temperature at which ξ i is defined. The critical value ξ * at the overproduced/underproduced boundary from Fig. 3 is ξ * ∼ 10 −2 with a weak dependence on the DM mass. Along the freeze-in trajectories in Fig. 2 a condition which is clearly required in order to have freeze-in of DM with mass greater than 1 MeV. Therefore, along the freeze-in curve seen in the left panel of Fig. 2, the product of couplings α D 2 is too small for the corresponding attractor solutions to raise these parameter points from the underproduced region to the overproduced region. In other words, along the freeze-in curve in α D -space, the initial DM production regime found in Fig. 3 is robust against the non-adiabatic evolution ofT , which gives us a simple way to understand the small-behavior of the curves corresponding to different ξ i in the left panel of Fig. 2. However, as either α D or increases, the temperature ratio given by the attractor solution eventually yields too much energy density in the hidden sector, necessitating a period of reannihilation to obtain the correct late-time relic abundance. Such a trajectory is illustrated by the yellow point (contour) in Fig. 2, left (right). This parameter point demonstrates that at large couplings, a hidden sector that would yield an underabundance of DM in the absence of thermalizing interactions in the hidden sector can develop an overabundance.
Finally, we also show with the orange dashed line in Fig. 2 the maximum value of for a given α D for which a traditional freeze-in solution is self-consistent, as given in Eq. (2.9). Meanwhile the brown point in the same figure is at the boundary of the region where rapid kinetic equilibration is a good approximation along the freeze-in line (α D 10 −8 , see Fig. 1). This leaves a notable portion of the α D -plane which can be handled self-consistently in either the non-interacting regime of Sec. 2 or the rapidly thermalizing regime of Sec. 3, depending on the presence or absence of a thermalized dark sector. This should not be a surprise: energetic dark particles produced from the SM plasma are underabundant compared to the thermal number abundance expected for the same ρ HS . Thus the rates for self-interactions of these frozen-in particles are small in comparison to the situation where an energetic DM particle with E χ ∼ T SM scatters off a colder thermal bath of dark particles, even when ρ HS is the same between the two scenarios. This can remain true even if the initial energy density in the dark sector is small, because the dark sector temperature will rapidly approach the attractor solution. Said another way, there are regions of -α D space where a minimal dark sector with zero initial abundance may undergo negligible self-scattering, but where a thermal initial population makes the approximation of rapid kinetic equilibration safe. This is yet another source of UV sensitivity that goes beyond the dependence on ξ i demonstrated here. For instance, the DM produced in the model represented by the brown dot could have very different predictions for its phase space distribution (as well as the number of relic dark mediators), depending on its cosmic history.
Depending on the initial temperature of the hidden sector, the couplings required to achieve the correct DM yield may be considerably smaller than those implied by traditional freeze-in. Indeed, if we drop the requirement of thermal contact with the SM, the kinetic mixing can vanish and the DM can still achieve the observed relic abundance in a decoupled hidden sector. However, if we take some nonzero amount of thermal contact to be a definition of glaciation, we can quantify the UV sensitivity of this scenario in terms of the initial temperature ratio ξ i . As shown in Fig. 2 (left), if we fine-tune ξ i to the critical value ξ * i exactly on the overproduced/underproduced boundary of Fig. 3, we end up with the correct relic density  by construction even for = α D = 0. 4 For all ξ i > ξ * i , decoupled hidden-sector freeze-out with = 0 is possible, and for all ξ i < ξ * i , sufficiently small α D will permit a traditional freeze-in solution. In this sense, glaciation is UV-sensitive for ξ i 3 × 10 −2 (MeV/m χ ) 1/3 . Interestingly, for ξ i < ξ * i , there are always points in the -α D plane below the traditional freeze-in curve, arising from a period of late-time leak-in supplemented by freeze-in (brown point and yield curve in Fig. 2). To account for this expanded parameter space, we propose that the freeze-in curve should be expanded to a "glaciation band" to account for this initial condition sensitivity of freeze-in; we explore the implications of this fact for direct detection experiments below.

Implications for direct detection
A key feature of the freeze-in scenario is the excellent discovery potential at terrestrial DMelectron scattering experiments, which can take advantage of the low velocity of the DM and the long-range nature of the light Z D mediator to make up for the small couplings required to match the observed relic abundance. Experimental results are typically expressed in terms of a fiducial DM-electron cross section, where µ χe is the DM-electron reduced mass and for simplicity we have assumed m Z D αm e in our choice of normalization. Since the dependence of σ e on the hidden sector couplings is given by 2 α D , we define the glaciation band for each DM mass as follows: 4 Due to the precise fine-tuning required, the ξ * i curve illustrated in Fig. 2 saturates at a finite value of αD due to accumulated rounding error in the numerical solutions to the Boltzmann equations.
• Upper boundary: 2 α D equal to the value at the intersection of the Ω χ /Ω DM = 1 contour with the thermalization contour (gray dashed in Fig. 2, left).
• Lower boundary: 2 α D equal to the minimum value achieved over all contours of Ω χ /Ω DM = 1 defined by ξ i < ξ * i .
By construction, the traditional freeze-in curve is enclosed in the glaciation band. The lower boundary of the glaciation band encompasses the region of parameter space where the relic abundance is dominated by freeze-in processes after accounting for a range of initial conditions. Meanwhile, the upper region of the glaciation band is UV-insensitive, as the attractor solution erases dependence on the initial temperature ratio. We show the glaciation band in Fig. 4 (left), along with constraints from SENSEI [46,51,54] which are the strongest for DM scattering through a light mediator in this region of parameter space. We see that direct detection has already ruled out large parts of the glaciation parameter space (see also [64,85]). Indeed, this can be visualized in the -α D plane as follows. For a given value of m χ , direct detection sets an upper bound on 2 α D , which is a line in the -α D plane (dashed lines in Fig. 4, right). The point at which this upper bound intersects the Ω χ /Ω DM = 1 contour represents the boundary of the equivalent exclusion region: any larger values of and α D are ruled out, and thus the remainder of the relic density contour for that m χ is ruled out. We show the result of this procedure in Fig. 4 (right), and see that direct detection constraints already rule out significant portions of the leak-in scenario (independent of ξ i ). The projected reach of Oscura [86] will cover the entire glaciation band for m χ > 1 MeV; if a positive signal is found at σ e below the traditional freeze-in line, that could either indicate a subdominant component of DM, or in the most optimistic case would offer the tantalizing possibility of directly probing the thermal history of a dark sector with a light mediator.

Conclusions
Models where DM freezes in through out-of-equilibrium production from the SM have emerged as important targets for developing terrestrial tests of (sub)-GeV-scale dark sectors. Carefully considering the predictions of freeze-in models is thus vital for understanding the information about the early universe that current and upcoming experiments will provide. Since DM never reaches thermal equilibrium in freeze-in models, there is necessarily some residual dependence on initial conditions in their predictions. In the case of "traditional" freeze-in, where DM does not interact after its production, this sensitivity is relatively minimal provided the DM-SM interaction is renormalizable, amounting to a constant and generically small offset of the total DM yield given specific couplings of the DM to the SM.
Another, richer scenario, of high experimental interest, is the case where DM interacts with the SM via a light mediator. Our results here have shown that models where DM freezes in through a kinetically-mixed light mediator can have a much more dramatic dependence on the initial conditions specified for the dark sector than do more traditional freeze-in scenarios.
Using the common and minimal reference model of Dirac fermion dark matter interacting with the SM via a kinetically-mixed dark photon, we have demonstrated a nontrivial dependence of the final DM yield on the initial conditions for the dark photon as well as the dark matter. We have shown that it is self-consistent to take the dark sector to be in internal kinetic equilibrium throughout the formation of the DM relic abundance in a large region of interest, and we parameterize the initial conditions for the dark sector in terms of ξ i , the initial ratio of dark to SM temperatures.
For sufficiently large values of the dark gauge coupling α D and the kinetic mixing parameter , the energy injection from the SM is large enough to overwhelm variations in the initial population density, meaning that the DM relic abundance is insensitive to variations in initial conditions.
However, for smaller values of α D and , the DM evolution within the hidden sector depends in detail on the initial population. In this region, the final DM relic abundance depends on the initial conditions, with different possible outcomes: if the temperature ratio is larger than a critical value ξ * i , the evolution of the number density is set by freeze-out in the hidden sector, but for smaller initial temperatures the final number density is determined by late freeze-in-like processes from the SM. Therefore, the initial population as well as the values of α D and determine the late-time abundance. In this region the predicted DM relic abundance exhibits a qualitatively new form of UV sensitivity.
We have pointed out that the freeze-in curve stops being a self-consistent experimental target for sufficiently large values of 2 α D and have clarified what happens to hidden sectors with couplings in this regime. We have shown that a sizeable portion of the resulting "glaciation band" is UV-insensitive, in the sense that variations in the initial conditions do not impact the final relic abundance obtained for a given parameter point. However, for the parameter space near the traditional freeze-in target, predictions for the final relic abundance do depend on the initial population of the hidden sector. Thus we are able to identify and quantify the residual UV dependence of the freeze-in scenario with light mediators, and clarify its consequences for experiments. We define the bottom of the glaciation band as the smallest SM-DM cross section that gives rise to DM through freeze-in processes from the SM, rather than through hidden sector freeze-out, and provide a simple prescription to compute this quantity. This glaciation band constitutes a robust and well-motivated target for near future DM-electron direct detection experiments such as Oscura.
Finally, we have provided a simple demonstration that the UV sensitivity of freeze-in with a light mediator goes beyond the dependence on a finite initial temperature. Since a frozen-in DM particle will scatter much more rapidly off of a cold particle from a pre-existing thermal population than off of another energetic frozen-in DM particle, there are regions of parameter space where both a non-interacting freeze-in solution and a kinetically-equilibrated glaciation solution can be self-consistent. In this region the DM phase space distribution will depend on initial conditions even if the DM yield does not.
In the limit of small ξ i , one may also start to ask whether the hidden sector would have time, in a given cosmological scenario, to approach internal kinetic equilibrium. The approach to internal thermal equilibrium can take an appreciable amount of time, even for dark sectors containing parametrically light mediators [64,87,88]. Such questions are particularly acute for the small values of α D needed to evade constraints on DM self-interactions for sub-MeV DM. For DM with mass below an MeV, after imposing constraints on DM self-interactions, the approximation of rapid kinetic equilibrium used here is applicable for a limited range of relatively large ξ i . However, the presence of a pre-existing dark sector population, whether equilibrated or not, will generically affect the DM phase space distribution in this mass range as well. Understanding the impact of scattering in this low-mass region is of particular interest, as the detailed shape of the phase-space distribution of light dark matter can be important for cosmological observables [77,89,90].
Determining the evolution of the DM phase-space distribution in the general out-ofequilibrium case requires solving the full Boltzmann hierarchy. Some work in this direction was recently done in [91] for a model with a heavy mediator and a constant matrix element. We expect that this task will be substantially harder for hidden sector with a light mediator, owing to the additional species that needs to be tracked and the need to carefully treat small momentum-transfer scatterings. However, freeze-in through a kinetically-mixed light dark photon is one of a very small number of cosmologically-viable models for sub-MeV DM, and thus this result is well worth pursuing.

A Collision terms for species at different temperatures
In this appendix, we derive the collision terms for the number density, energy transfer and momentum transfer rates. Unlike the standard case [92,93], where the species share the same temperature, we generalize the argument and work out the rates for the cases when the initial state particles have different temperatures. Therefore, our general focus will be on processes of the type 1 + 2 → 3 + 4, where particles 1 and 2 have different bath temperatures T andT respectively, i.e. T 1 = T = T 2 =T .
We will work under the Maxwell-Boltzmann approximation. Therefore, let us first review the relevant thermodynamic equations for particles in thermal equilibrium at temperature T that follow a Maxwell-Boltzmann distribution, namely f = e −(E−µ)/T , where E is the energy of the particle and µ its chemical potential. This leads to the expressions for the number density, energy density, and pressure, given by where g gives the internal degrees of freedom and K i are the modified Bessel functions.

A.1 Number density
We start with the derivation of the number density collision operator for particle 1, which reads 3 is the Lorentz-invariant phase space element. For simplicity, here we only consider the collision term governing the forward scattering and neglect the chemical potential, although including these effects is straightforward. The integral over two of the phase space differentials can be written in terms of the cross section σ (s) as where the two-body kinematic function λ(s, m 1 , m 2 ) is Then, the collision term can be written as follows C n 1 2→3 4 (T,T ) = −2g 1 g 2 dΠ 1 dΠ 2 λ 1/2 (s, m 1 , m 2 ) σ(s) f eq 1 (T ) f eq 2 (T ) (A.5) and the remaining phase space differentials are, in terms of the lab energies and incident angle, Here it is convenient to switch the integration variables to where γ and β are the boost factor and velocity, respectively. The Jacobian for this transformation is 8) and the integration limits are Putting everything together, the collision operator is (A.10) Using the dimensionless variable x i = m i /T i and focusing only on the integration over γ 1 and γ 2 , we have (A.11) Then, with the help of the rapidity γ 1 = cosh(w 1 ), β 1 γ 1 = sinh(w 1 ), γ r = cosh(w r ) and β r γ r = sinh(w r ), we can perform the following integration over w 1 Using cosh(θ 1 ) = x 1 and cosh(θ 2 ) = x 2 , the arguments of the exponential can be written as cosh(θ 1 ) cosh(w 1 ) + cosh(θ 2 ) cosh(w 1 ± w r ) =s cosh(w 1 ± φ) , (A.13) wheres = (x 2 1 + 2x 1 x 2 γ r + x 2 2 ) 1/2 and φ = sinh −1 x 2 sinh(w r ) s . (A.14) Shifting the integration variable t = w 1 ± φ we have Finally, we find that the collision operator is wheres min = x 1 + x 2 and σ(s) is evaluated at to perform the integral overs. Heres 2 plays a role reminiscent of the Mandelstam variable s but now with dependence on the bath temperatures. It is important to note that for elastic scattering processes, the collision term conserves particle number, i.e. C n 1 2↔1 2 (T,T ) = 0.

A.2 Energy transfer
The same calculation can be done for the collision operator describing the energy transfer rate for particle 1 with energy E = E 1 as (A.19) Then after performing a similar calculation as in the number density case, we have that Unlike the number density operator, the energy transfer for elastic scattering processes does not vanish, i.e., C ρ 1 2↔1 2 (T,T ) = 0.

A.3 Momentum transfer
Finally, similar to the energy transfer rate, we can define the momentum loss rate of an injected particle with temperature T and momentum p 1 through scattering off of a second particle with temperatureT and momentum p 2 by considering the average of the quantity where p 3 the momentum of the injected particle after the collision and in the second equality p 1 and the scattering angle are given in the center-of-mass frame. This expression for the momentum transfer-squared is just the Mandelstam variable −t in the center of mass frame. Then, the collision operator describing the momentum loss rate can be written as C p 1 2→3 4 (T,T ) = n 1eq (T )n 2eq (T ) σv∆p 2 (A.21) where σ T is the transfer cross section defined by σ T = dσ dΩ (1 − cos θ)dΩ [94,95]. Thus, using the previous results, the thermally-averaged momentum loss rate can be obtained by where p 2 1 , the average momentum-squared of an injected DM particle, is

A.4 Rapid kinetic equilibration
As explained in the text, we can now estimate when rapid kinetic equilibrium holds by requiring Γ p loss H. Our goal here is to verify that our choice of initial conditions for solving the Boltzmann equations is robust: if rapid kinetic equilibrium is obtained at some point while the DM is relativistic, it is maintained throughout all of the evolution of the DM number density. We are interested in the elastic scattering processes -namely Compton, Bhabha, and Møller -explicit cross-sections for which are given in Appendix B. Fig. 5 shows the momentum loss rate from each of these processes separately as well as the total, together with the Hubble parameter. For this example, we have chosen parameter values corresponding to the yellow point in Fig. 2 for which the value of 2 α D is small, meaning the energy transfer is small too and ξ i = 10 −3 . In this case, as can been seen in the right panel of Fig. 2, the hidden temperature starts evolving non-adiabatically on the attractor solution. For this parameter point, over the range of hidden sector temperatures for which we solve the Boltzmann equations (starting at x i = 10 −2 ξ i ), the total Γ p loss is always larger than the Hubble rate. Therefore, there is always self-consistency when the equations are solved numerically. To gain some intuition, let us consider the limiting case of T m χ , where we can obtain an approximate analytic expression to the momentum loss rate due to Bhabha scattering, as and for Compton, as where γ E is Euler's constant. If the hidden sector evolves adiabatically, ξ is a constant, in contrast to the non-adiabatic case, where the attractor solution is well-approximated by Eq. (3.3). While in the former case the rates scale as T , as generically expected, in the latter case the scaling is T 1/2 . As can be seen in Fig. 5 (black dotted lines), the semi-analytical estimates track the numerical solution perfectly. Thus, the approximate minimum value x min that ensures kinetic equilibrium will satisfy Γ Bhabha p loss (x min ) ≈ H (x min ), which gives x min ≈ 2 × 10 −5 4.3 × 10 −7 α D 5/3 2 × 10 −8 2/3 m χ 10MeV . (A.26) Moreover, we also need to ensure that kinetic equilibrium is maintained until the final DM number density has been achieved. For the parameter space considered here, we find that the Compton rate always preserves the kinetic equilibrium conditions for late times as shown in Fig. 5. The reason is that the Compton rate is not Boltzmann-suppressed, and furthermore at late times when the energy injection from the SM is negligible, the Compton rate has the same temperature scaling as the Hubble parameter, H ∝ T 2 . Therefore, once Hubble crosses the Compton rate from above, Compton dominates for all late times and kinetic equilibrium is maintained. where s min = max{4m 2 f , 4m 2 χ }. This recovers the well-known results from Ref. [92]. Similarly, using Eq. (A.20) the energy transfer rate reads where E = E 1 + E 2 and we have used σvE fi = 2 σvE 1 fi . Finally, we provide explicit formulae for the number and energy density rates for Z decays into DM, 29) where g Z = 3 gives the degrees of freedom of the Z boson.

B Cross sections
For reference, we present all of the 2 → 2 cross section and decay formulas we require in our Boltzmann equations. All cross sections here are summed, rather than averaged, over the final and initial states.

Decay of Z to DM, Z → χχ
The total decay width is where θ W is the weak mixing angle.
SM fermion annihilations to DM only through the dark photon, ff → χχ SM fermion annihilations to DM with Z D − Z contribution, ff → χχ Here we show the full annihilation cross section including the Z boson contribution, where the vector and axial couplings for a fermion f are C V = T 3 f − 2Q f sin 2 θ W , C A = T 3 f , g Z = e cos θ W sin θ W and Γ Z is the decay width of the Z boson.