Matter Power Spectrum of Light Freeze-in Dark Matter: With or without Self-Interaction

TL;DR

This paper investigates how light freeze-in dark matter affects the matter power spectrum, considering both scenarios with negligible and strong self-interactions, and finds that self-scattering can protect small-scale structures from free-streaming suppression.

Contribution

It introduces a numerical analysis of self-scattering effects on the matter power spectrum in light freeze-in dark matter models, highlighting the importance of decoupling time.

Findings

A 24 keV freeze-in dark matter is excluded by current WDM bounds.

Strong self-scattering can protect the power spectrum from free-streaming suppression.

Self-scattering decoupling time is a key parameter in structure formation.

Abstract

We study the free-streaming effect in a light freeze-in dark matter model. Naturally in the dark sector one can find dark matter related coupling, and such coupling may induce dark matter self-scattering. In case that such scattering is subdominant, the dark matter partition function is not thermal but determined by the freeze-in process, yet its high momentum side is generally also Boltzmann suppressed. We show that the matter power spectrum is very similar to a warm dark matter one in shape. When matched to the current WDM bound, a $24$~keV freeze-in dark matter is ruled out at $2\sigma$ confidence level. In case that the dark matter self-scattering is strong and decouples at a very late time, by a new numerical calculation we show that the early stage Brownian motion indeed protects the power spectrum against free-streaming suppression. However, such an effect cannot be characterized by a free-streaming length alone; we find that the self-scattering decoupling time is another necessary parameter. The currently interested dark matter self-interaction cross section $\sim\text{cm}^2/\text{g}$ is just marginal for such protection to be effective.

Equations

$$\mathcal{L}\supset \frac{1}{2}({\partial }_{\mu }\Phi {)}^2+\frac{1}{2}{\mu }^2{\Phi }^2-\frac{1}{4}{\lambda }_{\varphi }{\Phi }^4-\frac{1}{2}{\lambda }_{h\varphi }{H}^{†}H{\Phi }^2.$$

$$\frac{d{f}_{\varphi }}{d{E}_{\varphi }}=2{\int }_{E_{h\text{min}}(E_{\varphi })}^{E_{h\text{max}}(E_{\varphi })}\frac{d{E}_h}{\sqrt{(E_h^2-m_h^2)(1-4r^2)}}\frac{d{f}_h}{d{E}_h},$$

$$\frac{d}{dx}\frac{f_h}{s}=-\frac{{\Gamma }_{h\to \varphi \varphi }}{xH(x)}\frac{m_h}{E_h}\frac{f_h}{s},$$

$$\frac{d^{2}(f_{\phi}/s)}{dxdp_{\phi}}=\frac{2(p_{\phi}/E_{\phi})\Gamma_{h\to\phi\phi}}{xH(x)\sqrt{1-4r^{2}}}\int_{p_{h\text{min}}}^{p_{h\text{max}}}\frac{dp_{h}}{E_{h}^{2}}\Big{(}\frac{m_{h}}{E_{h}}+\frac{m_{h}}{T}\Big{)}\frac{e^{-\frac{E_{h}}{T}}}{s},$$

$$\frac{d(f_{\varphi }/s)}{d{y}_{\varphi }}=$$

$$\times$$

$${\partial }_{\tau }{\Psi }_0=$$

$${\partial }_{\tau }{\Psi }_1=$$

$${\partial }_{\tau }{\Psi }_2=$$

$${\partial }_{\tau }{\Psi }_{\ell }=$$

$$\lambda ={\int }_0^{t_{\text{eq}}}v(t)\frac{dt}{a(t)}={\int }_0^{a_{\text{eq}}}v(a)\frac{da}{a^{2}H(a)}.$$

$$\frac{T}{⟨p⟩}⟨v_0⟩=$$

$$=$$

$$m_{\text{FIDM}}\geq 24.3\Big{(}\frac{106.75}{g_{\ast s\text{FIDM}}}\Big{)}^{\frac{1}{3}}\Big{(}\frac{m_{\text{WDM}}}{6.3~{}\text{keV}}\Big{)}^{\frac{4}{3}}~{}\text{keV}.$$

$$\frac{{\rho }_c{\Omega }_{\chi }}{a_{\text{B}}^3}\frac{\sigma }{m}\frac{⟨p⟩}{T}\frac{T_0}{a_{\text{B}}m}=H_0\sqrt{\frac{{\Omega }_{\chi +b}}{a_{\text{B}}^3}+\frac{{\Omega }_{\gamma +\nu }}{a_{\text{B}}^4}+{\Omega }_{\Lambda }}$$

Full text

**I. Introduction. ** In the past, Weakly Interacting Massive Particle has played a central role in dark matter (DM) model building. For a particle with weak scale mass and more importantly weak scale interaction with the standard model (SM) particle ($\mathcal{O}(1)$ pb annihilation cross section), the relic abundance $\Omega_{\chi}h^{2}=0.1186$ Aghanim et al. (2018) as the most explicit measurement by far for DM arises naturally by the so-called thermal freeze-out mechanism Kolb and Turner (1990). However, current direct Aprile et al. (2018); Cui et al. (2017), indirect Ackermann et al. (2017) and collider Boveia and Doglioni (2018) searches for such a particle have substantially probed its parameter space, and the null results have reduced the motivation for this mechanism. Alternatively, freeze-in mechanism Ellis et al. (1984); McDonald (2002); Choi and Roszkowski (2005); McDonald and Sahu (2009); Hall et al. (2010) can also parametrically give correct relic abundance, in which DM is not present at the very beginning, but produced through scattering with the visible sector particles (see Bernal et al. (2017) and references therein).

Such DM particles naturally inherit the thermal motion of the visible sector particles. For light DM candidate which is becoming more popular either from direct detection perspective Knapen et al. (2017); Ibe et al. (2018) or motivated by potential DM indirect detection signal Babu and Mohapatra (2014); Biswas et al. (2018); Heeba et al. (2018) or even more astrophysically McDonald (2002); de Vega et al. (2014); An et al. (2018), a significant free-streaming (FS) effect is expected. However the freeze-in DM (FIDM. In literature it is often called “Feebly Interacting Massive Particle”, but here it is neither massive nor necessarily feeble in self-interaction.) is usually taken as nonthermal with its partition function determined in the freeze-in process Roland and Shakya (2017); Konig et al. (2016); Biswas and Gupta (2017); Bae et al. (2018), so for its matter power spectrum (MPS) strictly one cannot directly use a warm DM (WDM) one which is easily available. On the other hand, in DM model building it is natural to find (even order one) couplings in the dark sector (e.g. see McDonald (2002)) which induces DM self-scattering Spergel and Steinhardt (2000); Tulin and Yu (2018), and such scattering tends to bring the dark sector back into thermal equilibrium. Whether such self-scattering can thermalize the FIDM is up to parameters.

In this Letter, we will consider two opposite limits in detail in which the MPS of FIDM is determined differently. In the first case that the self-scattering is ineffective, we will calculate the partition function from the first principle and show that it is also Boltzmann suppressed on the high momentum side, then the MPS difference Konig et al. (2016) caused by the difference of the shape of the partition distribution compared to a thermal distribution is in general negligible. In the opposite limit that the self-scattering is very fast, so the FS gives way to Brownian motion early on, we will show that a total FS length is insufficient to characterize the MPS suppression, and a second parameter of the self-scattering decoupling time/scale factor will be necessary.

**II. The Weak Self-Scattering Limit. ** While the generic feature revealed by our calculation does not depend on specific model realization, we still use a certain freeze-in model for convenience. Our example model Yaguna (2011); Babu and Mohapatra (2014); Heeba et al. (2018) contains a real scalar gauge singlet extension of the SM Higgs sector

[TABLE]

Here $H$ is the SM Higgs doublet, and $\Phi$ can be decomposed into $\phi+u$, the particle and the vacuum expectation value. After both electroweak symmetry breaking and the $u$ turning on, the mixing angle between the SM Higgs and the new $\Phi$ sector $\theta\approx\frac{\lambda_{h\phi}uv}{m_{h}^{2}-m_{\phi}^{2}}$ is suppressed by the “feeble” interaction of $\lambda_{h\phi}$.

This simple model is motivated by the $3.55$ keV X-ray extragalactic line excess, since if $m_{\phi}=7.1~{}$keV then decay of $\phi$ through the SM Higgs portal is only kinetically open to photon pairs in the SM. In Heeba et al. (2018) the freeze-in process is calculated in detail, the post-electroweak phase transition freeze-in production is found to be dominant, and the self-scattering is insufficient to bring the DM sector back into thermal equilibrium through the $2\to 3$ process. The simplification here compared to models in Konig et al. (2016); Biswas and Gupta (2017) is that we do not need a new portal, as well as its independent Boltzmann equation.

We start the calculation of partition function with focus on the energy of the outgoing $\phi$. Let $r=\frac{m_{\phi}}{m_{h}}$, boosting the 4-momentum of the outgoing $\phi$ from the Higgs rest frame (in which $p_{\phi}^{\mu}=(\frac{m_{h}}{2},~{}\frac{m_{h}}{2}(1-4r^{2})^{\frac{1}{2}}\hat{r})$ and $\hat{r}=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$ is the unit vector in spherical coordinate) into a frame in which the Higgs energy is $E_{h}$, we get $E_{\phi}=\frac{1}{2}(E_{h}+((E_{h}^{2}-m_{h}^{2})(1-4r^{2}))^{\frac{1}{2}}\cos\theta)$. In this frame the spectrum distribution with $E_{\phi}$ is uniform since $d(\text{probability})/dE_{\phi}\propto\sin\theta d\theta d\phi/d\cos\theta$ is $\theta$ independent. So for differential partition function in the freeze-in process,

[TABLE]

where the factor $2$ counts the two $\phi$s produced in a single decay event, and $E_{h\text{min}}(E_{\phi})$ ($E_{h\text{max}}(E_{\phi})$) is the minimal (maximal) $E_{h}$ value to still give a certain $E_{\phi}$ energy.

In the expanding universe the Boltzmann equation is \frac{1}{a^{3}}\frac{d}{dt}\big{(}a^{3}f_{h}\big{)}=-\frac{m_{h}}{E_{h}}\Gamma_{h\to\phi\phi}f_{h}, where $a$ is the scale factor with normalization $a_{0}=1$, and we have simplified by using $\Gamma_{h\to\phi\phi}=\frac{1}{2m_{h}}\int\frac{d^{3}p_{\phi_{1}}d^{3}p_{\phi_{2}}(2\pi)^{4}\delta^{(4)}(p_{h}-p_{\phi_{1}}-p_{\phi_{2}})}{(2\pi)^{6}2E_{\phi_{1}}2E_{\phi_{2}}}\sum|\mathcal{M}|^{2}$, neglecting the inverse process of freeze-in. Let $x=\frac{m_{h}}{T}$, then the $x$-based Boltzmann equation reads

[TABLE]

where $s$ is the entropy density. Plugging Eq. 3 into the right hand side of Eq. 2 we will get the partition function of $\phi$ determined by the freeze-in process

[TABLE]

where we have used $pdp=EdE$ and Boltzmann distribution for Higgs.

Next let $x_{\phi}=\frac{p_{\phi}}{m_{\phi}}$, $x_{h}=\frac{p_{h}}{m_{h}}$ and $y_{\phi}=\frac{p_{\phi}}{T}=\frac{p_{\phi}}{m_{\phi}}\frac{m_{\phi}}{m_{h}}\frac{m_{h}}{T}=rx_{\phi}x$. Apparently $y_{\phi}\sim 1$, $x_{h}\sim 1$ and $x_{\phi}\sim r^{-1}\gg 1$. Note that $y_{\phi}$ also has the advantage that it is Hubble expansion invariant, it is the right variable for the FS partition function. Then integrating over $x$ we get

[TABLE]

where $x_{hl}=\frac{1}{2r}|\sqrt{(1+(\frac{y_{\phi}}{rx})^{2})(1-4r^{2})}-\frac{y_{\phi}}{rx}|$ is the integration lower bound, and the upper bound has been relaxed to infinity with help of Boltzmann suppression of Higgs.

Then we can calculate the differential freeze-in partition function numerically by Eq. Matter Power Spectrum of Light Freeze-in Dark Matter: With or without Self-Interaction. Here the reference temperature $T$ can be determined by the entropy conservation $g_{\ast s}(aT)^{3}=\text{constant}$, with reference values today $g_{\ast s,0}=2+\frac{7}{8}*2*3.046*\frac{4}{11}=3.93$ and $T_{0}=T_{\text{CMB}}=2.7255$ K. For arbitrary new physics, $g_{\ast s}$ should be taken as a free parameter. The SM gives $g_{\ast s}=g_{\ast}=106.75$ for universe early enough, and in our case, right after the electroweak phase transition it is slightly smaller Husdal (2016). In Fig. 1, we show the numerical result and compare it to the thermal differential partition functions. The overall normalization of the freeze-in differential partition function is determined by reproducing the relic abundance; here we choose a normalization similar in scale to the thermal ones for convenience of comparison.

The most useful feature for us is that on the large $y_{\phi}=\frac{p_{\phi}}{T}$ side the partition function is also suppressed exponentially, closely following the thermal distributions. All freeze-in partition function calculations found similar behavior due to the inevitable Boltzmann suppression of either parent particle or the center of mass energy, for the freeze-in process of both UV type and IR type Roland and Shakya (2017), and regardless of whether the portal achieves thermal equilibrium Konig et al. (2016); Biswas and Gupta (2017). In small $y_{\phi}$ regime, we can see a deviation from thermal distributions, and different models can give quite different shapes Konig et al. (2016); Biswas and Gupta (2017). However the bulk of FIDM particles resides on the $\frac{p_{\phi}}{T}\gtrsim 1$ region and dictates the MPS shape, so in general, the low momentum side is not important and can be approximately replaced by a thermal one such as the Fermi-Dirac one in the WDM model. We will see how well it works soon. A partition function that cannot be approximated by a single thermal distribution should be possible, in case that the multiple components with high and low momenta are comparable. However, such case is somewhat tuned in parameters, and the multiple components should be able to get approximated by multiple components of WDM. So in the following we will use this differential partition functions with $m_{h}=125$ GeV and $m_{\phi}=7$ keV for illustration, and the small shape difference for different $m_{\phi}$ value is ignored.

Further, we perform numerical analysis by solving the perturbation Boltzmann equation set for the evolution of the FS DM using camb Lewis et al. (2000). The general formulation follows Ma and Bertschinger (1995) with the choice of synchronous gauge. The DM follows the standard massive neutrino perturbative Boltzmann equation

[TABLE]

Here $\Psi_{\ell}$ is the fractional perturbation from the unperturbed partition function $f(\tilde{p})$ and the subscript labels the harmonic mode, $\tau$ is conformal time and $\tilde{p}=ap$ and $\tilde{E}=\sqrt{\tilde{p}^{2}+a^{2}m^{2}}$ are comoving momentum and energy respectively, $k$ is the interested Fourier mode, and $h$ and $\eta$ are the metric perturbation in synchronous gauge. Different from the redshifted Fermi-Dirac distribution of benchmark WDM model, the afore calculated differential partition function is inputted in the $\frac{d\ln f}{d\ln\tilde{p}}$ terms in the $\ell=0$ and $\ell=2$ equations.

In Fig. 2 we show the MPS ($\Delta(k)=\frac{k^{3}}{2\pi^{2}}P(k)$) of FIDM, in comparison with benchmark WDM ones. While tuning the characteristic FS velocities to be the same, we can see that the MPS of light DM is indeed not sensitive to the detailed shape of the unperturbed differential partition function such as the ones shown in Fig. 1 (at least in the case that the new partition function also contains a Boltzmann suppression on the large comoving momentum side)111In the context of small scale MPS suppression due to FS, another known class of unperturbed partition function arises in the late decaying DM model Huo (2011). However, the MPS can also be very similar to the thermal WDM one in shape while tuning the FS velocity. Note that an arbitrary unperturbed partition function may not match to a realistic particle physics model. On the other hand, if the small scale MPS suppression is not due to FS but due to acoustic oscillation in the dark sector, then the MPS shape can be quite different Cyr-Racine et al. (2016); Huo et al. (2018)., and rather indistinguishable from that of a thermal distribution. The more suppressed MPS of the $7$ keV FIDM model means that it is excluded. This exclusion can be understood by simply considering the FS velocity, which is the factor $\frac{\tilde{p}}{\tilde{E}}$ in Eq. 6. Viewing the Boltzmann equation set, when one ignores the small difference caused by non-thermal partition function such as that in our FIDM example, the FS velocity is the only factor which controls the MPS, and the more frequently used FS length Kolb and Turner (1990) is just its integration under certain fixed Hubble behavior from the bigbang till matter-radiation equality (or later)

[TABLE]

The convenient FS velocity to compare is the FS velocity extrapolated till today assuming no late virialization in structure formation,

[TABLE]

where for massless fermion $\frac{T}{\langle p\rangle}=\frac{180\zeta(3)}{7\pi^{4}}\approx 0.32$ and here for simplicity for all light DM we use this value and ignore the difference caused by partition function, and $g_{\ast s\text{FIDM}}$ is the effective $g_{\ast s}$ when the FIDM is produced. For the reason seen in the next section, we advocate it rather than the FS length. Then the matching above suggests that a WDM mass bound can be generally recast into an FIDM mass bound

[TABLE]

Here we have set the reference WDM bound to be the most recent one from the stellar stream observation Banik et al. (2019).

**III. The Strong Self-Scattering Limit. ** The random motion of a particle due to frequent (self) scattering is known as Brownian motion Brown (1828). The early DM Brownian motion in the universe will reduce the total FS length, by disabling FS in a specific direction from the bigbang to the decoupling of self-scattering. Here we show that such a non-free-streaming Brownian early stage will not work very effectively in protecting the small scale MPS from suppression.

As a simple approximation, we will set the FS velocity $\frac{\tilde{p}}{\tilde{E}}$ in Eq. 6 to be zero before the Brownian stage decoupling, which results in the same behavior with the Cold DM (CDM). Moreover, since elastic scattering will not statistically change the DM velocity with the latter only subject to redshift, after decoupling the FS velocity will retake the corresponding value as if there is no self-scattering. This is different from the old self-interacting DM MPS calculation de Laix et al. (1995), since here we do not have heating in the $3\to 2$ process. Such working approximation is easily implemented in camb Lewis et al. (2000), by turning off the DM velocity before a specific decoupling time.

In Fig. 3 in the above approximation we calculate the MPS for the models with an early Brownian stage, which shares the same late-stage FS velocity with the $5.3$ keV thermal WDM model Irsic et al. (2017). Similar to the previous analysis, although the equilibrium partition may not be a Fermi-Dirac one but be a bosonic one or distorted by chemical potential, we can still use the WDM one for convenience, and the difference should be negligible. We can see that the MPS with early Brownian stages are almost indistinguishable from the benchmark thermal WDM model without an early Brownian stage, unless the decoupling happens very late till when DM becomes very non-relativistic. If the FS length is calculated from Eq. 11 with the bigbang as the lower integration bound replaced by the Brownian decoupling, then actually till as high as a $50\%$ (the green curve) reduction is achieved, one can only see a mild protection against FS suppression. Moreover, the reduction of FS length in this way will not correspond to a unique MPS with the same total FS length but no Brownian early stage, for example matching the FS length corresponding to a Brownian decoupling DM velocity $v=10^{-3}~{}c$ we get $m_{\text{WDM}}=12.8$ keV for the thermal WDM model, but the MPS suppression for the latter is seen at a significantly smaller scale.

Apparently, the total FS length is not an excellent way to parameterize this effect. It is because the perturbation growth is not linear with the accumulation of FS length but the late time dominates. Even if an early Brownian stage prevents FS and makes the perturbation grow as effectively as CDM, such growth is driven by the dominant radiation-gravitation couple, and the magnitude is small at early times. And after the Brownian decoupling, the perturbation growth which is fast and dominant for CDM will still be erased by FS for light DM. On the other hand, if the self-scattering decoupling is pushed to the late side all the way till the matter radiation equilibrium or even later, we can indeed see the reduction of the FS effect.

To better quantify the effect of the Brownian stage for MPS, we treat the Brownian motion decoupling scale factor $a_{\text{B}}$ as a convenient parameter, and scan the two parameters and compare with the current Lyman-$\alpha$ bound of thermal WDM. Given the DM self-scattering cross section which is usually expressed as $\frac{\sigma}{m}$, $a_{\text{B}}$ can be solved through the $\Gamma\sim H$ decoupling criteria, namely the equation

[TABLE]

For a DM self scattering cross section $\frac{\sigma}{m}\sim\text{cm}^{2}/\text{g}$ Tulin and Yu (2018) and an extrapolated FS velocity $\frac{T_{0}}{m}\sim 10^{-8}$, interestingly $a_{\text{B}}$ is around the matter radiation equilibrium and the reduction of FS is visible. The result is shown in Fig. 4. Constraints are made based on the $\delta A$ criteria Murgia et al. (2017): we first calculate the one-dimensional power spectrum of $P_{\text{1d}}(k)=\frac{1}{2\pi}\int_{k}^{\infty}k^{\prime}dk^{\prime}P(k^{\prime})$ from three-dimensional MPS $P(k)$, then calculate \delta A=1-\int_{k_{\text{min}}}^{k_{\text{max}}}dk\big{(}P_{\text{1d}}(k)/P_{\text{1d}}^{\text{CDM}}(k)\big{)}/(k_{\text{max}}-k_{\text{min}}) with $k_{\text{min}}=0.5h/\text{Mpc}$ and $k_{\text{max}}=20h/\text{Mpc}$ and compare it to the $3.5$ keV ($5.3$ keV) thermal WDM model result. Such criteria should not be directly applied to the more recent WDM bound of $6.3$ keV based on stellar stream $+$ dwarf satellite count Banik et al. (2019), and similar recast work is beyond the current work.

**IV. Summary. ** In this letter we have explored several aspects of the light freeze-in dark matter model, of its partition function determined by the freeze-in process without sufficient self-scattering, of the numerically calculated matter power spectrum shape and its similarity to a warm dark matter one at the same free-streaming velocity, and oppositely that whether a fast self-scattering or Brownian early stage will result in a protection from free-streaming suppression, and how such protection can be parameterized. We can see a generic similarity to the warm dark matter free-streaming effect, characterized by the extrapolated free streaming velocity, as well as the need of introducing Brownian decoupling time as the second parameter in the latter case. The light freeze-in dark matter model is also severely constrained by recast WDM bound, currently with a bound of about $24$ keV.

Acknowledgments: We are grateful to useful discussions with Haipeng An and Hai-Bo Yu.