The Initial Spin Probability Distribution of Primordial Black Holes

TL;DR

This paper models the initial spin distribution of primordial black holes formed from early universe inhomogeneities, revealing that spins are typically at percent levels due to first-order tidal effects and shape asymmetries.

Contribution

It provides a theoretical framework using peak theory to predict the probability distribution of primordial black hole spins at formation, incorporating tidal torques and asphericity effects.

Findings

Typical spin parameter $a_s$ is about a few percent.

Spin distribution depends on the variance of overdensity and power spectrum width.

Narrow spectra suppress spin due to velocity shear alignment.

Abstract

We study the spin of primordial black holes produced by the collapse of large inhomogeneities in the early universe. Since such primordial black holes originate from peaks, that is, from maxima of the local overdensity, we resort to peak theory to obtain the probability distribution of the spin at formation. We show that the spin is a first-order effect in perturbation theory: it results from the action of first-order tidal gravitational fields generating first-order torques upon horizon-crossing, and from the asphericity of the collapsing object. Assuming an ellipsoidal shape, the typical value of the dimensionless parameter $a_{\rm s}=S/G_N M^2$, where $S$ is the spin and $M$ is the mass of the primordial black hole, is about $\sigma_\delta\sqrt{1-\gamma^2}/2\pi$. Here, $\sigma^2_\delta$ is the variance of the overdensity at horizon crossing and the parameter $\gamma$ is a measure of the width of the power spectrum giving rise to primordial black holes. One has $\gamma=1$ for monochromatic spectra. For these narrow spectra, the suppression arises because the velocity shear, which is strongly correlated with the inertia tensor, tends to align with the principal axis frame of the collapsing object. Typical values of $a_{\rm s}$ are at the percent level.

Equations

$${\chi }_{\mathrm{eff}}=\frac{M_1{\vec{S}}_1+M_2{\vec{S}}_2}{M_1+M_2}\cdot \hat{L}.$$

$$a_f=\frac{S_f}{G_N{M}_f^2}\simeq 0.69-0.56{\left(\frac{M_1-M_2}{M_1+M_2}\right)}^2,$$

$$a_{\mathrm{s}}=\frac{1}{2\pi }{\sigma }_{\delta }\sqrt{1-{\gamma }^2}\sim 10^{-2}\sqrt{1-{\gamma }^2},$$

$$S(\Sigma )=\frac{1}{16\pi G_N}{\int }_{\partial \Sigma }\mathrm{d}{S}_{\mu \nu }{D}^{\mu }{\varphi }^{\nu },$$

$$S(\Sigma )={\int }_{\Sigma }\mathrm{d}{S}_{\mu }{J}^{\mu }(\varphi ),$$

$$S(\Sigma )={\int }_{\Sigma }\mathrm{d}V{T}_{\mu }^0{\varphi }^{\mu }=\frac{4}{3}{\int }_{\Sigma }\mathrm{d}V\rho \vec{v}\cdot \vec{\varphi },$$

$$S_i=\frac{4}{3}{a}^4(\eta ){ϵ}_{ijk}\int {\mathrm{d}}^{3}x\sqrt{\gamma }\rho (\vec{x},\eta )(x-x_{\text{pk}}{)}^j(v-v_{\text{pk}}{)}^k,$$

$$\delta (\vec{x})=\frac{\delta \rho (\vec{x})}{\bar{\rho }}\simeq {\delta }_{\text{pk}}+\frac{1}{2}{\zeta }_{ij}(x-x_{\text{pk}}{)}^i(x-x_{\text{pk}}{)}^j>f{\delta }_{\text{pk}},{\zeta }_{ij}={\frac{{\partial }^2\delta }{\partial x^i\partial x^j}|}_{\text{pk}},$$

$$\delta \simeq {\delta }_{\text{pk}}-\frac{1}{2}{\sigma }_{\zeta }\sum _1^3{\lambda }_i(x^i-x_{\text{pk}}^i{)}^2.$$

$$\frac{2(\delta -{\delta }_{\text{pk}})}{{\sigma }_{\zeta }}=-\sum _1^3{\lambda }_i(x^i-x_{\text{pk}}^i{)}^2.$$

$$2\frac{{\sigma }_{\delta }}{{\sigma }_{\zeta }}(1-f)\nu =\sum _1^3{\lambda }_i(x^i-x_{\text{pk}}^i{)}^2,$$

$$a_i^2=2\frac{{\sigma }_{\delta }}{{\sigma }_{\zeta }}\frac{(1-f)}{{\lambda }_i}\nu .$$

$${\lambda }_i=\frac{\gamma \nu }{3}\left(1+{ϵ}_i\right),{ϵ}_i=\mathcal{O}\left(\frac{1}{\gamma \nu }\right)\mathrm{and}{a}_i^2\sim 6\frac{{\sigma }_{\delta }^2}{{\sigma }_{\times }^2},$$

$$\gamma =\frac{{\sigma }_{\times }^2}{{\sigma }_{\delta }{\sigma }_{\zeta }}$$

$$\displaystyle(v-v_{\text{\tiny pk}})^{k}=v_{l}^{k}(x-x_{\text{\tiny pk}})^{l},\qquad v_{l}^{k}=\frac{\partial v^{k}}{\partial x^{l}}\bigg{|}_{\text{\tiny pk}},$$

$$S_i$$

$$S_i=\frac{4}{3}{a}^4(\eta ){ϵ}_{ijk}{\overline{\rho }}_{\mathrm{rad}}(\eta )v_l^k{\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x{x}^j{x}^l=\frac{4}{3}{a}^4(\eta ){ϵ}_{ijk}{\overline{\rho }}_{\mathrm{rad}}(\eta )g_v(\eta ){\tilde{v}}_{kl}{\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x{x}^j{x}^l,$$

$$v_l^k(\eta )=g_v(\eta ){\tilde{v}}_l^k.$$

$${\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x=a_1{a}_2{a}_3{\int }_0^1{r}^2\mathrm{d}r{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\pi }\mathrm{d}\theta \sin \theta$$

$$x_1=a_{1}r\cos \varphi \sin \theta ,x_2=a_{2}r\sin \varphi \sin \theta ,x_3=a_{3}r\cos \theta .$$

$${\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x{x}^j{x}^l=\frac{4\pi }{15}{a}_1{a}_2{a}_3{\left[\begin{array}{ccc}{a}_1^2& 0& 0\\ 0& a_2^2& 0\\ 0& 0& a_3^2\end{array}\right]}_{jl}.$$

$${\vec{S}}^{(1)}=\frac{16\pi }{45}{a}^4(\eta ){\overline{\rho }}_{\mathrm{rad}}(\eta )g_v(\eta )a_1{a}_2{a}_3([a_2^2-a_3^2]{\tilde{v}}_{23},[a_3^2-a_1^2]{\tilde{v}}_{13},[a_1^2-a_2^2]{\tilde{v}}_{12}).$$

$${\vec{S}}^{(1)}$$

$${\alpha }_1=\frac{1}{{\lambda }_3}-\frac{1}{{\lambda }_2},{\alpha }_2=\frac{1}{{\lambda }_3}-\frac{1}{{\lambda }_1},{\alpha }_3=\frac{1}{{\lambda }_2}-\frac{1}{{\lambda }_1},$$

$$\Lambda ={\lambda }_1{\lambda }_2{\lambda }_3,R_{\ast }=\sqrt{3}\frac{{\sigma }_{\times }}{{\sigma }_{\zeta }}.$$

$$S_{\mathrm{ref}}(\eta )=\frac{4}{3}{a}^4(\eta ){\overline{\rho }}_{\mathrm{rad}}(\eta )g_v(\eta )R_{\ast }^5(1-f{)}^{5/2}.$$

$$S^{(1)}(\eta )=S_{\mathrm{ref}}(\eta )s_{\mathrm{e}}^{(1)}.$$

$$\rho (\vec{x},\eta )={\overline{\rho }}_{\mathrm{rad}}(\eta )+\delta \rho ={\overline{\rho }}_{\mathrm{rad}}(\eta )\left(1+\delta \right)={\overline{\rho }}_{\mathrm{rad}}(\eta )\left[1+{\delta }_{\text{pk}}+\frac{1}{2}{\zeta }_{ij}(x-x_{\text{pk}}{)}^i(x-x_{\text{pk}}{)}^j\right]$$

$$={\overline{\rho }}_{\mathrm{rad}}(\eta )\left(1+{\delta }_{\text{pk}}\right)+{\overline{\rho }}_{\mathrm{rad}}(\eta )\left[\frac{1}{2}{\zeta }_{ij}(x-x_{\text{pk}}{)}^i(x-x_{\text{pk}}{)}^j\right].$$

$$S_i=\frac{4}{3}a(\eta {)}^4{ϵ}_{ijk}{\overline{\rho }}_{\mathrm{rad}}(\eta )v_l^k\left[\left(1+{\delta }_{\text{pk}}\right){\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x{x}^j{x}^l+\frac{1}{2}{\zeta }_{mn}{\int }_{V_{\mathrm{e}}}{\mathrm{d}}^{3}x{x}^j{x}^l{x}^m{x}^n\right],$$

$$I^{1111}$$

Full text

1 Introduction

The physics of black holes has attracted a lot of attention since the recent discovery of the gravitational waves generated by the merger of rather massive black holes [1]. This interest is also motivated by the fact that PBHs might compose a fraction of the dark matter of the universe [2] (for a review, see Ref. [3]).

One particularly interesting property of the black holes is their spin. Larger integrated fluxes of gravitational waves are obtained for larger spins as they lead to smaller ISCO separations in binaries and to longer inspiral phases. The observations of gravitational waves sourced emitted during the merging of two massive black holes may be exploited to measure the spin of the final state, as well as the orbital projection of the effective spin of the black hole binary

[TABLE]

Here, $\vec{S}_{1}$ and $\vec{S}_{2}$ are the spins of the initial merging black holes in the binary and $\vec{L}$ their initial angular momentum. Currently, the handful number of gravitational wave measurements indicates that $\chi_{\rm eff}$ is compatible with zero [4].

There are currently two major astrophysical models of black hole spin distributions. One possibility is that the final spin of the binary system and the angular momentum are aligned, as expected for black hole binaries formed in a shared envelope evolution within galactic fields [5]. Alternatively, the initial spins have an isotropic distribution as expected in binaries originated in globular or stellar clusters in the proximity of active galactic nuclei [6]. There can be another mechanism which produces small effective spins. This is likely if the black holes are of primordial origin111Massive black holes of stellar origin may also be born with small natal spins [7]., hence the name Primordial Black Holes (PBHs), and they are produced by the collapse of sizeable inhomogeneities in the early universe. This case is particularly interesting because numerical relativity simulations provide evidence that two black holes with initially small spin merge into a bigger black hole of mass $M_{f}$ possessing a final spin $S_{f}$ of the order of [8]

[TABLE]

where $G_{N}$ is Newton’s constant. This result is robust to the details of the evolution, and to the choice of the initial binary orbital angular momentum in the case of quasi-circular orbits. This prediction for $a_{f}$ is compatible with current measurements.

PBHs can be generated in the early Universe through the collapse of large enough density perturbations. In particular, inflation can lead to PBH production when the comoving curvature perturbation ${\cal R}$ is enhanced on small scales (much smaller than the scales accessible to the CMB). The process of reheating after inflation transfers these sizeable inhomogeneities to radiation, and PBHs form when these extreme perturbations re-enter the horizon [3].

While the PBH spins are isotropically distributed at the PBH formation epoch, the probability distribution of the spin amplitude at PBH formation time is not known. This will be the focus of the present work. We shall study the PBH spin probability distribution at formation time. Such a discussion is, to the best of our knowledge, absent in the literature. To proceed, we will use standard results in peak theory [9] which is the appropriate framework since PBHs originate from peaks or, more precisely, from local maxima of the radiation density distribution.

We will show that the zeroth-order anisotropy of the collapsing object giving rise to a PBH is coupled to the first-order tidal gravitational field and generate a first-order torque. This happens if one takes into account the fact that the density profile around peaks generally deviates from spherical symmetry. This result bears a strong similarity with the well-known “tidal-torque theory” in large-scale structure [10, 12, 13, 14, 11, 15]. For instance, as demonstrated in Ref. [12], the classical result of Peebles [16] that the spin of halos is generated at second-order is a consequence of choosing a Lagrangian sphere to describe the collapsing object. In our case, the shape of the Lagrangian region which collapses to form a PBH is imprinted when the characteristic scale of the perturbation is superhorizon.

Our findings indicate that the spin does not vanish at first-order in perturbation theory because, in the homogeneous ellipsoid approximation adopted throughout this work (the gravitational collapse of such ellipsoidal perturbations has been investigated extensively in the large scale structure literature [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]), the lengths of the principal axes of the inertia tensor are different. However, another necessary condition to generate a non-zero spin is the non-vanishing of the off-diagonal components of the velocity shear. In other words, the ellipsoidal perturbation tracing the peak from which the PBH originates must have its inertia tensor misaligned with that of the velocity shear. This does not happen when the perturbations are on superhorizon scales. Notwithstanding, a torque is generated once the perturbation re-enters the horizon until it decouples from the background, i.e. until turnaround. This generates a small, albeit non-zero spin. After turnaround, we speculate that the angular momentum remains constant, as seen for the collapse of non-relativistic dust [16, 12, 14]. However, this issue should be explored further with numerical simulations.

We will show that the typical value of the dimensionless Kerr spin parameter $a_{\rm s}=S/G_{N}M^{2}$ (where $M$ is the mass of the PBH) is

[TABLE]

where $\sigma^{2}_{\delta}$ is the variance of the overdensity at horizon crossing and $\gamma$ measures how narrow is the power spectrum giving rise to the PBHs. Very narrow power spectra have $\gamma\simeq 1$ (a monochromatic power spectrum, or Dirac delta, has exactly $\gamma=1$). Hence, the smaller the PBH spin, the narrower the spectrum. As we shall see, the suppression factor due to the dependence on $\gamma$ arises because, as $\gamma$ approaches unity, the velocity shear tends to be more strongly aligned with the inertia tensor.

The conditional spin probability $a_{\rm s}P(a_{\rm s}|\nu)$ (where $\nu=\delta/\sigma_{\delta}$ parametrises the height of the peak) exhibits some interesting features: for higher peaks, the PBH spin shifts towards smaller values and, once the height of the peak is held fixed, slightly smaller values of the spins are obtained the steeper the power spectrum is.

Our paper is organised as follows. In Section 2, we start from peak theory to characterise the spin at first- and second-order in perturbation theory. In Section 3, we study the spin on superhorizon scales, while Section 4 is dedicated to the subhorizon treatment. Section 5 describes the PBH dimensionless Kerr parameter at formation time. Section 6 contains technical details of the computation of the spin probability distribution. Section 7 illustrates how the PBH spin is correlated with the shape of the power spectrum. Section 8 discusses the impact of the spin onto the PBH abundance. Finally, Section 9 summarizes our conclusions.

2 The spin of PBHs as local density maxima

The definition of spin or angular momentum in general relativity is a delicate issue, see for instance Ref. [28] for a detailed discussion and more recently Ref. [29]. The angular momentum represents a conserved quantity originated from rotational invariance. Naively, one might define it through the (3+1)-formalism as the conserved quantity associated with rotations at spatial infinity and to asymptotically flat spacetime observers. This would be the generalisation of the ADM momentum obtained by replacing the translation Killing vectors at spatial infinity by the corresponding rotational Killing vectors. The problem with this definition is that the resulting quantity does not transform as a vector under the change of the coordinates preserving the asymptotic properties at infinity. The issue arises from the existence of the so-called supertranslations at infinity [30]. However, if the spacetime possesses some symmetries, global quantities may be defined which are coordinate-independent through the technique introduced by Komar [31] which amounts to taking flux integrals of the derivative of the Killing vector associated with the symmetry over closed two-surfaces surrounding the matter sources. These quantities are conserved in the sense that they do not depend on the choice of the surface provided that all the matter is included in it. One can therefore define the amplitude of the Komar angular momentum on a given time-slicing $\Sigma$ through the Komar integral [28]

[TABLE]

where $\phi^{\nu}$ is the rotational Killing vector of the asymptotic flat metric. Using Gauss’ law and Einstein equations, one can rewrite this expression as

[TABLE]

where $J^{\mu}(\phi)=T^{\mu}_{\,\,\nu}\phi^{\nu}-T\phi^{\mu}/2$ is expressed in terms of the energy-momentum tensor. For a relativistic perfect fluid one gets

[TABLE]

where $\vec{v}$ is the velocity field and $\rho$ is the density field. In the following we will be concerned with the spin of perturbations collapsing around a local maxima of the overdensity. Therefore, expanding around the peak located at $\vec{x}_{\text{\tiny pk}}$, one obtains

[TABLE]

where $\vec{x}$ is a comoving coordinate and we have gone to conformal time. The presence of the velocity, which is a first-order quantity in cosmological perturbation theory, makes the spin at least a first-order quantity. Of course, by expanding the energy density $\rho(\vec{x},\eta)$ in perturbation theory one would obtain higher-orders contribution to the spin. However, in the rest of the paper we will be concerned with the spin of PBHs only at first-order.

2.1 First-order description of the PBH spin

The 3-dimensional volume $V_{\rm e}$ over which we perform the integral Eq. (2.4) could be significantly aspherical. To proceed further, we associate the PBHs to high peaks of the overdensity and use the triaxial ellipsoid approximation to prescribe the volume $V_{\rm e}$. We then use the standard results of peak theory and study the statistics of local maxima [9]. Expanding the overdensity around the peak up to second-order (in the derivatives), we obtain222We choose to work with the density contrast instead of the curvature perturbation to avoid un-physical IR effects arising from super-horizon modes, see for example Ref. [32].

[TABLE]

where the first gradient term is zero as we are dealing with maxima of the overdensity. The parameter $f$ is set by the criterion to form a PBH, that is that the matter nearby the peak be above a certain overdensity threshold. Performing a rotation of the coordinate axes to be aligned with the principal axes of length $\lambda_{i}$ of the constant-overdensity ellipsoids gives

[TABLE]

Here, $\sigma_{\zeta}$ is the characteristic rms variance of the components of $\zeta_{ij}$. We purposely avoid the spectral moment notation of Ref. [9] to remain as general as possible. As we will see later, the velocity shear and the density perturbations which contribute to the linear tidal-torque effect can expressed through the familiar spectral moments of the density power spectrum.

Eq. (2.6) can be rearranged as

[TABLE]

Focusing on the approximately ellipsoidal surface $\delta=f\delta_{\text{\tiny pk}}$ and defining the height of the peak $\nu$ in units of the rms overdensity $\sigma_{\delta}$ as $\nu=\delta_{\text{\tiny pk}}/\sigma_{\delta}$, leads to

[TABLE]

whose solutions define the boundary of the integration volume in the spin integral. They are the principal semi-axes of such an ellipsoid and are given by333Notice that $a^{2}_{i}$, being ratios of first-order quantities, are not small quantities and this is the ultimate reason why tidal-torque theory, when applied to dark matter halos, predicts a much larger spin parameter although one starts from a density fluctuation at recombination which is of order $(10^{-4}\div 10^{-3})$.

[TABLE]

Notice that in the limit of large $\nu$ one obtains the following useful expressions [9]

[TABLE]

where

[TABLE]

and $\sigma_{\times}^{2}$ is the characteristic cross-correlation between $\delta$ and $\zeta_{ij}$. Eq. (2.9) shows that the principal semi-axes behaves like $a_{i}^{2}\sim\sigma_{\delta}\nu/\sigma_{\zeta}|\zeta_{ij}|\sim R_{*}^{2}$, where $R_{*}$ is a characteristic scale to be defined below. Furthermore, the difference $|\lambda_{i}-\lambda_{j}|$ is of order unity because the “ellipticity” $\epsilon_{i}$ scales like $1/(\gamma\nu)$.

Next, we first expand the velocity field around the peak

[TABLE]

to get

[TABLE]

In this expression we have expanded at zeroth-order the density field $\rho(\vec{x},t)$, and defined the average density of the universe at a given time $\eta$, $\overline{\rho}_{\rm rad}(\eta)$, which we suppose to be dominated by radiation.

The volume $V_{\rm e}$ over which the integral is performed generally is a function of time. Internal and external tidal forces will gradually deform it until the whole perturbation collapses to a black hole. However, since the characteristic (comoving) size of the perturbation is $k_{\text{\tiny H}}^{-1}\sim M_{\text{\tiny PBH}}^{1/2}$, the deformation becomes significant only once the perturbation has re-entered the horizon and decoupled from the background, i.e. after turnaround. Before turnaround, we can approximate the perturbation as a rigid body of volume $V_{\rm e}$ whose boundary is the isodensity surface $\delta=f\delta_{\text{\tiny pk}}$. The value of $f<1$ is chosen such that $f\delta_{\text{\tiny pk}}$ equates the critical threshold $\delta^{\rm c}$ for PBH collapse. This critical density depends on the choice of gauge (see below).

The spin can be found by performing the integration over the ellipsoidal volume $V_{\rm e}$ appearing in its definition as

[TABLE]

where, without loss of generality, we have chosen $\vec{x}_{\text{\tiny pk}}=0$ and we have factorised out the time-dependence of $v_{l}^{k}$ as

[TABLE]

Eq. (2.15) defines the normalised velocity shear $\widetilde{v}_{kl}$ (where the scale factor can be normalised to unity today). We perform the integration by doing a change of coordinates to ellipsoidal coordinates of the form

[TABLE]

and identifying

[TABLE]

Thus we get

[TABLE]

Finally, Eq. (2.14) can be written as [11]

$\displaystyle\vec{S}^{(1)}=\frac{16\pi}{45}a^{4}(\eta)\overline{\rho}_{\rm rad}(\eta)g_{v}(\eta)a_{1}a_{2}a_{3}([a_{2}^{2}-a_{3}^{2}]\widetilde{v}_{23},[a_{3}^{2}-a_{1}^{2}]\widetilde{v}_{13},[a_{1}^{2}-a_{2}^{2}]\widetilde{v}_{12}).$

(2.19)

This shows that the spin does not vanish at first-order only if i) the lengths of the semi-axis are different and ii) the off-diagonal components of the velocity shear are non-zero, and misaligned with the inertia tensor. For later use, using Eq. (2.9) we can recast the spin in the form

[TABLE]

where

[TABLE]

are defined such that $\alpha_{i}\geq 0$, with $\alpha_{2}\geq\alpha_{1},\alpha_{3}$ due to the labelling choice $\lambda_{1}>\lambda_{2}>\lambda_{3}$. Here,

[TABLE]

The term inside the square brackets in the last line of Eq. (2.1) identifies the time-dependent part of the spin, that we will label as the reference spin

[TABLE]

It is the same for all peaks, while the information on the shape and the height of the peak is carried by the remaining term. The magnitude of the dimensionless spin can therefore be written as

[TABLE]

While $S_{\rm ref}(\eta)$ is common to all peaks, $s^{(1)}_{\rm e}$ depends on the shape and the height of the peaks. This parametrisation isolates the time-dependence in such a way that its time derivative may be thought of as the torque acting on matter in the neighbourhood of local density maxima.

As we can see, since $\alpha_{i}\sim|\lambda_{j}-\lambda_{k}|/\lambda_{j}\lambda_{k}\sim 1/(\gamma\nu)^{2}$, $s^{(1)}_{\rm e}$ genuinely is a quantity of order ${\cal O}(1)$ (up to multiplicative factors of $\gamma$ and $\nu$). The only first order quantity appearing in this calculation thus is the velocity shear $v_{l}^{k}$. Therefore, the spin $S^{(1)}(\eta)$ is truly first-order.

2.2 Second-order description of the PBH spin

For the sake of completeness, we present here the second-order expansion obtained taking one more term when perturbing the energy density

[TABLE]

This corresponds to consider the second-order effects in the spin

[TABLE]

where $V_{\rm e}$ is the ellipsoidal volume. We perform the integrations by doing a change of coordinates to ellipsoidal coordinates as in Eq. (2.16) and, defining $I^{jlmn}=\int_{V_{\rm e}}{\rm d}^{3}x\,x^{j}x^{l}x^{m}x^{n}$, one finds

[TABLE]

where all the other components are zero. The spin up to second-order is therefore given by

[TABLE]

where the components of $\vec{S}^{(2)}$ are

[TABLE]

We will now discuss the two relevant regimes, super- and sub-horizon.

3 The PBH spin before horizon crossing

In the early radiation-dominated Universe, PBHs are generated when highly overdense regions gravitationally collapse directly into a black hole. Before collapse, the comoving sizes of such regions are larger than the horizon length and the separate universe approach can be applied [33]. We can therefore expand at leading order in spatial gradients of the various observables, e.g. the overdensity. At this stage, one should fix the slicing and the threading of the spacetime manifold. The so-called constant mean curvature slicing (CMC, sometimes dubbed the uniform Hubble slicing) seems appropriate as it has been adopted to perform numerical relativity simulations to describe the formation of PBHs and to calculate the threshold for PBH formation [34]. In this slicing the equations for the lapse function and the gravitational potentials are similar to those one encounters in the maximal slicing gauge (in which the total curvature vanishes) in the asymptotically flat spacetime. We will provide further considerations about this gauge in section 4.

In the CMC slicing, the overdensity turns out to be [33]

[TABLE]

where $\mathcal{R}(\vec{x})$ is the comoving curvature perturbation. This expression is fully non-perturbative, being the only expansion in gradients. Notice that the coefficient is 3/2 times larger than the one in the comoving slicing (CG), as it is the corresponding threshold for PBH formation, i.e. $\delta^{\rm c}_{\text{\tiny cmc}}=3\delta^{\rm c}_{\text{\tiny CG}}/2$. We are going to consider, as a reference value, $\delta^{\rm c}_{\text{\tiny CG}}=0.45$, and therefore $\delta^{\rm c}_{\text{\tiny cmc}}=0.675$ [33], neglecting the shape dependence of the threshold [35, 36, 37].

As the universe expands, the overdensity grows. Regions where it becomes of order unity eventually stops expanding and collapse. This happens when the comoving scale of such a region becomes of the order of the horizon scale. Even though the gradient expansion approximation breaks down, it has been used to obtain an acceptable criterion for the PBH formation, as confirmed by nonlinear numerical studies [3], and we will follow the same strategy here. The spin grows until the system decouples from the background expansion and torques are reduced.

In the threading where the pertubed 0i component of the metric vanishes, the velocity in the long wavelength regime is given by [33]

[TABLE]

The velocity therefore scales like $a^{3}$ on superhorizon scales.

Since the velocity is proportional to the gradient of the overdensity, from Eq. (2.6) one immediately concludes that the off-diagonal entries of the matrix $\widetilde{v}_{ij}$ are vanishing. Therefore, before horizon-crossing the first-order PBH spin $\vec{S}^{(1)}$ is zero. The physical reason is that, in the coordinate frame aligned with the principal axis frame of the inertia tensor, the velocity shear $v_{l}^{k}$ is aligned with the inertia tensor and no spin can be generated independently of the deviation from sphericity of the collapsing region and at any order in perturbation theory in the comoving curvature perturbation $\mathcal{R}(\vec{x})$.

We emphasize that this is true for any particular realization of a density peak, regardless of the probability density of the random variables $\delta$, $v^{i}$ etc. As soon as $v^{i}$ is proportional to $\partial^{i}\delta$, diagonalizing $\zeta_{ij}(\vec{x}_{\text{\tiny pk}})$ will automatically diagonalize $\widetilde{v}_{kl}(\vec{x}_{\text{\tiny pk}})$.

As we can see from Eq. (2.29), the spin is also zero at second-order and, presumably, at any higher order since the relation $v^{i}\propto\partial^{i}\delta$ holds at any order in perturbation theory as long as $k\ll{\cal H}$. Furthermore, it is also true in the more familiar Newtonian longitudinal gauge where $\delta=4{\cal R}/3$ and the velocity is proportional to $\partial^{i}{\cal R}$ and therefore to $\partial^{i}\delta$.

4 The PBH spin after horizon crossing

From horizon crossing onwards, i.e. when the characteristic wavelength of the perturbations becomes smaller than the Hubble radius, the relation between the velocity and the overdensity changes. As a result, the spin can grow briefly owing to the linear tidal torque until the perturbation decouples from the background. Equating the free-fall and sound crossing timescale, the Jeans criterion tells us that gravitational instability occurs for perturbation with comoving wavenumber

[TABLE]

up to a factor of order unity. We have assumed $\delta_{\text{\tiny pk}}\sim 1$ in this estimate. Since the sound speed satisfies $c_{\text{\tiny S}}\simeq 1/\sqrt{3}$ deep in radiation domination, the perturbations must decouple from the background around horizon crossing in order to form PBHs. Otherwise, radiation pressure would quickly stabilise them. Assuming that turnaround occurs at horizon crossing, we can estimate the amount of angular momentum acquired by the perturbation through linear tidal-torque.

To be consistent with the previous section, we will work in the CMC slicing. Since the reader might be more familiar with the Newtonian longitudinal gauge we first summarise the basic steps to go from the Newtonian longitudinal gauge to the CMC gauge.

4.1 Relating Newtonian longitudinal gauge and the CMC gauge

To relate the CMC gauge and the Newtonian slicing and longitudinal gauge we follow the standard reference [38]. The subset of metric perturbations we will be concerned with are

[TABLE]

where $Q^{(0)}(x^{\mu})$ is the scalar harmonic and $Q_{i}^{(0)}(x^{\mu})$ and $Q_{ij}^{(0)}(x^{\mu})$ are the corresponding vector and tensor. Taking a coordinate transformation of the form

[TABLE]

one gets the following transformations

[TABLE]

We now impose the CMC gauge to be characterised by uniform-Hubble-constant hypersurfaces, that is

[TABLE]

and we impose this condition going from the Newtonian longitudinal gauge to the CMC gauge. In the Newtonian longitudinal gauge we use the notation of Ref. [39], for which $A=\Phi$, $B=H_{T}=0$ and $H_{L}=-\Psi=-\Phi$. We omit labels to designate longitudinal gauge quantities for the purpose of avoiding clutter. Furthermore, we have assumed zero shear or, equivalently, zero anisotropic stress (which is an excellent approximation at PBH formation since the mean free path is much shorter than any relevant scale.)

In a radiation-dominated universe the corresponding amplitude of change in the time coordinate is

[TABLE]

where the amplitude of the velocity is $v^{i}=-i(k^{i}/k)v$. This fixes the time-slicing. The corresponding change of the overdensity is then

[TABLE]

On superhorizon scales where $\Phi(\vec{k},0)=-2{\cal R}/3$ does not depend on time, this gives

[TABLE]

Since $\delta_{\text{\tiny cmc}}(\vec{k},0)=(2/3)(k/{\cal H})^{2}{\cal R}$, one correctly reproduces the condition $\delta(\vec{k},0)=-2\Phi(\vec{k},0)$ expected in the limit $k\ll{\cal H}$ for the Newtonian longitudinal gauge. Physically, this reflects the fact that a constant shift in the gravitational potential does not yield any observable effect. In other words, PBHs must trace peaks of the familiar (i.e. CMC) density field rather than peaks of the gravitational potential.

While we have fixed the CMC slicing, we have not determined the threading yet. We can do so by imposing $B^{\text{\tiny cmc}}=0$

[TABLE]

where we have inserted the zero to remind the reader that $B=0$ in the longitudinal gauge. This fixes

[TABLE]

Given the fact that the velocity transforms as

[TABLE]

the gauge-invariant velocity is $(v-H_{T}^{\prime}/k)$. This leads to

[TABLE]

This means that the amplitude of the velocity in the CMC gauge is

[TABLE]

where we have highlighted the fact that in the Newtonian longitudinal gauge $H_{T}=0$. Since

[TABLE]

we finally get [38]

[TABLE]

Since on superhorizon scales $v(\vec{k},0)\simeq-(1/3)(k/{\cal H}){\cal R}$, we get, in the same limit,

[TABLE]

which nicely reproduces the result mentioned in the previous Section. On subhorizon scales the CMC velocity coincides of course with the Newtonian longitudinal velocity.

4.2 Subhorizon description in the CMC gauge

The expressions (4.7) and (4.15) allow us to compute the overdensity and the velocity shear in the CMC gauge once the perturbations re-enter the horizon, in terms of quantities known in the Newtonian longitudinal gauge.

From now on we consider a perturbation with characteristic comoving scale entering the horizon at $k_{\text{\tiny H}}\gg k_{\rm eq}$, which is the regime relevant to PBH formation deep in radiation domination. Here, since the mean free path is very short, one can again ignore any anisotropic stress and take the Bardeen potentials to be $\Phi=\Psi$. Near horizon crossing the radiation density perturbation reads [39]

[TABLE]

where

[TABLE]

Correspondingly, the Fourier modes of the radiation bulk velocity are given by

[TABLE]

This expression holds as long as $k\gtrsim k_{\rm eq}$. Since the comoving Hubble parameter satisfies ${\cal H}=\eta^{-1}$, then horizon crossing occurs at $\eta_{\text{\tiny H}}=k_{\text{\tiny H}}^{-1}$, with $k_{\text{\tiny H}}\propto M_{\text{\tiny PBH}}^{-1/2}$. The multiplicative sine factor in the expression of the radiation bulk velocity becomes

[TABLE]

For $k=k_{\text{\tiny H}}$, we have $\sin(k_{\text{\tiny H}}/k_{\text{\tiny S}})\simeq 0.546$. Furthermore, the (physical) velocity shear at horizon crossing is given by

[TABLE]

while the density perturbation is

[TABLE]

Using the expressions (4.7) and (4.15) we finally get

[TABLE]

and

[TABLE]

At this point, it is convenient to introduce a normalised density $\nu$, shear $\widetilde{v}_{ij}$ and Hessian $\zeta_{ij}$ to retain the analogy with [11]:

[TABLE]

Here, $W(k)$ is the Fourier transform of a spherically symmetric window function with characteristic wavenumber $k_{\text{\tiny H}}$, whereas the “transfer functions” $T_{\delta}(k,\eta_{\text{\tiny H}})$ and $T_{v}(k,\eta_{\text{\tiny H}})$ are

[TABLE]

and

[TABLE]

Notice also that Eq. (4.2) makes it explicit that in the sub-horizon regime the velocity shear is not proportional to the gradients of the density contrast. As we can see from Eq. (4.2), the rms variance $\sigma_{\delta_{\text{\tiny cmc}}}$ and $\sigma_{\zeta_{\text{\tiny cmc}}}$ can be constructed from the spectral moments

[TABLE]

This shows that

[TABLE]

Finally we define

[TABLE]

where we approximated the contribution of the transfer functions to a constant, since in all the relevant cases the power spectrum is peaked at $k=k_{\text{\tiny H}}$. The numerical value is found to be $T_{v}(k_{\text{\tiny H}},\eta_{\text{\tiny H}})/T_{\delta}(k_{\text{\tiny H}},\eta_{\text{\tiny H}})\sim 0.5$.

Since $\widetilde{v}_{\text{\tiny cmc}ij}$ does not vanish and is not aligned with the Hessian of $\delta_{\text{\tiny cmc}}$, some angular momentum is generated at linear order. With the assumption that turnaround occurs close to horizon crossing, the corresponding reference spin is

[TABLE]

This determines the amplitude of the Kerr parameter at turnaround as we shall see next.

We conclude this section by some comments on the computation in the Newtonian longitudinal gauge. To calculate the spin acquired through linear tidal-torque, one should carefully handle the superhorizon contributions. While the velocity shear $\partial_{i}v^{j}$ is physical, the superhorizon limit $\delta=-2\Phi({\vec{k}},0)$ is irrelevant as it cannot affect the collapse of the perturbation. Therefore, in the expression of the Hessian of the density field

[TABLE]

using the first term in the right-hand side to define the shape of the ellipsoid would be equivalent to select peaks of the gravitational potential or, equivalently, the curvature. However, one cannot impose any constraint on the value of the gravitational potential at $\vec{x}=\vec{x}_{\text{\tiny pk}}$ because this would not lead to any observable effect. To remedy this problem, one should subtract the superhorizon contribution to the density in the calculation of the covariance etc. and consider instead

[TABLE]

As a result, the density Hessian that defines the ellipsoidal shape of the perturbation becomes

[TABLE]

which is suppressed by an additional factor of $(k/k_{\text{\tiny S}})^{2}$. This would ensure that the high density fluctuations which collapse to form PBHs trace the peaks of the radiation density field rather than the peaks of the curvature.

5 An estimate of the PBH dimensionless Kerr parameter at formation time

From now on, we will work at first-order in perturbation theory and, therefore, we will remove the label (1) to avoid cluttering notation. One dimensionless parameter which is of interest to us is the Kerr parameter

[TABLE]

such that

[TABLE]

Substituting $a^{4}(\eta_{\text{\tiny H}})\overline{\rho}_{\rm rad}(\eta_{\text{\tiny H}})=\overline{\rho}_{\rm rad}(\eta_{0})$ gives

[TABLE]

where $\eta_{0}$ identifies the present time. This quantity has been evaluated at the time $\eta_{\text{\tiny H}}$ when the relevant modes of the perturbation re-enter the horizon to collapse into a PBH. In the subhorizon regime, the leading time dependent factor in the velocity definition $g_{v}(\eta_{\text{\tiny H}})$ scales like $g_{v}(\eta_{\text{\tiny H}})=(1/2){\cal H}(\eta_{\text{\tiny H}})\sigma_{\delta_{\text{\tiny cmc}}}(\eta_{\text{\tiny H}})$. In the radiation phase the known behaviour ${\cal H}(\eta)\propto a^{-1}(\eta)$ leads to the relation

[TABLE]

where the subscript eq denotes the radiation-matter equality, $\Omega_{\text{\rm dm}}$ represents the present DM abundance, and where we used the fact that in the matter-dominated era ${\cal H}\sim a^{-1/2}$. We can compute the Hubble rate at crossing time

[TABLE]

where we used the relation

[TABLE]

to express the Hubble rate at the horizon crossing as a function of the primordial black hole mass $M$ and Newton coupling constant $G_{N}$. We then find

[TABLE]

For a power spectrum peaked at the scale $R_{*}$ one has typically $(1-f)\sim 1/3$ (in the CMC gauge) and $R_{*}\sim\sqrt{3}k_{\text{\tiny H}}^{-1}\sim\sqrt{3}{\cal H}^{-1}(\eta_{\text{\tiny H}})$ [37], such that

[TABLE]

and substituting Eq. (5.5) in the previous equation gives

[TABLE]

An estimate for the Kerr parameter at the time of formation is, thus, provided by the simple relation

$\displaystyle a_{\rm s}=A(\eta_{\text{\tiny H}})s_{\rm e}\sim\left[\frac{1}{2\pi}\sigma_{\delta_{\text{\tiny cmc}}}(\eta_{\text{\tiny H}})\right]s_{\rm e}.$

(5.10)

Notice that $a_{\rm s}$ scales correctly as a first-order quantity. To get the feeling of the numbers, we take $\sigma_{\delta_{\text{\tiny cmc}}}(\eta_{\text{\tiny H}})=\delta^{\rm c}_{\text{\tiny cmc}}/\nu\sim 0.08$ to get

[TABLE]

were we used the indicative value $\nu=8$. The details are presented in Sec. 7 and may change upon the shape of the power spectrum. More importantly, the determination of the probability distribution of $a_{\rm s}$ at the time of formation is directly related that of $s_{\rm e}$. Therefore, our next step is to study in detail the distribution of $s_{\rm e}$. From Eq. (2.1) one expects $s_{\rm e}$ to scale like (for some $i\neq j$)

[TABLE]

where we have taken into account that the velocity shear scales like $\sqrt{1-\gamma^{2}}$ (this point will become more clear in the next section). As a result, we expect $a_{\rm s}$ to be of the order of

[TABLE]

This estimate will be confirmed by the investigation of the probability distribution of $a_{\rm s}$ which will indeed be peaked around the value (5.13). In particular notice that in the limit of very narrow power spectra, for which $\gamma$ tends to unity, the spin vanishes. This has a physical reason: as we will show in the next section, see Eq. (6.6), the velocity shear has a strong tendency to align itself along the principal axis of the mass ellipsoid due to the strong correlation with the inertia tensor. In the limit $\gamma=1$ this alignment is total and no spin can be generated.

6 The PBH spin and the statistics of local maxima

Following [11], the starting point for the calculation of the probability of the spin $s_{\rm e}$ is a joint distribution which involves the sixteen variables

[TABLE]

where the last two matrices have six components, being both symmetric. We decide to drop the label CMC from all quantities to simplify the subsequent equations. We can create a vector $V$ of sixteen components with the joint distribution

[TABLE]

where the covariance matrix $\bf M$ is given by (in the following we will use the fact that $\langle V_{i}\rangle=0$)

[TABLE]

In the subhorizon regime, the sixteen variables are correlated as follows

[TABLE]

where the ellipsis stands for the other components and the remaining correlators are vanishing. We can rearrange the variables in the dimensionless forms

[TABLE]

which are correlated as follows

[TABLE]

where $\gamma=\sigma_{\times}^{2}/\sigma_{\delta}\sigma_{\zeta}$. Again, the zero correlators are not reported here. Since $\langle v_{A}\nu\rangle^{2}=\langle v_{A}^{2}\rangle\langle\nu^{2}\rangle=1$, then the variables $v_{A}$ and $\nu$ are correlated, so we can drop $v_{A}$ reducing to only fifteen independent variables. Notice that in the limit $\gamma=1$ the velocity shear is totally aligned with the inertia tensor.

Focusing now on the principal axes of the matrix $-\zeta_{ij}/\sigma_{\zeta}$, we can use its three eigenvalues plus the three Euler angles, describing the orientation of the principal axes, as new variables in the place of the $\zeta_{ij}$. In particular

[TABLE]

Being the system independent on these angles, we can integrate over them, leaving twelve independent variables. The corresponding distribution is

[TABLE]

where we have defined, starting from

[TABLE]

[TABLE]

and

[TABLE]

In the region around a peak we can make a Taylor expansion of the kind

[TABLE]

such that

[TABLE]

Integrating the distribution of Eq. (6.8) over space eliminates three variables more, the $\widetilde{\zeta}_{i}$, and leaves only nine variables. The distribution obtained in such a way describes the comoving number density of peaks in the element ${\rm d}\nu{\rm d}v_{B}{\rm d}v_{C}{\rm d}^{3}w_{i}$. Furthermore, we can integrate over $v_{B}$ and $v_{C}$, as these are unconstrained, leaving finally

[TABLE]

where

[TABLE]

and

[TABLE]

with $x_{*}=\gamma\nu$ and $w^{2}=w_{1}^{2}+w_{2}^{2}+w_{3}^{2}$. This last quantity identifies the squared radius of a polar coordinate system, with angles $\theta={\rm arccos}(u)$ and $\phi$, for the variables $w_{i}$, that we will use from now on.

The next step is to find the probability for the spin $s_{\rm e}$. To do so, we can rewrite Eq. (6.14) in the following way

[TABLE]

as a function of the previously defined dimensionless spin through Eq. (2.1)

[TABLE]

where

[TABLE]

In the limit $\gamma$ very close to unity or, equivalently, $\Gamma\ggg 1$ (that is for a monochromatic spectrum) the distribution (6.14) for $w$ approaches a product of Dirac deltas

[TABLE]

As a result, the value $\vec{w}=0$ is selected and the off-diagonal terms of the velocity shear vanish and therefore the velocity shear is aligned with the inertia tensor. This implies that the spin parameter of PBH collapsing from a monochromatic feature is zero. Furthermore, the expression (6.14) shows us that the typical value of the velocity shear scales like $1/\sqrt{\Gamma}=\sqrt{1-\gamma^{2}}$, a scaling we have been using to estimate $s_{\rm e}$ in Eq. (5.12).

Integrating now over $w$, $\phi$ and $\lambda_{i}$ (the integration over $\phi$ is replaced then by the integration over $\beta$) gives the distribution of peaks of given height $\nu$ and spin $s_{\rm e}$

$N_{\text{\tiny pk}}(\nu,s_{\rm e})=\frac{4Cs_{\rm e}}{R_{*}^{3}\nu^{5}}\int_{0}^{\infty}{\rm d}\lambda_{1}\int_{0}^{\lambda_{1}}{\rm d}\lambda_{2}\int_{0}^{\lambda_{2}}{\rm d}\lambda_{3}\int_{\alpha_{1}}^{\alpha_{2}}{\rm d}\beta\frac{e^{-Q_{5}}F(\lambda_{i})\Lambda T(\alpha_{3},\beta,s_{\rm e},\nu)}{\sqrt{|(\alpha_{1}^{2}-\beta^{2})(\alpha_{2}^{2}-\beta^{2})(\alpha_{3}^{2}-\beta^{2})|}},$

(6.21)

where

[TABLE]

and

[TABLE]

Here the argument of the Dawson’s integral $D(X)$ is given by

[TABLE]

The conditional differential probability for $s_{\rm e}$, given that the peak has a height $\nu$, is the given by [11]

[TABLE]

where the comoving differential peak density is expressed as

[TABLE]

in terms of the function

[TABLE]

The function $f(x)$ that appears in the previous formula is provided by the expression

[TABLE]

Fig. 1 shows the behaviour of the comoving peak density of Eq. (6.26) for several choices of $\gamma$. Fig. 2 shows the behaviour of the distribution $P(s_{\rm e}|\nu)$ for relevant choices of $\nu$. Finally, we can use the relation shown in Eq. (5.10) to find the distribution for the Kerr parameter $a_{\rm s}$, plotted in Fig. 3.

We note that the conditional probability distribution for the spin of a material in the vicinity of peaks of different heights shows a systematic shift towards smaller values for higher peaks if the parameter $\gamma$ is close to unity. Furthermore, for a given value of the height peak $\nu$, higher values of $\gamma$ provide slightly smaller values of spins, and one can check that the scaling with $\gamma$ is indeed like $\sqrt{1-\gamma^{2}}$.

6.1 The analytical approximation in the limit of high peaks

The formation of PBHs requires high thresholds and it is therefore interesting to investigate analytically the large $\nu$ limit of the probability $P(s_{\rm e}|\nu)$. To do so, we define a normalised dimensionless spin $h$ as444We note that Eq. (C.7) of Ref. [11] has the incorrect scaling with $\Gamma$. The factor $\Gamma^{1/2}$ should be $\Gamma^{-1/2}$.

[TABLE]

The scaling with $\sqrt{1-\gamma^{2}}$ is dictated by the scaling of the velocity shear, see Eq. (5.12). The distribution of the parameter $h$ can be analytically approximated as

[TABLE]

Such a distribution is the result of a best-fit that was performed for the values $\gamma=0.9$ and $\nu=8$. We checked numerically that it holds for the relevant parameter space related to the physics of PBH formation. In Fig. 4 we compared the numerical result and analytical approximation for the probability distribution $P(s_{\rm e}|\nu)$.

The analytical expression for the distribution of the Kerr parameter $a_{\rm s}$ is found to be in the large $\nu$ limit

$\displaystyle P(a_{\rm s}|\nu){\rm d}a_{\rm s}$ $\displaystyle=$ $\displaystyle\left(\frac{5\gamma^{6}\nu}{2^{9/2}\pi}\frac{\Gamma^{1/2}}{A(\eta_{\text{\tiny H}})}\right)\exp\left[-2.37-4.12\ {\ln\left(\frac{5\gamma^{6}\nu}{2^{9/2}\pi}\frac{\Gamma^{1/2}}{A(\eta_{\text{\tiny H}})}a_{\rm s}\right)}\right.$

(6.31)

$\displaystyle-$ $\displaystyle\left.1.53\ {\ln^{2}\left(\frac{5\gamma^{6}\nu}{2^{9/2}\pi}\frac{\Gamma^{1/2}}{A(\eta_{\text{\tiny H}})}a_{\rm s}\right)}-0.13\ {\ln^{3}\left(\frac{5\gamma^{6}\nu}{2^{9/2}\pi}\frac{\Gamma^{1/2}}{A(\eta_{\text{\tiny H}})}a_{\rm s}\right)}\right]{\rm d}a_{\rm s},$

where $A(\eta_{\text{\tiny H}})$ is defined in Eq. (5.10). This is the main result of our paper.

7 The PBH spin and the shape of the power spectrum

In the previous section we have seen that the spin distribution is characterised only by two parameters, namely $\nu$ and $\gamma$. The parameter $\nu$ describes the height of the peaks in terms of the variance and is determined by the required abundance of PBHs. In particular, assuming a Gaussian statistics (for its non-Gaussian extension, see Ref. [40]), the PBH mass fraction can be expressed as [3]

[TABLE]

where $\nu$ is defined as $\nu=\delta^{\rm c}_{\text{\tiny cmc}}/\sigma_{\delta_{\text{\tiny cmc}}}$, which is gauge independent. The approximate relation between $\beta$ and $\nu$ is shown in Fig. 5.

Assuming a monochromatic spectrum of masses and imposing the PBH to be the dark matter, one may deduce the relation between the mass fraction and the PBH mass to be

[TABLE]

For models forming PBH in the physically relevant mass range, the parameter $\nu$ lies in the range

[TABLE]

The power spectra of the density perturbation are directly calculable in terms of the power spectra of the comoving curvature perturbation $\zeta$ as

[TABLE]

The quantity ${\cal P}_{\zeta,{\rm tot}}(k)$ is the sum of two pieces, the smooth power spectrum giving the correct amplitudes for perturbations at the CMB scales and the term responsible for the formation of PBH at small scales, that is

[TABLE]

where $k_{\rm p}$ denotes a pivot scale, $A_{s}\sim 2\cdot 10^{-9}$ the corresponding amplitude and $n_{\zeta}$ the spectral index of scalar perturbations. To compute the momenta of these distributions, we are going to use the volume normalised Gaussian window function of the form $W(k)=\exp[-(kR_{\text{\tiny H}})^{2}/2]$ smoothing out perturbations on scales different from the characteristic scale $R\text{\tiny H}$ corresponding to the cosmological horizon at formation.

These momenta enter in the computation of the parameter $\gamma=\sigma_{\times\text{\tiny cmc}}^{2}/\sigma_{\delta_{\text{\tiny cmc}}}\sigma_{\zeta\text{\tiny cmc}}$, which contains all the relevant information on the shape of the power spectrum of the density perturbations. In the following subsections we consider few examples of typical power spectra and compute their characteristic parameter $\gamma$.

7.1 Log-normal power spectrum

We start by considering a log-normal power spectrum of the form

[TABLE]

In general, the factor $\gamma$ depends only on the parameter $\sigma$ as shown in Fig. 6.

We can define the parameter $Q=k_{*}/\Delta k$ to indicate how peaked is the power spectrum, where $\Delta k$ stands for the full width at half of the maximum value (FWHM). The value $Q$ for various shapes is plotted in Fig. 6. For narrow power spectra, the parameter $\gamma$ approaches unity, while for wider profiles it quickly reaches the asymptotic value $\gamma\sim 0.82$. For example, a choice of parameters consistent with the totality of dark matter composed by PBH of mass $M=10^{-12}M_{\odot}$ would be $A=0.066$, $\sigma=0.5$ and $k_{*}=3\cdot 10^{12}\,{\rm Mpc}^{-1}$ [41].

In the limit of small enough $\sigma$, the log-normal power spectrum can also serve as a good approximation of a Dirac power spectrum, usually used in the literature. For example, let us analyse the case of $\sigma=0.042$ corresponding to a highly peaked power spectrum. From the explicit expression of the power spectrum we can compute its momenta

[TABLE]

and

[TABLE]

Being independent on $A$, our choice of $\sigma$ determines the value $Q=10$.

7.2 Nearly flat power spectrum

For a power spectrum of the form [42]

[TABLE]

we may fix the parameters such that PBHs with masses $\sim M_{\odot}$ form the dark matter. This gives $n_{\rm p}\sim 0.96$, $A_{f}=0.0308$ and $k_{\rm min}=10^{3}{\rm Mpc}^{-1}$. For such a case we find $\gamma\simeq 0.88$.

7.3 Broken power law power spectrum

One can also consider the case of a power spectrum of the form [43]

[TABLE]

where an example of a parameter set giving the PBHs to form the dark matter is $m=3$, $n=0.5$, $A_{b}=0.0413$ and $k_{\rm p}=2\cdot 10^{6}\,{\rm Mpc^{-1}}$, corresponding to a population of PBH peaked at $M\sim M_{\odot}$. This set gives $\gamma=0.85$.

8 The impact of the spin onto the PBH abundance and the spin distribution function

Once we have calculated the spin distribution of PBH at formation time, we may study the impact of the spin on the abundance of PBHs. This exercise has been already provided in Ref. [44] where, however, the spin distribution has been assumed to be flat. In order to compare with their results, we are going to assume the same parameter dependence

[TABLE]

where \gamma_{\text{\tinyM}}=0.3558 [45], \gamma_{S}=(5/2)\gamma_{\text{\tinyM}}=0.8895 and $K=0.005685$ [46]. The parameter $q$ describes the rotation and is related to the Kerr parameter $a_{\rm s}=S/G_{N}M^{2}$ by the relation

[TABLE]

The parameter $\delta_{\rm c}(q)$ is in turn related to the density fluctuations $\delta$ by the relation

[TABLE]

such that we can easily go from the variables $\delta$ and $q$ to $M$ and $a_{\rm s}$, respectively.

The PBH density parameter at formation epoch is given by

[TABLE]

where, for demonstrative purposes, we take $P(\delta)$ to be a Gaussian distribution with variance $\sigma$. Since the integration measure is transformed as

[TABLE]

the density abundance becomes

[TABLE]

where $P(a_{\rm s})$ is the distribution of the Kerr parameter $a_{\rm s}$ we have previously calculated. The PBH spin distribution is therefore finally obtained by computing

[TABLE]

where we have to perform the integration over the masses $M$. The amount of PBHs formed at a given epoch by the collapse of a region is described by555 As we have seen, the generation of a first-order non-zero spin during the collapse of the overdensity is due to the small deviations from spherical symmetry. One should recall that this has an impact also on the threshold value at which the collapse takes place [47].

[TABLE]

In Fig. 7 we plot ${\rm d}\Omega/\beta{\rm d}a_{\rm s}$ and ${\rm d}\Omega/\beta{\rm d}\ln M$ for the distribution $P(a_{\rm s})$ of the Kerr parameter obtained numerically.

We observe that the Kerr parameter is distributed with a peak which shifts towards smaller values for higher values of $\gamma$, while the impact on the mass distribution is rather negligible.

Our results assume that PBH spin after the formation does not evolve appreciably. For a rough estimate we report here the considerations put forward in Ref. [44]. After formation, the spin of PBH evolves until the present day under the effect of the torque generated by the background radiation fluid. The interaction with the background radiation has the effect of reducing the spin. In particular, one can assess the loss of $S$ as (passing to cosmic time and setting $8\pi G_{N}$ to unity)

[TABLE]

where we have estimated the force generating the torque as $F_{\rm rad}\sim\overline{\rho}_{\rm rad}M^{2}a_{\rm s}$ [48, 49, 44] and used the fact that in a radiation dominated universe $H^{2}\sim\overline{\rho}_{\rm rad}$. On the other hand, the PBH mass evolves by accretion from the surrounding radiation fluid as

[TABLE]

Solving the previous equation tells us that the mass evolution is described by [50]

[TABLE]

where $t_{\text{\tiny H}}$ and $M_{\text{\tiny H}}$ are the initial time and mass, respectively. One can conclude that for PBH with initial masses $M_{\text{\tiny H}}$ smaller than the horizon mass $t_{\text{\tiny H}}$, their accretion is small. Finally, the rate of change of the spin parameter is estimated to be

[TABLE]

resulting in an evolution given by [44]

[TABLE]

Therefore, one might expect that the spin parameter evolution can be safely neglected after formation when M_{\text{\tiny H}}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip 0.5pt\sim}\hss}\raise 1.0pt\hbox{<}}t_{\text{\tiny H}}.

9 Conclusions

In this paper, we have computed the spin probability distribution of PBHs at formation time. We have shown that two ingredients are crucial to generate a spin at first-order in perturbation theory: the collapsing object must be non-spherical, and the velocity shear be misaligned with the inertia tensor of the object. If these conditions are satisfied, a spin of order $10^{-2}$ is generated owing to the action of a first-order tidal torque as the perturbation re-enters the horizon. The spin probability distribution of PBHs is such that smaller values of the spin are achieved for a collapsing perturbation with higher peak significance and narrower power spectrum.

Our results could be improved in several ways. On the one hand, we have assumed that the initial distributions of the overdensity and velocity are Gaussian. It would be worth assessing what is the impact of a primordial non-Gaussianity by extending our results from peak theory to non-Gaussian distributions. Since PBHs are created from the tails of the probability distribution, we expect this may have an impact on the spin as well.

On the other hand, we have not carefully followed the evolution of the spin after formation. Even though our rough estimates indicate that changes will not be appreciable, it has been argued that spinning PBHs suffer from super-radiant instabilities in the radiation phase (for masses $M\lesssim 0.1a_{\rm s}M_{\odot}$) and this may lead to a sizeable reduction of the spin $\delta a_{\rm s}\sim a_{\rm s}$ [51]. We leave these interesting points for future investigation.

Let us conclude with a brief comparison to the recent literature on the same topic. In Ref. [44], the PBH spin distribution has been derived by integrating the probability $P(s_{\rm e},\delta)$ (or better their proxies $a_{\rm s}$ and $M_{\text{\tiny PBH}}$) over $\delta$ under the assumptions that there is no correlation between the initial overdensity and $s_{\rm e}$, and the probability density for $s_{\rm e}$ is flat. This represents a limiting case in which the initial overdensity is allowed to have a large angular momentum. Our computations show that the assumption of a flat distribution for $s_{\rm e}$ is not a good starting point.

Furthermore, another piece of work, Ref. [29], appeared on the same subject of PBH spin at formation time while we were completing our draft. Even though the spin probability distribution is not discussed there, our results coincide with theirs as far as the percent level of PBH spin is concerned. We also agree on the statement that the PBH spin should vanish in the limit of very narrow power spectra, even though for a different reason. Namely, we find that, in this limit, the probability for the off-diagonal elements of the velocity shear is peaked around zero, thus delivering zero spin. This follows from the tendency of the velocity shear to align with the inertia tensor for very narrow power spectra, which inhibits the generation of spin. By contrast, Ref. [29] have calculated the spin at second-order in perturbation theory, $a_{\rm s}\sim{\cal R}^{2}$, in analogy with the argument of [16] that the total angular momentum contained in a spherical proto-galaxy starts at second-order in perturbation theory. As we explained in the introduction, this result is a consequence of the choice of a Lagrangian sphere to describe the collapsing object. When the ellipsoidal shape of the density profile around the peaks is taken into account, the spin can grow at first order for generic power spectra, in analogy with the argument of Refs. [10, 12].

Acknowledgments

V. DL., G.F. and A.R. are supported by the Swiss National Science Foundation (SNSF), project The Non-Gaussian Universe and Cosmological Symmetries, project number: 200020-178787. V.D. acknowledges support by the Israel Science Foundation (grant no. 1395/16).

Appendix A Covariance matrix

In this Appendix we compute the covariant matrix ${\bf M}$, defined in Eq. (6.3). All quantities are defined in the CMC gauge and we decide to drop the label to simplify the notation. Let us start with the computation of the correlator $\langle{\delta^{2}}\rangle$. We know from the definitions that

[TABLE]

thus

[TABLE]

Now we can consider all the others. Indeed

• •

$\langle{\delta}\,\zeta_{ij}\rangle$:

[TABLE]

• •

$\langle\zeta_{i}\zeta_{j}\rangle$:

[TABLE]

• •

$\langle\zeta_{ij}\zeta_{mn}\rangle$:

[TABLE]

therefore one finds $\langle\zeta_{11}\zeta_{11}\rangle=3\langle\zeta_{11}\zeta_{22}\rangle=3\langle\zeta_{12}^{2}\rangle=...=\sigma_{\zeta}^{2}/5$;

• •

$\langle\delta\,\widetilde{v}_{ij}\rangle$:

[TABLE]

• •

$\langle\zeta_{ij}\widetilde{v}_{mn}\rangle$:

[TABLE]

therefore one finds $\langle\zeta_{11}\widetilde{v}_{11}\rangle=3\langle\zeta_{11}\widetilde{v}_{22}\rangle=3\langle\zeta_{12}\widetilde{v}_{12}\rangle=...=\sigma_{\times}^{2}/5\sigma_{\delta}$;

• •

$\langle\widetilde{v}_{ij}\widetilde{v}_{mn}\rangle$:

[TABLE]

therefore one finds $\langle\widetilde{v}_{11}^{2}\rangle=3\langle\widetilde{v}_{11}\widetilde{v}_{22}\rangle=3\langle\widetilde{v}_{12}^{2}\rangle=...=1/5$.