A Preferred Mass Range for Primordial Black Hole Formation and Black Holes as Dark Matter Revisited

TL;DR

This paper revisits the idea that primordial black holes within a specific mass range could account for dark matter, emphasizing the role of non-thermal histories and spectral index in their formation and abundance.

Contribution

It demonstrates that a preferred mass range for primordial black holes as dark matter is robust under certain spectral conditions and discusses formation likelihood and relic density constraints.

Findings

Black holes in a specific mass range could constitute all dark matter.

Non-thermal histories favor formation of black holes in this mass range.

Observational constraints limit the relic density of primordial black holes.

Abstract

Bird, et. al. and Sasaki, et. al. have recently proposed the intriguing possibility that the black holes detected by LIGO could be all or part of the cosmological dark matter. This offers an alternative to WIMPs and axions, where dark matter could be comprised solely of Standard Model particles. The mass range lies within an observationally viable window and the predicted merger rate can be tested by future LIGO observations. In this paper, we argue that non-thermal histories favor production of black holes near this mass range -- with heavier ones unlikely to form in the early universe and lighter black holes being diluted through late-time entropy production. We discuss how this prediction depends on the primordial power spectrum, the likelihood of black hole formation, and the underlying model parameters. We find the prediction for the preferred mass range to be rather robust assuming a blue spectral index less than two. We consider the resulting relic density in black holes, and using recent observational constraints, establish whether they could account for all of the dark matter today.

Figures

Equations

$$t_{\text{mbox}PBH}\le t_{\text{mbox}caustic},$$

$$W_{\text{mbox}caustic}\simeq {\delta }_{\text{mbox}M}^{3/2},$$

$${\delta }_{\text{mbox}M}\equiv \frac{\delta M}{M}\sim M^{\frac{(1-n)}{6}},$$

$$W_s\simeq C{\delta }_{\text{mbox}M}^p.$$

$$\beta (M)=\frac{{\rho }_{\text{mbox}PBH}(M)}{{\rho }_{\text{mbox}tot}}\simeq W_s{W}_{\text{mbox}caustic}\simeq C{\delta }_{\text{mbox}M}^p{\delta }_{\text{mbox}M}^{3/2}.$$

$${\dot{\rho }}_m+3H{\rho }_m$$

$${\dot{\rho }}_r+3H{\rho }_r$$

$$3H^2{m}_p^2$$

$${\rho }_m={\rho }_m^{(0)}\left(\frac{e^{-\Gamma t}}{a^3}\right).$$

$${\rho }_r$$

$$\Delta$$

$$\beta (t_r)$$

$$\beta (t_r)=\frac{4M}{3T_r}{\left(\frac{n}{s}\right)|}_{t=t_r},$$

$$Y_{\text{mbox}PBH}\equiv \frac{n}{s}=\left(\frac{3T_r}{4M}\right)\beta (t_r),$$

$${\Omega }_{\text{mbox}PBH}=\frac{{\rho }_{\text{mbox}PBH}}{{\rho }_c}=\frac{M{Y}_{\text{mbox}PBH}s}{{\rho }_c},$$

$${\Omega }_{\text{mbox}PBH}{h}^2$$

$$M_{\text{mbox}min}=\frac{3H_{\text{mbox}osc}^2{m}_p}{H_{\text{mbox}osc}^3}=3\frac{m_p}{H_{\text{mbox}osc}}\simeq 3\frac{m_p}{m_{\sigma }}.$$

$$M_{\text{mbox}max}={\alpha }^{\frac{1}{n+3}}{\left(\frac{M_{\text{mbox}C}}{m_p}\right)}^{\frac{n-1}{n+3}}{\left(\frac{m_p}{m_{\sigma }}\right)}^{\frac{12}{n+3}}{m}_p,$$

Full text

I Introduction

The question “Did LIGO detect dark matter?” was raised recently by the authors of Bird et al. (2016). Following the first LIGO observation of gravitational waves GW150914 Abbott and et. al. (2016), the authors examined the possibility that all of dark matter could consist of 30 solar mass black holes. Making this assumption, they calculated the expected merger rate from binaries forming in dark matter halos. Instead, in Sasaki et al. (2016) the authors examine a similar scenario (see also Orlofsky et al. (2016); Carr et al. (2016)), but allow the fraction of dark matter that is comprised of primordial black holes (PBHs) to be a free parameter. In both papers, it is noted that there are considerable astrophysical uncertainties associated with merger rates and the dark matter distribution, but given the uncertainties, their results can be consistent with the value measured by LIGO. Curiously, existing observational constraints leave a window in precisely this mass range (near a solar mass) to allow black holes to comprise all of the dark matter Bird et al. (2016). However, one important theoretical piece to the story is lacking: why are black holes in this mass range favored over other masses? We will argue that non-thermal histories offer a possible explanation.

A generic prediction of beyond the Standard Model (BSM) physics is the existence of additional scalar fields (moduli), which frequently predict an epoch of early matter domination following inflation. In many approaches the mass and decay rate of these scalars is determined by symmetry breaking processes near the TeV scale – the physics of which has been well-studied, particularly for its implications for WIMPs and axions, collider physics, and fundamental questions such as the electroweak hierarchy and the Cosmological Moduli Problem (see Kane et al. (2015) for a review and guide to the literature).

In this paper, we are interested in another prediction of non-thermal histories during the post-inflationary epoch. Unlike in a thermal history, during a matter phase, PBHs can form on all sub-Hubble scales, resulting in a continuum of PBH masses Georg et al. (2016); Cotner and Kusenko (2016); Khlopov (2010); Polnarev and Khlopov (1981); Khlopov and Polnarev (1980); Polnarev and Khlopov (1981); Harada et al. (2016). Depending on the primordial power spectrum, this can lead to strong constraints on these non-thermal epochs and the underlying particle theory Georg et al. (2016). Another aspect of these models is that the phase ends when the scalar decays, again ‘reheating’ the universe. The associated entropy production dilutes any pre-existing relics (such as baryons or dark matter), but the decay can also be a replenishing source of baryons Allahverdi et al. (2011); Kane et al. (2011) or dark matter Moroi and Randall (2000); Acharya et al. (2008). Here, we examine the effect of the non-thermal epoch and scalar decay on the formation and survival of PBHs as dark matter.

The rest of the paper is organized as follows. In the next section, we review the formation of PBHs in an early matter phase and the resulting prediction for the dark matter relic density taking into consideration entropy production. In Section III, we establish an argument for a preferred mass range for the PBHs, and calculate the corresponding PBH abundance for each allowed mass. We then establish constraints on the allowed values of the spectral index of the primordial power spectrum, which is important for establishing both the mass and abundance of the most dominant PBHs today. We present our conclusions in Section IV.

II Black Hole Formation in Early Matter Phases

What sets a non-thermal universe apart from a thermal universe is that PBHs can continue to form after perturbations cross the Hubble radius. Matter (scalar) perturbations will continue to grow because the Jeans pressure is negligible during epochs of scalar oscillations. This means that the usual radiation domination calculations for PBH formation are not applicable.

In matter domination, the probability of forming PBHs stems from two different considerations. Firstly, the over-dense region must collapse to its Schwarzschild radius before a caustic can form at its center (from which particles falling inwards might be dispersed and prevent the formation of a PBH). This is also related to the need for the over-density to be sufficiently homogeneous once the perturbations break away from the background expansion (when they become non-linear) Khlopov (2010). This implies

[TABLE]

where $t_{\mbox{\tiny{PBH}}}$ is the PBH formation time and $t_{\mbox{\tiny{caustic}}}$ is the time it takes a caustic to form. In Polnarev and Khlopov (1981) it is shown that the probability that (1) is obeyed can be approximated by

[TABLE]

where the density contrast is given by

[TABLE]

where $n$ is the spectral index.

The second consideration takes into account if the gravitational collapse is sufficiently spherical and therefore results in a point like singularity and not a 1-dimensional (cigar) or 2-dimensional (Zel’dovich pancake) singularity. This probability was calculated in Doroshkevich (1970) and is

[TABLE]

This result comes from analytically evolving the initial perturbations to the non-linear regime and then tracking their behavior through modeling the collapse as a Tolman-like solution. One can introduce a deformation tensor with random values (reflecting the nearly Gaussian initial spectrum) and then use this to find the probability of spherical collapse. After diagonalizing this tensor and finding the corresponding eigenvalues one arrives at the likelihood for PBH formation Polnarev and Khlopov (1981). Following Doroshkevich (1970) the authors of Polnarev and Khlopov (1981) find $p=5$ and $C=2\times 10^{-2}$. In what follows we will allow for a range of values for $C$ and $p$ and see that this uncertainty will not alter our main conclusions.

The mass fraction in PBHs of mass $M$ is then given by the total probability for their formation by

[TABLE]

In a radiation dominated universe $\rho_{\mbox{\tiny tot}}\sim 1/a^{4}$, whereas for PBHs which have not evaporated ($M>10^{16}$ g) we have $\rho_{\mbox{\tiny PBH}}\sim 1/a^{3}$. Thus, the mass fraction grows as $\beta\sim a(t)$. However, in a matter dominated phase $\beta$ remains constant until the time of reheating as long as there is not significant entropy production. Notably, in the scalar dominated phases considered here, entropy production plays an important role, which we now discuss.

II.1 Entropy Production in a Non-thermal History

In this section, we establish the effect of entropy production on PBH formation. To be specific, we will focus on models where the non-thermal phase results from the displacement of scalar moduli following inflationary reheating (see Kane et al. (2015) for a review). Following reheating, the universe is radiation dominated. However, the scalar oscillations will quickly come to dominate the energy density since their energy density dilutes more slowly than radiation. The non-thermal phase begins around the time when the scalar oscillations become comparable to the radiation density and we will denote the radiation density at this time by $\rho_{r}^{(0)}$. Then, the radiation density during the non-thermal epoch will be given by $\rho_{r}=\rho_{r}^{(0)}+\rho_{\mbox{\tiny gen}}$, where the additional contribution represents the radiation generated by scalar decays111We follow closely Gorbunov and Rubakov (2011) and refer to that textbook for a more detailed discussion..

The relevant equations are

[TABLE]

We can solve (6)

[TABLE]

We are interested in establishing when entropy production first becomes important and so we consider times $t_{\ast}<t\ll\Gamma^{-1}$. Thus, we have $\Gamma t\ll 1$ and so $\rho_{m}\sim a^{-3}$ or $\rho_{m}=4m_{p}^{2}/(3t^{2})$. This corresponds to a matter dominated universe222This is not strictly correct since the radiation gives an equivalent contribution at $t_{\ast}$. However, as is the case with the textbook treatment of standard radiation/matter equality it is a very good approximation and adequate for our purposes here (see e.g. Chapter 2 of Padmanabhan (1993)). $H\simeq 2/(3t)$. Using this solution in (7) and solving, we find

[TABLE]

where $\rho_{\mbox{\tiny init}}$ and $\rho_{\mbox{\tiny gen}}$ are the initial radiation and radiation generated from scalar decays, respectively. In arriving at (10) we have used that $\rho_{m}(t_{\ast})=\rho_{r}(t_{\ast})$ at the beginning of the non-thermal phase and $\rho_{\mbox{\tiny gen}}(t_{\ast})=0$.

The production of entropy then becomes important when the amount of radiation generated is comparable to the initial radiation. Setting $\rho_{\mbox{\tiny init}}(t)=\rho_{\mbox{\tiny gen}}(t)$ we find this occurs at a time $t=t_{\ast}^{2/5}\Gamma^{-3/5}$ and so entropy production is important long before $t=\Gamma^{-1}$.

Prior to this, any relics that have ‘frozen-out’ (such as baryons or dark matter) will have constant co-moving number density $Y\equiv n_{x}/s$ where $n_{x}$ is their number density and $s$ the entropy density. However, once entropy production becomes important, any existing relics will be diluted until entropy production ceases near $\Gamma^{-1}$. The resulting abundance is then given by $Y_{x}\rightarrow Y_{x}\Delta^{-1}$. With $\Delta$ given by

[TABLE]

where we used $s\sim\rho_{r}^{3/4}$ and (10). If the co-moving number density instead became constant at a time $t_{f}>t_{\ast}$ then $t_{\ast}\rightarrow t_{f}$ in (11).

II.2 Entropy Production and the PBH Abundance

The previous section showed that one effect of entropy production is to dilute any existing relics. The aim of this section is to establish the effect on the abundances of PBHs.

We can account for the effect of entropy production on the mass fraction of PBHs using (11). As discussed above, for the PBHs forming during the matter phase, the initial mass fraction at the time of PBH formation $\beta(t_{f})$ will remain constant until the time of reheating – except for the effect of entropy production. Thus, using (11), for PBHs forming at $t_{f}$ the corresponding mass fraction at the time of reheating is given by

[TABLE]

where $T_{r}\sim g^{-1/4}(\Gamma m_{p})^{1/2}$ is the reheat temperature and $t_{f}$ is the time of PBH formation with $\beta(t_{f})$ the corresponding mass fraction given by (5). We have also introduced the horizon crossing time $t_{H}$. It is important to emphasize that in the matter phase, PBHs can form after Hubble radius crossing and so, in general, for a given mode $t_{f}>t_{H}$. We have left this dependence in Eq. (12) to account for this effect (in the previous literature it is commonly assumed $t_{f}=t_{H}$, but this is only an accurate approximation when $w=p/\rho\neq 0$). We have chosen fiducial values that are representative of low reheat scenarios and normalized the PBH mass to the solar mass $M_{\odot}=2.0\times 10^{33}$ g. We will justify these choices shortly.

Following reheating due to the scalar decay, the co-moving density in PBHs, $Y_{\mbox{\tiny PBH}}$, will remain fixed. We can relate this to the mass fraction at reheating

[TABLE]

or

[TABLE]

where we used $\rho\sim sT$ after reheating. We can then find the critical density in PBHs (for a given mass $M$) today

[TABLE]

where $\rho_{c}$ is the critical density. Utilizing (12) and (14) we then have

[TABLE]

where $\beta(t_{f})$ is the undiluted mass fraction given by (5), and $\rho_{c}/s_{0}=3.6\times 10^{-9}h^{2}$ GeV, with $h$ the Hubble parameter in units of $100$ km$/$s$/$Mpc (we take $h=0.7$).

II.3 Expected Mass Range

Having established the effect of entropy production on the mass fraction and abundance of PBHs, we now consider the possible ranges of masses. The minimum PBH mass is determined by the size of the Hubble volume at the beginning of the matter dominated phase. Matter phases resulting from the oscillations of a scalar field Georg et al. (2016) start with the onset of oscillations when $H_{\mbox{\tiny osc}}\simeq m_{\sigma}$. The energy density at that moment within the Hubble volume is $\rho=M_{\mathrm{min}}H_{\mbox{\tiny{osc}}}^{3}$, and we have

[TABLE]

The maximum PBH mass is given by Georg et al. (2016)

[TABLE]

where $\alpha=3.6\times 10^{-22}$, $M_{\mbox{\tiny C}}=10^{57}h^{-1}$ g ($1$ GeV = $1.8\times 10^{-24}$ g). As explained in Georg et al. (2016), this result comes from looking at the last PBHs to form at the time of scalar decay taking into account the evolution of perturbations on sub-Hubble scales as well. This is because the later in the universe a PBH is produced, the more massive it typically is, due to the increased size of the Hubble volume and the prolonged growth of the perturbations. The decay width of the scalar (moduli) is typically given by $\Gamma=c\,m_{\sigma}^{3}/m_{p}^{2}$, where $c$ is slightly model-dependent, but typically an order one number. The lighter the scalar is, the later it decays and the larger the maximum PBH mass. Therefore, we establish that the scalar mass driving the non-thermal phase determines both the minimal and maximal mass and thus the entire range of expected PBH masses.

III Non-thermal Histories and a Preferred Mass Range for PBH Dark Matter

We have seen that the scalar mass $m_{\sigma}$ controls the range of PBH masses and also the reheat temperature $T_{r}\sim(\Gamma m_{p})^{1/2}\sim m_{\sigma}^{3/2}/m_{p}^{1/2}$. Therefore, the mass fraction (12) and the PBH abundance (16) basically depend on two a priori free parameters – the mass of the scalar $m_{\sigma}$ and the spectral index $n$. In this section, we argue that the mass is in fact not a free parameter in fundamental approaches to BSM physics that predict non-thermal histories. We then use the preferred mass to find the expected relic abundance of PBHs and use existing observations to place restrictions on the spectral index.

III.1 The PBH / Dark Matter Relic Density

We have seen that PBH formation in a matter phase is vastly different than that in a thermal universe. We instead expect PBHs to form over a range of masses $M_{\mathrm{min}}\lesssim M\lesssim M_{\mathrm{max}}$ determined by the scalar mass $m_{\sigma}$. In fundamental approaches to BSM, this mass is typically set by a scale of symmetry breaking (which must be at or above the Electro-weak scale) and, as such, is connected to other aspects of the BSM theory. As a well-studied example, we briefly consider supersymmetry (SUSY), but emphasize this expectation occurs in other approaches to BSM physics.

If SUSY is broken by gravity or anomaly mediation, the mass of the scalar is tied to the gravitino mass and the scale of SUSY breaking, as $m_{\sigma}\sim m_{3/2}\sim\Lambda^{2}/m_{p}$, where $\Lambda$ is the scale of SUSY breaking. This same scale controls other aspects of the theory, such as the squark masses and the level of flavor changing neutral currents (FCNC). Thus, the scalar mass is really a consequence of SUSY breaking, and so is the duration of the non-thermal phase Kane et al. (2015). Observationally, existing bounds from the Large Hadron Collider (LHC) push squark masses to be considerably higher than the electro-weak scale, implying that the scalar (moduli) masses would also have to be higher for consistency with collider data. However, such a ‘split-spectrum’ (with TeV-scale masses) Wells (2005); Arkani-Hamed et al. (2005); Arvanitaki et al. (2013) was actually first proposed by Wells in 2004 Wells (2005). This was long before LHC constraints on SUSY, as it could improve SUSY’s phenomenological successes (e.g. improved gauge coupling unification), while also addressing challenging issues for SUSY like FCNC and CP violation. What was less appreciated at the time is that pushing the symmetry breaking scale slightly higher could also resolve the so-called Cosmological Moduli Problem Kane et al. (2015). Allowing for this modest hierarchy in scales, raising the SUSY breaking to the TeV range then predicts scalar (moduli) masses around $m_{\sigma}\simeq 10-100$ TeV in a large range of BSM models Kane et al. (2015). Although SUSY represents the best studied example, there are also string theory and extra-dimensional approaches to the hierarchy problem that make similar predictions for TeV-scale moduli. Independently, one could also simply ask the question, what was the lowest possible reheat temperature in models of inflation? Giudice et al. (2001). The answer to this question leads to a similar conclusion as the BSM considerations above. The universe must be thermalized by the time of BBN ($T_{\mathrm{r}}\simeq 5$ MeV), which corresponds to $m_{\sigma}\simeq 50$ TeV.

Given this motivation, we consider a $m_{\sigma}\simeq 50$ TeV scalar and find the corresponding prediction for the range of possible PBH masses. The lightest PBHs are given by (17) and correspond to $M_{\mbox{\tiny min}}\simeq 10^{9}$ g. The maximal mass is given by (18). The exact value depends on the spectral index of the primordial power spectrum, which is the only remaining free parameter (although it is observationally restricted as we discuss shortly). For a slightly blue spectral index ($1<n<2$), we find a maximal PBH mass near the solar mass range, or slightly above. This implies that in an early matter phase driven by scalar(s), prior to BBN, that PBHs will form in the mass range $10^{9}\;\mbox{g}\lesssim M\lesssim 1000M_{\odot}$.

Our main finding is that, given this mass range, entropy production will dilute the lighter PBHs, since they form earlier, whereas the most massive PBHs form last, near the end of reheating. That is, from (12) we see the amount of entropy dilution depends on the formation time as compared to the reheat time ($t_{f}/\Gamma^{-1})$. This suggests that the last PBHs to form will be the most abundant. Which corresponds to the maximal allowed PBH mass given by (18).

It is important to note that this maximal mass does depend on the spectral index. We find that for a blue spectrum (at pre-BBN scales) with $n=1.6$ that the maximal mass is then $24\;M_{\odot}$ – near the mass detected by LIGO and suggested by Bird, et. al. in Bird et al. (2016). However, for lower (higher) values the maximal mass will also be lower (higher). The corresponding abundance today is then given by (16), also using (5).

Our results are presented in Fig. 1, where we consider a range of values for the spectral index. We also allow for uncertainties in the initial formation probability. As mentioned above, this probability relies on estimating the evolution of perturbations in the non-linear regime, a task that is familiar from studies of large scale structure and is typically done through N-body simulations. We are not aware of such studies prior to BBN (as a thermal history is typically assumed), and instead we rely on analytic estimates and allow for a range of uncertainty, where the classic literature Khlopov (2010); Polnarev and Khlopov (1981); Khlopov and Polnarev (1980); Polnarev and Khlopov (1981) found $\beta(t_{f})\sim\delta_{M}^{13/2}$, and we instead find $\beta(t_{f})\sim\delta_{M}^{7/2}$. Regardless of this uncertainty, we see from Fig. 1 that for $n=1.6$, these PBHs will be a small fraction of the dark matter. Most optimistically ($p=2,C=10^{-1}$ in (4)), we find $\Omega_{\mbox{\tiny PBH}}=2.8\times 10^{-7}$. The PBH contribution to dark matter increases with the spectral index, but as we discuss next, this also implies stronger tension with observations due to the abundance of lighter PBHs.

III.2 Constraints on the Spectral Index

The spectral index on the scales we are considering is unknown. In this sense, we have traded the question of what PBH mass is preferred for the unknown value of the spectral index. However, we find that although entropy dilution can significantly reduce the number of PBHs with $M<M_{\mbox{\tiny max}}$, increasing the spectral index will increase the formation probability of lighter PBHs as seen from (3) and (5). Thus, we can use existing observations to place an upper bound on the spectral index. In Figs. 3 and 3 we focus on $M_{\mbox{\tiny PBH}}=10^{16}$ g PBHs, since these are within the mass range $10^{15}$ g – $10^{17}$ g most strongly restricted by CMB and $\gamma$-ray observations Carr et al. (2010). In Fig. 3 we present the formation probability for differing values of the spectral index reflecting the theoretical uncertainties. We present the same predictions in Fig. 3, but taking into account the effect of entropy dilution and we include the most recent observational constraints Clark et al. (2016). We find that even with entropy dilution and the uncertainties in the formation probability, that $n>1.6$ would be in conflict with existing observations. Instead we find that $n=1.6$ is marginally compatible with existing data given the theoretical uncertainties in the formation process. It is somewhat intriguing that this corresponds exactly to the maximal mass range corresponding to $M_{\mbox{\tiny PBH}}\simeq 30\;M_{\odot}$. However, we realize this is both speculative, and we also emphasize again that even in this case we find that PBHs would not make up all the cosmological dark matter.

IV Conclusions

In this paper we have considered the process of PBH formation in a non-thermal history. Assuming the non-thermal phase results from the oscillations of a $10-100$ TeV scalar field (which is well-motivated by BSM physics) we have a prediction that PBHs will form in a mass range $10^{9}\;\mbox{g}\lesssim M\lesssim 1000M_{\odot}$, although the maximal mass depends on the spectral index of the primordial power spectrum. If we then account for the effect of entropy production as the scalar decays, we find the most abundant PBHs will be the heaviest. Other PBHs will form, and this leads to a restriction on the spectral index of around $n<1.6$. We find that for a spectral index near $n=1.6$ that non-thermal histories predict that the most abundant PBHs will be around $M_{\mbox{\tiny PBH}}\simeq 30\;M_{\odot}$ – motivating the proposal of Bird et al. (2016) – however, we also find they would only be a small fraction of the total amount of dark matter. Even with such a small relic density, it should still be possible to further establish model building constraints using Pulsar timing as discussed recently in Orlofsky et al. (2016). We leave this investigation to future work.

Acknowledgements

We thank Simeon Bird, Steven Clark, Bhaskar Dutta, Yu Gao, Nicholas Orlofsky, and Louis Strigari for useful discussions. S.W. thanks the Michigan Center for Theoretical Physics for hospitality. This work was supported in part by NASA Astrophysics Theory Grant NNH12ZDA001N and DOE grant DE-FG02-85ER40237.

• Bird et al. (2016) S. Bird, I. Cholis, J. B. Munoz, Y. Ali-Haimoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli,  and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
• Sasaki et al. (2016) M. Sasaki, T. Suyama, T. Tanaka,  and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016), arXiv:1603.08338 [astro-ph.CO] .
• Abbott and et. al. (2016) B. P. Abbott and et. al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016).
• Orlofsky et al. (2016) N. Orlofsky, A. Pierce,  and J. D. Wells,   (2016), arXiv:1612.05279 [astro-ph.CO] .
• Carr et al. (2016) B. Carr, F. Kuhnel,  and M. Sandstad, Phys. Rev. D94, 083504 (2016), arXiv:1607.06077 [astro-ph.CO] .
• Kane et al. (2015) G. Kane, K. Sinha,  and S. Watson, Int. J. Mod. Phys. D24, 1530022 (2015), arXiv:1502.07746 [hep-th] .
• Georg et al. (2016) J. Georg, G. Sengor,  and S. Watson, Phys. Rev. D93, 123523 (2016), arXiv:1603.00023 [hep-ph] .
• Cotner and Kusenko (2016) E. Cotner and A. Kusenko,   (2016), arXiv:1612.02529 [astro-ph.CO] .
• Khlopov (2010) M. Yu. Khlopov, Res. Astron. Astrophys. 10, 495 (2010), arXiv:0801.0116 [astro-ph] .
• Polnarev and Khlopov (1981) A. G. Polnarev and M. Y. Khlopov, Soviet Journal of Astrophysics 25, 406 (1981).
• Khlopov and Polnarev (1980) M. Y. Khlopov and A. G. Polnarev, Phys. Lett. B97, 383 (1980).
• Polnarev and Khlopov (1981) A. G. Polnarev and M. Y. Khlopov, Astronomicheskii Zhurnal 58, 706 (1981).
• Harada et al. (2016) T. Harada, C.-M. Yoo, K. Kohri, K.-i. Nakao,  and S. Jhingan, Astrophys. J. 833, 61 (2016), arXiv:1609.01588 [astro-ph.CO] .
• Allahverdi et al. (2011) R. Allahverdi, B. Dutta,  and K. Sinha, Phys. Rev. D83, 083502 (2011), arXiv:1011.1286 [hep-ph] .
• Kane et al. (2011) G. Kane, J. Shao, S. Watson,  and H.-B. Yu, JCAP 1111, 012 (2011), arXiv:1108.5178 [hep-ph] .
• Moroi and Randall (2000) T. Moroi and L. Randall, Nucl.Phys. B570, 455 (2000), arXiv:hep-ph/9906527 [hep-ph] .
• Acharya et al. (2008) B. Acharya, P. Kumar, K. Bobkov, G. Kane, J. Shao,  and S. Watson, JHEP 0806, 064 (2008), arXiv:0804.0863 [hep-ph] .
• Doroshkevich (1970) A. G. Doroshkevich, Astrophysics 6, 320 (1970).
• Gorbunov and Rubakov (2011) D. S. Gorbunov and V. A. Rubakov, Introduction to the theory of the early universe (World Scientific, Hackensack, 2011).
• Padmanabhan (1993) T. Padmanabhan, Structure formation in the universe (Cambridge university press, 1993).
• Wells (2005) J. D. Wells, Phys. Rev. D71, 015013 (2005), arXiv:hep-ph/0411041 [hep-ph] .
• Arkani-Hamed et al. (2005) N. Arkani-Hamed, S. Dimopoulos, G. F. Giudice,  and A. Romanino, Nucl. Phys. B709, 3 (2005), arXiv:hep-ph/0409232 [hep-ph] .
• Arvanitaki et al. (2013) A. Arvanitaki, N. Craig, S. Dimopoulos,  and G. Villadoro, JHEP 02, 126 (2013), arXiv:1210.0555 [hep-ph] .
• Giudice et al. (2001) G. F. Giudice, E. W. Kolb,  and A. Riotto, Phys. Rev. D64, 023508 (2001), arXiv:hep-ph/0005123 [hep-ph] .
• Carr et al. (2010) B. J. Carr, K. Kohri, Y. Sendouda,  and J. Yokoyama, Phys. Rev. D81, 104019 (2010), arXiv:0912.5297 [astro-ph.CO] .
• Clark et al. (2016) S. Clark, B. Dutta, Y. Gao, L. E. Strigari,  and S. Watson,   (2016), arXiv:1612.07738 [astro-ph.CO] .