A possible mass distribution of primordial black holes implied by LIGO-Virgo
AA possible mass distribution of primordial black holes implied byLIGO-Virgo
Heling Deng ∗ Physics Department, Arizona State University, Tempe, AZ 85287, USA
Abstract
The LIGO-Virgo Collaboration has so far detected around 90 black holes, some of which havemasses larger than what were expected from the collapse of stars. The mass distribution of LIGO-Virgo black holes appears to have a peak at ∼ M (cid:12) and two tails on the ends. By assuming thatthey all have a primordial origin, we analyze the GWTC-1 (O1&O2) and GWTC-2 (O3a) datasetsby performing maximum likelihood estimation on a broken power law mass function f ( m ), withthe result f ∝ m . for m < M (cid:12) and f ∝ m − for m > M (cid:12) . This appears to behave betterthan the popular log-normal mass function. Surprisingly, such a simple and unique distributioncan be realized in our previously proposed mechanism of PBH formation, where the black holesare formed by vacuum bubbles that nucleate during inflation via quantum tunneling. Moreover,this mass distribution can also provide an explanation to supermassive black holes formed at highredshifts. ∗ Electronic address: [email protected] a r X i v : . [ a s t r o - ph . C O ] F e b . INTRODUCTION Over the past few years, the LIGO-Virgo Collaboration has detected gravitational wavesemitted by about 50 inspiraling and merging black hole binaries [1, 2], and many more eventsare anticipated in the near future. It is intriguing that most of the black holes have masses ∼ M (cid:12) , which certainly provides indicative information of the mass distribution of blackholes in the universe.The origin of these black holes is so far unknown. While the possibility that they areordinary astrophysical black holes from stellar collapse (possibly from different channels) areunder active investigations [3–6], a fascinating speculation is that LIGO-Virgo has detectedprimordial black holes (PBHs) [7–9]. PBHs are hypothetical black holes formed by deviationsfrom density homogeneity in the early universe, before any large scale structures and galaxies.Their masses can range from the Planck mass ( ∼ − g) to many orders of magnitude largerthan the solar mass ( M (cid:12) ∼ g). After the release of LIGO-Virgo results, much effort hasbeen dedicated to constrain PBHs and the role they play in dark matter [9–13]. It is nowgenerally recognized that they are unlikely to constitute all dark matter due to the smallmerger rate inferred by LIGO-Virgo.It is also of great interest to constrain particular PBH mass distributions with the LIGO-Virgo results [10, 14–17]. Among them the log-normal mass function [18] is the most popularbecause it is a good approximation for various PBH mechanisms [19–21]. However, if allLIGO-Virgo black holes were PBHs, the peak of the log-normal function should be at ∼ M (cid:12) , which is incompatible with what was observed. The possibility of LIGO-Virgo blackholes being a mixture of two populations of black holes was recently considered in ref. [22],which concludes that astrophysical black holes that dominate the mass range m < ∼ M (cid:12) together with PBHs given by a log-normal (or critical collapse) mass function are much moreprobable than PBHs only.In this paper, inspired by ref. [22], we perform a maximum likelihood estimation tofit the LIGO-Virgo datasets with a broken power law only. More sophisticated techniques(e.g., Bayesian analysis [23, 24]) are available, but we believe likelihood analysis shouldsuffice for our purposes, given the limited number of detected events and uncertainties inPBH formation. We shall also neglect some details in analyzing the formation of PBHs andbinaries, as long as they are not expected to affect the results by orders of magnitude. For2nstance, following ref. [22], we shall not consider the impact from black hole’s spin andmass accretion.Apart from the single peak at ∼ M (cid:12) in the distribution of LIGO-Virgo black holes,our investigation of the broken power law mass spectrum is mainly motivated by a PBHmechanism we proposed in recent years. In a series of works [25–27], we studied PBHs formedby vacuum bubbles that possibly nucleate during inflation. Depending on its size afterinflation ends, a bubble will turn into a black hole in the either subcritical or supercriticalregime, and the mass distributions of black holes in these two regimes could obey differentpower laws. We would like to know if the best-fit parameters from likelihood analysis ofLIGO-Virgo data is compatible with this mechanism. Surprisingly, the broken power lawmass function implied by the LIGO-Virgo black holes can indeed be realized in our model.The rest of the paper is organized as follows. In section II, we will first discuss themerger rate of PBHs and the detection probability of LIGO-Virgo, and then apply maximumlikelihood estimation to find the best-fit parameters for the broken power law mass function.In section III, we will study the mechanism of PBH formation from primordial bubblesand how could produce the LIGO-Virgo black holes. Conclusions will be summarized anddiscussed in section IV. We set c = (cid:126) = G = 1 throughout the paper. II. PBHS AND LIGO-VIRGO EVENTS
It will be assumed that all black holes detected by LIGO-Virgo are primordial. They areformed during the radiation era and behave like dark matter, with their number and massdensity diluted by Hubble expansion. Two neighboring black holes may collide as theirgravitational attraction defeats the Hubble stretching before the radiation-dust equality.Disturbance from the surrounding environment, a typical example being a third nearby blackhole exerting a tidal torque, may impede the head-on collision, leading to the formation ofan inspiraling binary. If the coalescence time of the binary is comparable to the age of theuniverse, gravitational waves from the merger could possibly be detected when they reachthe earth. Roughly speaking, whether a merger event can be heard depends on the detector’ssensitivity, the masses of the two black holes (source masses), the time when the merger tookplace (source redshift), and the sky location and orientation of the binary system. In orderto estimate how often a merger event can be recorded, we also need the merger rate of the3inaries, which is determined by the underlying mechanism of PBH formation.
A. PBH mass function and merger rate
We characterize the mass distribution of PBHs by mass function ψ ( m ), defined by ψ ( m ) = mρ CDM d n d m . (1)Here m is the black hole mass, d n is the number density of black holes within the massrange ( m, m + d m ), and ρ CDM is the energy density of cold dark matter. Since black holesand dark matter are diluted by cosmic expansion in the same way, ψ ( m ) is a constant overtime. Integrating ψ ( m ) gives the total fraction of dark matter in PBHs: f P BH ≡ ρ P BH ρ CDM = (cid:90) ψ ( m )d m, (2)where ρ P BH = (cid:82) m d n is the energy density of PBHs. If all dark matter is made of PBHs,we have f P BH = 1.Another function often used to describe the PBH mass distribution is f ( m ) ≡ mψ ( m ) , (3)which can be interpreted as the fraction of dark matter in PBHs at m within the mass range∆ m ∼ m . This function is particularly useful if the spectrum is relatively broad.The PBH merger rate in the early universe is well studied in the literature [9–13]. In thispaper we take the formula of differential merger rate from ref. [14]:d R ≈ . × Gpc yr f − / P BH (cid:18) tt (cid:19) − / (cid:18) MM (cid:12) (cid:19) − / η − / S ( ψ, f P BH , M ) ψ ( m ) ψ ( m )d m d m . (4)Integrating d R over all masses gives the number of merger per Gpc per year. Here m and m are the two masses in the binary, M ≡ m + m and η ≡ m m /M ; t is the time whenthe merger occurs, and t is the present time; S is a suppression factor accounting for theeffects from other matter components, including nearby black holes, and S = O (1) for anot particularly wide spectrum. In our calculations we shall approximate S by a simplifiedexpression given in ref. [22], but setting S = 1 would not lead to much difference.To roughly estimate the largest f P BH we can have from the LIGO-Virgo results, we canassume a monochromatic mass function with ψ ( m ) = f P BH δ ( m − M (cid:12) ), i.e., all PBHs have4ass 30 M (cid:12) . Setting t = t and S = 1, we have (cid:90) d R ∼ f / P BH
Gpc − yr − . (5)The merger rate implied by LIGO-Virgo is O (10 − − yr − , which gives f P BH ∼ − .Therefore, if all PBHs are of masses around 30 M (cid:12) , they can at most contribute to 0 .
1% ofthe dark matter.The distribution of LIGO-Virgo black holes is clearly not monochromatic. In this paperwe shall consider the 10 black hole merger events from the GWTC-1 catalog [1] and 34 eventsfrom GWTC-2 (discarding GW190719 and GW190909) [2], so there are 88 black holes inour dataset. In fig. 1, we depict the black hole numbers in different mass ranges. Thedistribution appears to have a peak at m ∗ ∼ M (cid:12) and two tails on the ends, with moreblack holes having masses less than m ∗ than those on the other end. The simplest massfunction that can be conceived is a broken power law, ψ ( m ) = f P BH m ∗ (cid:0) α − − α − (cid:1) − ( m/m ∗ ) α − , m < m ∗ ( m/m ∗ ) α − , m > m ∗ , (6)which satisfies (cid:82) ψ d m = f P BH . We further assume α > α < f ( m ) = mψ ( m ) has a peak at m ∗ , which means most contribution to the PBH energy density comesfrom black holes with masses around m ∗ . In the next subsection, we will analyze the LIGO-Virgo results and find the best-fit parameters for (6).
B. Maximum likelihood estimation
Given the differential merger rate found by eq. (4), the expected number of merger signalsreaching the earth per unit time within the ranges ( m , m + d m ), ( m , m + d m ), and( z, z + d z ) can be evaluated as [28]d N ( m , m , z ) = 11 + z d V c d z d R d m d m d m d m d z, (7)where (1 + z ) − accounts for the time difference between the source frame and the detectorframe, and V c is the comoving Hubble volume. For later convenience, we defineΛ( m , m , z ) ≡
11 + z d V c d z d R d m d m (8)5
20 40 60 80 m / M d ln( N )/ dm FIG. 1: Numbers of detected LIGO-Virgo black holes in different mass ranges, rescaled such thatthe areas add up to 1. There are 88 black holes in total. The black hole masses are simply takento be the median values reported in refs. [1, 2].
In addition, we ought to take into account the fact that not all merger events can beobserved by LIGO-Virgo. The detection probability of an event depends on the sensitivityof the instruments, the waveform of the signal as well as the extrinsic parameters of thebinary system, i.e., its sky location and orientation. Integrating out the extrinsic parametersfor a given detector gives its detection probability p det ( r ) , with r ≡ ρ c /ρ ( m , m , z ). Here ρ c is the threshold signal-to-noise ratio above which a signal can be detected, usually taken as ρ c = 8; and ρ is the signal-to-noise ratio for a merger located directly above the detector. ρ is defined by ρ ( m , m , z ) = 2 (cid:115)(cid:90) | ˜ h ( f ) | S n ( f ) d f , (9)where ˜ h ( f ) is the Fourier transform of the signal, and S n ( f ) is the power spectrum of thedetector’s strain noise. In this work, ˜ h ( f ) and S n ( f ) are taken from refs. [22, 29] and refs.[1, 2] respectively. Then the detection probability of a particular event can be found by theparameter fit of p det ( r ) given in the appendix of ref. [28]. Since S n ( f ) is unique to a certainrun, p det ( r ) takes different values in, e.g., O1&O2 and O3a of LIGO-Virgo.Now the total expected number of detection during a certain run can be found by N e = T o (cid:90) p det ( m , m , z )d N, (10)6here T o is the observation time interval (166.6 days for O1&O2 and 183.3 days for O3a).Therefore, once we are given the PBH mass function and the detector’s sensitivities ( S n ( f )),we are able to obtain from eq. (10) the expected number of observable events.We assume that (observable) mergers occur randomly with a constant rate, so they reachthe detectors following a Poisson process. Let N o be the actual number of observed events.Then the probability of this happening during the observation interval is p P oisson ∝ N N o e e − N e . (11)The likelihood of a single detected event with source masses m ,o , m ,o and redshift z o is p o ∝ p det ( m ,o , m ,o , z o )Λ( m ,o , m ,o , z o ) (cid:82) p det ( m , m , z )d N . (12)Here we have made an assumption that the source masses and redshift can be perfectlydetermined by detection. In practice, we take m ,o , m ,o and z o to be the median valuesof the posterior samples reported by LIGO-Virgo [1, 2]. In reality, their values are givenwith probability distributions resulting from measurement uncertainties. Note that thedenominator in eq. (12) is ∝ N e . Then the likelihood of all black hole merger eventsdetected by LIGO-Virgo is L = p P oisson N o (cid:89) i =1 p io ∝ e − N e N o (cid:89) i =1 p det ( m i ,o , m i ,o , z io )Λ( m i ,o , m i ,o , z io ) , (13)where each event is denoted by a superscript i . By maximizing the log-likelihood ln L withthe GWTC-1 and GWTC-2 datasets, we will be able to attain the constraints on a PBHmass function.There are four parameters in a broken power law (6): α , α , m ∗ and f P BH . Scanningthrough the relevant part of the parameter space, we found the largest log-likelihood ln L max at α ≈ . , α ≈ − , m ∗ ≈ M (cid:12) , f P BH ≈ . . (14)These results are also partly shown in fig. 2. The power law mass function from differentPBH mechanisms was discussed in, e.g., refs.[26, 31, 32], which suggested f ( m ) ∝ m − / (for large masses). These models are obviously disfavored by LIGO-Virgo. In order to account for the measurement uncertainties, the standard approach to calculate p o is to replacethe numerator in eq. (12) by (cid:104) Λ( m ,o , m ,o , z o ) /π ( m ,o , m ,o , z o ) (cid:105) , where m ,o , m ,o and z o are the poste-rior samples, π is the corresponding source prior, and the brackets denote an average over all samples ofthat event [23, 30]. .2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 m * m * FIG. 2: Three 2-D slices of the 4-D parameter space with f P BH = 0 . L max , and the shaded regions are where the parameter valuesgive likelihoods L with ln ( L max / L ) ≤ . Left: α ≈ . Middle: α ≈ − Right: m ∗ ≈ M (cid:12) . The unit for m ∗ is M (cid:12) . As a comparison, we did the same analysis for the log-normal mass function ψ ( m ) = f P BH m √ πσ exp (cid:20) − ln ( m/m c )2 σ (cid:21) . (15)The most probable parameters are m c ≈ M (cid:12) , σ ≈ . , f P BH ≈ . , (16)which are consistent with the results found in refs. [10, 14, 22, 24]. The difference betweenthe maximum log-likelihood of broken power law and that of log-normal is ln( L BP L / L LN ) ≈ . However, considering that we are evaluating the likelihood ratio with two different models,and that the broken power law has one more free parameter, we cannot say with muchconfidence that the broken power law is a better model at the moment in explaining theLIGO-Virgo results.Nevertheless, using the Metropolis–Hastings algorithm, we have drawn 50000 randomsamples from the probability distribution p det ( m , m , z )d N with the best-fit parameters,both for the broken power law and the log-normal model, and compared the resulting massdistributions with that from LIGO-Virgo black holes. As we can see from fig. 3, although itis difficult to compare models under likelihood analysis, the broken power law mass functionappears to be a better fit than the log-normal.In the next section we will discuss a physical mechanism where PBHs could form withsuch a simple distribution. 8
20 40 60 80 m / M d ln( N )/ dm FIG. 3: Mass distributions of black holes from LIGO-Virgo (blue), the best-fit broken power law(orange) and log-normal (green).
III. PBHS FROM PRIMORDIAL BUBBLES
In a series of works [25–27], we studied in detail a mechanism of PBH formation whereblack holes are formed by vacuum bubbles that nucleate during inflation. After inflationends, a bubble will eventually turn into a black hole in the either subcritical or supercriticalregime. PBHs formed in the two regimes might obey different power laws, and the transitionregion could just be the peak at ∼ M (cid:12) seen by LIGO-Virgo. We shall first briefly describehow these black holes are formed, and then discuss how they are constrained by observationsincluding LIGO-Virgo. A. PBH formation
Vacuum bubbles could nucleate during inflation via quantum tunneling if there is a (pos-itive) “true” vacuum near the inflationary (quasi-de Sitter) vacuum in a multidimensionalfield potential. The bubble interior has an energy density ρ b smaller than the inflationaryenergy density ρ i . A typical bubble expands rapidly, almost at the speed of light. Afterinflation ends, inflatons outside the bubble turn into radiation with density ρ r ≈ ρ i , and thebubble will run into the surrounding radiation fluid with a huge Lorentz factor. The bub-ble continues to grow, but will eventually come to a halt with respect to the Hubble flow,9ecause all the forces acting on the bubble wall, including the interior vacuum pressure,exterior radiation pressure, and the wall’s surface tension, point inwards. A bubble is called“subcritical” if it simply collapses into a black hole after reentering the Hubble horizon.A supercritical bubble, on the other hand, will inflate without bound due to the repulsivevacuum inside. However, instead of consuming our universe, the bubble grows into a babyuniverse, which is connected to us by a wormhole. For an exterior observer, the wormholeis a spherical object that will eventually turn into a black hole as the “throat” pinches off.Once the relation between the black hole mass in these two regimes and the bubble’s size atthe end of inflation is determined, we are able to obtain the mass spectrum of these blackholes from the bubbles’ size distribution.In ref. [26], we considered an ideal scenario where radiation outside the bubble canbe completely reflected by the bubble wall, which implies strong interactions between thebubble field and the standard model particles. In this setting, the mass function was foundto be f ∝ m − / in the supercritical regime. As we know from the previous section, such adistribution is disfavored by the LIGO-Virgo results.In ref. [27], we studied the other extreme possibility, where the bubble wall is transparent,and radiation can freely flow inside. In this case, the trajectory of the bubble wall before itceases to grow with respect to the Hubble flow can be estimated by assuming an FRW metricdominated by radiation outside the bubble and matching the spacetimes on two sides of thewall. Let r be the the bubble’s comoving radius. The equation of motion of the bubble wallcan be found to be¨ r + (cid:0) − a ˙ r (cid:1) H ˙ r + 2 a r (cid:0) − a ˙ r (cid:1) + (cid:16) ρ b σ + 6 πσ (cid:17) (1 − a ˙ r ) / a = 0 , (17)where the overdot represents the derivative with respect to cosmic time t , H ≡ ˙ a/a = (2 t ) − is the Hubble parameter in the exterior FRW universe, and σ is the surface tension of thebubble wall. Let t i be the time when inflation ends, then the scale factor is defined by a = ( t/t i ) / . The “initial” conditions of eq. (17) is r ( t i ) = r i and ˙ r ( t i ) = (cid:0) − γ − i (cid:1) / ,where r i can be smaller or much larger than the Hubble horizon H − i at t i , and γ i (regardedas a free parameter ) is the Lorentz factor of the bubble wall with respect to the Hubbleflow at t i . See discussion in section IV. γ i , the trajectory of the wall can be approximated by a ˙ r ≈
1. Let t s be the time when the wall comes to a stop with respect to the Hubble flow. The bubble’scomoving radius at t s can then be estimated as r s ≡ r ( t s ) ∼ r i + ( a s − H − i , (18)where a s ≡ ( t s /t i ) / . We are particularly interested in the case where γ i is so large that r s (cid:29) r i . In this case, r s ∼ a s H − i , and hence the bubble’s physical radius at t s can beestimated as R s ≡ a s r s ∼ a s H − i ∼ t s , (19)which means the bubble size is comparable to the Hubble horizon at t s .
1. subcritical bubble
A bubble with sufficiently small r i would continue to expand after t s , until it reachesa maximum physical size R max . It will then shrink and eventually collapse into a blackhole. By eq. (19), the bubble reenters the horizon at ∼ t s , then its evolution will nolonger be affected by Hubble expansion significantly, which means the mass of the bubbleitself (excluding radiation inside) will almost conserve after then. When the bubble radiusreaches R max , the kinetic energy of the bubble wall vanishes. Furthermore, since the bubbleis expected to shrink upon horizon reentry, we have R max ∼ R s . Hence, the mass of theresulting black hole from a subcritical bubble is m ∼ πρ b R s + 4 πσR s − π σ R s . (20)Here the three terms can be understood as the volume energy, the surface energy and thesurface binding energy, respectively.As r i is increased to a critical value, R max can no longer be reached, because a sufficientlylarge bubble dominated by its interior vacuum would inflate. When this happens, we enterthe supercritical regime.
2. supercritical bubble
Due to the third term on the left hand side in eq. (17), a bubble with a larger r i experiences a smaller “friction” as it grows, and therefore has a larger t s . By the time t s ,11he bubble’s vacuum density ρ b could become even larger than the exterior radiation density ρ r ( t s ) (note that ρ r ( t i ) = ρ i > ρ b ), i.e., ρ r ( t s ) = ρ i a − s < ∼ ρ b . (21)Then by eq. (19), R s ∼ a s H − i = (cid:18) π ρ i a − s (cid:19) − / > ∼ (cid:18) π ρ b (cid:19) − / = H − b . (22)Here H − b is the de Sitter horizon associated with the bubble vacuum. R s > ∼ H − b thenimplies that the bubble would inflate before t s and thus create a wormhole, which laterturns into a black hole. By the estimate in ref. [27], the black hole mass in this supercriticalregime is m ∼ H i (cid:0) r s − c s a s H − i (cid:1) , (23)where c s is the speed of sound in the radiation fluid, and the value of a s is determined bysolving eq. (17). By the estimate in eq. (19), we have m ∼ H i r s ∼ a s H − i ≈ t s , which isthe horizon mass at t s .Therefore, by numerically integrating eq. (17) till ˙ r = 0, we are able to find out the blackhole masses in the subcritical and the supercritical regimes. These two regimes are expectedto be connected by a region where the black hole mass transits smoothly from one regimeto the other. In practice, we evaluate both (20) and (23), and take the smaller one to bethe black hole mass for a certain bubble. B. Size distribution and mass function
Now that we have the masses of black holes from bubbles with various r i , we need thedistribution of r i in order to obtain the mass function. Bubbles formed earlier expand tolarger sizes, but they are rarer due to cosmic expansion. By assuming that bubbles areformed with a vanishing size and that the bubble worldsheet is the future light cone of thenucleation point, one can find the number density of bubbles having radius in the interval( r i , r i + d r i ) at t ( t > t i ) to be [25]d n ( t ) ≈ λ d r i a ( t ) (cid:0) r i + H − i (cid:1) . (24)12here λ is the dimensionless bubble nucleation rate per Hubble volume per Hubble time.Then by the definition of f ( m ) (eqs. (1) and (3)), we have f ( m ) ∝ λm H / i (cid:0) r i + H − i (cid:1) d r i d m . (25)The relation between r i and m can be found by numerically solving eq. (17) and usingeqs. (20) and (23). Therefore, the mass function is completely determined by the followingfive parameters: the Lorentz factor of the bubble wall at the end of inflation γ i , the bubblenucleation rate λ , the inflationary density ρ i , the vacuum density of the bubble interior ρ b ,and the bubble wall tension σ . In the following we shall use η i ∼ ρ / i , η b ∼ ρ / b , η σ ∼ σ / , (26)to characterize ρ i , ρ b and σ . They represent the energy scales of inflation, bubble interiorand bubble wall, respectively.Depending on the parameter values, the mass function (25) can have very different shapesand can be wide or relatively narrow. Several examples of f ( m ) are shown in fig. 4, wherewe have fixed all other parameters except for the Lorentz factor γ i . An intriguing featurerelevant to our discussion is that when γ i is sufficiently large, black holes in the supercriticalregime near the transition approximately follow a power law f ( m ) ∝ m − (by a numericalfit to the curves) , which is what LIGO-Virgo implies for PBHs with m > m ∗ ∼ M (cid:12) . Moreover, there is a peak in the transition region for some values of γ i , and the massfunction can be approximated by different power laws on two sides of the peak. Now thequestion is whether suitable bubble parameters can be found such that subcritical blackholes obeys f ( m ) ∝ m . . C. Observational constraints
The answer is yes. In fact, there are many sets of bubble parameters that can give thedesired mass function. In fig. 5 we show an example of f ( m ) with γ i ≈ , λ ≈ − ,η i = O (1) TeV, η b = O (0 .
1) GeV and η σ = O (10) TeV (black, solid curve). It agrees wellwith the broken power law mass function with the parameters in eq. (14) (orange, dashed A semi-analytic discussion of this relation can be found in the appendix. For supercritical bubbles witha huge γ i , the mass function near the transition satisfies f ∝ m − . . IG. 4: Some examples of the mass function f ( m ) given by (25). We have fixed parameters exceptfor the Lorentz factor γ i . For small γ i the maximum of f ( m ) appears at the transition region. As γ i increases, another peak develops in the subcritical regime. For a sufficiently large γ i , supercriticalblack holes near the transition region approximately follow a power law f ∝ m − (the red straightline is ∝ m − ). For large masses, all spectra approach f ∝ m − / . The mass function in black (thethird one from the right) agrees well with the broken power law suggested by LIGO in the relevantmass range, as shown in fig. 5. lines). As a comparison, we also show the best-fit log-normal function (green, dotted curve).The light grey and colored areas are regions excluded by different astrophysical observations(see ref. [33] for a recent review). Strictly speaking, these upper bounds are valid only forPBHs with a very narrow spectrum, and so are not supposed to be used to compare with anextended mass function f ( m ). (The method of applying these bounds on a broad spectrumwas discussed in ref. [34].) Loosely speaking, however, as long as ψ ( m ) = f ( m ) /m does nothave a plateau over a relatively big range, we can still constrain a model by placing f ( m )near the bounds while avoiding hitting them.The most stringent bound for us in fig. 5, which covers the mass range around 30–300 M (cid:12) ,comes from Planck. Ref. [35] studied how disk or spherical accretion of a halo around PBHscould affect CMB, which strongly constrains the population of PBHs in the range around1–10 M (cid:12) . The broken power law function implied by LIGO-Virgo seems in tension withthe model of disk accretion (light blue), but is free from the bound imposed by spherical14 IG. 5: Observational constraints on our mechanism. The light grey, light blue and light purpleregions are excluded for PBHs with a monochromatic spectrum [35, 38–42] (adapted from fig. 10in ref. [33]). The black solid curve is a possible distribution from our PBH mechanism. The orangedashed lines are the best-fit broken power law, with f ∝ m . for m < M (cid:12) and f ∝ m − for m > M (cid:12) . The green dotted curve is the best-fit log-normal mass function. accretion (light purple). However, we note that the effect of PBH accretion was not takeninto account in our discussion. It is pointed out in refs. [36, 37] that efficient PBH accretionbefore the reionization epoch tends to relax the tension with these bounds.Besides being able to account for the LIGO-Virgo black holes, another interesting featureof our mass function is that it may provide an explanation to supermassive black holes(SMBHs) at the center of most galaxies [43, 44]. These black holes have masses from ∼ –10 M (cid:12) , which cannot be explained by standard accretion models [45]. Moreover,observations of quasars indicate that many of them were already in place at high redshifts(see ref. [46] for a review). Recently, a black hole as large as ∼ M (cid:12) was discovered at z ≈ .
6, which greatly challenges the conventional formation theory of astrophysical blackholes [47]. One is then led to speculate that SMBHs were seeded by PBHs, which could havelarge masses by birth [19, 48–51]. At the present time ( t ), the number density of PBHs ofmass ∼ m is approximately given by n ( m ) ∼ ρ CDM ( t ) f ( m ) m ≈ × (cid:18) M (cid:12) m (cid:19) f ( m ) Mpc − . (27)It was shown in refs. [35] that PBHs at z ∼ m ∼ –10 M (cid:12) can attain sufficient15ccretion and grow to SMBHs even if f ( m ) is as small as 10 − . We can see from fig.5 that the best-fit broken power law mass function gives f (10 M (cid:12) ) > ∼ − , which gives n (10 M (cid:12) ) ∼ . − . This is larger than the present density of galaxies ∼ . − .Therefore, SMBHs could indeed have been seeded by these PBHs. IV. CONCLUSIONS AND DISCUSSION
In this work we have used maximum likelihood estimation to analyze the GWTC-1 andGWTC-2 datasets released by the LIGO-Virgo Collaboration, assuming that all of the blackholes in the 44 confident merger events are PBHs with a mass spectrum f ( m ) satisfying abroken power law. We found that the best-fit f ( m ) has a peak at m ∗ ≈ M (cid:12) , consistentwith the observations, and that f ∝ m . for m < m ∗ and f ∝ m − for m > m ∗ . Theseblack holes can constitute about 0 .
1% of the dark matter.Such a unique mass function can be realized in a mechanism of PBH formation we pro-posed in recent years, where the black holes are formed by vacuum bubbles that couldnucleate during inflation. In general, the inflaton field rolls in a multidimensional potentialcontaining a number of minima. Bubble nucleation can occur during inflation if a minimumhas energy scale lower than the inflationary scale. After inflation ends, depending on theirsize, these bubbles will turn into black holes either by simple collapse (subcritical), or bycreating a wormhole (supercritical). The resulting mass function of the black holes mayobey different power laws in the two regimes, which are connect by a transition region thatcould possibly be at m ∗ . Surprisingly, if the bubble walls have a sufficiently large Lorentzfactor at the end of inflation, which is typically the case, the mass function in the super-critical regime near m ∗ is indeed approximately given by f ∝ m − . With properly chosenparameters (including the Lorentz factor, the energy scales of bubble vacuum, bubble walland inflation), we can have f ∝ m . for m < m ∗ .With this mass function, we can also have a sufficient number of seeds ( m ∼ M (cid:12) )that could grow into SMBHs observed at the galactic centers, which is difficult to explainwith conventional accretion models. In addition, as we can see in fig. 4, it is also possibleto estimate the distribution of stupendously large black holes [52] from the large mass end,where f ∝ m − / . This however will not be discussed further in the present work.Another noticeable feature of our mass function is that it suggests a relatively low in-16ationary scale. During the slow roll, the ultimate Lorentz factor of the bubble wall for acomoving observer outside the bubble is γ ∼ η i M P l /η ( i )3 σ , where η ( i ) σ is the energy scale ofthe wall tension σ during inflation [27, 53]. We assumed γ to be the Lorentz factor of thewall as it runs into the ambient radiation after inflation, i.e., γ i = γ . However, at the endof inflation, the inflaton quickly rolls down into our vacuum, which could cause a drasticchange in σ since the shape of the barrier in the field potential could change significantly.Therefore, we regarded γ i and σ as two free parameters. With the parameter values thatlead to the black solid mass function in fig. 5, where the inflationary scale is η i = O (1) TeV,we have η σ /η ( i ) σ = O (100–1000). If we require a larger η i , the ratio η σ /η ( i ) σ is even larger,which seems less likely. Of course, this should be determined by the configuration of themultidimensional potential, which is beyond the scope of the present work.An inflationary model at the TeV scale was constructed in ref. [54]. In this model ofhybrid inflation, the inflaton is directly coupled to the Higgs field such that the symmetryis restored even at a temperatures lower than the electroweak scale. As the Higgs poten-tial becomes unstable later, the fields roll down in random directions, leading to non-trivialHiggs configurations. In the presence of CP-violation, this might produce a baryon asym-metry, which is referred to as cold electroweak baryogenesis [55–58], and has been studiedextensively with simulations [59–62]. Such a process might also be a source of primordialmagnetic fields [63, 64] and stochastic gravitational waves [65].With many more merger events to be detected by LIGO-Virgo in the near future, themass distribution of black holes can be determined with increasing certainty. Excited as weare, that our mechanism is the only factory of LIGO-Virgo black holes can be ruled out ifthe mass function turns out to have a very different shape, such as one with a slope (in thelog-log plot) much smaller than -4 for large black holes. Appendix: Mass distribution of supercritical black holes
In this appendix, we will find a semi-analytic solution to the bubble wall’s equation ofmotion for relatively large bubbles, and thus determine the mass function of supercriticalblack holes, which was found numerically in section III, and is shown in fig. 4.17he bubble wall’s equation of motion after inflation is (eq. (17))¨ r + (cid:0) − a ˙ r (cid:1) H ˙ r + 2 a r (cid:0) − a ˙ r (cid:1) + K (1 − a ˙ r ) / a = 0 , (28)where we have defined K ≡ ρ b /σ + 6 πσ . This equation can be rewritten as˙ u + 32 t u + 2 γar + K = 0 , (29)where u ( t ) ≡ (cid:112) γ − γ ≡ (1 − a ˙ r ) − / . As discussed in section III, if we assume ahuge γ i (or u ( t i )), the trajectory of the wall can be approximated by a ˙ r ≈
1, which gives r ( t ) ≈ r i + 2 √ t i ( √ t − √ t i ) . (30)Then when the wall comes to a stop with respect to the Hubble flow, we have r s ∼ r i + a s H − i , (31)where we have used r i (cid:29) H − i , which is typically the case for supercritical bubbles in ourdiscussion. In what follows we shall discuss two limits: (a) r i (cid:29) a s H − i , and (b) r i (cid:28) a s H − i . (a) r i (cid:29) a s H − i In this limit, r s ∼ r i . Then by eq. (23), the resulting black hole mass is m ∝ r i . (32)Then it can easily be shown that the mass function (eq. (25)) becomes f ( m ) ∝ m r i d r i d m ∝ m − / . (33)This explains the behavior of f ( m ) on the end of large masses, as shown in fig. 4. (b) r i (cid:28) a s H − i In this limit, the bubbles are assumed to still be in the supercritical regime, but theirinitial radius r i is negligible compared to r s ∼ a s H − i . By eq. (30), we can approximate eq.(29) by ˙ u + 32 t u + 2 ur i √ t/ √ t i + 2 t + K = 0 . (34)18ow we can roughly divide the time interval between t i and t s by two stages: before andafter t ≈ r i / t i ≡ t m , which determines whether r i √ t/ √ t i or 2 t dominates the denominatorin the third term of the above equation.Before t m , we have ˙ u + 32 t u + 2 √ t i r i √ t u + K = 0 , (35)which has an analytic solution u ( t ) = e − √ t i t/r i t / (cid:34) γ i t / i e − t i /r i + K Γ(5 , − t i /r i ) − Γ(5 , − √ t i t/r i )512( √ t i /r i ) (cid:35) . ( t < t m ) (36)Here we have taken into account the initial condition u ( t i ) = γ i . Then at t m , we have u ( t m ) = e − t / m (cid:34) γ i t / i e − t i /r i + K Γ(5 , − t i /r i ) − Γ(5 , − √ t i /r i ) (cid:35) (37)After t m , eq. (34) becomes ˙ u + 52 t u + K = 0 , (38)with the solution u ( t ) = u ( t m ) (cid:18) t m t (cid:19) / + 27 Kt (cid:34)(cid:18) t m t (cid:19) / − (cid:35) . ( t > t m ) (39)At t s , we have u ( t s ) = 0. After some algebra, substituting u ( t m ) by eq. (37) in eq. (39)leads to t s ≈ . (cid:16) γ i K (cid:17) / t / i r / i . (40)In the computations we have assumed the relation γ i > − Kr i /t i , which is true in thecases we are interested in. This estimate for t s has also been verified numerically.Now by eq. (23), we have m ∼ t s ∝ r / i . (41)Then the mass function becomes f ( m ) ∝ m r i d r i d m ∝ m − . . (42)This explains the behavior of f ( m ) in the supercritical regime near the transition, as shownin fig. 4. Note that this relation is but an approximation. For a smaller γ i , it was numericallyfound that the power for m gets slightly larger, giving f ∝ m − , as shown in fig. 5.19 cknowledgments I would like to thank Tanmay Vachaspati, Ville Vaskonen and Alex Vilenkin for stimu-lating discussion and comments. This work is supported by the U.S. Department of Energy,Office of High Energy Physics, under Award No. de-sc0019470 at Arizona State University. [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X , 031040 (2019), 1811.12907.[2] R. Abbott et al. (LIGO Scientific, Virgo) (2020), 2010.14527.[3] K. W. K. Wong, K. Breivik, K. Kremer, and T. Callister (2020), 2011.03564.[4] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez,F. Antonini, D. E. Holz, and C. Pankow (2020), 2011.10057.[5] C. L. Rodriguez, K. Kremer, S. Chatterjee, G. Fragione, A. Loeb, F. A. Rasio, N. C. Weath-erford, and C. S. Ye (2021), 2101.07793.[6] M. Fishbach, Z. Doctor, T. Callister, B. Edelman, J. Ye, R. Essick, W. M. Farr, B. Farr, andD. E. Holz (2021), 2101.07699.[7] S. Bird, I. Cholis, J. B. Mu˜noz, Y. Ali-Ha¨ımoud, M. Kamionkowski, E. D. Kovetz, A. Rac-canelli, and A. G. Riess, Phys. Rev. Lett. , 201301 (2016), 1603.00464.[8] S. Clesse and J. Garc´ıa-Bellido, Phys. Dark Univ. , 142 (2017), 1603.05234.[9] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. , 061101 (2016),1603.08338.[10] M. Raidal, V. Vaskonen, and H. Veerm¨ae, JCAP , 037 (2017), 1707.01480.[11] Y. Ali-Ha¨ımoud, E. D. Kovetz, and M. Kamionkowski, Phys. Rev. D , 123523 (2017),1709.06576.[12] V. Vaskonen and H. Veerm¨ae, Phys. Rev. D , 043015 (2020), 1908.09752.[13] J. Garriga and N. Triantafyllou, JCAP , 043 (2019), 1907.01455.[14] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veerm¨ae, JCAP , 018 (2019), 1812.01930.[15] Z.-C. Chen and Q.-G. Huang, Astrophys. J. , 61 (2018), 1801.10327.[16] A. D. Gow, C. T. Byrnes, A. Hall, and J. A. Peacock, JCAP , 031 (2020), 1911.12685.[17] Y. Wu, Phys. Rev. D , 083008 (2020), 2001.03833.[18] A. Dolgov and J. Silk, Phys. Rev. D , 4244 (1993).
19] S. Clesse and J. Garc´ıa-Bellido, Phys. Rev. D , 023524 (2015), 1501.07565.[20] S. Blinnikov, A. Dolgov, N. K. Porayko, and K. Postnov, JCAP , 036 (2016), 1611.00541.[21] K. Kannike, L. Marzola, M. Raidal, and H. Veerm¨ae, JCAP , 020 (2017), 1705.06225.[22] G. H¨utsi, M. Raidal, V. Vaskonen, and H. Veerm¨ae (2020), 2012.02786.[23] A. Hall, A. D. Gow, and C. T. Byrnes, Phys. Rev. D , 123524 (2020), 2008.13704.[24] K. W. K. Wong, G. Franciolini, V. De Luca, V. Baibhav, E. Berti, P. Pani, and A. Riotto(2020), 2011.01865.[25] J. Garriga, A. Vilenkin, and J. Zhang, JCAP , 064 (2016), 1512.01819.[26] H. Deng and A. Vilenkin, JCAP , 044 (2017), 1710.02865.[27] H. Deng, JCAP , 023 (2020), 2006.11907.[28] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz,T. Bulik, and F. Pannarale, Astrophys. J. , 263 (2015), 1405.7016.[29] P. Ajith et al., Phys. Rev. D , 104017 (2008), [Erratum: Phys.Rev.D 79, 129901 (2009)],0710.2335.[30] I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. Roy. Astron. Soc. , 1086 (2019),1809.02063.[31] H. Deng, J. Garriga, and A. Vilenkin, JCAP , 050 (2017), 1612.03753.[32] V. De Luca, G. Franciolini, and A. Riotto, Phys. Lett. B , 135550 (2020), 2001.04371.[33] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama (2020), 2002.12778.[34] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen, and H. Veerm¨ae, Phys. Rev. D96 , 023514(2017), 1705.05567.[35] P. D. Serpico, V. Poulin, D. Inman, and K. Kohri, Phys. Rev. Res. , 023204 (2020),2002.10771.[36] V. De Luca, V. Desjacques, G. Franciolini, P. Pani, and A. Riotto (2020), 2009.01728.[37] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP , 044 (2020), 2005.05641.[38] C. Alcock et al. (MACHO), Astrophys. J. , 281 (2000), astro-ph/0001272.[39] P. Tisserand et al. (EROS-2), Astron. Astrophys. , 387 (2007), astro-ph/0607207.[40] M. Oguri, J. M. Diego, N. Kaiser, P. L. Kelly, and T. Broadhurst, Phys. Rev. D , 023518(2018), 1710.00148.[41] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. , 161102 (2019), 1904.08976.[42] Y. Inoue and A. Kusenko, JCAP , 034 (2017), 1705.00791.
43] D. Lynden-Bell, Nature , 690 (1969).[44] J. Kormendy and D. Richstone, Ann. Rev. Astron. Astrophys. , 581 (1995).[45] Z. Haiman, Astrophys. J. , 36 (2004), astro-ph/0404196.[46] J. Kormendy and L. C. Ho, Ann. Rev. Astron. Astrophys. , 511 (2013), 1304.7762.[47] F. Wang, J. Yang, X. Fan, J. F. Hennawi, A. J. Barth, E. Banados, F. Bian, K. Boutsia,T. Connor, F. B. Davies, et al., The Astrophysical Journal Letters , L1 (2021).[48] S. G. Rubin, A. S. Sakharov, and M. Yu. Khlopov, J. Exp. Theor. Phys. , 921 (2001), [J.Exp. Theor. Phys.92,921(2001)], hep-ph/0106187.[49] R. Bean and J. Magueijo, Phys. Rev. D66 , 063505 (2002), astro-ph/0204486.[50] N. Duechting, Phys. Rev.
D70 , 064015 (2004), astro-ph/0406260.[51] B. Carr and J. Silk, Mon. Not. Roy. Astron. Soc. , 3756 (2018), 1801.00672.[52] B. Carr, F. Kuhnel, and L. Visinelli, Mon. Not. Roy. Astron. Soc. , 2 (2021), 2008.08077.[53] V. Berezin, V. Kuzmin, and I. Tkachev, Phys. Rev. D , 2919 (1987).[54] E. J. Copeland, D. Lyth, A. Rajantie, and M. Trodden, Phys. Rev. D , 043506 (2001),hep-ph/0103231.[55] J. Garcia-Bellido, D. Y. Grigoriev, A. Kusenko, and M. E. Shaposhnikov, Phys. Rev. D ,123504 (1999), hep-ph/9902449.[56] L. M. Krauss and M. Trodden, Phys. Rev. Lett. , 1502 (1999), hep-ph/9902420.[57] J. Smit and A. Tranberg, in (2003), pp. 509–513, hep-ph/0210348.[58] T. Konstandin and G. Servant, JCAP , 024 (2011), 1104.4793.[59] A. Tranberg and J. Smit, JHEP , 012 (2006), hep-ph/0604263.[60] A. Tranberg, J. Smit, and M. Hindmarsh, JHEP , 034 (2007), hep-ph/0610096.[61] Z.-G. Mou, P. M. Saffin, and A. Tranberg, JHEP , 010 (2017), 1703.01781.[62] Z.-G. Mou, P. M. Saffin, and A. Tranberg, JHEP , 197 (2018), 1803.07346.[63] A. Diaz-Gil, J. Garcia-Bellido, M. Garcia Perez, and A. Gonzalez-Arroyo, JHEP , 043(2008), 0805.4159.[64] Z.-G. Mou, P. M. Saffin, and A. Tranberg, JHEP , 075 (2017), 1704.08888.[65] J. Garcia-Bellido, D. G. Figueroa, and A. Sastre, Phys. Rev. D , 043517 (2008), 0707.0839., 043517 (2008), 0707.0839.