A universal threshold for primordial black hole formation

TL;DR

This paper demonstrates that the threshold for primordial black hole formation is approximately universal when expressed via the compaction function, with a specific critical value supported by numerical simulations, and provides an analytic formula for this threshold.

Contribution

It introduces a universal threshold for primordial black hole formation based on the compaction function and derives an analytic expression for it, supported by numerical evidence.

Findings

Threshold value for black hole formation is approximately 0.40.

The threshold depends mainly on the curvature of the compaction function at its maximum.

Numerical simulations support the universality of the threshold.

Abstract

In this letter, we argue and show numerically that the threshold to form primordial black holes from an initial spherically symmetric perturbation is, to an excellent approximation, universal, whenever given in terms of the compaction function averaged over a sphere of radius $r_m$, where $r_m$ is the scale on which the compaction function is maximum. This can be understood as the requirement that, for a black hole to form, each shell of the averaged compaction function should have an amplitude exceeding the so-called Harada-Yoo-Kohri limit. For a radiation dominated universe we argued, supported by the numerical simulations, that this limit is $\delta_c = 0.40$, which is slightly below the one quoted in the literature. Additionally, we show that the profile dependence of the threshold for the compaction function is only sensitive to its curvature at the maximum. We use these results to provide an analytic formula for the threshold amplitude of the compaction function at its maximum in terms of the normalised compaction function curvature at $r_m$.

Equations

$$d{s}^2\simeq -d{t}^2+a(t{)}^2\left[\frac{d{r}^2}{1-K(r)r^2}+r^{2}d{\Omega }^2\right].$$

$$\mathcal{C}(r)=f(\omega )K(r)r^2,$$

$$q\equiv -\frac{{\mathcal{C}}^{\prime \prime }(r_m)r_m^2}{4\mathcal{C}(r_m)}$$

$$K_b(r)=\frac{\mathcal{C}(r_m)}{f(\omega )r_m^2}{e}^{\frac{1}{q}\left(1-{\left[\frac{r}{r_m}\right]}^{2q}\right)}.$$

$${\bar{\mathcal{C}}}_{\mathrm{c}}\equiv \frac{3}{r_m^3}{\int }_0^{r_m}{\mathcal{C}}_c(x)x^{2}dx,$$

$${\bar{\mathcal{C}}}_{\mathrm{c}}=\frac{3}{2}{e}^{\frac{1}{q}}{q}^{-1+\frac{5}{2q}}\left[\Gamma \left(\frac{5}{2q}\right)-\Gamma \left(\frac{5}{2q},\frac{1}{q}\right)\right]{\delta }_c,$$

$${\bar{\mathcal{C}}}_{\mathrm{c}}=\frac{2}{5}$$

$${\delta }_c=\frac{4}{15}{e}^{-\frac{1}{q}}\frac{q^{1-\frac{5}{2q}}}{\Gamma \left(\frac{5}{2q}\right)-\Gamma \left(\frac{5}{2q},\frac{1}{q}\right)},$$

$$K_1$$

$$K_2$$

$$K_3$$

$$g(n,k_p,r)$$

$$+i\left\{-E_{4+n}(i{k}_{p}r)+E_{4+n}(-i{k}_{p}r)\right\}],$$

$$\frac{\delta \rho }{\rho }\propto {\delta }_{\mathrm{D}}(\frac{r}{r_m}-1).$$

$$\mathcal{C}(r)\propto \frac{1}{r}{\int }_0^r\frac{\delta \rho (x)}{\rho }{x}^{2}dx$$

$$\mathcal{C}(r)=\frac{r_m}{r}\mathcal{C}(r_m)\theta (r-r_m).$$

$${\bar{\mathcal{C}}}_{\mathrm{c}}={\mathcal{C}}_c^{\mathrm{hom}}(\omega ),$$

Full text

I Introduction

In a Friedmann-Robertson-Walker (FRW) universe filled with a single fluid component having equation of state $p=\omega\rho$, a spherically symmetric local perturbation can be approximated, to leading order in gradient expansion (super-horizon scales), as111The metric is gauge dependent. Nevertheless, at leading order in gradient expansion, several gauges give the same result. To fix ideas we considered the Kodama-Sasaki gauge sasaki_gauge .

[TABLE]

Here, the local “gravitational potential” ($K(r)r^{2}$) parameterizes the initial curvature perturbation. We assume in the rest of the paper that the relevant fluctuations are all at super-horizon scales.

In yokoyama , numerical simulations were used to argue that the threshold for the amplitude of an over-density peak forming a spherically symmetric black hole in a FRW universe, only depends upon two master parameters: the integral of the initial $K(r)$ and the edge of the over-density distribution. Musco musco recently refined the arguments of yokoyama by showing that the threshold may be more conveniently given in terms of the amplitude of the gravitational potential at its maximum ($r=r_{m}$), as already noticed in shibata , and that it mainly depends upon the functional form (shape) of the gravitational potential up to $r_{m}$. More precisely, in musco , the threshold was given in terms of the “compaction function” shibata at super-horizon scales ${\cal C}(r)$ (here and after we shall simply call it compaction function)222The definition in shibata differs by a factor $2$ with respect to the one used here.: The compaction function, which closely resembles the Schwarzschild gravitational potential, is defined as twice the local excess-mass over the co-moving areal radius. At super-horizon scales it is (in units $G_{N}=1$)

[TABLE]

where $f(\omega)\equiv 3\left(1+\omega\right)/\left(5+3\omega\right)$. From this, one finds $r_{m}$ as the first root of ${\cal C}^{\prime}(r)=0$. For a radiation dominated universe, $\omega=1/3$ and so $f(\omega)\equiv 2/3$.

Regularity – the gravitational potential within a vanishingly small volume must be zero – ensures that $K(r)r^{2}\rightarrow 0$ for $r\rightarrow 0$. Thus, the behavior of $K(r)$ around the origin plays little role in black hole formation. In addition, the threshold for primordial black hole formation should be quite insensitive to the behaviour beyond $r_{m}$ as already numerically noticed in musco . The reason is simple: the threshold is the amplitude above which a “virtual” black hole of zero mass is formed. Therefore, all the over-density beyond $r_{m}$ will be diffused away while that just in the vicinity of $r_{m}$ will hinder collapse. Hence, we also expect the threshold to be very weakly dependent on the exact scale of $r_{m}$. However, above threshold, a larger part of the initial profile would be involved in the collapse. This is why the mass of the black hole with peak density above threshold also depends non-trivially on the scale $r_{m}$ on which the compaction function peaks.

Summarizing, while we expect the threshold to depend on profile shape, this dependence should come mainly from ${\cal C}$ around $r_{m}$. Since $r_{m}$ is the scale on which ${\cal C}^{\prime}=0$, we expect the threshold to depend primarily on ${\cal C}^{\prime\prime}$ at $r_{m}$. We explore this further in the next section.

II Approximating the curvature

We have checked numerically that the family of centrally peaked exponentials used in musco , for the purpose of obtaining the threshold for black hole formation, is an efficient basis with which to approximate any compaction function around its maximum.

By defining the parameter

[TABLE]

we consider the basis

[TABLE]

Note that $K_{b}(r)\to{\cal C}(r_{m})/f(\omega)r_{m}^{2}\theta(r_{m}-r)$ as $q\to\infty$. This ‘homogeneous sphere’ limit will be useful below. In contrast, $K_{b}\to{\cal C}(r_{m})/f(\omega)r^{2}$ and so $C(r)\to C(r_{m})$ as $q\to 0$. In addition, note that $r_{m}$ only defines the units of length for the scaling of $K_{b}$. Thus, as already argued before, although the mass depends on $r_{m}$, the threshold value, defined at super-horizon scales, does not and, it is constant in time333 Of course this is only strictly correct at leading order in gradient expansion. At the full non-linear level the threshold would depend upon time. At super-horizon scales, where we put the initial conditions for the gravitational collapse, however, this dependence is negligible..

We have tested our basis by considering a representative class of curvatures for the case of a radiation dominated universe ($\omega=1/3$). We have found that the threshold for black hole formation, obtained by the use of our basis, only differs by $(1\div 2)\%$ from the one obtained by the exact curvature profiles considered.

In the next section we use the basis (4) to provide an analytical formula for the thresholds. We demonstrate that our formula accurately reproduces the numerical results obtained from the publicly available code for black hole formation of albert . In turn, this will also show numerically our claim that the basis (4) well-approximates any realistic desired curvature for the calculation of the threshold.

III Universal threshold

As noticed by musco , the threshold for ${\cal C}(r_{m})$ is not universal: it depends upon the shape of the curvature profile. This implies that, if initial conditions for primordial black hole (PBH) formation are generated during inflation (see e.g. vicente ; inf ), then the threshold for PBH formation strongly depends on the form of the inflationary power spectrum germani . For example, for an almost gaussianly distributed over-density field, the mean profile around a rare peak is proportional to the two-point correlation function of the field and hence, the power spectrum bbks ; germani 444 In the estimates of PBH abundances by the use of peak theory in over-densities (e.g. germani ), one needs a threshold for the central amplitude of the un-smoothed over-density, rather than ${\cal C}(r_{m})$ (see however sg for a non-linear approach). The relation between these two quantities has been worked out in vicente and only depends on the power spectrum of curvature perturbations calculated with linear analysis.. Similarly, the parameter $q$, which will give us information about the non-linear evolution of the perturbation, only depends upon the details of the inflationary power spectrum and thus is solely related to the linear analysis of inflation.

What we show below is the remarkable fact that, nevertheless, the threshold for the average compaction function is, within $2\%$ with respect to the simulations, universal.

Let us define

[TABLE]

where ${\cal C}_{c}(r)$ is the critical compaction function for generating a black hole with zero mass. By using the basis (4), we have

[TABLE]

where $\Gamma(x)$ is the gamma function, $\Gamma(x,y)$ the incomplete gamma function and we have defined555Note that in musco $\delta_{c}$ is $\delta_{m}$. $\delta_{c}\equiv{\cal C}_{c}(\rm r_{m})$. Note that $\delta_{c}=3\frac{\delta\rho}{\rho}\frac{r_{m}^{2}}{a^{2}H^{2}}$ at super-horizon scales musco , where $\frac{\delta\rho}{\rho}$ is the over-density and $H$ is the Hubble expansion.

Following musco , if the initial perturbation is not already a black hole, the compaction function is bounded by $f(\omega)$. When radiation dominates then $\delta_{c}\to f(\omega)=2/3$ as $q\rightarrow\infty$, as shown numerically in musco . At large $q$ we have $\bar{\cal C}_{\rm c}\sim(3/5)\,\delta_{c}$ so $\bar{\cal C}_{\rm c}=2/5$. Our assumption, that we will prove both numerically and argue in the following, is that

[TABLE]

for any value of $q$. (E.g., this implies that when $q\to 0$ then $\delta_{c}\to 2/5$.) More generally, we find that the threshold for different curvature profiles $K(r)$ (when radiation dominates) is given by

[TABLE]

where $q$ is given by (3).

III.1 Numerical checks

In this section, we use the publicly available code developed in albert , to check the accuracy of (7) and (8).

First, Fig. 1 shows that (8) (solid curve) provides a good description of the threshold $\delta_{c}^{N}$ measured in the simulations (symbols) for profiles parametrized by (4). The subpanel shows the relative errors $d\equiv|\delta_{c}^{N}-\delta_{c}^{A}|/\delta_{c}^{N}$: the agreement is better than $\sim 98\%$ for all $q$.

To show that our results are more general than (4), we have also considered the following families for $K(r)$:

[TABLE]

where

[TABLE]

with $E_{n}(x)\equiv\int_{1}^{\infty}e^{-xt}\,dt/t^{n}$.

Fig. 2 illustrates that $K_{b}$ is able to provide a good approximation to (9), (10) and (11) around $r_{m}$, for a few representative parameter choices. The (oscillating) profile, $K_{3}$, is related to specific templates for inflationary power spectrum vicente , as explained in albert , and represents a non-trivial test of our claims. Indeed, this family of curves differs from the others as, generically, $k_{p}\neq r_{m}^{-1}$. Hence, $r_{m}$ (the scale on which $C^{\prime}=0$) does not define a characteristic scale for $K$. Nevertheless, as we have checked numerically, the associated threshold value still mainly depends on the behaviour of $K$ around its own $r_{m}$.

Figs. 3, 4 and 5 are similar in format to Fig. 1: they show that (7), using $q$ calibrated by fitting to $K$, provides a good description of $\delta_{c}^{N}$. The upper inner plot of Fig. 3 shows that $\delta_{c}^{N}$ reaches values that are slightly smaller than the 0.41 limit quoted in harada and later in musco , but they do not drop below the 0.4 limit of our Eq. (7).

III.2 Analytical argument

The threshold for the average compaction function ($\bar{\cal C}_{\rm c}=2/5=0.4$) is very close to the so-called Harada-Yoo-Kohri (HYK) limit which was analytically found to be $\sim 0.41$. The value of the second significant digit is related to assumptions about the Jeans length of the perturbation harada . As already mentioned, in Fig. 3, we show that this limit is actually closer to our theoretical value of $0.40$. Nevertheless, we shall still call the minimal threshold the HYK limit, as the interpretation of it will not change.

The HYK limit is the threshold for which a very sharply peaked over-density profile would collapse into a zero mass black hole, as discussed in musco . Let us then approximate an initial over-density to be a Dirac delta function

[TABLE]

One can always find an initial time where the linear approximation is good enough and find germani (see also sg )

[TABLE]

and therefore

[TABLE]

Because, as discussed above, what happens at $r>r_{m}$ is not crucial for the calculation of the threshold, we can approximate (15) as a very thin shell with finite amplitude ${\cal C}(r_{m})$ positioned at $r=r_{m}$. In other words we shall cut-off the tail in (15). The HYK limit indicates that such a shell would collapse and form a zero mass black hole if ${\cal C}(r_{m})\sim 0.4$.

Now suppose we have a continuum of concentric shells forming a homogeneous ball. This ball would then collapse to a black hole of zero mass if each shell had the same amplitude equal to the HYK threshold. Our averaging relates the problem of a generic compaction function shape to this homogeneous one.

We then conjecture that the same would happen for any $\omega$ and so, for a generic fluid matter, the threshold would be obtained for

[TABLE]

where ${\cal C}_{c}^{\rm hom}$ is the threshold for a homogeneous ball. As a first approximation, one then may be tempted to consider the functional form harada $\bar{\cal C}_{\rm c}\sim f(\omega)\sin^{2}\left[\frac{\pi\sqrt{\omega}}{1+3\omega}\right]$, however, we have numerically checked that, for $\omega\neq 1/3$, the HYK formula for ${\cal C}_{c}^{\rm hom}(\omega)$ leads to large errors.

While for radiation we could exactly fix $\bar{\cal C}_{\rm c}$ by using the limit of very peaked compaction function, we cannot do the same for other $\omega$. We then leave for future work the extensive proof of our conjecture.

IV Conclusions

Primordial black holes can account for the majority of dark matter if they are in the range of $\left[10^{-16},10^{-12}\right]\,M_{\odot}$ (see e.g. last ). The seeds for primordial black hole formation might be generated by large statistical fluctuations during inflation. The abundance of these statistical fluctuation, and in turn of the generated PBHs, is extremely sensitive to the threshold required to form a PBH sg ; germani ; vicente ; jaume . This threshold, depending upon the full non-linear evolution of the system, is typically given at super-horizon scales where the leading order in gradient expansion is an excellent approximation (see e.g. albert ).

To date, analytical estimates of it (see for example carr ; harada ) are insufficiently accurate, so numerical analyses have been employed (for the latest results see musco ; albert ).

In this paper, we have shown that although the threshold to form a PBH is initial curvature profile dependent, as noticed by musco , the threshold for the mean (i.e. volume averaged) compaction function within a sphere of radius $r=r_{m}$, is, to a very good approximation, universal and equal to the one obtained in the Harada-Yoo-Kohri limit. We used this remarkable result to provide an analytical formula for the threshold that only depends upon the normalised second derivative of the compaction function at its maximum. Specifically, for a radiation dominated universe, the threshold for a compaction function ${\cal C}(r)$ is given by Eq. 8.

Acknowledgements.

AE and CG would like to thank Vicente Atal, Jaume Garriga and Ilia Musco for many illuminating discussions. CG is supported by the Ramon y Cajal program and partially supported by the Unidad de Excelencia Maria de Maeztu Grant No. MDM-2014-0369. AE and CG are partially supported by the national FPA2016-76005-C2-2-P grants. AE is supported by the Spanish MECD fellowship FPU15/03583.