---
author:
- 'Siri Chongchitnan$^{a,b}$'
- 'Teeraparb Chantavat$^{c,d}$'
- 'and Jenna Zunder$^{e}$'
bibliography:
- 'pbh\_ref.bib'
title: Extreme Primordial Black Holes
---

Introduction
============

Primordial black holes (PBHs) originate from large inflationary perturbations that subsequently collapse into black holes in the early Universe (for reviews, see [@sasaki; @garcia; @carr; @kashlinsky]). LIGO gravitational wave events[^1] over the past few years have given rise to the resurgence of PBHs not only as a viable dark matter candidate, but also as potential pregenitors of massive black holes ($\gtrsim30M_\odot$) that can typically give rise to the observed amplitude of gravitational waves.

Given an inflationary scenario, it would be useful to predict the mass of the most massive PBHs expected within a given observational volume. Such a calculation would serve as an additional observational test of competing inflationary theories. The primary aim of this work is to present such a framework, whilst demonstrating the method for a particular model of PBH formation.

The framework discussed is based on previous work by one of us [@me_void1; @me_void2] in the context of extreme cosmic voids, as well as previous work by Harrison and Cole [@harrison; @harrison2] on extreme galaxy clusters. Our main result will be the probability density function (pdf) for the most massive PBHs expected in an observational volume. We will apply the framework to a simple model of PBH formation and demonstrate the soundness of the calculations.

Throughout this work we will use the cosmological parameters for the LCDM model from Planck [@planck18].

Quantifying PBH abundance
=========================

Fraction of the Universe in PBHs
--------------------------------

In this section we will derive an expression for the cosmological abundance of PBHs, namely =,where $\rho\sub{PBH}$ is the mean cosmic density in PBHs, and $\rho\sub{crit}$ is the critical density. Typically we will be interested in the abundance of PBHs within a certain mass range (say, $\Omega\sub{pbh}(>M)$,  the fraction of the Universe in PBHs of mass greater than $M$). The PBH abundance naturally depends on how primordial perturbations were generated (the shape of the primordial power spectrum of curvature perturbations), details of the PBH collapse mechanism (structure formation theory), and thermodynamical conditions during the radiation era when PBHs were formed. We will obtain an expression for $\Omega\sub{PBH}$ that depends on all these factors.

One viable approach to begin modelling the PBH abundance is to borrow and modify the structure-formation theory from the extended Press-Schechter (PS) theory [@press], but with a modification of the collapse threshold. This approach has been widely used in previous work to model PBH abundances ([@chong_pbh; @ballesteros; @wang; @byrnes; @young]). We present the key equations below.

In the PS formalism, the probability that a region within a window function of size $R$, containing mass $M$, has density contrast in the range $[\delta, \delta+d\delta]$ is given by the Gaussian distribution P()d=[1]{}[1]{}e\^[-\^2/2\^2]{},where $\sigma$ is the variance of the primordial density perturbations $\delta$ smoothed on scale $R$, \^2(R)=\_0\^W\^2(kR)\_(k) k.We choose the Fourier-space window function $W$ to be Gaussian[^2]: W(x)=e\^[-x\^2/2]{},and the power spectrum $\mc{P}_\delta(k)$ depends on the PBH formation mechanism.

Assuming that PBHs originate from Fourier modes that re-entered the Hubble radius shortly after inflation ends ( during radiation era, when $R$ becomes comparable to $k^{-1}=(aH)^{-1}$), $\sigma(k)$ can be expressed in terms of the primordial curvature power spectrum, $\mc{P}_{\mc{R}}(k)$, as [@liddle] \^2(k)=\_[-]{}\^q W\^2(qk\^[-1]{})(qk\^[-1]{})\^4 T\^2(q, k\^[-1]{}) \_(q),where the transfer function, $T(q,\tau)$, is given by: T(q, ) = [3y\^3]{}, y .

Numerical simulations suggest that the initial mass, $M$, of a PBH formed when density perturbation of wavenumber $k$ re-enters the Hubble radius, is known to be a fraction of the total mass, $M_H$, within the Hubble volume ($M_H$ is usually called the ‘horizon mass’). In this work, we follow recent literature in modelling $M$ as [@musco; @byrnes] M= K(-\_c)\^ M\_H.where we take $\delta_c=0.45$ (the threshold overdensity for collapse into a PBH during radiation era), with $K=3.3$ and $\gamma=0.36$.

For a given Hubble volume with horizon mass $M_H$, the corresponding temperature, $T$, satisfies the equation [@wang] M\_H = 12\^[3]{}\^[1/2]{}T\^[-2]{}, where $m\sub{Pl}$ is the Planck mass, and the effective degree of freedom $g_{*,\rho}(T)$, corresponding to energy density $\rho$, can be numerically obtained as described in [@saikawa]. The latter reference also gave the fitting function for the effective degree of freedom $g_{*,s}(T)$ corresponding to entropy $s$, which we will also need.

Using the extended PS formalism, one obtains the following expression for $\beta_{M_H}$, the fraction of PBHs within a Hubble volume containing mass $M_H$ [@niemeyer; @byrnes] \_[M\_H]{} &= 2\_[\_c]{}\^ P()\
&=\_[\_c]{}\^ B\_[M\_H]{}(M) M.The integrand $B_{M_H}(M)$ can be interpreted as the probability density function (pdf) for PBH masses on logarithmic scale at formation time. Using and , one finds B\_[M\_H]{}(M)&= [K(k\_H)]{}\^[1+1/]{},\
&,\
[k\_H\^[-1]{}]{}&=3.74510\^[6]{}\^[-1/2]{}\^[1/4]{}\^[-1/3]{}.

We next consider an important quantity $f(M)$, the fraction of dark matter in the form of PBHs of mass $M$. For our purposes, the expression for $f(M)$ can be obtained by integrating the pdf $B_{M_H}(M)$ over all logarithmic horizon masses, weighted by a thermodynamical factor $\tau(M_H)$ [@wang]: f(M) &\
&=[\_m]{}\_[-]{}\^ (M\_H) B\_[M\_H]{}(M) M\_H,\
(M\_H)&,where $\Omega\sub{CDM}$ and $\Omega_m$ are the cosmic density parameters for cold dark matter and total matter (CDM+baryons) respectively. $T\sub{eq}$ is the temperature at matter-radiation equality.

Once we have calculated the PBH fraction, $f(M)$, the total fraction of PBHs in dark matter can be calculated by integrating over all PBH masses, f=\_[M]{}\^f(M) M,(we will discuss $M\sub{min}$ in the next section). Finally, the fraction of the Universe in PBHs of mass $>M$ can simply be integrated as (&gt;M)=\_[M]{}\^f(M\^) M\^.

PBH number count
----------------

In analogy with the abundance of massive galaxy clusters (see [@mobook] for a pedagogical treatment), the differential number density of PBHs at present time ( the PBH ‘mass function’) can be expressed as: &= -[|M]{} =[|M]{}f(M), where $\bar\rho$ is the present-day mean cosmic density. In an observational volume covering the fraction $f\sub{sky}$ of the sky up to redshift $z$, we would find the total number of PBHs to be N(z)= f \^[z]{}\_[0]{} z\_[M(z)]{}\^M , where d$V/$d$z$ is the Hubble volume element given by &= [4H(z)]{}\^2,\
H(z)&H\_0\^[1/2]{}, where the cosmic densities $\Omega_i$ have their usual meaning. In this work, we will assume that $f\sub{sky}=1$.

It remains to discuss the redshift dependence in the integral . For the mass function, we introduce redshift dependence by the replacement $\sigma\to D(z)\sigma$ where $D(z)$ is the growth function normalized so that $D(0)=1$ (effectively this means that PBHs grow like dark matter). This technique has been used elsewhere ([@murgia]) to extrapolate to higher redshifts.

The lower bound in the mass integral in is the minimum PBH mass (at formation time) below which a PBH would have evaporated by redshift $z$. For $z=0$, it is well known that [@hawking], M(z=0)=5.110\^[14]{}2.610\^[-19]{} M\_.At higher redshifts, the minimum initial mass can be estimated assuming some basic properties of black holes. We outline the calculations in the Appendix. We found that $M\sub{min}$ remains within the same order of magnitude for a wide range of redshift (see Fig. \[Mmin\]). Therefore, for models which generate an observationally interesting abundance of PBHs, it is sufficient to make the approximation $M\sub{min}(z)\approx M\sub{min}(0)$ in Eq. \[Nto\]. We have checked that this makes no numerical difference for the models studied in this work.

PBH formation: a case study
===========================

The log-$\delta$ model
----------------------

It is well known that the simplest models of single-field slow-roll inflation cannot produce observable abundance of PBHs unless the primordial power spectrum is very blue (although this has firmly been ruled out by CMB constraints). Viable inflation models which generate an interesting density of PBHs are potentials that typically produce sharp features in the primordial power spectrum, so as to generate power at small scales [@kawasaki; @garcia2; @drees]. In this work, we will represent a generic primordial power spectrum with a sharp feature using a delta function spike in $\ln k$, (k)= A\_D(k - k\_0),where $\delta_D$ is the Dirac delta function. The constants $A$ and $k_0$ parametrize the amplitude and location of the spike in the resulting matter power spectrum. This log $\delta$-function model was previously studied in Wang [@wang] in the context of gravitational wave production by PBHs.

Observational constraints
-------------------------

A range of observational constraints, including CMB anisotropies and microlensing observations, have placed upper bounds on $f(M)$, the PBH fraction in CDM (see, for example, [@carr2]). Nevertheless, the published bounds assume that all PBHs have the same mass. These so-called constraints on $f(M)$ were traditionally the main quantity of interest in the literature, as there is a wide range of observational techniques that can place upper bounds on $f(M)$ over several decades of $M$.

If we now assume that PBHs are formed across a spectrum of masses, the monochromatic upper bounds, denoted $f\super{mono}\sub{max}(M)$, must be corrected using procedures such as those previously presented in [@carr2; @kuhnel; @azhar; @lehmann]. These studies have only relatively recently gained traction, but are nevertheless indispensable if PBHs were to be taken as a serious candidate for dark matter and GW sources.

![The colour gradient shows the theoretical values of the PBH to CDM ratio, $f\sub{PBH}$ (Eq. \[fint\]) as a function of the parameters $A$ and $M_H$ in the log-$\delta$ model. The white line indicates the observational upper bound $f\sub{PBH,max}$ (Eq. \[correct\]) converted from monochromatic constraints in the literature.[]{data-label="fpbh"}](fpbh_max_2.pdf){width="5in"}

The upshot from these studies is that the corrected upper bound for the total PBH fraction in CDM, $f\sub{PBH,max}$, is given by [@azhar] f= \^[-1]{}.The result from applying this correction to monochromatic constraints on the log-$\delta$ model is shown in Fig. \[fpbh\]. The figure shows the colour-coded magnitude of the PBH fraction $f\sub{PBH}$ (Eq. \[fint\]), as a function of model parameter $A$ (vertical axis) and $k_0$ (horizontal axis, converted to the corresponding horizon mass through Eq. ). The white line shows the corrected upper bound $f\sub{PBH,max}$. In other words, the region below the white line is the allowed parameter space for the log-$\delta$ model given current observations[^3].

The upper bound is increasing in the domain shown, until $M_H\sim10^2M_\odot$, where the dip corresponds to the more stringent constraint from the CMB anisotropies, since black-hole accretion effects can significantly alter the ionization and thermal history of the Universe [@ricotti].

Another interesting observation from the figure is the values of $f\sub{pbh}$ along the white line. The maximum occurs when the spike is at $M_H=10^{-8}M_\odot$, with $f\sub{PBH}\approx0.46$, and the minimum at $M_H=10^3M_\odot$, with $f\sub{PBH}\approx1.6\times10^{-3}$. This means that present constraints allows the log-$\delta$ model to consolidate almost half of all dark matter into PBHs. However, this comes from imposing a spike at very small scales where nonlinear effects (as seen in numerical simulations [@hutsi; @inman]) have not been incorporated into our calculations.

Extreme PBHs
============

Having established a method to calculate the PBH number count and mass function, we now set out to derive the probability distribution of the most massive PBHs expected in an observational volume. Our calculation is based on the exact extreme-value formalism previously used in the context of massive galaxy clusters [@harrison; @harrison2] and cosmic voids [@me_void1; @me_void2]. We summarise the key concepts and equations below.

Exact extreme-value formalism
-----------------------------

From the PBH number count , we can construct the probability density function (pdf) for the mass distribution of PBHs with mass in the interval $[\log M,\log M+\D\log M]$ within the redshift range $[0,z]$ as f\_[&lt;z]{}(M)= [fN]{} \^[z]{}\_[0]{} z . To verify that this function behaves like a pdf, one can see that clearly $f_{<z}(M)\in[0,1]$ for $M\in[M\sub{min},\infty)$, and, by comparing and , we also have the required normalization $$\int_{-\infty}^{\infty} f_{<z}(M)\,\D\log M=1.$$

The cumulative probability distribution (cdf), $F(M)$, can then be constructed by integrating the pdf as usual: F(M)=\_[M]{}\^[M]{} f\_[&lt;z]{}(m) m.This gives the probability that an observed PBH has mass $\leq M$.

Now consider $N$ observations of PBHs drawn from a probability distribution with cdf $F(M)$. We can ask: what is the probability that the PBH observed will have mass $M^*$? The required probability, $\Phi$, is simply the product of the cdfs: (M\^\*,N)=\_[i=1]{}\^N F\_i(MM\^\*)=F\^N(M\^\*),assuming that PBH masses are independent, identically distributed variables. As $\Phi$ is another cdf, the pdf of PBH can be obtained by differentiation: (M\^\*,N)=F\^N(M\^\*)=Nf(M\^\*)\[F(M\^\*)\]\^[N-1]{}.It is also useful to note that the peak of the extreme-value pdf (the turning point of $\phi$) is attained at the zero of the function X(M) = (N-1)f\^2 +F ,as can be seen by setting $\text{d}\phi/\text{d}\log M^*=0$.

In summary, starting with the PBH mass function, one can derive the extreme-value pdf for PBHs using Eq. \[eev\].

Application to the log-$\delta$ model
-------------------------------------

![The extreme-value probability density function for PBHs assuming the log-$\delta$ model with spike at $M_H=10^2 M_\odot$, assuming that $N=10$, $10^2$, $10^3$ observations. []{data-label="fpbh"}](phi.pdf){width="4.5in"}

Figure \[fpbh\] shows the pdfs of extreme-mass PBHs given for $N=10^2 ,10^3$ and $10^4$ observations (this figure summarises the key results of this work.) We assume the log-$\delta$ model with the power-spectrum spike at $M_H=10^2M_\odot$. The pdfs are not symmetric but have a positive skewness, consistent with previous derivations of extreme-value pdfs [@me_void1; @me_void2]. As $N$ increases, the peaks of the pdf naturally shifts towards higher values of $M^*$, with increasing kurtosis ( more sharply peaked).

![Profile of the extreme-value pdf peaks (solid line) and the 5th/95th percentiles (dashed lines) for the log-$\delta$ model as the spike location is varied (horizontal axis), assuming 100 observations of black holes. The vertical axis shows the location of the peaks. The relationship is linear to a good approximation (Eq. \[lin\]).[]{data-label="band"}](band.pdf){width="3.5in"}

When we vary the location of the spike (whilst keeping $N$ fixed, and using values of $A$ that saturate the upper bound shown in Fig. \[fpbh\]), we obtain an almost linear variation as shown in Fig. \[band\] (in which $N=100$). Each vertical slice of this figure can be regarded as the profile of the extreme-value pdf, with the peak of the pdf being along the solid line, whilst the 5th and 95th percentiles are shown in dashed lines. The band is linear to a good approximation, with the peak $M^*\sub{peak}$ satisfying the relation M\^\*2.3M\_H.Thus, given the assumptions of the log-$\delta$ model, the current observational constraints on $f\sub{pbh}$ do not rule out massive PBHs associated with LIGO events. Early-universe mechanisms that produce a spike in the primordial power spectrum at $M_H\sim10M_\odot$ can, in principle, produce PBHs of mass within the same order of magnitude.

It is also interesting to consider how tightening observation bounds will affect the extreme-value pdfs. Fig. \[shift\] shows what happens in this situation in the model with $M_H=10^3M_\odot$ (with $N=100$), supposing that the upper bound on $f\sub{pbh}$ is tightened to 50% of the current values (a realistic prospects from CMB experiments such as Euclid [@habouzit]). We see that, in line with expectation, the pdf shifts to smaller masses by $\sim20\%$, whilst the distance between the 5% and 95% percentiles shrinks by $\sim30\%$.

![Extreme-value pdfs for the log-$\delta$ model with $M_H=10^3M_\odot$ given current constraints (thick line) and futuristic constraints. The pair of vertical dotted lines on each pdf indicate the 5th and 95th percentiles. If observational constraints on $f\sub{pbh}$ were tightened to 50% the current values, the peak of the extreme-value pdf would shift downwards by $20\%$, whilst the inter-percentile distance would shrink by $\sim30\%$. []{data-label="shift"}](fig_shift.pdf){width="4in"}

Conclusions and discussion
==========================

In this work, we have established a framework to calculate the mass distribution of most massive PBHs expected within a given observational volume. The calculations were based mainly on four main ingredients:

the PBH formation mechanism (details of inflation or the shape of $\mc{P_R}(k)$)

the abundance of massive objects (modified Press-Schechter theory)

the exact extreme-value formalism

the observational constraints on $f\sub{pbh}$ (the PBH fraction in CDM).

We applied our formalism to the log-$\delta$ model, a prototype of models with a spike in the power spectrum. They are generically associated with inflationary models that produce interesting densities of PBH (via a phase transition in the early Universe). Our main results are the extreme-value pdf shown in Fig. \[fpbh\]. The fact that the location of the power-spectrum spike is close to the peak of the resulting extreme-value pdf gives assurance that our calculations are sound, and can thus be applied to many inflationary models known to produce PBHs. In future work, we will present a survey of extreme-value pdfs for a range of inflationary scenarios.

Some avenues for further investigation include studying the effect of changing the mass function (for example, extending the Sheth-Tormen mass function to PBHs [@me_pbh]), as well as understanding the role of PBH clustering and merger [@kohri; @young2; @tada] which will serve to strengthen the validity of the extreme-value formalism presented here.

Appendix: The minimum initial mass of an unevaporated black hole at redshift $z$. {#appendix-the-minimum-initial-mass-of-an-unevaporated-black-hole-at-redshift-z. .unnumbered}
=================================================================================

Consider a Schwarzchild black hole. Its decay rate depends on three variables, namely, 1) the spin ($s$) of the particles it decays into, 2) the energy ($E$) of those particles, and 3) the instantaneous mass $(M)$ of the black hole. By summing over all the emitted particles, the decay rate of a black hole can be expressed as $$\label{pintegral}
    \frac{{\rm d}M}{{\rm d}t} = -\frac{1}{2\pi \hslash c^2} \sum_{j} \Gamma_j \int {\rm d}E \frac{E}{\exp\big(8\pi G E M/\hslash c^3 \big) - (-1)^{2s_j}},\tag{A1}$$ where the sum is taken over all emitted particle species. The integral is taken over $(0,\infty)$ for massless particles, or $(\mu_j,\infty)$ for massive particles with rest energy $\mu_j$. $\Gamma_j$ is the dimensionless absorption probability, and $s_j$ is the spin of the $j$th species.

MacGibbon [@MacGibbon1991] showed that Eq. (\[pintegral\]) can be written as $$\label{macgibbon1}
    \frac{{\rm d}M}{{\rm d}t} = -5.34\times10^{22} f\sub{emit}(M) M^{-2} \text{ kg s}^{-1}.\tag{A2}$$ The function $f\sub{emit}(M)$ is given in a rather complicated piecewise form in Eq. (7) in [@MacGibbon1991]. Eq. (\[macgibbon1\]) can be inverted and integrated to yield the evaporation timescale, $\tau_{\rm evap}$, as $$\tau_{\rm evap} = \big(1.87266\times10^{-23}\ \text{s kg}^{-1} \big) \int_{M_f}^{M_i} {\rm d}M f(M)^{-1} M^2,\tag{A3}$$ where $M_i$ is the initial mass of the black hole and $M_{f}$ is the final mass. Letting $M_f = 0$ and $M_i = M_*$, we can obtain $M_*$ as a function of $z$ by solving the non-linear equation: $$\tau_{\rm evap}(z)\bigg|_{M_f = 0, M_i = M_*} = t_{\rm univ}(z),\tag{A4}$$ where $t_{\rm univ}$ is the age of the universe at redshift $z$. The minimum initial black hole mass as a function of redshift is shown in Figure \[fig:mmin\]. The figure closely resembles Fig. 1 of [@MacGibbon1991], although we believe that the labelling of the two curves in that figure should be exchanged.

![\[fig:mmin\] Minimum initial mass of a black hole which have not evaporated, observed at redshift $z$. The solid line is the result assuming the LCDM universe with Planck 2018 parameters. For comparison, the dashed line assumes the Einstein-de Sitter (EdS) cosmology $(\Omega_m=1, \Omega_\Lambda=0)$.[]{data-label="Mmin"}](Mmin.pdf){width="4in"}

[^1]: <https://www.ligo.org/detections.php>

[^2]: See [@ando] for an interesting study of how window functions affect the inferred PBH abundances.

[^3]: We use observational constraints summarised in Fig 3 of [@carr], and interpolated the upper bounds to obtain an approximate functional form for $f\super{mono}\sub{max}(M)$.
