Primordial Black Holes and Co-Decaying Dark Matter

TL;DR

This paper explores how co-decaying dark matter models can lead to primordial black holes forming in the early universe, potentially accounting for a significant dark matter fraction and explaining LIGO-detected black holes.

Contribution

It derives the mass fraction and distribution of primordial black holes from co-decaying dark matter models, linking early universe physics to observable black hole phenomena.

Findings

Primordial black holes can constitute a significant fraction of dark matter.

Black hole masses can be near the solar-mass range, relevant for LIGO detections.

Models predict a non-monochromatic mass distribution of black holes.

Abstract

Models of Co-Decaying dark matter lead to an early matter dominated epoch -- prior to BBN -- which results in an enhancement of the growth of dark matter substructure. If these primordial structures collapse further they can form primordial black holes providing an additional dark matter candidate. We derive the mass fraction in these black holes (which is not monochromatic) and consider observational constraints on how much of the dark matter could be comprised in these relics. We find that in many cases they can be a significant fraction of the dark matter. Interestingly, the masses of these black holes can be near the solar-mass range providing a new mechanism for producing black holes like those recently detected by LIGO.

Figures

Equations

$$\xi \equiv \frac{s_{\text{mbox}DS}}{s_{\text{mbox}SM}}\sim \mathcal{O}(10^{-1}).$$

$${\rho }_{\text{mbox}DS}=m{s}_{\text{mbox}DS}=m\xi s_{\text{mbox}SM}=m\xi \frac{{\rho }_{\text{mbox}SM}}{T_{\text{mbox}SM}},$$

$$N_{\mbox{\tiny dur}}=\log a_{\Gamma}-\log a_{\mbox{\tiny{MD}}}=\frac{1}{3}\log\left(\frac{s_{\mbox{\tiny{SM,MD}}}}{s_{\mbox{\tiny{SM,$\Gamma$}}}}\right),$$

$$s_{\mbox{\tiny{SM,$\Gamma$}}}\simeq\frac{\Gamma_{\mbox{\tiny{B}}}^{2}m_{p}^{2}}{m\xi},\quad s_{\mbox{\tiny{SM,MD}}}\simeq(m\xi)^{3}\Rightarrow N_{\mbox{\tiny dur}}\simeq\frac{2}{3}\log\left(\frac{m_{\mbox{\tiny{A}}}^{2}\xi^{2}}{m_{p}\Gamma_{\mbox{\tiny{B}}}}\right),$$

$$f_{A}\equiv\frac{\Omega_{A}}{\Omega_{DM}}=\left(\frac{1\;$\mathrm{p}\mathrm{b}$}{\sigma}\right)\left(\frac{m_{\mbox{\tiny{A}}}}{1\;$\mathrm{G}\mathrm{e}\mathrm{V}$}\right)^{2}\left(\frac{10^{-18}\;$\mathrm{G}\mathrm{e}\mathrm{V}$}{\Gamma_{\mbox{\tiny{B}}}}\right).$$

$$\vec{r}(t,\vec{q})=a(t)\left[\vec{q}-b(t){\nabla }_q\Phi (\vec{q})\right],$$

$$\ddot{b}+2H\dot{b}=4\pi G\rho b,$$

$$\delta \vec{v}\equiv \frac{d\vec{r}}{dt}-H\vec{r}=a\frac{d\vec{x}}{dt}=-a\dot{b}{\nabla }_q\Phi ,$$

$$\rho (t,\vec{r})=\frac{⟨\rho (t)⟩}{\det \left(D\right)},$$

$$J_{ij}=\frac{\partial r_i}{\partial q_j}=a(t){\delta }_{ij}-a(t)b(t)\frac{{\partial }^2\Phi }{\partial q_i\partial q_j},$$

$$\rho (t,\vec{r})=\frac{⟨\rho (t)⟩}{\left(1+b(t)\alpha \right)\left(1+b(t)\beta \right)\left(1+b(t)\gamma \right)},$$

$$\delta (t,\vec{x})=\frac{\delta \rho }{\bar{\rho }}\simeq -\left(\alpha +\beta +\gamma \right)b(t),$$

$$w({\phi }_1\dots {\phi }_6)\text{differential}{\phi }_{1\dots 6}=\frac{|C_{pq}^{-1}{|}^{\frac{1}{2}}\exp \left[-\frac{1}{2}{\sum }_{p,q}{\phi }_p{C}_{pq}^{-1}{\phi }_q\right]}{(2\pi {)}^3}\text{differential}{\phi }_{1\dots 6}$$

$$C_{pq}\equiv \text{expectationvalue}{\phi }_p{\phi }_q=\text{expectationvalue}{D}_{ij}{D}_{kl}=\frac{{\sigma }^2}{15}({\delta }_{ij}{\delta }_{kl}+{\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}),$$

$$C=\frac{{\sigma }^2}{15}\left(\begin{array}{cc}A& 0\\ 0& \text{mathbbm}1\end{array}\right),\mathrm{with}A=\left(\begin{array}{ccc}3& 1& 1\\ 1& 3& 1\\ 1& 1& 3\end{array}\right).$$

$$w(D_{ij})\text{differential}[6]D_{ij}=\frac{675\sqrt{5}}{16{\pi }^3{\sigma }^6}\exp [N_D]\text{differential}[6]D_{ij},$$

$$N_D=-\frac{3}{{\sigma }^2}\left[(D_{11}^2+D_{22}^2+D_{33}^2)-\frac{1}{2}(D_{11}{D}_{22}+D_{11}{D}_{33}+D_{22}{D}_{33})+\frac{5}{2}(D_{12}^2+D_{13}^2+D_{23}^2)\right].$$

$$w(D_{ij})\text{differential}[6]D_{ij}=w(\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi )\left|\frac{{\partial }^6{D}_{ij}}{\partial (\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi )}\right|\text{differential}[6](\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi ).$$

$$N_D\to N_{\alpha ,\beta ,\gamma }=\frac{-3}{{\sigma }^2}\left({\alpha }^2+{\beta }^2+{\gamma }^2-\frac{1}{2}(\beta \gamma +\alpha \beta +\alpha \gamma )\right),$$

$$\left|\frac{{\partial }^6{D}_{ij}}{\partial (\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi )}\right|=(\alpha -\beta )(\alpha -\gamma )(\beta -\gamma )\sin \theta ,$$

$$\begin{array}{rl}\underset{0}{\overset{\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\underset{0}{\overset{\pi }{\int }}& w(\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi )\left|\frac{{\partial }^6{f}_{ij}}{\partial (\alpha ,\beta ,\gamma ,\theta ,\psi ,\phi )}\right|\text{differential}\theta \text{differential}\psi \text{differential}\phi \text{differential}[3](\alpha ,\beta ,\gamma )\\ & =\frac{675\sqrt{5}}{8\pi {\sigma }^6}\exp (N_{\alpha ,\beta ,\gamma })(\alpha -\beta )(\alpha -\gamma )(\beta -\gamma )\text{differential}[3](\alpha ,\beta ,\gamma ).\end{array}$$

$$\underset{-\infty }{\overset{\infty }{\int }}d\alpha \underset{\alpha }{\overset{\infty }{\int }}d\beta \underset{\beta }{\overset{\infty }{\int }}d\gamma \left(\frac{675\sqrt{5}}{8\pi {\sigma }^6}\right)\exp (N_{\alpha ,\beta ,\gamma })(\alpha -\beta )(\alpha -\gamma )(\beta -\gamma )=1.$$

$${\beta }_{\text{mbox}collapse}=\underset{-\infty }{\overset{0}{\int }}d\alpha \underset{\alpha }{\overset{\infty }{\int }}d\beta \underset{\beta }{\overset{\infty }{\int }}d\gamma \Theta \left[S(\alpha ,\beta ,\gamma )\right]W(\alpha ,\beta ,\gamma ),$$

$$S=\frac{\gamma +\beta }{2\alpha }-(1-x),$$

$$u=\frac{{\rho }_0-{\overline{\rho }}_1}{{\overline{\rho }}_1},x=\frac{r_g}{r_1},$$

$${\beta }_{\text{mbox}inhom}\sim \frac{1}{\sigma }{\int }_0^{x^{3/2}}{e}^{-(u/\sigma {)}^2}\text{differential}u\sim x^{3/2},$$

$${\beta }_{\text{mbox}inhom}\simeq 3.7x^{3/2},x\ll 1.$$

$$\beta ={\beta }_{\text{mbox}inhom}{\beta }_{\text{mbox}collapse}\simeq \left(3.7x^{3/2}\right)\times \left(5.6\times 10^{-2}{x}^5\right)\simeq 0.2{\delta }_m^{13/2},$$

$$M_{\text{mbox}min}=\frac{3m_p^2}{H_{eq}},$$

$$3H_{eq}^2{m}_p^2={\rho }_{\text{mbox}DS}+{\rho }_{\text{mbox}SM}=2{\rho }_{\text{mbox}SM}=\frac{{\pi }^2}{15}{g}_{\ast }{T}_{\text{mbox}SM}^4,$$

Full text

I Introduction

Conventional ideas for the cosmic origin and microscopic nature of dark matter (DM) are in growing tension with observations Tanabashi et al. (2018). A simple and elegant DM candidate is the thermally produced weakly interacting massive particle (WIMP). However, after years of searches utilizing direct and indirect detection techniques, and colliders, there is still no sign of the WIMP. Moreover, little is known about the evolution of the universe prior to Big Bang Nucleosynthesis (BBN). Thermal production of dark matter requires the early universe to be radiation dominated and that the DM is in thermal equilibrium with the Standard Model (SM) – but this need not be the case. Indeed, fundamental theories and top-down approaches to inflationary model building suggest that other histories are possible (and often more likely) with one example being an early matter dominated epoch (EMDE) prior to BBN Kane et al. (2015). In addition, the requirement to achieve successful inflation implies that the inflaton must couple very weakly to other sectors. In many cases, this implies that the SM and any hidden sectors would be gravitationally coupled at best. This suggests that hidden sectors may decouple very early from the SM and we need not expect them to be at the same temperature (see however Adshead et al. (2016)). Establishing the expected temperature of hidden sectors is very important for restricting model building as future experiments, such as CMB-S4, will significantly improve constraints on $N_{eff}$, and such constraints rely on a knowledge of the hidden sector temperature as compared with the SM Abazajian et al. (2016).

Co-Decaying DM (Co-Decay) is one possible alternative to the standard WIMP paradigm. As we review below, Co-Decay Dror et al. (2016, 2018); Dery et al. (2019) (and also ‘cannibalistic’ DM Pappadopulo et al. (2016); Farina et al. (2016)) posits that DM decoupled in the very early universe from the SM while relativistic. The dark sector then evolves to become non-relativistic and has a temperature differing from the SM. Another interesting property is that, upon becoming non-relativistic, the dark sector particles can dominate the energy density, leading to an EMDE until one (or more) of the particles decay to the SM ensuring a radiation dominated universe prior to BBN. As shown in Dror et al. (2018), this EMDE can lead to enhanced DM substructure due to rapid density perturbation growth Erickcek and Sigurdson (2011); Fan et al. (2014); Erickcek et al. (2016). The question we ask in this paper is: can these structures further collapse to form primordial black holes (PBHs), and could PBHs be a significant fraction of the DM?

PBHs, and whether they could be all or part of the DM, have recently received renewed interest both in model building Georg et al. (2016); Georg and Watson (2017); Cotner and Kusenko (2017a, b); Allahverdi et al. (2018); Byrnes et al. (2018); Kawasaki and Takhistov (2018); Cai et al. (2018); Kohri and Terada (2018); Cotner et al. (2018); Allahverdi et al. (2018); Gregory et al. (2017); Cotner and Kusenko (2017a); Carr et al. (2017a); Dolgov (2018); Domcke et al. (2017); Soni et al. (2017); Quintin and Brandenberger (2016); Belotsky et al. (2018), and for their observational implications Clark et al. (2017); Dalianis (2018); Clark et al. (2018); Takhistov (2019); Guo et al. (2019); Takhistov (2018); Cole and Byrnes (2018); Emami and Smoot (2018); Inoue and Kusenko (2017); Gong and Kitajima (2017); Clark et al. (2018); Lehmann et al. (2018); Dror et al. (2019). In Georg et al. (2016); Georg and Watson (2017), we considered whether PBHs could form in the presence of string moduli that lead to an EMDE, and whether they could be a significant component of the DM. An important difference of PBH formation in these models is that, unlike many approaches to PBH formation that predict a monochromatic spectrum, the extended EMDE leads to continual production of PBHs over a range of masses. For the EMDE of Georg and Watson (2017), CMB constraints on the lightest PBHs then forced the total abundance of PBH DM to be negligible.

In this paper, we consider PBH formation in the Co-Decay scenario, where the duration of the EMDE is expected to be much shorter lived than that of Georg and Watson (2017). This will weaken some of the constraints on the PBH formation process leading to the possibility of PBHs making up a significant fraction of the DM. The rest of the paper is as follows. In the next section, we review Co-Decay focusing on the aspects relevant for PBH and DM production. In Section III, we provide a somewhat detailed discussion of PBH production in an EMDE closely following the early work of Novikov (1975); Zeldovich (1970); Zeldovich and Novikov (1983); Doroshkevich (1970); Doroshkevich and Shandarin (1978); Polnarev and Khlopov (1981). In Section IV, we present our main results giving the extended mass function of PBHs from Co-Decay and placing constraints on their DM abundance using observations. We provide our conclusions in the final section.

II Co-Decaying Dark Matter

Co-Decaying Dark Matter Dror et al. (2016, 2018) posits the existence of a dark sector with at least two particles that decoupled from the SM in the very early universe. In the simplest model with only two particles, one is stable – providing a possible dark matter candidate – whereas the second can decay to the SM. For simplicity we can take these dark sector particles to be nearly mass degenerate and we label them $A$ (which plays the role of DM) and $B$ (which decays to the SM), and denote their mass as $m$.

At some early point in the cosmological evolution, the dark sector and SM fall out of equilibrium with each other. Until the $B$ particle decays with rate $\Gamma_{\mbox{\tiny{B}}}$, there will be negligible entropy transfer between the dark sector and SM. Therefore, each sector’s entropy is approximately conserved until that point. It will be useful below to consider the ratio of the entropy densities of the dark sector and SM at decoupling,

[TABLE]

After decoupling, $AA\leftrightarrow BB$ interactions111Dror et al. (2016) also considered the possibility of number changing processes between the $A$ and $B$ particles. While this does affect the relic abundance for the DM, it does not affect the duration of the matter domination and this duration is the main relevant quantity for the consideration of PBH formation, as we will see below. For this reason, we neglect the effect of ‘cannibalism’ for the rest of the paper. (with threshold s-wave annihilation rate denoted $\sigma$) maintain equilibrium between the $A$ and $B$ number density. When the temperature of the dark sector becomes order of the mass of the dark sector particles $T_{\mbox{\tiny{DS}}}\sim m$, they become non-relativistic and scale cosmologically as a pressure-less gas ($p=0$). In many cases the dark sector particles have an initial abundance that after becoming non-relativistic will naturally lead to an EMDE prior to BBN Dror et al. (2018). Once the dark sector particles become non-relativistic $T_{{\mbox{\tiny{DS}}}}<m$ we have

[TABLE]

which implies that dark sector-SM equality corresponds to $T_{{\mbox{\tiny{SM}}}}=T_{eq}=\xi m$. If an EMDE begins, it will last until the time at which the $B$ particles begin decaying to the SM, $H^{-1}>\Gamma_{\mbox{\tiny{B}}}^{-1}$. The duration of the EMDE is conveniently expressed in terms of e-folds, $N=\log a$. Assuming instantaneous decay, the duration is given by

[TABLE]

where MD denotes the onset of the EMDE, and the final equality follows from conservation of co-moving entropy222In practice, entropy production will occur continuously throughout the EMDE as discussed in Dror et al. (2018), however using the instantaneous approximation here will not change our main conclusions.. At the end of the EMDE, we have that $H\sim\Gamma_{\mbox{\tiny{B}}}$ and we find the number of e-folds

[TABLE]

where we also used $s_{\mbox{\tiny{SM}}}\sim T^{3}$. For reasonable choices of parameters, the EMDE will lead to $0\leq N_{\mbox{\tiny dur}}\leq 20$ Dror et al. (2018).

Due to the dark sector interactions, the $A$ particle number density will also decrease with the decaying particle’s number density until dark sector interactions become negligible, which follows from their velocity-averaged cross-section. To obtain the present day relic density, one must calculate the number density of the $A$ particles at this freeze-out point $H\approx n_{\mbox{\tiny{A}}}\langle\sigma v\rangle$, and then evolve it forward to the current time. The relic abundance predicted by Co-Decay is then Dror et al. (2016)

[TABLE]

Co-Decay abundance can match the observed DM relic abundance for a large portion of parameter space and is plotted in Fig. 2.

We have seen that Co-Decay can induce an EMDE. In a matter dominated universe, perturbations grow linearly with the scale factor, and the same result was shown to hold in the Co-Decay context Dror et al. (2018). In that paper, it was found that sub-structures could grow during the EMDE. Given that the dark matter decoupled from the SM far before the EMDE, this can lead to very concentrated substructures that may survive until today and lead to boosted signals for indirect detection experiments. Instead, here we want to see whether non-linear growth of these structures can lead to the formation of PBHs. In the next section, we review the formation of PBHs in a matter epoch, and then we estimate the size and amount of PBHs we expect to form.

III PBH Formation in an Early Matter Phase

The density contrast on sub-Hubble scales grows as the scale factor $\delta\equiv\delta\rho/\rho\sim a(t)$ during the EMDE driven by Co-Decay Dror et al. (2018). Once these fluctuations grow to $\delta\sim{\cal O}(1)$ linear perturbation theory fails and one must turn to other methods. In this paper, we will use the Zel’dovich Approximation Zeldovich (1970); Zeldovich and Novikov (1983) to estimate the abundance of PBHs. This method has been used by a number of authors to estimate PBH production in an EMDE Novikov (1975); Zeldovich (1970); Zeldovich and Novikov (1983); Doroshkevich (1970); Doroshkevich and Shandarin (1978); Polnarev and Khlopov (1981) – and more recently in Kokubu et al. (2018); Harada et al. (2016); Georg and Watson (2017); Georg et al. (2016); Malec and Xie (2015). The approximation has been shown to be in good agreement with full N-body simulations until the time of shell-crossing (caustic formation) Coles et al. (1993) and we leave a full N-body simulation to future work. The majority of this section can be understood by piecing together various parts of the references mentioned above, but here we want to present a self-contained, but brief review, and establish our notation for the next section.

III.1 The Zel’dovich Approximation

The Zel’dovich Approximation makes investigating the non-linear regime of density perturbations feasible by doing perturbation theory in the distance a particle moves from its initial value as opposed to requiring $\delta\ll 1$. Given that the perturbations have their initial values set by inflation, we expect the particles to be nearly homogeneous and isotropic (also in good agreement with observations). Then, as the perturbations enter the Hubble radius and eventually the non-linear regime, one can estimate the probability of PBH formation by establishing whether the collapse preserves enough sphericity to collapse to a PBH. To make this quantitative, we can write the physical location of the particles separated into their background value and the separation resulting from evolution of the perturbations,

[TABLE]

where $a(t)$ is the scale factor so that $\vec{r}=a(t)\vec{x}$ are the physical coordinates, and the co-moving coordinate $\vec{x}$ is split into its initial location $\vec{q}$ (Lagrangian coordinate) and the change resulting from the perturbations – where we again emphasize we are taking the particle separation to be small, not the density perturbation.

Initially, when a mode enters the Hubble radius we have $b(t_{H})=0$, and then the time dependence of the fluctuation leads to particle separation. This will evolve in the linear regime as

[TABLE]

which is the same as the equation for the linearized density perturbation in a pressure-less universe. Therefore, we have $b(t)\sim a(t)\sim t^{2/3}$ for the time dependent growth in the linear regime. The other term in (6) is related to the velocity perturbation as

[TABLE]

implying that the velocity perturbation can be written as a gradient, i.e. the fluid is irrotational. Using the Poisson equation (see e.g. Brandenberger (2004)) one can then show $\delta(t,\vec{r})=b(t)\nabla^{2}_{q}\Phi(q)$.

We now want to consider the non-linear regime ($\delta\sim{\cal O}(1)$). For small particle displacements in (6) we can define a linear map between $\vec{q}$ and $\vec{r}$ allowing us to relate the energy density at later times to the initial average energy density as

[TABLE]

where

[TABLE]

where $J$ is the Jacobian of the coordinate transformation and the shift in particles trajectories, $D_{ij}=J_{ij}-a(t)\delta_{ij}$, is the deformation tensor. Given the irrotational flow, the Jacobian is symmetric and diagonalizable, and we can express (9) in terms of the eigenvalues as

[TABLE]

where $\alpha,\beta,$ and $\gamma$ are the eigenvalues of the deformation tensor. It follows from (11) that in the linear regime

[TABLE]

as expected. In general, negative eigenvalues correspond to a direction collapsing. If there is more than one negative eigenvalue, the more negative one will imply the corresponding direction collapses faster. To get PHB formation, we need the eigenvalues to be roughly the same (symmetric, spherical collapse), and we will see that this requirement leads to a suppression in the formation rate since it corresponds to a less likely configuration. Also, if we consider one negative eigenvalue (say $\alpha$) in (11), at the time $b(t_{c})=-\alpha^{-1}$ the energy density is infinite. This is when a caustic forms (particle locations intersect), and then (10) can not be diagonalized. We will next consider the distribution of the eigenvalues in the deformation tensor, which will be important for addressing caustics and sufficiently spherical collapse.

III.2 Probability Distribution Function of the Eigenvalues

As discussed above, the deformation tensor, $D_{ij}$, is a symmetric matrix with six entries contained in the Jacobian (10). The distribution of its entries were calculated some time ago by Doroshkevich Doroshkevich (1970) assuming the fluctuations to be random, be Gaussian-distributed, and have zero mean.

The multivariate Gaussian distribution function of the six independent entries can be written

[TABLE]

where we have defined $\varphi_{1,2,3}=D_{ii}$ for $i=1,2,3$ (diagonal components) and $\varphi_{4,5,6}=D_{ij}=D_{ji}$ for $i\neq j$ (off-diagonal components). The covariance matrix $C_{pq}$ gives the correlation between the deformation tensor elements

[TABLE]

where $\sigma^{2}$ is the variance and $C_{pq}$ is a $6\times 6$ and block-diagonal matrix:

[TABLE]

Using this expression for $C_{pq}$ in (13) gives the probability distribution function (PDF) for the values of the deformation tensor

[TABLE]

with

[TABLE]

Again, recalling that the Jacobian (10) is diagonalizable, and therefore so is the deformation tensor, we can diagonalize it through an $SO(3)$ rotation $D=R^{T}\lambda R$ to find the six unknown elements in terms of the eigenvalues $\lambda=\mbox{diag}(\alpha,\beta,\gamma)$. The elements of the matrix $R_{ij}$ are the directional cosines of the new coordinate system with respect to the old one. The angles specifying the rotation are the Euler angles $\theta$, $\psi$, $\varphi$ (see, e.g. Goldstein et al. (2002)) and in this new coordinate system the PDF (16) becomes

[TABLE]

This expression can then be simplified by using the explicit components $D_{ij}=R\indices{{}_{i}^{k}}\lambda_{kl}R\indices{{}^{l}_{j}}$ in (17) and one finds

[TABLE]

does not depend on the Euler angles of the rotation, and the Jacobian to transform (16) is

[TABLE]

which, surprisingly, only depends on one Euler angle $\theta$.

After integration over the Euler angles the PDF becomes

[TABLE]

It is important to know the range of validity of this PDF. We want it to be positive semi-definite, so we notice that when $\gamma$ becomes greater than $\beta$, our PDF changes sign (and similarly for $\alpha$ and $\beta$). Therefore (without loss of generality), we impose that $-\infty\leq\alpha\leq\beta\leq\gamma\leq\infty$, so that our PDF will be properly normalized333We note that this PDF essentially recovers the same result of Doroshkevich (1970); Doroshkevich and Shandarin (1978), however our answer differs by a factor of five in the argument of the exponential, as Doroshkevich multiplied (14) by a factor of five to make the value of the first three diagonal components equal to the variance $\sigma^{2}$.

[TABLE]

Given this PDF we can now calculate the probability for spherical collapse, $\beta_{\mbox{\tiny collapse}}$. As discussed in Section III.1, we need at least one negative eigenvalue for collapse to occur, and for spherical collapse we need all eigenvalues negative and approximately equal. We can see the PDF is suppressed relative to a pure Gaussian distribution due to the factors $(\alpha-\beta)(\alpha-\gamma)(\beta-\gamma)$ reflecting that it is rare to have perfectly spherical collapse, and forming objects like disks (pancakes) are more likely since typically one expects the eigenvalues to differ significantly. Recalling we have ordered the eigenvalues so that $\alpha$ is the most negative (and must be negative for collapse) we have for the probability of collapse to a PBH

[TABLE]

where $W$ is the integrand of (22), $\Theta$ is the Heaviside theta function, and $S$ is a “shape function”. As mentioned above, it is not enough that $\alpha<0$ for PBH formation, and the shape function places constraints on the values of $\beta$ and $\gamma$ for a given $\alpha$. In Polnarev and Khlopov (1981), the authors supposed that a black hole could not form unless the matter clump was nearly spherical. In this case, one obtains the following shape function

[TABLE]

where $x\equiv r_{g}/r_{1}$ is the ratio of the Schwarzschild radius to the radius of the collapsing mass at the time it goes non-linear $\delta\sim{\cal O}(1)$. One can determine that this limits both $\beta$ and $\gamma$ to be greater than $\alpha(1-x)$, yet also less than $\alpha$ (i.e. the initial deformation is nearly symmetric and the eigenvalues are approximately equal). When $x\ll 1$ one can replace the argument of the exponent in the PDF with ${-9}\alpha^{2}/2$, and this recovers the result of444Performing the same integral with the full exponential argument leads to a complicated polynomial in $\arctan(x)$, but we have shown that this numeric result recovers that of Polnarev and Khlopov (1981) for $x\lesssim 0.2$. Polnarev and Khlopov (1981), namely $\beta_{\mbox{\tiny collapse}}=0.02\times x^{5}$. However, other shape functions have been considered. In particular, the authors of Harada et al. (2016) derived a different shape function. There, instead of relying on a nearly spherical collapse, they made use of the “hoop conjecture” Malec and Xie (2015); Misner et al. (1973). This conjecture states that a black hole will form from a distribution of mass if and only if the maximum circumference of the distribution is less than the Schwarzschild radius associated with its total mass. Again using the $x\ll 1$ approximation, an analytic estimate can be done and it was found $\beta_{\mbox{\tiny collapse}}=5.6\times 10^{-2}\,x^{5}$ – we will use this value in what follows.

III.3 Obtaining the Mass Fraction

To calculate the mass fraction of the universe contained within PBHs, we must consider two probabilities that determine whether an over-dense region will collapse and form a PBH. The probability for sufficient spherical collapse was found above to be $\beta_{\mbox{\tiny collapse}}=5.6\times 10^{-2}\,x^{5}$. However, we also need to consider whether or not the initial inhomogeneity is conducive to collapse (which is also equivalent to determining whether caustics prevent PBH formation) – we will call this probability $\beta_{\mbox{\tiny inhom}}$. This probability describes the difference between the background and perturbed matter densities. Using the Lemaitre-Tolman-Bondi (LTB) model Lemaitre (1997); Tolman (1934); Bondi (1947) (an exact solution of the Einstein equations) it is possible to capture this difference. Using this one can characterize the formation of two important objects for the black hole: 1) the singularity, and 2) the apparent horizon. The universe abhors naked singularities, so black hole formation can only occur when the apparent horizon forms before the singularity. In the LTB formalism, we can define quantities that describe the degree of inhomogeneity for a region of space and a characteristic size of a dark matter clump in that region when the clump becomes non-linear in its density perturbations

[TABLE]

where $\rho_{0}$ is the matter density at the center of the clump and $\overline{\rho}_{1}$ is the mean matter density inside $r_{1}$. It was shown in Polnarev and Khlopov (1981) that the formation of an apparent horizon precedes the formation of a singularity only if $u\lesssim x^{3/2}$. Assuming that the inhomogeneity is a Gaussian and random variable one finds

[TABLE]

where it is again assumed that $x\ll 1$.

This calculation has recently been refined in Kokubu et al. (2018). There it was argued that the result of Polnarev and Khlopov (1981) slightly underestimates $\beta_{\mbox{\tiny inhom}}$ due to the fact that the presence of a naked singularity can only reach the outside universe if a null signal propagating from the singularity can escape past an apparent horizon. In other words, the time between singularity and apparent horizon formation must be less than the time it takes a null vector to traverse the distance between these objects. An analytic result can be found in the limit of small initial perturbations

[TABLE]

Given this result, the total mass fraction in PBHs is then the product of the two probabilities

[TABLE]

where in the last step we have replaced $x$ by the mass density fluctuation $\delta_{m}$ at Hubble radius crossing – which follows from scaling arguments of the clump and perturbations Polnarev and Khlopov (1981). This result is about an order of magnitude larger than that used in previous PBH studies (e.g. Georg and Watson (2017); Georg et al. (2016)). When performing calculations within the context of a specific model (in this paper the model will be Co-Decaying dark matter), we use this mass fraction estimate to obtain concrete results for the mass fraction as a function of black hole mass.

IV PBHs from Co-decay

IV.1 Minimum PBH Mass

As discussed above, the EMDE begins when the energy density of the dark sector surpasses the SM radiation, corresponding to temperatures $T_{{\mbox{\tiny{SM}}}}<T_{eq}=\xi m$. A very conservative estimate of the minimum PBH mass is then found by assuming that the mass contained in the horizon volume at this time collapses to form a PBH Georg et al. (2016). As discussed above, lacking special features in the primordial power spectrum, prior to this time radiation pressure will prevent the growth of smaller PBHs. The minimal PBH mass is then given by the horizon as Georg et al. (2016)

[TABLE]

where $H_{eq}$ is the Hubble parameter when the dark sector becomes equal to the SM radiation. We emphasize that this is a very conservative estimate and most likely underestimates the mass of the PBH given that at this time half of the energy density is still in radiation and the pressure could have a substantial effect. As we will see, establishing the minimal mass is important for avoiding the strongest constraints on PBHs, and so the more massive the first PBHs, the better for avoiding the stringent bounds near $M_{\mbox{\tiny{PBH}}}\approx 10^{15}$ g. That is, increasing the mass would only lead to weaker constraints.

From the Friedmann equation at this time

[TABLE]

using $T_{{\mbox{\tiny{SM}}}}=T_{eq}=\xi m$ we find the Hubble parameter

[TABLE]

Using this result (29) leads to the minimum PBH mass

[TABLE]

where $M_{\odot}=2.0\times 10^{33}$ g is the solar-mass.

IV.2 Maximum PBH Mass

The maximum mass of PBHs formed during a matter-dominated phase is determined by those perturbations collapsing just before reheating. Since the pressure is negligible in a matter-dominated phase, we must account for sub-horizon growth – the most massive PBHs need not correspond to the horizon scale (see Georg et al. (2016); Georg and Watson (2017)). Therefore, the largest PBHs will form when $\delta_{M}(t_{r})\sim\mathcal{O}(1)$ where

[TABLE]

where we normalize the density perturbation and mass to the CMB scale, $\delta_{\mbox{\tiny{C}}}\simeq 3.8\times 10^{-6}$ and $M_{\mbox{\tiny{C}}}=10^{57}h^{-1}\;$\mathrm{g}, $h$ is the Hubble parameter in units of 100 km/s \mathrm{M}\mathrm{p}\mathrm{c}$^{-1}$, $n$ is the spectral tilt of the primordial power spectrum, and the labels $r$ and $H$ denote reheating (decay of B particles) and horizon crossing, respectively. Note that we have also used the fact that if the perturbation magnitudes are seeded by the primordial power spectrum, we can write $\delta_{M}\sim M^{(1-n)/6}$ Georg et al. (2016). Later, this identification will allow us to write the black hole formation probability as a function of mass. Since we are in an EMDE ($a(t)\propto t^{2/3}$), we can furthermore substitute expressions for the scale factors at the various times. The Hubble time at the time of decay is $H^{-1}=\Gamma_{\mbox{\tiny{B}}}^{-1}$, and this along with the relation between the horizon mass and scale factor at the time of horizon crossing of the modes resulting in the most massive PBH (which in an EMDE does not correspond to the horizon at reheating) implies that (33) is given by

[TABLE]

where $\alpha\simeq 6.9\times 10^{-23}$. Considering different values of the tilt of the power spectrum we find

[TABLE]

where our choice of the fiducial decay rate corresponds to a reheat temperature of $T_{r}=12.6$ MeV.

At this point, it is important to note that, for certain ranges of $\Gamma_{\mbox{\tiny{B}}}$ and $n$ (given $m_{\mbox{\tiny{A}}}$ and $\xi$), one could find that $M_{min}>M_{max}$! Via inspection of (32) and (IV.2), one learns that this occurs at low values of $n$ and large values of $\Gamma_{\mbox{\tiny{B}}}$. Smaller values of $n$ imply that matter perturbations require more time to become non-linear and collapse to a black hole. The duration of the matter dominated phase decreases with increasing decay rate (for a given dark matter mass). Since the perturbations require more time to grow and yet have less time to do so, there are combinations of the $n$ and $\Gamma_{\mbox{\tiny{B}}}$ such that, in order to become non-linear, a perturbation would have to enter the horizon and begin growing before the onset of the matter dominated phase. This implies that no black holes can form in a matter dominated phase characterized by these values of $n$ and $\Gamma_{\mbox{\tiny{B}}}$. Therefore, we make sure to exclude these portions of parameter space: for example, the $f_{\mbox{\tiny{PBH}}}/f_{\mbox{\tiny{A}}}$ curves are cut-off when they encounter these parameters (see Fig. 4).

IV.3 PBH Abundance

The mass fraction in PBHs during a EMDE phase is given by the probability that a PBH forms. As discussed in Section III.3, this probability for a given PBH mass $M$ is given by

[TABLE]

However, this only tells us the mass fraction at the time of their formation. In order to compare to $\Omega_{DM}$, we must calculate the present day mass fraction. The decays of the dark sector particles will induce a dilution of PBHs after they form. Including this dilution factor, the mass fraction at the time of reheating is given by Georg and Watson (2017)

[TABLE]

where $t_{f}$ is the time at which the black hole forms. To further expand this formula, we use the fact that, during a matter dominated era, we can identify

[TABLE]

where $t_{H}$ is the time at which a given perturbation mode enters the horizon and $M$ is the mass contained within a sphere of radius $H^{-1}$. Therefore, we can substitute $t_{f}=t_{H}\tau$, where $\tau=t_{f}/t_{H}>1$. The expression we obtain for $\beta(t_{r})$ represents a lower bound on the mass fraction of dark matter contained within black holes.

To determine the total portion of dark matter residing within black holes, we cannot simply integrate this mass fraction. To find the total mass fraction, we must account for the enhancement of the black hole density during the radiation dominated era that follows reheating. The black hole mass fraction will then remain approximately constant after matter-radiation equality to the present day. To this end, we calculate the function $\psi(M)$:

[TABLE]

Denoting the scale factor today as $a_{0}$ and using that $a_{eq}/a_{0}=(z_{eq}+1)^{-1}$, $a_{r}/a_{0}\simeq T_{0}/T_{r}$, $z_{eq}=3365$, and $T_{0}\simeq 230\,$\mu\mathrm{e}\mathrm{V}$$, we can determine $\psi$ for a few fiducial values of the spectral tilt:

[TABLE]

In the past, all PBH constraints were cast in terms of “monochromatic” mass functions (see, for instance, Carr et al. (2016) for a discussion of these constraints). This means that most constraints cited in various papers only apply when all black holes are characterized by the same mass. More general cosmological scenarios can lead to mass distribution functions that are extended in mass, like $\psi$ above. Therefore, some recent work has been devoted to adapting the monochromatic constraints to extended mass functions Carr et al. (2017b); Bellomo et al. (2018). This usually means that the constraints become more restrictive. However, a recent work Lehmann et al. (2018) has shown that a maximum allowed PBH mass fraction can be derived independently of a given mass distribution shape. When we restrict $\psi(M)$ for various parameters of the Co-Decay model, we can thus determine the models viability in two ways: (1) find whether $\psi(M)$ falls into an area excluded by monochromatic constraints, and (2) determining whether the total mass fraction (see below) is larger than that allowed by analysis similar to Lehmann et al. (2018). If either of these scenarios occur, Co-Decay with those parameters is excluded (see Fig. 3 for comparisons between the constraints and $\psi$ for various parameters).

The quantity $\psi(M)dM$ represents the PBH dark matter fraction in the interval $(M,M+dM)$, which means that Carr et al. (2017c)

[TABLE]

Upon integration, we find that we can restrict the total fraction in terms of the $\chi$ parameter defined above (see Fig. 4). However, in this case, we must still pick a value for $m_{\mbox{\tiny{A}}}$, which in turn implies that $\chi$ here plays the role of $\Gamma_{\mbox{\tiny{B}}}$.

Our constraints on the Co-Decay PBHs appear in Figure 4. We see from Fig. 4 that as $n_{s}$ increases, so does the portion of dark matter residing in black holes. This makes sense because the larger the spectral tilt, the larger the density contrast a given mass scale has upon horizon entry. The difference between the solid ($\sigma=10^{6}$ pb) and dashed lines ($\sigma=10^{9}$ pb) resides in the fact that for a given mass and decay rate a higher annihilation threshold implies a smaller A particle abundance (see also Fig. 2). The black solid and dashed lines represent the maximum ratio as derived from Lehmann et al. (2018). We use the most restrictive value for the maximum allowed total, which is determined from the constraint set ĀBC from Lehmann et al. (2018).

There are two interesting features in the ratio dependence of $\chi$. The most obvious of these is the sharp downturn to zero at certain values of low $\chi$ (large $\Gamma_{\mbox{\tiny{B}}}$). For each value of $n_{s}$, this occurs at the value of the decay rate for which the PBH maximum mass approaches the minimum mass. As discussed above, this implies that there is no black hole formation throughout the matter dominated era. The more surprising of the features is the gradual decrease in mass fraction for increasing $\chi$ (decreasing $\Gamma_{\mbox{\tiny{B}}}$). Decreasing decay rate implies both a lengthening of the matter domination and an increase in the maximum mass for the black holes. This suggests that for smaller decay rates, there would be more time to form black holes and a larger range of masses that can be collapsed to black holes.

Two factors contribute to the decreasing nature of the total mass fraction. As can be seen from (5), the A particle abundance is inversely proportional to the decay rate. As the EMDE lengthens, the amount of relic A particles increase as well, which would cause a decrease in the ratio of the mass fractions. Secondly, the mass distribution function is determined by the probability of collapse at the time of reheating. Due to the microphysics of Co-Decay, the black holes formed early in the matter dominated phase are diluted due to the entropy transfers between the dark and visible sectors. For this reason, the longer the matter domination, the more dilution of the relics that occurs.

We would now like to consider whether there are portions of the allowed parameter space in which PBHs constitute an appreciable fraction of the universe’s DM budget. To this end, we recall from Fig. 2 that $f_{A}=1$ at $\chi=10^{24}$ GeV for $\sigma=10^{6}$ pb, and at $\chi=10^{27}$ GeV for $\sigma=10^{9}$ pb. For $n_{s}<2$, PBHs never make up more than a negligible portion of DM. We can see this from the fact that for the larger cross-section values, the maximum fraction of DM in black holes is $10^{4}f_{A}$, and occurs at $\chi\sim 10^{21}$. The value of $f_{A}$ at this point is $\sim 10^{-8}$. At a spectral tilt of 1.8, Co-Decay can only account for one ten-thousandth of all DM.

When we push the spectral tilt up to 2 and have a dark sector annihilation rate of $10^{9}$ pb, we see that if $\chi\in(10^{24},10^{27})$ GeV, PBHs can account for all of DM, while still avoiding the constraints in Fig. 3: $f_{\mbox{\tiny{PBH}}}=10^{4}f_{A}=1$ at $\chi=10^{24}$ GeV. When $\chi\sim\mathcal{O}(10^{25})$ GeV, the A particles and PBHs have equal abundances, but their total is only $\mathcal{O}(.01)$. Moving to higher $\chi$ (lower $\Gamma_{\mbox{\tiny{B}}}$) brings us closer to the point where $f_{A}=1$. When that is the case, Fig. 4 still shows that PBHs can account for $\sim\mathcal{O}(.01)-\mathcal{O}(.1)$ of DM.

The above plot and reasoning leads us to suspect that, although it is unlikely that PBHs are all of DM, early-forming, solar-mass black holes can be non-negligible fractions of the total DM budget for appreciable ranges of parameter space. Thus, Co-Decay DM offers an interesting explanation for the presence of an abundance of near solar-mass black holes.

V Conclusions

In this paper, we have considered the production of PBHs in models of Co-Decay. We found the mass fraction and relic abundance for varying values of the primordial power spectrum. We have seen that observational constraints place important restrictions on the parameter space for Co-Decay and the amount of PBHs that can comprise DM. We find that given the range of masses that are produced that the mass range near or around a solar-mass is the most promising. Depending on the further evolution of these PBHs, this could provide a new mechanism for accounting for near solar-mass black holes, like those detected by LIGO. An important challenge is then to see if these PBHs can evolve to form binary pairs and determine their merger rate. We leave this to future work.

Acknowledgements

We thank Jeff Dror, Bhaskar Dutta, Adrienne Erickcek, Benjamin Lehmann, Eric Kufick, Alex Kusenko, and Stefano Profumo for useful discussions. This research was supported in part by NASA Astrophysics Theory Grant NNH12ZDA001N and DOE grant DE-FG02-85ER40237.

• Tanabashi et al. (2018) M. Tanabashi et al. (Particle Data Group), Phys. Rev. D98, 030001 (2018).
• Kane et al. (2015) G. Kane, K. Sinha,  and S. Watson, Int. J. Mod. Phys. D24, 1530022 (2015), arXiv:1502.07746 [hep-th] .
• Adshead et al. (2016) P. Adshead, Y. Cui,  and J. Shelton, JHEP 06, 016 (2016), arXiv:1604.02458 [hep-ph] .
• Abazajian et al. (2016) K. N. Abazajian et al. (CMB-S4),   (2016), arXiv:1610.02743 [astro-ph.CO] .
• Dror et al. (2016) J. A. Dror, E. Kuflik,  and W. H. Ng, Phys. Rev. Lett. 117, 211801 (2016), arXiv:1607.03110 [hep-ph] .
• Dror et al. (2018) J. A. Dror, E. Kuflik, B. Melcher,  and S. Watson, Phys. Rev. D97, 063524 (2018), arXiv:1711.04773 [hep-ph] .
• Dery et al. (2019) A. Dery, J. A. Dror, L. Stephenson Haskins, Y. Hochberg,  and E. Kuflik,   (2019), arXiv:1901.02018 [hep-ph] .
• Pappadopulo et al. (2016) D. Pappadopulo, J. T. Ruderman,  and G. Trevisan, Phys. Rev. D94, 035005 (2016), arXiv:1602.04219 [hep-ph] .
• Farina et al. (2016) M. Farina, D. Pappadopulo, J. T. Ruderman,  and G. Trevisan, JHEP 12, 039 (2016), arXiv:1607.03108 [hep-ph] .
• Erickcek and Sigurdson (2011) A. L. Erickcek and K. Sigurdson, Phys. Rev. D84, 083503 (2011), arXiv:1106.0536 [astro-ph.CO] .
• Fan et al. (2014) J. Fan, O. Ozsoy,  and S. Watson, Phys. Rev. D90, 043536 (2014), arXiv:1405.7373 [hep-ph] .
• Erickcek et al. (2016) A. L. Erickcek, K. Sinha,  and S. Watson, Phys. Rev. D94, 063502 (2016), arXiv:1510.04291 [hep-ph] .
• Georg et al. (2016) J. Georg, G. Sengor,  and S. Watson, Phys. Rev. D93, 123523 (2016), arXiv:1603.00023 [hep-ph] .
• Georg and Watson (2017) J. Georg and S. Watson, JHEP 09, 138 (2017), [JHEP09,138(2017)], arXiv:1703.04825 [astro-ph.CO] .
• Cotner and Kusenko (2017a) E. Cotner and A. Kusenko, Phys. Rev. D96, 103002 (2017a), arXiv:1706.09003 [astro-ph.CO] .
• Cotner and Kusenko (2017b) E. Cotner and A. Kusenko, Phys. Rev. Lett. 119, 031103 (2017b), arXiv:1612.02529 [astro-ph.CO] .
• Allahverdi et al. (2018) R. Allahverdi, J. Dent,  and J. Osinski, Phys. Rev. D97, 055013 (2018), arXiv:1711.10511 [astro-ph.CO] .
• Byrnes et al. (2018) C. T. Byrnes, P. S. Cole,  and S. P. Patil,   (2018), arXiv:1811.11158 [astro-ph.CO] .
• Kawasaki and Takhistov (2018) M. Kawasaki and V. Takhistov, Phys. Rev. D98, 123514 (2018), arXiv:1810.02547 [hep-th] .
• Cai et al. (2018) Y.-F. Cai, X. Tong, D.-G. Wang,  and S.-F. Yan, Phys. Rev. Lett. 121, 081306 (2018), arXiv:1805.03639 [astro-ph.CO] .
• Kohri and Terada (2018) K. Kohri and T. Terada, Class. Quant. Grav. 35, 235017 (2018), arXiv:1802.06785 [astro-ph.CO] .
• Cotner et al. (2018) E. Cotner, A. Kusenko,  and V. Takhistov, Phys. Rev. D98, 083513 (2018), arXiv:1801.03321 [astro-ph.CO] .
• Gregory et al. (2017) R. Gregory, D. Kastor,  and J. Traschen, JHEP 10, 118 (2017), arXiv:1707.06586 [hep-th] .
• Carr et al. (2017a) B. Carr, T. Tenkanen,  and V. Vaskonen, Phys. Rev. D96, 063507 (2017a), arXiv:1706.03746 [astro-ph.CO] .
• Dolgov (2018) A. D. Dolgov, Usp. Fiz. Nauk 188, 121 (2018), [Phys. Usp.61,no.2,115(2018)], arXiv:1705.06859 [astro-ph.CO] .
• Domcke et al. (2017) V. Domcke, F. Muia, M. Pieroni,  and L. T. Witkowski, JCAP 1707, 048 (2017), arXiv:1704.03464 [astro-ph.CO] .
• Soni et al. (2017) A. Soni, H. Xiao,  and Y. Zhang, Phys. Rev. D96, 083514 (2017), arXiv:1704.02347 [hep-ph] .
• Quintin and Brandenberger (2016) J. Quintin and R. H. Brandenberger, JCAP 1611, 029 (2016), arXiv:1609.02556 [astro-ph.CO] .
• Belotsky et al. (2018) K. M. Belotsky, V. I. Dokuchaev, Y. N. Eroshenko, E. A. Esipova, M. Yu. Khlopov, L. A. Khromykh, A. A. Kirillov, V. V. Nikulin, S. G. Rubin,  and I. V. Svadkovsky,   (2018), arXiv:1807.06590 [astro-ph.CO] .
• Clark et al. (2017) S. Clark, B. Dutta, Y. Gao, L. E. Strigari,  and S. Watson, Phys. Rev. D95, 083006 (2017), arXiv:1612.07738 [astro-ph.CO] .
• Dalianis (2018) I. Dalianis,   (2018), arXiv:1812.09807 [astro-ph.CO] .
• Clark et al. (2018) S. Clark, B. Dutta, Y. Gao, Y.-Z. Ma,  and L. E. Strigari, Phys. Rev. D98, 043006 (2018), arXiv:1803.09390 [astro-ph.HE] .
• Takhistov (2019) V. Takhistov, Phys. Lett. B789, 538 (2019), arXiv:1710.09458 [astro-ph.HE] .
• Guo et al. (2019) H.-K. Guo, J. Shu,  and Y. Zhao, Phys. Rev. D99, 023001 (2019), arXiv:1709.03500 [astro-ph.CO] .
• Takhistov (2018) V. Takhistov, Phys. Lett. B782, 77 (2018), arXiv:1707.05849 [astro-ph.CO] .
• Cole and Byrnes (2018) P. S. Cole and C. T. Byrnes, JCAP 1802, 019 (2018), arXiv:1706.10288 [astro-ph.CO] .
• Emami and Smoot (2018) R. Emami and G. Smoot, JCAP 1801, 007 (2018), arXiv:1705.09924 [astro-ph.CO] .
• Inoue and Kusenko (2017) Y. Inoue and A. Kusenko, JCAP 1710, 034 (2017), arXiv:1705.00791 [astro-ph.CO] .
• Gong and Kitajima (2017) J.-O. Gong and N. Kitajima, JCAP 1708, 017 (2017), arXiv:1704.04132 [astro-ph.CO] .
• Lehmann et al. (2018) B. V. Lehmann, S. Profumo,  and J. Yant, JCAP 1804, 007 (2018), arXiv:1801.00808 [astro-ph.CO] .
• Dror et al. (2019) J. A. Dror, H. Ramani, T. Trickle,  and K. M. Zurek,   (2019), arXiv:1901.04490 [astro-ph.CO] .
• Novikov (1975) I. D. Novikov, Astronomicheskii Zhurnal 52, 1038 (1975).
• Zeldovich (1970) Ya. B. Zeldovich, Astron. Astrophys. 5, 84 (1970).
• Zeldovich and Novikov (1983) I. B. Zeldovich and I. D. Novikov, Relativistic astrophysics. Volume 2. (Chicago, IL, University of Chicago Press, 1983, 751 p. Translation., 1983).
• Doroshkevich (1970) A. G. Doroshkevich, Astrofizika 6, 581 (1970).
• Doroshkevich and Shandarin (1978) A. G. Doroshkevich and S. F. Shandarin, Astronomicheskii Zhurnal 22, 653 (1978).
• Polnarev and Khlopov (1981) A. G. Polnarev and M. Y. Khlopov, Astronomicheskii Zhurnal 25, 406 (1981).
• Kokubu et al. (2018) T. Kokubu, K. Kyutoku, K. Kohri,  and T. Harada,  (2018), arXiv:1810.03490 [astro-ph.CO] .
• 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] .
• Malec and Xie (2015) E. Malec and N. Xie, Phys. Rev. D91, 081501 (2015), arXiv:1503.01354 [gr-qc] .
• Coles et al. (1993) P. Coles, A. L. Melott,  and S. F. Shandarin, Monthly Notices of the Royal Astronomical Society 260, 765 (1993).
• Brandenberger (2004) R. H. Brandenberger, Lect. Notes Phys. 646, 127 (2004), arXiv:hep-th/0306071 [hep-th] .
• Goldstein et al. (2002) H. Goldstein, C. Poole,  and J. Safko, Classical Mechanics (Addison Wesley, 2002).
• Misner et al. (1973) C. W. Misner, K. S. Thorne,  and J. A. Wheeler, Gravitation (W. H. Freeman, San Francisco, 1973).
• Lemaitre (1997) G. Lemaitre, Gen. Rel. Grav. 29, 641 (1997), [Annales Soc. Sci. Bruxelles A53,51(1933)].
• Tolman (1934) R. C. Tolman, Proc. Nat. Acad. Sci. 20, 169 (1934), [Gen. Rel. Grav.29,935(1997)].
• Bondi (1947) H. Bondi, Mon. Not. Roy. Astron. Soc. 107, 410 (1947).
• 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] .
• Barnacka et al. (2012) A. Barnacka, J. F. Glicenstein,  and R. Moderski, Phys. Rev. D86, 043001 (2012), arXiv:1204.2056 [astro-ph.CO] .
• Graham et al. (2015) P. W. Graham, S. Rajendran,  and J. Varela, Phys. Rev. D92, 063007 (2015), arXiv:1505.04444 [hep-ph] .
• Niikura et al. (2017) H. Niikura et al.,   (2017), arXiv:1701.02151 [astro-ph.CO] .
• Griest et al. (2014) K. Griest, A. M. Cieplak,  and M. J. Lehner, Astrophys. J. 786, 158 (2014), arXiv:1307.5798 [astro-ph.CO] .
• Tisserand et al. (2007) P. Tisserand et al. (EROS-2), Astron. Astrophys. 469, 387 (2007), arXiv:astro-ph/0607207 [astro-ph] .
• Zumalacarregui and Seljak (2018) M. Zumalacarregui and U. Seljak, Phys. Rev. Lett. 121, 141101 (2018), arXiv:1712.02240 [astro-ph.CO] .
• Allsman et al. (2001) R. A. Allsman et al. (Macho), Astrophys. J. 550, L169 (2001), arXiv:astro-ph/0011506 [astro-ph] .
• Koushiappas and Loeb (2017) S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. 119, 041102 (2017), arXiv:1704.01668 [astro-ph.GA] .
• Brandt (2016) T. D. Brandt, Astrophys. J. 824, L31 (2016), arXiv:1605.03665 [astro-ph.GA] .
• Monroy-Rodriguez and Allen (2014) M. A. Monroy-Rodriguez and C. Allen, Astrophys. J. 790, 159 (2014), arXiv:1406.5169 [astro-ph.GA] .
• Ali-Haimoud and Kamionkowski (2017) Y. Ali-Haimoud and M. Kamionkowski, Phys. Rev. D95, 043534 (2017), arXiv:1612.05644 [astro-ph.CO] .
• Carr et al. (2017b) B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen,  and H. Veermae, Phys. Rev. D96, 023514 (2017b), arXiv:1705.05567 [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] .
• Bellomo et al. (2018) N. Bellomo, J. L. Bernal, A. Raccanelli,  and L. Verde, JCAP 1801, 004 (2018), arXiv:1709.07467 [astro-ph.CO] .
• Carr et al. (2017c) B. Carr, T. Tenkanen,  and V. Vaskonen, Phys. Rev. D96, 063507 (2017c), arXiv:1706.03746 [astro-ph.CO] .