Relativistic Freeze-in

We study production of scalar dark matter via the freeze--in mechanism in the relativistic regime, focussing on the simplest Higgs portal model. We derive the corresponding relativistic reaction rates based on the Bose--Einstein statistics taking into account the thermal mass effects as well as the change in the Higgs degrees of freedom at the electroweak phase transition. The consequent constraints on the Higgs portal coupling are obtained.


Introduction
The freeze-in dark matter (DM) production mechanism [1], [2] is an attractive possibility for generating the observed DM abundance. It makes use of a feeble coupling between the Standard Model (SM) and dark matter, which leads to slow production of DM by the SM thermal bath. In this case, dark matter never thermalizes and its abundance accumulates over time. DM production can take place both in relativistic and non-relativistic regimes, depending on the dominant mode for a given dark matter mass.
The existing analyses typically employ the non-relativistic approximation replacing the Bose-Einstein distribution function with the Maxwell-Boltzmann one. While this can often be justified, it can also significantly underestimate the reaction rates at high temperatures. Depending on the model parameters, the rates can differ by orders of magnitude [3]. In this work, we derive the relativistic reaction rates based on the Bose-Einstein statistics in the framework of Higgs portal dark matter [4]- [6]. These calculations parallel our recent work on self-interacting scalar dark matter [3]. We also take into account the effects of thermal mass corrections and the electroweak phase transition, both of which make a significant impact on the result.
Recent work on the Higgs portal dark matter freeze-in which employs various degrees of simplification can be found in Refs. [7]- [18]. 1

The set-up
We consider the simplest real scalar dark matter model with the Higgs portal interaction where s is a real scalar with the potential Here, the symmetry s → −s stabilizes the scalar which plays the role of dark matter (DM). The freeze-in DM production takes place in the Standard Model thermal bath through the Higgs coupling when λ hs 1. For the coupling values in the ballpark of λ hs ∼ 10 −11 , one obtains the correct DM relic abundance, while DM never thermalizes (assuming also that λ s is small enough).
There are two reactions in which DM is produced: Higgs annihilation hh → ss and Higgs decay h → ss, if allowed kinematically. The latter is only possible at temperatures below the electroweak scale, so we distinguish (a) symmetric phase, T T EW and H = 0 (b) broken phase, T < T EW and H = 0 Here T EW ∼ v is the EW phase transition temperature. In the symmetric phase, the gauge bosons are massless and there are 4 massive Higgs degrees of freedom h i such that In this case, dark matter is produced through annihilation, In the broken phase, one can use the unitary gauge H = (0, (h + v)/ √ 2) T , which is singular at H = 0, such that This allows for both annihilation and decay, if allowed kinematically, Note that only one Higgs degree of freedom contributes in this case. 2 The total DM yield is given by the sum of yields in the two regimes. Here we omit complications having to do with the transition period at T ∼ T EW and approximate the reaction rates using the step function θ(T − T EW ). While in the broken phase relativistic effects are unimportant, at high temperature they are significant. In this case, it is important to take into account the thermal mass corrections which represent the leading thermal effects. For the Higgs field, the corrected mass-squared is with m h0 being the zero temperature Higgs mass; g 1,2 , y t are the gauge and top quark Yukawa couplings, and λ h is the Higgs self-coupling. In the symmetric phase, On the other hand, we assume that DM is not thermalized and m s0 does not receive significant thermal corrections (which would not be suppressed by λ hs ).
We also note that the DM mass changes during the phase transition and receives an extra contribution, For λ hs in the range of interest, this effect is negligible unless s has an MeV (or below) mass. However, for such light dark matter only the decay production mode in the broken phase is important, so only the total m 2 s matters. For heavier DM, we make no distinction between m s0 and m s .

Relativistic reaction rates
In this section, we compute the 2 → 2 and 1 → 2 relativistic reaction rates necessary for evaluation of the DM relic abundance. We follow closely our earlier work [3] where analogous computations for self-interacting scalar DM have been performed.
The a → b reaction rates per unit volume are (9) Here M a→b is the QFT transition amplitude, in which we also absorb the initial and final state symmetry factors; f (p) is the Bose-Einstein momentum distribution function. In thermal equilibrium, f (p) can be written in a covariant form as where u µ is the 4-velocity of our reference frame relative to the gas rest frame in which u = (1, 0, 0, 0) T . For freeze-in production of DM, the final state enhancement factors 1 + f (p j ) can be set to 1 since DM is not thermalized and its abundance is much lower than that in equilibrium.

hh → ss reaction
Consider the reaction hh → ss due to the quartic interaction term. Following Gelmini and Gondolo [19], we write it as where the Møller velocity for the incoming particles is given by Here F (p 1 , p 2 ) = (p 1 · p 2 ) 2 − m 4 h and the cross section is defined by where, in our convention, |M 2→2 | 2 includes the symmetry factors for identical particles in the final and initial states. The cross section is conveniently calculated in the center-of-mass (CM) frame. Following [3], we introduce and use the parametrization in terms of the half-the-CM energy E and rapidity η, where θ and φ are the polar and azimuthal angles, respectively. This corresponds to p = Λ(E, 0, 0, 0) T with Λ being a Lorentz transformation. As shown in [3], for any g(p 1 , p 2 ), where in the integrand one must set k 0 = 0, |k| = E 2 − m 2 h in k-dependent quantities, and Ω p,k are the solid angles in p-and k-spaces. In our case, the angular integrations can be performed explicitly. In the CM frame, the cross section is a function of E only and the angular dependence appears solely in the Bose-Einstein factors, where k 3 = |k| cos θ k . Computing the integrals, we get The QFT cross section for the process hh → ss is where, in our convention, we include the symmetry factor 1/(2!2!) associated with identical particles in the final and initial states. Inserting this in the reaction rate, we get where E is half the center-of-mass energy and we have factored out the symmetry factor 1/(2!2!). This rate is to be multiplied by 4 in the symmetric phase.

h → ss reaction
Let us now consider the decay mode h → ss. Again, it is convenient to go over to the rest frame of the decaying particle with momentum p, so that p = Λ(E, 0, 0, 0) T . Parametrizing p as we did in the previous subsection, we have d 3 p 2E = δ(E 2 −m 2 h ) sinh 2 η E 3 dE dη dΩ. Using the decay width Γ definition we find where we have used which includes a symmetry factor 1/2! due to identical particles in the final state. In the non-relativistic Maxwell-Boltzmann limit, this agrees with the result in [2].

Implications
The computed reaction rates are to be used in the Boltzmann equation in order to determine the DM density evolution. Compared to their Maxwell-Boltzmann counterparts, these rates are enhanced due to the Bose-Einstein distribution function peaking at low momenta. In general, the Bose-Einstein rates can exceed the Maxwell-Boltzmann ones by orders of magnitude [3], however the effect is sensitive to the thermal mass: for larger masses it is less pronounced. In the case at hand, the Higgs field receives a large thermal correction due the gauge and top quark couplings. The resulting enhancement is therefore modest as shown in Fig. 1. It reaches 50% for the annihilation mode and 20% for the decay.
It is important to note that the inclusion of the thermal mass regulates the high-T behaviour of the rates which is equivalent to curing the infrared divergence as m h,s → 0. Let us set m s = 0 and consider the limit m h → 0. In this case, we find that the 2 → 2 rate diverges as ln m h which is unphysical. Including the thermal mass, we get which also represents the high-T behaviour. Here c is a constant depending on the couplings. Therefore, the rate exhibits the expected scaling behaviour T 4 .

The Boltzmann equation
The evolution of the dark matter number density n(t) is governed by the Boltzmann equation. In our case, it takes the forṁ where the dot denotes a time derivative, H is the Hubble rate; the factor of 2 is due to production of 2 DM particle in each reaction, and the θ-functions take into account the change in Higgs degrees of freedom and the vertices at the electroweak phase transition. Here the inverse processes have been neglected given that the DM density is far below its value in equilibrium. For our purposes, we may omit the contributions like W W → h → ss (see e.g. [15]), which are only possible in the broken phase. We find that if m h > 2m s , the DM yield is dominated by the decay in the broken phase and the annihilation processes can be neglected altogether. For m h < 2m s , the production takes place predominantly in the relativistic regime and the above mode is irrelevant. In conclusion, although the gauge boson contribution is similar in magnitude to hh → ss in the broken phase, its effect is inconsequential.

Analytic solution in the relativistic regime
Qualitative behaviour of n(t) in the relativistic regime at T > T EW can be understood analytically. It is convenient to trade the time variable for T . Due to SM entropy conservation, T 3 a 3 = const, where a is the scale factor. This impliesṪ = −HT and the Boltzmann equation can be written as Here the Hubble rate is given by with g * being the number of SM degrees of freedom and M Pl being the reduced Planck mass. The last term in Eq. 25 scales as T 2 at high temperature, therefore it is convenient to define a constant Imposing the boundary condition that the DM abundance vanish at some initial temperature T 0 , the solution to the Boltzmann equation reads We thus see that the total number of produced DM quanta is proportional to 1/T at temperatures substantially below T 0 . The result is conveniently expressed in terms of the ratio of the DM number density to the SM entropy density, with g * s being the number of degrees of freedom contributing to the entropy. The correct DM relic abundance corresponds to The total DM yield from (28) depends on the temperature T * at which relativistic production stops. For m s > m h , we find that T * m s and subsequent non-relativistic production is negligible. In that case, Eq. 30 requires independently of m s . Here we have used g * 107 and also assumed that the reheating temperature has been high enough, T m s .
In this case, the annihilation mode makes a negligible contribution. The decay mode activates at T T EW and dominates the DM production. At T 20 GeV, the Higgs number density becomes too small and the process terminates.
The numerical values of the necessary couplings (31),(32) are close to those obtained in [7]. The correction factors derived in our work tend to compensate each other leaving the original estimate almost intact.

Non-thermalization constraint
An essential ingredient in our considerations is the assumption of non-thermalization of dark matter. Since thermal equilibrium implies detailed balance, we must require which is satisfied for where n eq (T ) is the equilibrium density at the temperature of the SM thermal bath. Our solution (28) for n(T ) then implies the constraint where we have taken T ∼ m s as the lowest temperature at which relativistic production can take place. The couplings in the range of interest satisfy this constraint. The thermalization bound on λ s has recently been obtained in [3]. The calculation is quite complicated since it involves a 2 → 4 process and depends on the DM number density, while for our purposes it is sufficient to take a conservative bound λ s < 10 −5 .

Conclusion
The aim of this work is to study the freeze-in production of dark matter in the relativistic regime, focussing on the simplest Higgs portal model. To this end, we have derived the relativistic reaction rates with the Bose-Einstein distribution function and solved analytically as well as numerically the corresponding Boltzmann equation. Obtaining a physically meaningful result requires inclusion of the Higgs thermal mass correction.
We have also made several other improvements over previous analyses. In particular, we take into account the electroweak phase transition. This implies, for example, that above the critical temperature the trilinear Higgs-DM vertices do not exist, which forbids Higgs decay and gauge boson annihilation into dark matter. The relevant DM production mode is instead annihilation of 4 massive Higgs degrees of freedom through the Higgs portal vertex. The required couplings leading to the correct relic abundance are given by Eqs. 31,32.
The overall numerical effect of our improvements happens to be modest compared to a straightforward (albeit unjustified) extrapolation of the non-relativistic result to high temperatures. Nevertheless, our systematic approach results in a more reliable evaluation of the dark matter yield and the relevant constraints on parameter space.