Relativistic Freeze-in

TL;DR

This paper investigates the production of scalar dark matter through the relativistic freeze-in process in the Higgs portal model, accounting for thermal effects and phase transition impacts to derive relevant constraints.

Contribution

It provides the first detailed derivation of relativistic reaction rates for scalar dark matter production, including thermal mass effects and electroweak phase transition considerations.

Findings

Derived relativistic reaction rates with Bose-Einstein statistics

Established constraints on Higgs portal coupling

Analyzed thermal mass effects and phase transition impacts

Abstract

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.

Figures

Equations

$$V_{hs}=\frac{1}{2}{\lambda }_{hs}{H}^{†}H{s}^2,$$

$$V_s=\frac{1}{4!}{\lambda }_s{s}^4+\frac{1}{2}{m}_{s0}^2{s}^2.$$

$$V_{hs}=\frac{1}{4}{\lambda }_{hs}\sum _i{h}_i^2{s}^2.$$

$$h_i{h}_i\to ss.$$

$$V_{hs}=\frac{1}{4}{\lambda }_{hs}{h}^2{s}^2+\frac{1}{2}{\lambda }_{hs}vh{s}^2+\frac{1}{4}{\lambda }_{hs}{v}^2{s}^2.$$

$$hh\to ss,h\to ss.$$

$$m_h^2\simeq m_{h0}^2+\left(\frac{3}{16}{g}_2^2+\frac{1}{16}{g}_1^2+\frac{1}{4}{y}_t^2+\frac{1}{2}{\lambda }_h\right)T^2,$$

$$m_s^2=m_{s0}^2+\frac{1}{2}{\lambda }_{hs}{v}^2.$$

$${\Gamma }_{a\to b}=\int \left(\prod _{i\in a}\frac{d^3{\mathbf{p}}_i}{(2\pi {)}^{3}2E_i}f(p_i)\right)\left(\prod _{j\in b}\frac{d^3{\mathbf{p}}_j}{(2\pi {)}^{3}2E_j}(1+f(p_j))\right)|{\mathcal{M}}_{a\to b}{|}^2(2\pi {)}^4{\delta }^4(p_a-p_b).$$

$$f(p)=\frac{1}{e^{\frac{u\cdot p}{T}}-1},$$

$${\Gamma }_{2\to 2}=(2\pi {)}^{-6}\int d^3{\mathbf{p}}_1{d}^3{\mathbf{p}}_{2}f(p_1)f(p_2)\sigma (p_1,p_2)v_{\mathrm{M}\text{o}\mathrm{l}},$$

$$v_{\mathrm{M}\text{o}\mathrm{l}}=\frac{F(p_1,p_2)}{E_1{E}_2}.$$

$$\sigma (p_1,p_2)=\frac{1}{4F(p_1,p_2)}\int |{\mathcal{M}}_{2\to 2}{|}^2(2\pi {)}^4{\delta }^4(p_1+p_2-k_1-k_2)\prod _i\frac{d^3{\mathbf{k}}_i}{(2\pi {)}^{3}2E_{k_i}},$$

$$p=\frac{p_1+p_2}{2},k=\frac{p_1-p_2}{2},$$

$$p^0=E\cosh \eta ,$$

$$p^1=E\sinh \eta \sin \theta \sin \varphi ,$$

$$p^2=E\sinh \eta \sin \theta \cos \varphi ,$$

$$p^3=E\sinh \eta \cos \theta ,$$

$$\int \frac{d^3{\mathbf{p}}_1}{2E_1}\frac{d^3{\mathbf{p}}_2}{2E_2}g(p_1,p_2)=2{\int }_{m_h}^{\infty }dE\sqrt{E^2-m_h^2}{E}^2{\int }_0^{\infty }d\eta {\sinh }^2\eta \int d{\Omega }_{p}d{\Omega }_{k}g(p_1,p_2),$$

$$u\cdot p_1=E\cosh \eta +\sqrt{E^2-m_h^2}\sinh \eta \cos {\theta }_k,$$

$$u\cdot p_2=E\cosh \eta -\sqrt{E^2-m_h^2}\sinh \eta \cos {\theta }_k,$$

$$\int d{\Omega }_{p}d{\Omega }_{k}f(p_1)f(p_2)=8{\pi }^2\frac{2T}{\sqrt{E^2-m_h^2}\sinh \eta (e^{\frac{2E\cosh \eta }{T}}-1)}\ln \frac{\sinh \frac{E\cosh \eta +\sqrt{E^2-m_h^2}\sinh \eta }{2T}}{\sinh \frac{E\cosh \eta -\sqrt{E^2-m_h^2}\sinh \eta }{2T}}.$$

$${\sigma }^{\mathrm{CM}}(E)=\frac{1}{2!2!}\frac{{\lambda }_{hs}^2}{64\pi E^2}\frac{\sqrt{E^2-m_s^2}}{\sqrt{E^2-m_h^2}},$$

$${\Gamma }_{2\to 2}$$

$$2m_h\Gamma =\int \left(\prod _{j\in f}\frac{d^3{\mathbf{p}}_j}{(2\pi {)}^{3}2E_j}\right)|{\mathcal{M}}_{1\to 2}{|}^2(2\pi {)}^4{\delta }^4(p-p_f),$$

$${\Gamma }_{1\to 2}=\frac{{\lambda }_{hs}^2{v}^2{m}_h^2}{64{\pi }^3}\sqrt{1-\frac{4m_s^2}{m_h^2}}{\int }_1^{\infty }dx\frac{\sqrt{x^2-1}}{e^{\frac{m_h}{T}x}-1},$$

$$\Gamma =\frac{{\lambda }_{hs}^2{v}^2}{32\pi m_h}\sqrt{1-\frac{4m_s^2}{m_h^2}},$$

$${\Gamma }_{2\to 2}\propto T^4\ln \frac{T}{m_h}\to c{T}^4,$$

$$\dot{n}+3Hn=(4-3\theta (T_{\mathrm{EW}}-T))\times 2{\Gamma }_{2\to 2}+\theta (T_{\mathrm{EW}}-T)\times 2{\Gamma }_{1\to 2},$$

$$T\frac{dn}{dT}-3n+8\frac{{\Gamma }_{2\to 2}}{H}=0.$$

$$H=\sqrt{\frac{{\pi }^2{g}_{\ast }}{90}}\frac{T^2}{M_{\mathrm{Pl}}},$$

$$\kappa \equiv 8\frac{{\Gamma }_{2\to 2}}{H{T}^2}.$$

Full text

1. 1 Introduction

2. 2 The set–up

3. 3 Relativistic reaction rates

4. 3.1 $hh\rightarrow ss$ reaction

5. 3.2 $h\rightarrow ss$ reaction

6. 3.3 Implications

7. 4 The Boltzmann equation

8. 4.1 Analytic solution in the relativistic regime

9. 4.2 Non–thermalization constraint

10. 5 Conclusion

1 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].111Ref. [13] obtains a different parametrization for the relativistic rates. However, it does not take into account the thermal mass corrections nor the effect of the electroweak phase transition.

2 The set–up

We consider the simplest real scalar dark matter model with the Higgs portal interaction

[TABLE]

where $s$ is a real scalar with the potential

[TABLE]

Here, the symmetry $s\rightarrow-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 $\lambda_{hs}\ll 1$. For the coupling values in the ballpark of $\lambda_{hs}\sim 10^{-11}$, one obtains the correct DM relic abundance, while DM never thermalizes (assuming also that $\lambda_{s}$ is small enough).

There are two reactions in which DM is produced: Higgs annihilation $hh\rightarrow ss$ and Higgs decay $h\rightarrow ss$, if allowed kinematically. The latter is only possible at temperatures below the electroweak scale, so we distinguish

(a) symmetric phase, $T\gtrsim T_{\rm EW}$ and $\langle H\rangle=0$

(b) broken phase, $T<T_{\rm EW}$ and $\langle H\rangle\not=0$

 

Here $T_{\rm EW}\sim 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

[TABLE]

In this case, dark matter is produced through annihilation,

[TABLE]

In the broken phase, one can use the unitary gauge $H=(0,(h+v)/\sqrt{2})^{T}$, which is singular at $\langle H\rangle=0$, such that

[TABLE]

This allows for both annihilation and decay, if allowed kinematically,

[TABLE]

Note that only one Higgs degree of freedom contributes in this case.222The would-be Goldstone boson contributions are recovered in the high energy limit through longitudinal gauge boson scattering. In this case, the Higgs propagator is cancelled by the gauge boson polarization vectors producing an effective quartic vertex.

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\sim T_{\rm EW}$ and approximate the reaction rates using the step function $\theta(T-T_{\rm 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

[TABLE]

with $m_{h0}$ being the zero temperature Higgs mass; $g_{1,2},y_{t}$ are the gauge and top quark Yukawa couplings, and $\lambda_{h}$ is the Higgs self–coupling. In the symmetric phase, $m_{h0}^{2}=-\lambda_{h}v^{2}$ in the convention $V_{h}=\lambda_{h}(h^{2}-v^{2})^{2}/4$. 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 $\lambda_{hs}$).

We also note that the DM mass changes during the phase transition and receives an extra contribution,

[TABLE]

For $\lambda_{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_{s}^{2}$ matters. For heavier DM, we make no distinction between $m_{s0}$ and $m_{s}$.

3 Relativistic reaction rates

In this section, we compute the $2\rightarrow 2$ and $1\rightarrow 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\rightarrow b}$ reaction rates per unit volume are

[TABLE]

Here ${\cal M}_{a\rightarrow 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

[TABLE]

where $u_{\mu}$ 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.

3.1 $hh\rightarrow ss$ reaction

Consider the reaction $hh\rightarrow ss$ due to the quartic interaction term. Following Gelmini and Gondolo [19], we write it as

[TABLE]

where the Møller velocity for the incoming particles is given by

[TABLE]

Here $F(p_{1},p_{2})=\sqrt{(p_{1}\cdot p_{2})^{2}-m_{h}^{4}}$ and the cross section is defined by

[TABLE]

where, in our convention, $|{\cal M}_{2\rightarrow 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

[TABLE]

and use the parametrization in terms of the half–the–CM energy $E$ and rapidity $\eta$,

[TABLE]

where $\theta$ and $\phi$ are the polar and azimuthal angles, respectively. This corresponds to $p=\Lambda(E,0,0,0)^{T}$ with $\Lambda$ being a Lorentz transformation. As shown in [3], for any $g(p_{1},p_{2})$,

[TABLE]

where in the integrand one must set $k_{0}=0,|{\bf k}|=\sqrt{E^{2}-m_{h}^{2}}$ in $k$-dependent quantities, and $\Omega_{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,

[TABLE]

where $k^{3}=|{\bf k}|\cos\theta_{k}$. Computing the integrals, we get

[TABLE]

The QFT cross section for the process $hh\rightarrow ss$ is

[TABLE]

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

[TABLE]

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.

3.2 $h\rightarrow ss$ reaction

Let us now consider the decay mode $h\rightarrow ss$. Again, it is convenient to go over to the rest frame of the decaying particle with momentum $p$, so that $p=\Lambda(E,0,0,0)^{T}$. Parametrizing $p$ as we did in the previous subsection, we have ${d^{3}{\bf p}\over 2E}=\delta(E^{2}-m_{h}^{2})\;\sinh^{2}\eta\;E^{3}$ $\;dE~{}d\eta~{}d\Omega$. Using the decay width $\Gamma$ definition

[TABLE]

we find

[TABLE]

where we have used

[TABLE]

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].

3.3 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}\rightarrow 0$. Let us set $m_{s}=0$ and consider the limit $m_{h}\rightarrow 0$. In this case, we find that the $2\rightarrow 2$ rate diverges as $\ln m_{h}$ which is unphysical. Including the thermal mass, we get

[TABLE]

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}$.

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 form

[TABLE]

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 $\theta$–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 $WW\rightarrow h\rightarrow 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\rightarrow ss$ in the broken phase, its effect is inconsequential.

4.1 Analytic solution in the relativistic regime

Qualitative behaviour of $n(t)$ in the relativistic regime at $T>T_{\rm 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 $\dot{T}=-HT$ and the Boltzmann equation can be written as

[TABLE]

Here the Hubble rate is given by

[TABLE]

with $g_{*}$ being the number of SM degrees of freedom and $M_{\rm 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

[TABLE]

Imposing the boundary condition that the DM abundance vanish at some initial temperature $T_{0}$, the solution to the Boltzmann equation reads

[TABLE]

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,

[TABLE]

with $g_{*s}$ being the number of degrees of freedom contributing to the entropy. The correct DM relic abundance corresponds to

[TABLE]

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_{*}\simeq m_{s}$ and subsequent non–relativistic production is negligible. In that case, Eq. 30 requires

[TABLE]

independently of $m_{s}$. Here we have used $g_{*}\simeq 107$ and also assumed that the reheating temperature has been high enough, $T\gg m_{s}$.

Our numerical solution to the full Boltzmann equation (24) supports this result. The evolution of $Y(T)$ for representative values of $m_{s}$ is shown in Fig. 2. It grows as $1/T$ at high temperature and freezes–in at $T_{*}\sim m_{s}$ if $m_{s}>m_{h}$. For $m_{s}\ll m_{h}$, we find that the required coupling is

[TABLE]

In this case, the annihilation mode makes a negligible contribution. The decay mode activates at $T\lesssim T_{\rm EW}$ and dominates the DM production. At $T\lesssim 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.

4.2 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

[TABLE]

which is satisfied for

[TABLE]

where $n_{\rm eq}(T)$ is the equilibrium density at the temperature of the SM thermal bath. Our solution (28) for $n(T)$ then implies the constraint

[TABLE]

where we have taken $T\sim 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 $\lambda_{s}$ has recently been obtained in [3]. The calculation is quite complicated since it involves a $2\rightarrow 4$ process and depends on the DM number density, while for our purposes it is sufficient to take a conservative bound $\lambda_{s}<10^{-5}$.

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.