Dynamic Freeze-In: Impact of Thermal Masses and Cosmological Phase Transitions on Dark Matter Production

The cosmological abundance of dark matter can be significantly influenced by the temperature dependence of particle masses and vacuum expectation values. We illustrate this point in three simple freeze-in models. The first one, which we call kinematically induced freeze-in, is based on the observation that the effective mass of a scalar temporarily becomes very small as the scalar potential undergoes a second order phase transition. This opens dark matter production channels that are otherwise forbidden. The second model we consider, dubbed vev-induced freeze-in, is a fermionic Higgs portal scenario. Its scalar sector is augmented compared to the Standard Model by an additional scalar singlet, $S$, which couples to dark matter and temporarily acquires a vacuum expectation value (a two-step phase transition or `vev flip-flop'). While $\langle S \rangle \neq 0$, the modified coupling structure in the scalar sector implies that dark matter production is significantly enhanced compared to the $\langle S \rangle = 0$ phases realised at very early times and again today. The third model, which we call mixing-induced freeze-in, is similar in spirit, but here it is the mixing of dark sector fermions, induced by non-zero $\langle S \rangle$, that temporarily boosts the dark matter production rate. For all three scenarios, we carefully dissect the evolution of the dark sector in the early Universe. We compute the DM relic abundance as a function of the model parameters, emphasising the importance of thermal corrections and the proper treatment of phase transitions in the calculation.


I. INTRODUCTION
"In order for the light to shine so brightly, the darkness must be present." This quote, attributed to Sir Francis Bacon, subsumes much of modern day cosmology. The Universe as we know it, with its abundance of bright galaxies, could not have formed without the presence of large amounts of dark matter (DM). DM drives the formation of structure; the gravitational collapse of primordial density fluctuations leads to dense objects like galaxy clusters, galaxies, and stars, with at least one of the latter harbouring life. Even though one of the lifeforms likes to describe itself as intelligent, it is still very much in the dark about dark matter, its origin and its nature.
For a long time, the best-motivated scenario to understand DM has been the Weakly Interacting Massive Particle (WIMP) scenario: DM particles are hypothesised to be heavy (m DM 100 GeV) and to have weak, but non-negligible interactions with Standard Model (SM) particles. In the very early Universe, these interactions keep the DM and SM sectors in kinetic and chemical equilibrium, until eventually Hubble expansion dilutes the DM density to the extent that DM annihilation into SM particles ceases. This freeze-out typically occurs at temperatures T where x ≡ m DM /T 25 [1].
While freeze-out is arguably still the leading scenario for explaining the DM abundance in the Universe, the lack of experimental evidence for DM particles, in spite of a vigorous, multi-pronged search program [2][3][4][5][6][7][8], motivates the study of alternative mechanisms [9][10][11][12][13][14][15][16][17][18][19][20]. Freeze-in models assume that the initial DM abundance after inflation was zero and that DM particles couple so weakly to the particles in the thermal bath that they never reach thermal equilibrium [17][18][19][20][21]. Consequently, DM annihilation does not determine the relic abundance. Instead, the observed abundance is the result of processes with DM in the final state, typically at x 1-5. The small coupling between DM and SM particles also implies a significantly greater challenge for all experimental probes of the nature of DM.
In this paper, we consider the impact of finite temperature effects on the DM abundance in freeze-in models. Such effects are manifold: first, particle masses receive corrections from thermal loops, implying that the kinematics of DM production is in general T -dependent. Certain production channels may be open during some epochs of cosmological history, but kinematically closed during others. Closely related to this effect is the T dependence of the effective scalar potential V eff , which implies that not only scalar masses (the second derivatives of V eff ), but also their vacuum expectation values (vevs) change with time in the early Universe. These vevs, in turn, will affect gauge boson and fermion masses (of course, gauge boson and fermion masses also receive direct corrections from self-energy diagrams evaluated at T = 0, although the contribution to fermion masses will be unimportant in the scenarios we discuss). The most interesting case is where the scalar potential develops several disjoint minima and transitions from one to another in a phase transition. The phase transition, for which the scalar vev is an order parameter, can be first order, second order, or a mere cross-over. The latter (and perhaps least interesting case) is realised in the SM [22][23][24][25][26].
In standard freeze-out scenarios, thermal effects are usually negligible since they are small at x 1, i.e., at temperatures much lower than the masses of the involved particles. For freeze-in at x 1-5, on the other hand, they can be large and have a decisive impact on the DM abundance. Demonstrating this with several examples is the main topic of this paper.
To focus on the important effects and avoid unnecessary complications, we consider a toy model consisting of the standard model, a new gauge singlet scalar, and one or two gauge singlet dark sector fermions. Since we focus on freeze-in, we imagine some couplings to be 1. As discussed below, these small couplings will often be technically natural (i.e., protected from large radiative corrections). In other cases, we assume a small coupling to simplify the analysis, but in these cases we do not expect a larger coupling to spoil the overall picture. The scenarios we consider could be motivated from a wide range of UV theories, including SUSY and extra dimensional scenarios [19]. Scenarios with gauge singlets at the weak scale are notoriously difficult to test at colliders and will not be ruled out at the LHC. The scenarios discussed in sections III and IV will be best probed through their O(1) Higgs portal couplings, but the full parameter space will not be probed until a 100 TeV collider has collected 3 ab −1 of integrated luminosity [27]. The scenario discussed in section II is even harder to exclude, due to the small Higgs portal coupling.
The impact of thermal effects in the early Universe on the DM abundance has been previously discussed for instance in [28,29], where it was argued that DM could be temporarily unstable in the early Universe, so that its abundance would be controlled by its decay rates and by the temperature of the phase transition that stabilises it. Similar thermal effects have also been considered in, e.g., [30,31].
The outline of the paper is as follows. In section II we discuss a scenario, dubbed "kinematically induced freeze-in", in which the kinematics of DM production is controlled by the T -dependent masses of a new scalar S and a new dark sector fermion ψ. DM particles χ freeze in through their couplings to S and ψ, but for most of cosmological history DM production via ψ → Sχ is kinematically forbidden. However, as S transitions from a phase with S = 0 to a phase with S = 0, its mass drops close to zero and DM production becomes kinematically allowed for a short period of time. We emphasise that, in a more economical version of this model, S could be replaced by the SM Higgs itself. In section III we consider an alternative freeze-in model -a variant of the fermionic Higgs portal scenario -in which the DM production rate in the dominant channels is proportional to S . We call this scenario "vev-induced freeze-in". There is a large region of parameter space where its scalar potential undergoes a two-step phase transition ("vev flip-flop"), i.e., the Universe starts out in a S = 0 phase, followed by an epoch with S = 0. The electroweak phase transition ends this epoch and reverts the Universe to S = 0 (but a nonzero vev for the SM Higgs). It is thus the two phase transitions that control the DM abundance today. In a third model, which will be the topic of section IV, S controls mixing between the DM particle and a second new fermion. This mixing, in turn, opens up production channels that are otherwise inaccessible, thus boosting freeze-in production. Hence we call this scenario "mixing induced freeze-in". We summarise and conclude in section V.

II. KINEMATICALLY INDUCED FREEZE-IN: TEMPERATURE-DEPENDENT MASSES AND THRESHOLDS
In thermal quantum field theory, particle masses can receive temperature-dependent corrections from self-energy diagrams and thus become functions of T themselves. For instance, the SM Higgs mass is m h (T = 0) = 125 GeV today, but was much larger in the very early Universe and close to zero around the time of the electroweak cross-over at T EW 165 GeV. In this section we discuss a scenario where the kinematics of the DM freeze-in rate are controlled by the mass of a new real scalar S.

II.1. Toy Model
We consider a simple toy model -a variant of the fermionic Higgs portal scenario -whose particle content is given in table I. Besides the real scalar S and the Dirac fermion χ, which is the DM candidate, we introduce a second new Dirac fermion ψ. All new particles are SM singlets, and ψ and χ are charged under a Z 2 symmetry. We remark already here that one could imagine a variation of the model in which S is replaced by the SM Higgs field itself. The relevant terms in the Lagrangian at dimension four are The first line of eq. (1) contains the standard kinetic terms and the fermion mass terms. In the second line of eq. (1), we identify the Yukawa couplings between S, χ and ψ. We assume y χ and y ψχ to be tiny to avoid full thermalisation of χ and thus allow for DM production via freeze-in rather than freeze-out. The smallness of y χ and y ψχ could be motivated, for instance, by extra-dimensional scenarios where χ could be localised far away from ψ and S along the fifth dimension. The coupling y ψ on the other hand, is assumed to be sizeable so that ψ and S remain in thermal equilibrium until T m ψ , m S . Alternatively, one can also introduce an extra particle -for instance a second new scalar S with negligible couplings to the DM particle χ, but appreciable couplings to ψ and to the SM sector -to achieve the same goal. In fact, in the numerical results shown below we will assume this second possibility because it simplifies the dynamics of the temperature-dependent effective scalar potential and opens up larger regions of parameter space than the vanilla scenario from eq. (1).
The first line of the scalar potential V (H, S) in eq. (2) contains the mass terms and quartic couplings for S and H. We assume µ 2 H , µ 2 S > 0, so that not only H, but also S, may obtain a vev at tree level. The second line of eq. (2) consists of a cubic coupling for S proportional to λ S3 and of the cubic (∝ λ p3 ) and quartic (∝ λ p4 ) Higgs portal couplings. We will assume that λ S3 , λ p3 , and λ p4 are 1. For λ S3 this is just a simplifying assumption that could be relaxed at the expense of unnecessarily complicating our analysis. The smallness of λ p3 and λ p4 could again be motivated in extra-dimensional scenarios by localising S and H far from each other along the fifth dimension. We hypothesise, however, that λ p4 is still large enough to keep S in thermal contact with the SM particles at temperatures T m S , when DM freeze-in happens. We find that these conditions are satisfied for λ p4 ∼ 10 −3 .
Where necessary, we will decompose H into its components according to H = G + , (h + iG 0 / √ 2), where h is the neutral CP even SM-like Higgs boson and G ± , G 0 are the Goldstone modes. Moreover, we will often use the definitions v S ≡ S and v H ≡ h for the vacuum expectation values of S and h.
Freeze-in of χ can proceed through the decays ψ → Sχ or S → χψ, (depending on the relative magnitude of the χ, S, and ψ masses), and through the 2 → 2 reactions ψS → χS and SS → χψ. The Feynman diagrams for these four processes are depicted in fig. 1. We will focus on masses such that m ψ + m χ > m S (T = 0) and m ψ − m χ < m S (T = 0). This implies that the decays ψ → Sχ and S → χψ are kinematically forbidden today. In the very early Universe, however, m S receives large thermal corrections ∝ T , which can lift its value above m ψ + m χ and thus open up the channel S → χψ. Later, around the time when S develops a non-zero vev, m S (T ) approaches zero and the channel ψ → Sχ becomes temporarily available. The decay rates and annihilation cross sections for the processes in fig. 1 are In the last two expressions, we have taken the limit m ψ m χ and m S 0. Moreover, we have set the width of ψ to Γ ψ = 0. In eq. (6), we have also set v S = 0 because the full expression is fairly lengthy, and we will see below that the channel SS → χψ is very subdominant when v S = 0. In our numerical analysis below, we of course use the full expressions.
It is important that, thanks to the non-zero S at low temperatures, S can decay through its Higgs portal coupling λ p4 (or also through λ p3 ). If this decay was not present, S would have a relic abundance that would be too large. For m S below the W and Z thresholds, the decay rate is Here, the sum runs over light fermions, f = e, µ, τ, u, d, s, c, b, and y f = √ 2m f /v H are the corresponding Yukawa couplings. The h-S mixing angle is In these expressions v S (T ) and v H (T ) are the S and Higgs vevs, respectively. S also couples to SM particles via annihilations. After the electroweak phase transition, the main annihilation process is SS → h * → ff , with cross section in the non-relativistic limit. In this expression, C f is a colour factor. We will choose λ p4 10 −3 , which makes S decays, inverse decays, and annihilations fast enough to keep S in thermal equilibrium with the SM at all T m χ , where DM freeze-in dominantly occurs. We have verified that this is always possible for m χ > m S 3 GeV. If S were not in thermal equilibrium during DM freeze-in, the model would not be invalidated, but the dark and visible sector temperatures would differ, complicating the analysis.
We will also assume that, while DM freezes in, ψ remains in thermal equilibrium with the SM sector, either through ψψ ↔ SS or through interactions with a second new scalar S as explained below eq. (2). The cross section for ψψ → SS in the non-relativistic limit is where v rel is the relative velocity of the annihilating ψ particles. Note the v 2 rel suppression of the annihilation cross section. Eventually, ψ will freeze out, and we ensure that this happens late enough for its relic abundance to make up a subdominant contribution to the DM density (less than 1%). Note that ψ is not absolutely stable, but can decay even at T = 0 through ψ → χ(S * /H * → ff ). The rate of this decay is (for m ψ − m χ > m S , (m ψ − m χ )/m χ 1, and m f = 0) given by Here, the sum runs over all kinematically accessible SM fermion species. Since we treat fermions as massless, eq. (11) will not be accurate near any of the fermion mass thresholds. Even though ψ freezes out with only a subdominant relic abundance, its decays may violate constraints if they happen around the time of recombination [32][33][34]. We therefore demand τ ψ ≡ 1/Γ ψ 10 11 sec or τ ψ 10 5 sec [33]. We have verified that for the parameter region we will discuss in the following, λ p4 and thus θ hS (T ) can indeed be adjusted such that τ ψ 10 11 sec. Even in this case, care must be taken that the residual abundance of ψ at the time of decay is tiny ( 10 −10 × (τ ψ /10 15 sec) times the DM abundance) to avoid anomalous reionization [35]. We have checked that this can be automatically achieved for DM masses 1 GeV. In this case, τ ψ is large because the only decay channels available to ψ are suppressed by the small Yukawa couplings of light quarks and leptons. At larger DM mass, the easiest way of ensuring compatibility of the model with reionization constraints is to introduce an auxiliary fermion χ , with m χ m χ , and with a Yukawa coupling to S and ψ chosen such that the decay ψ → χ ff is much slower than freeze-in of χ, but still occurs significantly before recombination (τ ψ 10 5 sec [33]).

II.2. The Effective Potential
To correctly describe the evolution of the scalar sector of the model from eqs. (1) and (2) in the hot early Universe, we must go beyond the tree level potential and consider the finite temperature effective potential V eff which includes radiative corrections and thermal effects. Since we assume the portal couplings λ p3 and λ p4 to be small, we can treat the evolution of the dark sector potential as independent from that of the visible sector potential (we will consider the case of large portal couplings in sections III and IV). We begin with the approximate tree level potential in the dark sector, The effective potential V eff is defined in the usual way [36]: one first rewrites the generating functional, , as a functional of v S (x) instead of the external source field J(x).

(Here, Z[J] is the partition function.) This is achieved by relating
/δJ. Note that, in the presence of an x-dependent external source, v S (x) becomes a function of x as well. The effective action Γ eff , which is in turn the spacetime integral of the effective potential V eff , is given by a Legendre transform: . We see that Γ eff has the property that δΓ eff [v S ]/δv S = 0, that is, the vacuum configuration v S , including all quantum corrections, is obtained from Γ eff using a variational principle. Of course, Γ eff itself needs to be computed perturbatively.
As we outline in more detail in appendix A, the leading corrections that distinguish V eff from V tree are the Coleman-Weinberg term V CW that corresponds to resummed 1-loop diagrams at T = 0 [37], the thermal one-loop contribution [38], V T , and the resummed series of higher order "daisy" diagrams, V daisy [39][40][41][42]. With our assumption that y ψχ , λ p3 , and λ p4 are tiny, loops involving H and χ are negligible. Loops involving ψ could be relevant at temperatures not too far below m ψ if y ψ 0.01. In the following we will assume y ψ 0.01 to simplify the analysis. As explained below eq. (2) this will require a different mechanism for keeping ψ in thermal equilibrium throughout DM freeze-in. We have verified that our toy model can be phenomenologically successful also for larger y ψ . For y ψ 0.01, the only relevant diagrams contributing to V eff are those involving the quartic coupling λ S4 . In other words, the sums in eqs. (A2), (A7) and (A8) run only over S. The coefficient n i , which can be interpreted as counting degrees of freedom (although see [42]), is n S = 1. As a function of the field value, the mass of S is given by V daisy also depends on the thermal, or Debye, mass, which is given by the 1-loop self energy at non-zero T . The Debye mass of S is given by [39] With the effective potential in hand, we can now consider the evolution of m S and v S ≡ S as a function of T . This allows us to describe the phase transition, which is analogous to the electroweak phase transition in the SM. We use the program CosmoTransitions [43][44][45][46] to track the minimum of the effective potential, to find S as a function of T , and to determine the mass of S as the second derivative of V eff (S). Although in our particular toy model it would be easy to do the computation without invoking CosmoTransitions, we still use it for consistency with sections III and IV.
The effective potential at several temperatures is shown in fig. 2 (left), while the behaviour of S and m S (T ) is shown in fig. 2 (right). In the left panel we see the well known behaviour of a second order phase transition or cross-over: at high temperatures, T T c , the effective potential has its minimum at S = 0. At the critical temperature, T c , the minimum begins to move away from S = 0 and a non-zero vev begins to develop. At the present time, near T = 0 GeV, the effective potential has its minimum at S 9 GeV. In fig. 2 (right), we similarly see that at high temperatures, S has no vev, and its effective mass is large thanks to thermal corrections ∝ T . As the temperature drops, the mass of S drops as well and approaches zero at the phase transition. This behaviour can be understood also from the green curve in fig. 2 (left): at the transition temperature, its second derivative at the minimum is zero. After the phase transition, the vev and the mass of S grow and quickly approach their present day values, m S (T = 0) = 5 GeV and S 9 GeV. Figure 2. The new scalar effective potential at a relatively high temperature, at T = T c and at T ∼ 0 GeV, (left). A constant term has been subtracted at each T so that V eff (0) = 0. The evolution of the new scalar mass and vev with temperature, (right), with the temperature of the phase transition indicated. Note that we neglect contributions from ψ loops to V eff , assuming y ψ 0.01.
As well as this temperature dependence of the S mass, the fermions may also have temperature dependent mass contributions. Since we take the tree level fermion masses to be much larger than the S mass, self-energy diagrams evaluated at T = 0 will only give a small contribution near the phase transition (fermions have no zero Matsubara mode, so the self-energy contributions for fermions are smaller than those for bosons). However, if y ψ is large, then there can be significant corrections to the ψ mass from the Lagrangian term y ψψ Sψ when S obtains its vev. The T dependent ψ mass is then The crucial point for us is that, thanks to the behaviour of m S (T ) and m ψ (T ), the freeze-in channel ψ → χS is kinematically closed long before and long after the dark sector phase transition, while around the transition temperature, it is open and DM freeze-in can proceed efficiently.

II.3. Results and Discussion
In fig. 3 (top-left) we show the instantaneous DM yield from each freeze-in channel in fig. 1, including (solid lines) and ignoring (dashed lines) finite temperature contributions to the scalar masses and vevs and to the fermion masses. If we ignore the finite temperature corrections, only two channels (Sψ → χS and SS → χψ) contribute to the χ abundance. These contributions are largest at high temperatures (or small x, where x = m χ /T ), where the abundance of ψ is not Boltzmann suppressed and where S have enough energy to produce the heavier states ψ and χ. The instantaneous freeze-in yield reduces smoothly as the temperature reduces, except for small steps where SM particles freeze-out and the effective number of relativistic degrees of freedom in the Universe, g eff , changes.
If finite temperature effects are included, then two new channels contribute at different times. At very high temperatures, S has a sufficiently large mass that it can decay to χψ. As χ and ψ are The instantaneous change in the yield (top-left) and the resulting freeze-in relic abundance extrapolated to the present day (top-right) for our benchmark point, including (solid lines) and ignoring (dashed lines) finite temperature effects. Also shown is the relative final abundance (normalised to the abundance when finite temperature effects are included) through each channel as a function of m χ for m S (T = 0) = 5 GeV (bottom-left) and m S (T = 0) = 10 GeV (bottom-right). much heavier than S at T = 0 GeV, this channel is only open at very high temperatures (x 0.1), and it no longer contributes at lower temperatures. As the Universe approaches the phase transition in the dark sector, the mass of S reduces until it becomes smaller than m ψ (T ) − m χ at x 1.7 (the mass of ψ is constant before S obtains its vev). At this point, the decay ψ → χS becomes kinematically possible and χ is produced. This happens at a rate much larger than that via the Sψ → χS and SS → χψ channels because the latter channels are suppressed by the off-shellness of the intermediate ψ propagator. The ψ → χS channel reaches its maximum rate around the dark phase-transition, where m S (T ) goes to zero. As the temperature further reduces, the mass of S increases and the mass of ψ reduces until the channel closes at x 4.5. Comparing the rates with and without including finite temperature effects, we see that these effects are relevant in all channels. The rate of Sψ → χS shows a peak at the dark phase-transition because the intermediate s-channel ψ propagator in the third diagram of fig. 1 can go nearly on-shell when m S (T ) is small around the transition temperature. This is essentially a manifestation of the infrared divergence of the corresponding amplitude.
The resulting relic abundance of χ extrapolated to zero redshift is shown in fig. 3 (top-right). The extrapolated abundance at a given x ≡ m χ /T is obtained by rescaling the number density at this time by the subsequent expansion of the Universe and normalising to the critical density today. Here we clearly see that the dominant contribution to the relic abundance comes from the ψ → χS channel. We emphasise that if finite temperature effects were not included, the calculated relic abundance would be incorrect by a factor O(10 4 ). For the benchmark parameters chosen, the resulting χ abundance matches the observed relic abundance for y ψχ = 2.57 × 10 −12 .
In fig. 3 (bottom-left) we show the abundance produced through each channel as a function of m χ , normalised to the abundance when finite temperature effects are included. We keep the mass difference between ψ at T = 0 and χ fixed at 4 GeV. The fraction of χ produced through channels other than ψ → χS is always small, but it is smallest for m χ 50 GeV. At lower values of m χ , S no longer requires significant thermal energy to produce χ via Sψ → χS and SS → χψ, so the ψ → χS channel produces a smaller fraction of χ. As the value of m χ increases, processes which occur at lower temperatures receive more Boltzmann suppression than those occurring at higher temperatures. This means that the amount of χ produced through Sψ → χS and SS → χψ (which is important at high temperatures) is mildly reduced whereas the amount produced via ψ → χS (which is important around the phase transition) has greater Boltzmann suppression, reducing its relative importance.
Finally, fig. 3 (bottom-right) shows the freeze-in abundance of χ through the different channels for m S = 10 GeV, where now the phase transition occurs at T = 36 GeV. We see that the picture is qualitatively similar to fig. 3 (bottom-left). For m s = 10 GeV, there is a milder reduction in the relative importance of the ψ → χS channel at large m χ , due to a milder Boltzmann suppression of the ψ abundance at the phase transition. In both fig. 3 (bottom-left) and (bottom-right), the yields due to Sψ → χS are similar including or ignoring finite temperature corrections. This channel predominantly produces χ at high temperatures, where the particle momentum dominates over the particle masses. For the small Yukawa coupling y ψ = 0.01 and large quartic coupling λ S4 = 1, the SS → χψ production rate is much larger when S has a vev, which explains the difference seen between the curves that include or ignore the finite temperature corrections in this channel.
For definiteness in our numerical calculations, we fix a reheating temperature T R = 1 TeV which is sufficiently large so that any freeze-in at higher temperatures will produce negligible abundance of χ.
We finish this section by noting a particularly simple model which shares many features with the toy model discussed above. If the SM is extended by two dark sector fermions, a SM gauge singlet, χ, and an su(2) L doublet with hypercharge, ψ, then the SM Higgs can play the role of S above. The Lagrangian termψHχ + h.c. can lead to processes which produce χ via freeze-in. If m ψ − m χ < m h (T = 0) and m χ < m ψ , then the channel ψ → χH will be open only when the mass of the SM Higgs is reduced during the second order electroweak phase transition. For ψ lighter than 1.1 TeV, it freezes-out as a subdominant component of the dark matter abundance [47]. The remaining relic abundance can then be provided by the freeze-in of χ. The calculation of the abundance is somewhat complicated and we defer this to later work.

III. VEV-INDUCED PRODUCTION WITH A VEV FLIP-FLOP
In the previous section we have highlighted the potential importance of including finite temperature corrections to particle masses in calculations of DM freeze-in. In this section, we consider a model with a fermionic DM candidate χ and with a scalar sector identical to the one in eq. (1), but  focussing on a different region of parameter space. Namely, we now take the Higgs portal couplings so large that a two-step phase transition (or "vev flip-flop" [28]) is realised [27,[48][49][50][51][52][53][54] (see also [55]). In other words, the Universe goes through a phase where the new scalar S obtains a non-zero vev, but this vev jumps back to zero in a first order electroweak phase transition. The value of S will control the DM freeze-in rate, so DM can only be efficiently produced during a relatively short time interval. The final DM abundance is determined not only by the relevant coupling constants, but also by the length of this time interval. We dub this mechanism "vev-induced production."

III.1. Toy Model
The field content of the dark sector in this toy model is shown in table II. As in section II, our dark matter candidate is a Dirac fermion, χ, which is a SM gauge singlet. We assume that it is stabilised by a Z 2 symmetry. The DM mass m χ and the scalar mass parameter |µ S | are taken to be 100 GeV. The relevant terms in the Lagrangian are with V (H, S) again given by eq. (2). DM will freeze-in via the Yukawa coupling y χ ; consequently, this coupling needs to be tiny. As in section III.1, this could be motivated in extra-dimensional scenarios by localising χ far away from the other fields along a fifth dimension. λ S3 and λ p3 are assumed to be small as well to simplify the analysis. Note that an extra global Z 2 symmetry is restored if y χ , λ S3 , and λ p3 are set to zero. Therefore, small values for these couplings are natural in the 't Hooft sense [56]. Setting λ p3 1 means we can ignore mixing between the SM Higgs boson and S at T = 0. In this limit, the mass of S is given by at tree level and T = 0. We see that for µ S v H (T = 0) and λ p4 1, a situation can be realised where λ p4 v 2 H (T = 0)/2 > µ 2 S . In this case, v S (T ) ≡ S = 0 as long as v H (T ) = 0 (barring thermal corrections for the moment), but when v H (T ) becomes significantly different from zero, the term quadratic in S in the scalar potential experiences a sign flip, making v S (T ) = 0 energetically favourable in the broken phase of electroweak symmetry. v S (T ) = 0 is also realised at very early time thanks to thermal corrections to V eff . These corrections are large especially when λ S4 1. This behaviour is the gist of the vev flip-flop, and it defines the parameter region we will be interested in in the following: µ S v H (T = 0), λ p4 λ S4 1.

III.2. The Effective Potential & The Vev Flip-Flop
To quantitatively compute the effective potential for the model defined in eq. (16), the same methods as in section II.2 can be applied, but since the Higgs portal coupling λ p4 is no longer neg-ligible, we need to consider the joint evolution of the visible and dark scalar sector. In other words, we need to treat V eff as a function of both S and H. As explained in section II.2 and appendix A, the one-loop contribution to the effective potential, V CW (h, S) + V T (h, S), depends on the field dependent masses of all particles with couplings to the scalars. In the model from eq. (16), the field dependent masses of W i , B and t are the same as in the SM. In particular the gauge boson mass eigenvalues are m 2 W ± = 1 4 g 2 h 2 , m 2 Z = 1 4 (g 2 + g 2 )h 2 , m 2 γ = 0, and m t = y t h/ √ 2. Here, g and g are the su(2) L and u(1) Y gauge couplings, respectively, and y t is the top quark Yukawa coupling. The mass matrix of the neutral CP-even Higgs bosons is while the mass of the neutral CP-odd and the charged component of H are given by We can see that in regions where both h and S are non-zero, there will be mixing between the associated particles. Note that with our simplifying assumption λ p3 0, we can neglect h-S mixing in the S = 0, h = 0 phase at T = 0. The sums in V CW (h, S) and V T (h, S) (see eqs. (A2) and (A7)) now run over i ∈ {h, S, G 0 , G + , W i , B, t}. For h and S, this is understood to mean summing over the neutral CP-even mass eigenstates, determined by diagonalising the mass matrix in eq. (18). The coefficients n i for the SM fields are n h = n s = n G 0 = 1, n G + = 2, n W i = n B = 3, and n t = −12 [42]. The Debye masses of h and S, relevant in the computation of V daisy (see eq. (A8)) are The Debye masses of W i , B and t are the same as in the SM (see appendix A) as these particles do not couple to S. The sum in eq. (A8) runs over i ∈ {h, G 0 , G + , S, W i , B}. The contribution of χ to the effective potential is negligible due to the smallness of y χ and is therefore dropped in our calculations. The behaviour of the effective potential is illustrated in fig. 4 for a parameter point featuring a two-step phase transition or vev flip-flop. At early times (top left panel of fig. 4), V eff is dominated by the finite temperature and therefore approximately parabolic. Consequently, the SM Higgs and S have zero vacuum expectation values. As the universe expands and cools down (top right panel of fig. 4), the finite temperature corrections become similar in magnitude to the tree-level terms and the effective potential develops minima at S = 0. There is typically a second order phase transition, so the Universe immediately enters a phase where S has a non-zero vev. After further cooling, new minima at h = 0 develop. These become the global minima at some critical temperature T c . However, there will now be a barrier between the S = 0 and the h = 0 minima, so the Universe cannot immediately transition into the global minimum, but undergoes a short period of supercooling. The subsequent phase transition is first order and proceeds via bubble nucleation, when at the nucleation temperature T n it becomes energetically favourable for bubbles of the new phase to expand and fill the entire universe. Typically, one finds T n T c , but in narrow regions of parameter space, T n may also be significantly below T c . As in section II,  we have used CosmoTransitions [43][44][45][46] to determine the nucleation temperature. Numerically, CosmoTransitions computes T n by determining the temperature at which S E (T )/T drops below a critical value of 140. Here S E (T ) is the minimum Euclidean action corresponding to a transition between the two potential minima [43]. The effective potential at T = T n is shown in the bottom left panel of fig. 4. The h = 0 minima then deepen as T goes to zero, and the universe remains in a phase where S = 0 and h = 0.
In fig. 5 we show the masses and vevs of the new scalar S and the SM Higgs doublet as a function of temperature. At high temperatures S is zero, but at lower temperatures it obtains a non-zero value. This situation persists until the first order electroweak phase transition at T n = 136 GeV. At this point, S goes to zero while the SM Higgs vev becomes non-zero. h then gradually increases until it attains its T = 0 value of 246 GeV. We can see that, as in fig. 4, the scalars receive large finite temperature corrections to their masses at high temperatures. For the parameters chosen here, both m S and m h become smaller than their T = 0 values between the phase transitions, similar to what we found for m S in section III.2.

III.3. Dark Matter Freeze-In and Relic Abundance
In the model defined in eq. (16), the coupling λ p4 between the new scalar field, S, and the SM Higgs doublet needs to be of order one for the vev flip-flop to occur (see discussion below eq. (17)). Sizeable λ p4 in turn means that at T m S , m H , the scalar S is in thermal equilibrium with the SM sector. χ, on the other hand, never comes into thermal equilibrium because of our assumption that y χ is tiny. Instead, a small abundance of χ is produced via freeze-in, facilitated by the processes S → χχ, H † H → χχ, and SS → χχ. The corresponding Feynman diagrams are shown in fig. 6, and the decay rates and annihilation cross sections are The last expression should be understood as the cross section for one of the two components of the doublet H to annihilate with its antiparticle into χχ. The first process, S → χχ, is kinematically forbidden when m S (T ) < 2m χ , but since thermal corrections may drive m S (T ) to large values at T m S (T = 0), it may be allowed at early times. Whether or not this production channel is important will thus depend on the reheating temperature. We will be particularly interested in T R not too far above the electroweak scale, as in this case the dynamics of the vev flip-flop are most important for DM physics. The other two freeze-in processes can be mediated either by the cubic scalar couplings λ S3 , λ p3 , or by the quartic couplings λ S4 , λ p4 if v S = 0. We will focus on the parameter region where the cubic couplings are small because this is the region where Figure 6. The main production modes for the DM particle χ in the vev-induced freeze-in scenario.
freeze-in depends most strongly on v S and thus on the dynamics of the vev flip-flop. Moreover, as explained in section III.1, λ S3 , λ p3 1 is technically natural in the 't Hooft sense. We note, however, that λ p3 should not be exactly zero. Otherwise, the small relic abundance of S could not decay away and would violate direct detection limits. The small vev that S has even at T = 0 when λ p3 = 0 would not affect our results. It is also important to note that the restrictions we impose here on the parameters of the scalar potential and on T R are not necessary to make the model phenomenologically viable. There are other large regions of parameter space where the DM relic abundance can be successfully generated, albeit without strong involvement of the vev flip-flop.
The DM production rate and the resulting abundance (extrapolated to zero redshift) are shown in fig. 7 for two illustrative parameter points of the vev-induced freeze-in scenario. The first parameter point shown in fig. 7 (top panels) is characterised by a low reheating temperature T R = 500 GeV. In this case, DM production is entirely dominated by 2 → 2 processes proportional to v S , so the dynamics of the vev flip-flop are crucial in this case. We see that the DM production rate dY /dx rises rapidly after S develops a vev at T 250 GeV. Between the two phase transitions, dY /dx follows the evolution of v S , and the contribution from the v S -dependent channels is 2-3 orders of magnitude larger than the contribution from vev-independent SS → χχ annihilation via λ S3 . Among the vev-dependent channels, SS → χχ dominates over H † H → χχ mainly because λ S4 > λ p4 . The small drop in dY /dx immediately before the electroweak phase transition is due to the onset of Boltzmann suppression of H and S. After the electroweak phase transition at T 136 GeV, the v S -dependent production processes cease. Before and after the two phase transitions, only the v S -independent channel SS → χχ via λ S3 is active, but its overall contribution is tiny. At x 0.5, this channel is suppressed as the s-channel mediator, S, is very heavy at these high temperatures. Beyond the electroweak phase transition this channel gradually reduces, due to the Boltzmann suppression of S.
At the second parameter point shown in fig. 7 (bottom panels) the behaviour of the vevdependent DM production channels and of vev-independent production via SS → χχ (mediated by λ 3 ) is similar to the top panels. However, as the second parameter point features a larger T R = 10 TeV, the decay channel S → χχ, which does not depend on v S and on the vev flip-flop, is kinematically allowed at x 0.2, when thermal corrections lift m S above 2m χ . In this case, this production channel dominates the final abundance.
We conclude that, for low T R , it is the dynamics of the vev flip-flop that determines the DM abundance today. For high T R , it is the thermal corrections to m S (T ), which in turn depend on the couplings in the scalar sector, especially λ S4 . In either case, the inclusion of thermal effects in the computation of the DM relic density is essential. To emphasise this point, we show in fig. 7 also the production rate and abundance that would be obtained if thermal corrections to V eff (and thus the two-step phase transition) were neglected in the calculation (dashed blue lines). In this case, only the processes SS → χχ and H † H → χχ with cross sections proportional to the small couplings λ 2 S3 and λ 2 p3 , respectively, would contribute. For our choice λ p3 = 0, only SS → χχ is open. The production rate in this channel is non-zero at T R (or even during preheating [57], which we neglect here assuming it is very rapid) and first rises slightly as there is more time to freeze-in  Figure 7. Evolution of the DM production rate (left) and the DM abundance extrapolated to z = 0 (right) as a function of x ≡ m χ /T for two different parameter points in the vev-induced freeze-in scenario defined by eq. (16). Solid curves correspond to the different production mechanisms shown in fig. 6. The dashed blue line indicates the result one would obtain if thermal corrections and the vev flip-flop were neglected. We see that for a reheating temperature T R just above the electroweak scale (top panels), the DM abundance is entirely dominated by the processes H † H → χχ and SS → χχ, whose rates are greatly enhanced when S = 0. For larger T R (bottom panels), the process S → χχ, which is independent of S , becomes allowed thanks to thermal corrections, even though at T = 0, m S < 2m χ .
at greater x. At x 0.5, the rate begins to drop as Boltzmann suppression becomes significant. Eventually, λ S3 -mediated production via SS → χχ freezes out.
We further explore the crucial importance of thermal corrections in fig. 8, which shows a cut through the model's parameter space in the plane spanned by the zero temperature mass of S, m S (T = 0), and the quartic Higgs portal coupling λ p4 . The pixelated region shows where the twostep phase transition occurs. The colour coding quantifies the ratio of the DM abundance obtained including thermal corrections to the abundance if these corrections were neglected for a low T R . We see that thermal effects dominate the abundance by up to four orders of magnitude. For fixed m χ , they are largest at small m S (T = 0), where thermal corrections to m S are most important, and at large λ p4 , where the S = 0 phase lasts longer. For a high T R ∼ 10 TeV, the freeze-in abundance is dominated by the S → χχ channel and thermal effects dominate the abundance by around four orders of magnitude over the whole parameter space shown. In the white area at the top of fig. 8, the global minimum of V eff at T = 0 would be the one with S = 0, i.e., electroweak symmetry would never be broken. In the white region at the bottom of the plot, S never acquires a non-zero vev.

IV. VEV-INDUCED MIXING WITH A VEV FLIP-FLOP
Let us now move to a third scenario illustrating the importance of thermal effects on the DM abundance in the Universe. The scenario discussed in the following, which we dub mixing induced freeze-in, is based on the same particle content as the kinematically induced freeze-in model from section II with an extra discrete symmetry. The model's Lagrangian is thus given by eqs. (1) and (2), with the two forbidden Yukawa terms removed. The most important term for DM freezein is again the Yukawa coupling y ψχψ Sχ. However, we now assume the DM candidate χ and the new scalar S to have masses around the electroweak scale, with m S > m χ . The auxiliary new fermion ψ is assumed to be much heavier (see table III). The idea is that the reheating temperature is low, T R < m ψ , so that ψ never comes into thermal equilibrium and DM production via ψ → Sχ (the channel we had focused on in section II) does not occur. Instead, the main DM production channels will be S → χχ, facilitated by χ-ψ mixing through the Sχψ coupling, and SS → χχ, mediated by a t-channel ψ (see fig. 9). The former process is of particular interest to us because χ-ψ mixing depends on the vev of S. For the parameters of the scalar potential, we consider values similar to the ones we chose (and motivated) in section III: negligible λ S3 , λ p3 , but sizeable Higgs portal and dark sector quartic couplings, λ p4 , λ S4 O(1), to induce a vev flip-flop. Thus, DM production via S → χχ will be open for a limited amount of time while v S = 0. In the following, we  will study the interplay of the two production processes, focusing in particular on the importance of v S and the vev flip-flop. The decay width for S → χχ is and the cross section for SS → χχ reads where we have taken the limit m ψ √ s, m S (T ), m χ , which is justified in view of our assumption m ψ T R . The mixing angle between χ and ψ is The dynamics of freeze-in via mixing are illustrated in fig. 10. In analogy to fig. 7 in section III, this figure shows the DM production rate and the extrapolated DM abundance as a function of x ≡ m χ /T . We again consider one parameter point with a low reheating temperature, T R = 150 GeV (top panels), and one parameter point with a higher reheating temperature, T R = 500 GeV (bottom panels). For low T R , the DM abundance today is dominated by the decay S → χχ, which is only possible for v S = 0. For the parameters shown in the top panels of fig. 10, this phase is already realised at T R . The rate for S → χχ increases along with v S (T ) and then drops sharply to zero at x 0.7 as the electroweak phase transition switches v S off in favour of non-zero v H . Similar behaviour is also seen for higher reheating temperature (bottom panels of fig. 10), but in this case, the overall importance of S → χχ does not dominate that of SS → χχ. As the latter process is independent of v S , it leads to DM production immediately after reheating (or already during preheating, which we neglect here), while S → χχ is only activated when v S becomes non-zero at  Figure 10. Evolution of the DM production rate (left) and the DM abundance extrapolated to z = 0 (right) as a function of x ≡ m χ /T for two different parameter points in the mixing induced scenario from section IV. Solid curves correspond to the different production mechanisms shown in fig. 9. The dashed curves indicates the result one would obtain if thermal corrections and the vev flip-flop were neglected. We see that for a reheating temperature T R just above the electroweak scale (top panels), the DM abundance is dominated by the mixing-induced decay S → χχ, which is facilitated by the vev flip-flop. For larger T R (bottom panels), the process SS → χχ, which is independent of S , becomes relevant.
x 0.2. Note from eq. (26) that for T R < m ψ (the case realised in fig. 10), the DM production rate via SS → χχ scales ∝ T 3 /m 2 ψ . In other words, freeze-in is ultraviolet-dominated [58]. The dependence of SS → χχ on thermal corrections (namely the temperature dependence of m S (T )) is relatively weak. It is reflected for instance in a jump in the rate at the electroweak phase transition. For comparison, we show also the production rate and DM yield in the absence of thermal corrections (dashed lines). We see that a calculation neglecting these corrections would fairly accurately predict freeze-in via SS → χχ, but would completely miss the mixing-induced channel S → χχ.
We further explore the dependence of thermal effects on the parameters of the mixing induced scenario in fig. 11. The two panels in this figure show different slices through the parameter space:  Figure 11. The ratio of the DM abundance including over ignoring finite temperature effects in the mixing induced scenario. We show a cut through the parameter space in the plane spanned by the zero temperature mass of S, m S (T = 0), and the quartic Higgs portal coupling λ p4 (left) and a slice of parameter space in the m S (T = 0)-λ S4 plane (right). We see that the impact of thermal effects is largest at large portal coupling λ p4 and small self-coupling λ S4 , where the S = 0 phase lasts longer. The black outline indicates the point considered in fig. 10. one in the m S (T = 0)-λ p4 plane (cf. fig. 8) and one in the m S (T = 0)-λ S4 plane. We observe that thermal effects -in this case the v S -dependent process S → χχ -are most relevant both when the Higgs portal coupling λ p4 is large and when the quartic self-coupling λ S4 or the T = 0 S mass is small. For larger λ p4 or smaller m S (T = 0), the v S = 0 vacuum becomes deeper and the phase during which S has a vev and S → χχ is open becomes longer. For small λ S4 , the v S = 0 phase begins earlier, again implying that there is more time for DM production via mixing. In any case, we see that thermal effects are at most O(1) in most of the parameter space, but in some regions can modify the DM abundance today by an order of magnitude.

V. SUMMARY AND CONCLUSIONS
In this work we have considered the impact of finite temperature corrections on the freeze-in of dark matter. We have highlighted several effects which can have a dramatic impact on dark matter production. We have illustrated the impact of these effects in three toy models, which demonstrate 'kinematically induced', 'vev induced' and 'mixing induced' freeze-in, respectively.
In 'kinematically induced freeze-in', the dominant production channel of dark matter may be closed at zero temperature, but may be open in the early universe as temperature-dependent particle masses vary. Although calculationally complex, a simple realisation of this is the SM Higgs coupling to two dark sector fermions. We have highlighted the analogous effect in a realistic toy model consisting of a new scalar which is weakly coupled to the SM, and two dark sector fermions. We show that a calculation ignoring the temperature-dependent scalar mass produces an estimate of the dark matter abundance which is incorrect by a factor of O(10).
If instead, the new scalar couples significantly to the SM Higgs, the Higgs portal coupling can induce a two-step phase transition (or "vev flip-flop"). In this case the new scalar may obtain a vev for some time, which then disappears when the SM Higgs obtains its vev. This can lead to 'vev induced' production, where dominant channels of dark matter production only open when the new scalar has a vev. We illustrate the effect in a phenomenologically viable toy model. For reheating temperatures around the electroweak scale, dark matter production occurs mainly via 'vev induced' production, whereas for higher reheating temperatures (T R a few TeV) 'kinematically induced freeze-in' dominates the production. In both cases, these finite temperature effects can easily change the relic abundance by several orders of magnitude.
Finally we consider a scenario where the temporary vev of a new scalar leads to mixing between fermionic dark matter and another fermion, which we call 'mixing induced freeze-in'. This mixing may then offer a dark matter production mode which is not available at zero-temperature. We again consider a realistic toy model and show that dark matter can be dominantly produced through this temperature-dependent channel. The error introduced by ignoring this effect can be as large as a factor of ∼ 10, but depends crucially on a reheating temperature around the electroweak scale and a long period of vev induced mixing. When the reheating temperature is much higher or the vev induced mixing is weak, the standard calculation provides a reliable estimate.

ACKNOWLEDGMENTS
It is a pleasure to thank Christophe Grojean, Matthias König, and Andrea Thamm for useful discussions. MJB would like to thank CERN for warm hospitality during part of this work. This  For a 1 → 2 decay of the form B 1 → χB 2 (where χ is the DM candidate and B 1 , B 2 are its interaction partners), the general Boltzmann equation iṡ Here, dΠ i ≡ d 3 p i /(2π) 3 denotes the three-dimensional phase space integral for particle i = χ, B 1 , B 2 ; p i are the particles' momenta and f i their momentum distribution functions in the primordial plasma. M B 1 →B 2 χ and M B 2 χ→B 1 are the transition amplitudes of the freeze-in reaction and its inverse. In freeze-in scenarios, the abundance of χ remains well below its equilibrium abundance, so f χ 0. Consequently, the second term in square brackets in eq. (B1) can be dropped. We also neglect Pauli blocking and stimulated emission, so that (1 ± f B i ) = 1. We can then writė (B2) We assume that B 1 is in thermal equilibrium and that E B 1 T so that f B 1 e −E B 1 /T . We expect this assumption to introduce an error of around 15% in our final abundances [20]. We obtaiṅ where K 1 is a modified Bessel function of the second kind and g B 1 is the number of degrees of freedom of B 1 . It is convenient to normalise the number densities n i to the entropy density s to factorise the trivial dilution of n i due to the expansion of the Universe. This leads to the yield Y i ≡ n i /s. Introducing the dimensionless evolution variable x ≡ m χ /T , the Boltzmann equation takes its final form where H(x) is the Hubble rate.
2. Freeze-In via Annihilation: B 1 B 2 → χB 3 Following similar steps as in appendix 1, and using the definition of the Møller velocity v Møl [60], the Boltzmann equation for a 2 → 2 process of the form B 1 B 2 → χB 3 readṡ n χ + 3Hn χ = g B 1 g B 2 dp 3 B 1 dp 3 We simplify this expression following ref. [61]. To this end, we define where E i and p i are the particles' energies and three-momenta, respectively, and θ is the angle between p B 1 and p B 2 . It is straightforward to show that Moreover, we have [61] where is the modulus of the momentum of B 1 and B 2 in the centre of mass frame. The integrand on the right hand side of the Boltzmann equation (B5) is independent of E − and depends on E + only through the exponential e −E + /T . Evaluating the integrals over E + and E − yields [61] n χ + 3Hn Note that, in counting the degrees of freedom g B 1 , g B 2 , care must be taken that each degree of freedom counted towards g B 1 should be able to annihilate with each degree of freedom counted towards g B 2 . This is usually true for spin and colour degrees of freedom (which the cross sections are typically averaged over). Particles and antiparticles, on the other hand, should be treated as different initial states, not as different degrees of freedom of the same initial state. The same is true for different components of an su(2) L multiplet. Expressed in terms of particle yields rather than number densities, eq. (B12) transforms into For the special case B 2 =B 1 , this simplifies to