Spin Distribution of Primordial Black Holes

TL;DR

This paper estimates the spin distribution of primordial black holes, concluding they are mostly slowly rotating, based on recent studies of gravitational collapse phenomena.

Contribution

It provides a novel estimation of primordial black hole spins using recent critical phenomena research in rotating gravitational collapse.

Findings

Primordial black holes are mostly slowly rotating.

The study connects gravitational collapse physics with black hole spin distribution.

Abstract

We estimate the spin distribution of primordial black holes based on the recent study of the critical phenomena in the gravitational collapse of a rotating radiation fluid. We find that primordial black holes are mostly slowly rotating.

Equations

$$M$$

$$J$$

$${\delta }_c(q)$$

$$P(\delta )=\frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\frac{{\delta }^2}{2{\sigma }^2}\right),$$

$$\beta (M_H)={\int }_0^{\infty }Q(q)dq{\int }_{{\delta }_c(q)}^{\infty }P(\delta )d\delta \simeq \mathcal{N}{\int }_0^{\infty }dq\frac{\sigma }{\sqrt{2\pi }{\delta }_c(q)}\exp \left(-\frac{{\delta }_c(q{)}^2}{2{\sigma }^2}\right),$$

$$\beta \simeq \frac{\mathcal{N}{\sigma }^2}{2{\delta }_{c0}\sqrt{2{\delta }_{c0}K}}\exp \left(-\frac{{\delta }_{c0}^2}{2{\sigma }^2}\right).$$

$${\Omega }_{PBH}=\frac{\mathcal{N}}{M_H}{\int }_0^{\infty }dq{\int }_{{\delta }_c(q)}M(\delta )P(\delta )d\delta .$$

$$a=\frac{J}{M^2}=\frac{C_J}{C_M^2}|\delta -{\delta }_c(q){|}^{{\gamma }_J-2{\gamma }_M}q=\frac{C_J}{C_M^2}{\left(\frac{M}{C_M}\right)}^{1/2}q,$$

$$q$$

$$\delta$$

$$dqd\delta =\frac{C_M}{{\gamma }_M{C}_J}{\left(\frac{M}{C_M}\right)}^{-3/2+1/{\gamma }_M}dadM,$$

$${\Omega }_{PBH}$$

$$\frac{d{\Omega }_{PBH}}{d\ln M}\simeq \frac{\mathcal{N}}{2{\gamma }_M\sqrt{2{\delta }_{c0}K}}\frac{C_M}{M_H}{\left(\frac{M}{C_M}\right)}^{1+1/{\gamma }_M}\exp \left[-\frac{1}{2{\sigma }^2}\left({\delta }_{c0}^2+2{\delta }_{c0}{\left(\frac{M}{C_M}\right)}^{1/{\gamma }_M}\right)\right].$$

$$\frac{d{\Omega }_{PBH}}{d\ln M}$$

$$\frac{d{\Omega }_{PBH}}{da}$$

$$M_{\ast }\simeq {\left(\left(1-\frac{{\gamma }_M}{2}\right)\frac{{\sigma }^2}{{\delta }_{c0}}\right)}^{{\gamma }_M}{C}_M.$$

$$\frac{d{\Omega }_{PBH}}{da}$$

$$F_{rad}$$

$$\dot{J}\simeq -M{F}_{rad}\simeq -H^2{M}^{3}a,$$

$$\dot{M}\simeq {\rho }_{rad}{M}^2\sim H^2{M}^2.$$

$$M\simeq \frac{t}{1+\frac{t}{t_i}(t_i/M_i-1)},$$

$$\dot{a}$$

$$a\simeq a_i\exp \left(\alpha M_i/t_i\left(\frac{t_i}{t}-1\right)\right),$$

$$C_M=({\delta }_{c0}/{\sigma }^2{)}^{{\gamma }_M}(1+{\gamma }_M{)}^{-{\gamma }_M}{M}_H\simeq 5.117M_H.$$

$$M=C_0^{{\gamma }_M}(\eta -{\eta }_{\ast 0}{)}^{{\gamma }_M},$$

$$\eta =A\delta ,\text{leavevmode}\mathrm{with}\text{leavevmode}A:={\eta }_{\ast 0}/{\delta }_{c0}\simeq \mathrm{3.05513.}$$

$$M={\left(C_{0}A\right)}^{{\gamma }_M}{\left(\delta -{\delta }_{c0}\right)}^{{\gamma }_M},$$

$$J={\left(\bar{\eta }-{\bar{\eta }}_{\ast }\right)}^{{\gamma }_J}\overline{{\Omega }_{\ast }},$$

$$J={\left(\bar{\eta }-{\bar{\eta }}_{\ast }\right)}^{{\gamma }_J}\overline{{\Omega }_{\ast }}$$

$$=$$

$$q:={\bar{\Omega }}_{\ast }\times {\left(C_{0}A\right)}^{{\gamma }_M/2},$$

$$J$$

Full text

I Introduction

LIGO-VIRGO Collaboration has finally detected gravitational waves and also revealed the existence of binary black holes (BHs) in our universe Abbott:2016blz . In particular, the source of the first gravitational event, named GW150914, was reported to be a merger of a binary system of BHs with mass $30{\rm M}_{\odot}$, which is larger than that expected for a standard stellar BH. The origins of such massive BH binaries are proposed, ranging from an isolated stellar binary system to dense stellar clusters ligo2 , or even to primordial black holes (PBHs) Clesse:2016vqa ; Sasaki:2016jop ; Bird:2016dcv (see macho ; macho2 for earlier discussions).

In the near future, we expect that a large number of binary BH systems will be detected as the sources of gravitational wave events, and we can discuss the origin of binary BH systems from a statistical perspective. Among such statistical quantities, the mass function and the spin distribution of the BHs should be useful knn ; Clesse:2016vqa . The mass function of PBHs is studied in nj ; jy ; cks by applying the critical phenomena near the threshold of the BH formation choptuik ; ec ; koike . While the spin distribution of massive BH binaries formed from Population III stars has been recently studied in knn , as far as we aware, little is known about the spin distribution of PBHs.

Recently, Baumgarte and Gundlach bg ; bg2 , motivated by analytical study by gundlach , performed numerical simulations of the collapse of a rotating radiation fluid and found the critical behavior of the angular momentum. It is of immediate interest to apply these findings to study rotating PBHs. In this paper, we investigate the mass function and the spin distribution of the PBHs based on the recent study of the critical phenomena in the formation of rotating BHs.

This paper is organized as follows. In Sec. II, first we formulate the density parameter of the PBHs, $\Omega_{\rm PBH}$, and then we derive the mass function and the spin distribution of PBHs analytically for a simplified situation. We also study the evolution of the spin of a PBH after formation. Sec. III is devoted to the conclusion. We use the units $G=c=1$.

II Mass function and spin distribution

Recently, the critical phenomena in the formation of rotating BHs from radiation fluids have been investigated bg ; bg2 , and it has been found that the angular momentum $J$ also obeys the scaling relation bg and the mass and the angular momentum of the BHs depend on the two parameter family of the initial data bg2 :

[TABLE]

where $\gamma_{M}\simeq 0.3558$ koike and $\gamma_{J}=(5/2)\gamma_{M}\simeq 0.8895$. The numerical value of $\gamma_{J}$ is predicted analytically from the study of non-spherical perturbations of the critical (self-similar) solution gundlach and is confirmed numerically bg . Here $q$ denotes the parameter that characterizes the rotation of the initial data and $\delta_{c0}$ and $K$ are constants. Although the calculations are performed in asymptotically flat spacetime, we expect that similar behavior exists in the collapse of radiation fluids in the expanding universe.

For demonstrative purposes, we limit ourselves to a Gaussian probability distribution function $P(\delta)$ for density fluctuations $\delta$ :

[TABLE]

where $\sigma$ is the root-mean-square fluctuation amplitude and we assume a flat distribution function $Q(q)$ for $q$. Then the volume fraction of the region collapsing into BHs at a given epoch is given by

[TABLE]

where ${\cal N}$ is a normalization constant and we have performed the integration with respect to $\delta$ by expanding $P(\delta)$ around $\delta=\delta_{c}(q)$. The integral can be performed by expanding the integrand around $q=0$:

[TABLE]

Here, we assume $\delta_{c0}/\sigma\gg 1$. The density parameter of PBHs at a given epoch is then given by

[TABLE]

Let us introduce a dimensionless specific angular momentum (spin) parameter $a=J/M^{2}$ that corresponds to the dimensionless Kerr parameter. From Eq. (1) and Eq. (2), $a$ can be rewritten as

[TABLE]

where $\gamma_{J}=(5/2)\gamma_{M}$ is used in the last equality. $q$ and $\delta$ are now written in terms of $M$ and $a$:

[TABLE]

We can then make a change of variables from $(\delta,q)$ to $(M,a)$ in Eq. (7) to calculate the mass-spin distribution function. Since the integration measure is transformed as

[TABLE]

we write $\Omega_{PBH}$ in terms of $M$ and $a$ as

[TABLE]

II.1 Mass Function

We can obtain the mass function when we perform the $a$ integration in Eq. (12). By expanding the integrand around $a=0$, the integral can be calculated to obtain the mass function

[TABLE]

The mass function has a peak at $M_{max}=\left((1+\gamma_{M})\sigma^{2}/\delta_{c0}\right)^{\gamma_{M}}C_{M}$ and drops steeply at $2M_{max}$. Since the horizon mass $M_{H}$ is the important mass scale in the problem and a typical black hole mass is around $M_{H}$, we equate $M_{max}$ with the horizon mass $M_{H}$ jy . Then we find

[TABLE]

where we have used Eq. (5). The mass function is the same as given in jy ; nj . The mass function is shown in Fig. 1.

II.2 Spin Distribution

Next, we perform the $M$ integration in Eq. (12) to obtain the spin distribution function:

[TABLE]

The integrand has a maximum at $M=M_{*}$, where $M_{*}$ for $a<1$ is given approximately by

[TABLE]

Here, we assume $\delta_{c0}>(M_{*}/C_{M})^{1/\gamma_{M}}$. Then, by making use of the saddle point approximation around $M=M_{*}$ and further assuming $\delta/\sigma\gg 1$, which may be the case for PBH formation nj ; gl , the integral can be calculated to obtain the spin distribution function:

[TABLE]

The distribution function is given approximately by a Gaussian function of $a$.

PBHs in cosmologically interesting numbers are formed in the early universe for $\sigma/\delta_{c0}\simeq 0.1-0.2$ gl ; cksy , which corresponds to $\beta(M_{H})\simeq 5\times 10^{-24}-4\times 10^{-7}$. The value of $\delta_{c0}$ (the threshold of the density perturbation in the comoving slice) is somewhat uncertain: Carr estimated $\delta_{c0}\simeq 1/3$ carr . Recent numerical and analytical works suggest slightly larger value $\delta_{c0}\simeq 0.4\sim 0.5$ shibata ; musco1 ; harada . For definiteness, we adopt $\delta_{c0}=1/3$, $\gamma_{M}=0.3558$, and $\sigma/\delta_{c0}=0.15$, and by using the numerical values given in bg2 , which corresponds to $C_{M}\simeq 5.118M_{H},C_{J}\simeq 26.19M_{H}^{2}$ and $K\simeq 0.005685$ in our notation,111The derivation of these numerical values in our notation is given in Appendix. we show the spin distribution function in Fig. 2. We find that PBHs formed in the early universe are mostly only slowly rotating, $a<0.4$. We note that the results do not change much for $\delta_{c0}=0.4$. Note also that the distribution function is not applicable near $a\simeq 1$ because the deviations from the scaling law are large bg2 .222Ref. bg observed that Eq. (3) is valid for $\Omega<0.3$, where $\Omega$ is the parameter controlling the angular momentum, which corresponds to $a<0.7$ in our notation. We have assumed a flat distribution for $q$ in deriving Eq. (17), for simplicity. For a decreasing function $Q(q)$, $d\Omega/da$ would be much narrowed near $a=0$.

II.3 Spin Evolution

In order to connect the spin of PBHs at formation to that of the present time, we calculate the evolution of the spin of PBHs due to torque by the background radiation fluid and mass accretion.

A rotating BH sweeps the background radiation fluid and thus receives momentum from the radiation hogan ; macho2 . The force is estimated as

[TABLE]

where $\rho_{rad}$ is the energy density of the radiation. Then, the loss of angular momentum due to the torque of the force is estimated by

[TABLE]

where we have used the Friedmann equation $H^{2}\sim\rho_{rad}$. The mass of PBH grows due to the accretion of background radiation on the PBH:

[TABLE]

The solution of Eq. (20) is given by ch

[TABLE]

where $M_{i}$ is the mass at the initial time $t_{i}$. Hence, for a PBH whose mass is lighter than the horizon mass, $M_{i}<t_{i}$, the mass of PBH grows little by accretion ch .

Since $M$ remains almost constant, the spin evolution is estimated as

[TABLE]

The equation can be integrated to obtain

[TABLE]

where $\alpha$ is a constant of $O(1)$. Hence, the spin evolution is negligible for $M_{i}<t_{i}$, $a\sim a_{i}$. 333We note that spinning PBHs may suffer from super-radiant instability in the radiation dominated era because photons interacting with a plasma (free electrons) acquire an effective mass equal to the plasma frequency. Such a super-radiant instability due to cosmic plasma is effective for $M\lesssim 0.1aM_{\odot}$pani , and a PBH satisfying this relation would lose its mass and angular momentum.

III Conclusion

Motivated by the recent study of the critical behavior of angular momentum in the collapse of a rotating radiation fluid bg ; bg2 , we calculate the distribution function of the spin of PBHs. We found that most PBHs are slowly rotating, $a<0.4$. This is basically because for larger $q$ the threshold density $\delta_{c}(q)$ becomes larger and hence the probability of PBH formation is suppressed. Note that the result should depend on the distribution function for $q$. Here, we have assumed a flat distribution for $q$ for simplicity. Naively, for PBHs the parameter $q$, which is related to the initial rotational mode, can be calculated in the cosmological perturbation theory peebles , and in the standard scenario the distribution for $q$ would be expected to have a peak around $q=0$. Hence, we expect that the spin distribution of PBHs is much narrowed near $a=0$ in more realistic situation. However, it is important and should be interesting to investigate the distribution function of the initial rotational mode, and we leave it as a future study. We also estimate the spin evolution after the formation and find that it is expected to be negligible.

The evolution of Population III star binaries, which could be the sources of $30\,{\rm M}_{\odot}$ BH binaries, has been recently studied by knn , and it is found that the spins of a large fraction of the resulting BHs are high: $a\sim 1$. Therefore, we expect that the spins of binary BHs can be the probe of the origin of binary BHs.

Acknowledgments

This work is supported by MEXT KAKENHI Grant Number 15H05894 (TC), 15H05888 (SY) and 16H01103 (SY) and by JSPS KAKENHI Grant Number 15K17659 (SY) and in part by Nihon University (TC). We thank Teruaki Suyama and Jun’ichi Yokoyama for useful comments and Paolo Pani for useful information.

Appendix: Conversion of Units

In order to plot $d\,\Omega_{\rm PBH}/da$ as a function of $a$, here we show the values of input parameters ($C_{M},C_{J},K,\delta_{c0},\sigma,\gamma_{M},$ and $\gamma_{J}$) which we use in the integration with respect to $M$. Except for the values of $\gamma_{J}$, $C_{J}$ and $K$, we can use the values of input parameters in the case $q=0$. Following Refs. koike ; nj ; jy ; bg ; bg2 , we use $\delta_{c0}=1/3$, $\sigma/\delta_{c0}=0.15$ and $\gamma_{M}=0.3558$. In Ref. jy , it was found that in the case $q=0$, the mass function has a peak at $M_{\rm max}:=C_{M}(\delta_{c0}/\sigma^{2})^{-\gamma_{M}}(1+\gamma_{M})^{\gamma_{M}}$ and $M_{\rm max}$ can be identified with $M_{H}$. Following this paper, we can estimate

[TABLE]

Next, let us consider the parameters related to angular momentum, following Ref. bg2 . Eq. (20) in bg2 , given by

[TABLE]

is for the $\Omega=0$ sequence, where $\Omega$ is a control parameter related to angular momentum and hence it corresponds to $q$ in our notation. As shown in Eq. (20), $\eta_{\ast 0}=1.0183772$ and $C_{0}=0.28$. Comparing our notation with the above formula, we can find that we have a relation between $\delta$ and $\eta$

[TABLE]

Substituting this relation into (25), we have

[TABLE]

and then $C_{M}$ is given by $\left(C_{0}A\right)^{\gamma_{M}}\simeq 0.946$. However, we do not know the overall normalization scale in Ref. bg2 , and here we denote it by $M_{\rm GB}$, i.e., $C_{M}=\left(C_{0}A\right)^{\gamma_{M}}\simeq 0.946\,M_{\rm GB}$. By comparing this expression with Eq. (24), $M_{\rm GB}\simeq 5.41\,M_{H}$.

Based on the above expression, let us evaluate $C_{J}$ and $K$ in our notation. Eq. (21b) in bg2 is given by

[TABLE]

where $\bar{\eta}=C_{0}(\eta-\eta_{\ast 0})$ (from Eq. (20) in bg2 ) and $\bar{\eta}_{\ast}=K_{\rm GB}\bar{\Omega}_{\ast}^{2}$ has been used. Let us rewrite the above equation in terms of $\delta$ in our notation by using Eq. (26).

[TABLE]

As shown above, $(C_{0}\,A)^{\gamma_{M}}=0.946\,M_{\rm GB}$, and hence $\left(C_{0}\,A\right)^{\gamma_{J}}=\left(0.946\,M_{\rm GB}\right)^{\gamma_{J}/\gamma_{M}}$. By using $\gamma_{J}/\gamma_{M}=5/2$, we have $\left(C_{0}\,A\right)^{\gamma_{J}}\simeq 0.8704\,M_{\rm GB}^{5/2}$. However, the dimension of $J$ is $[{\rm mass}^{2}]$, and hence $\bar{\Omega}_{\ast}$ should have $[{\rm mass}^{-1/2}]$. Introducing a dimension less parameter $q$ as

[TABLE]

the above equation for $J$ can be written as

[TABLE]

with $K\simeq 0.005685$ and $C_{J}=\left(C_{0}\,A\right)^{2\gamma_{M}}\simeq 0.8949\,M_{\rm GB}^{2}\simeq 26.19\,M_{H}^{2}$.