Identifying the most crucial parameters of the initial curvature profile for primordial black hole formation
Tomohiro Nakama, Tomohiro Harada, A. G. Polnarev, Jun'ichi Yokoyama
aa r X i v : . [ g r- q c ] J a n RESCEU-45/13, RUP-13-11
Identifying the most crucial parameters of the initial curvatureprofile for primordial black hole formation
Tomohiro Nakama,
1, 2
Tomohiro Harada, A. G. Polnarev, and Jun’ichi Yokoyama
2, 5 Department of Physics, Graduate School of Science,The University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan Research Center for the Early Universe (RESCEU),Graduate School of Science, The University of Tokyo,Bunkyo-ku, Tokyo 113-0033, Japan Department of Physics, Rikkyo University, Toshima, Tokyo 175-8501, Japan Astronomy Unit, School of Physics and Astronomy,Queen Mary University of London,Mile End Road, London E1 4NS, United Kingdom Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), WPI, TODIAS,The University of Tokyo, Kashiwa, Chiba 277-8568, Japan (Dated: July 12, 2018)
Abstract
Primordial black holes (PBHs) are an important tool in cosmology to probe the primordial spec-trum of small-scale curvature perturbations that reenter the cosmological horizon during radiationdomination epoch. We numerically solve the evolution of spherically symmetric highly perturbedconfigurations to clarify the criteria of PBHs formation using an extremely wide class of curvatureprofiles characterized by five parameters, (in contrast to only two parameters used in all previouspapers) which specify the curvature profiles not only at the central region but also at the outerboundary of configurations. It is shown that formation or non-formation of PBHs is determinedessentialy by only two master parameters one of which can be presented as an integral of curvatureover initial configurations and the other is presented in terms of the position of the boundary andthe edge of the core.
PACS numbers: . INTRODUCTION It is well known that a region with large amplitude curvature profile can collapse to a primordialblack hole (PBH) [1, 2]. The probability of existence of such regions and hence the abundanceof PBHs depends on the physical conditions in the very early universe. PBHs are formed soonafter the region enters the cosmological horizon during the radiation-dominated epoch. Even if thecontribution of the energy density of PBHs to the total cosmic energy density is extremely smallat formation, it drastically increases in the course of subsequent cosmic expansion during radiationdomination. For this reason PBHs can be used as a unique and powerful tool for the explorationof very small-scale structure of the early universe. This tool is really unique since there are noalternative techniques to study this small-scale structure.PBHs with different masses have different cosmological and astrophysical implications. ThePBHs with mass smaller than ∼ g are especially interesting because by now they would haveevaporated through Hawking radiation [3] and their abundance is constrained by the effects ofemitted high-energy particles and photons on the Big-Bang Nucleosynthesis [4–9], the gamma-raybackground [10–12], galactic and extragalactic antiprotons [13] as well as on the cosmic microwavebackground radiation (CMB).PBHs with larger masses would still be present and are constrained by dynamical and lensingeffects [14] and by the stochastic gravitational wave background [15, 16]. All these constraints areupdated and summarized in [17]. For the mass range between 6 × g and 4 × g the moststringent constraint is given by the cosmic dark matter density [18].Even if PBHs would never have been detected these constraints (being reformulated in terms ofthe constraints on the amplitude of small-scale perturbations) will provide valuable information oninflationary cosmological models [19–22], which predict formation of super-horizon curvature per-turbations [23–26]. So far their large-scale components have been precisely tested by observationsof CMB [27, 28] and large-scale structure. It is important to probe the perturbation spectrum onsignificantly smaller scales as well in order to obtain more helpful information to single out thecorrect inflation model. Indeed there exist a number of inflation models that predict features atsome small scales [29–42] which may lead to PBHs abundances not compatible with observationalconstraints. So we can gain insights into primordial perturbations on small scales even from theabsence of observed PBHs. Or if PHBs are detected in the future it will extend our insight intoinflationary models drastically.In order to consistently use PBHs as a cosmological tool to probe the spectrum of small-scalecurvature perturbations , the relation between their mass function and the primordial perturbationspectrum must be clarified. This enterprise involves several steps. First we must clarify the forma-tion condition of PBHs as a functional of the primordial curvature profiles of perturbed regions.Second we must calculate the resultant mass of a PBH for each given profile. Finally we mustcalculate the probability distribution functional of curvature perturbation spectrum which givesrealization probability of curvature profiles leading to PBH formation. This paper is predominantlydevoted to the first issue and touches the second one in a preliminary manner. We will considerthe second issue in more detail in another paper.Originally the problem of PBH formation was studied analytically by considering the balancebetween gravity and pressure gradient which hampers contraction, yielding a simple analytic cri-terion of PBH formation [43, 44] 13 . ¯ δ hc , (1.1) here ¯ δ hc is the energy density perturbation averaged over the overdense region evaluated at thetime of horizon crossing. This criterion has long been used in papers on theoretical prediction ofPBH abundance (but has recently been refined in [60]). In this simple picture, since ¯ δ hc is thedensity perturbation averaged over the overdense region, the dependence on the profile or shapeof perturbed regions has not been taken into account.However, recent numerical analyses have shown that the condition for PBH formation doesdepend on the profile of perturbation [45, 46] (see also [47, 48]). Both [45] and [46] used two-parameter families of the initial profile (one corresponding to the amplitude of overdensity andthe other to the width of the transition region) and obtained two parametric conditions of PBHformation. It was clear from the above publications that one parametric description was notsufficient. However it was not clear whether the two-parametric description is good enough (say,in comparison with three-parametric description). In the present paper, we extend these precedinganalyses by making much more numerical computations of PBH formation based on the initialcurvature profile including many more parameters, adopting the five-parameter family of profiles.We show that the criterion of PBH formation can still be expressed in terms of two crucial (master)parameters which represent the averaged amplitude of over density in the central region and thewidth of transition region at outer boundary, even though the considered profiles belong to the five-parametric family. Thereby we have obtained extra evidence that the two-parametric descriptionmentioned above is self-consistent and provided a reliable physical interpretation of the criterionwe have obtained.In our previous paper [49] (hereafter PNY), we investigated the time evolution of spherical per-turbed regions embedded in a flat Friedmann-Lemaitre-Robertson-Walker (FLRW) universe whilethe region is outside the cosmological horizon, using an asymptotic expansion. In order to calculatethe growth of the perturbed region after the horizon reentry, the Einstein equations need to besolved numerically. In this paper, we first present the results of the numerical computation in theslicing used in our previous work, with the initial data generated using the asymptotic expansion.In this slicing, at spatial infinity the time coordinate coincides with that in the background FLRWuniverse , so this slicing is sometimes referred to as cosmic time slicing in the literature.During the numerical computation a singularity appears at the centre after the apparent horizonis formed and as a result the computation has to be stopped and the subsequent accretion onto thePBH cannot be followed. That is, the eventual mass of the PBH cannot be obtained. In order toavoid this problem, we employ what is called the null slicing, which is also discussed in this paper.Finally, note that PBHs are formed only from extremely high peaks of perturbation, correspond-ing to the tail of the probability distribution of primordial perturbation. The shape of these peakshas been calculated (see, for example, [50, 51]) and has turned out to be nearly spherically sym-metric and monotonic near the peak. So in this paper we continue to consider only the sphericallysymmetric profiles which are monotonic near the centre.The rest of the paper is organized as follows. In § II, we briefly discuss the initial condition interms of the asymptotic expansion of the Einstein equations. In § III we present the basic equationsused in the numerical computations. Then in § IV we review the time evolution of the perturbationin the cosmic time slicing . In § V we present the two parameters crucial for PBH formaiton. § VI is devoted to the determination of the PBH masses using the null slicing and § VII to conclusionand discussion. I. SETTING UP THE INITIAL CONDITION
Assuming spherical symmetry, it is convenient to divide the collapsing matter into a system ofconcentric spherical shells and to label each shell with a Lagrangian comoving radial coordinate r .Then the metric can be written in the form used by Misner and Sharp [52]: ds = − a dt + b dr + R ( dθ + sin θdφ ) , (2.1)where R , a and b are functions of r and the time coordinate t . We consider a perfect fluidwith the energy density ρ ( r, t ) and pressure P ( r, t ) and a constant equation-of-state parameter γ , P ( r, t ) = γρ ( r, t ). We express the proper time derivative of R as U ≡ ˙ Ra , (2.2)with a dot denoting a derivative with respect to t .We define the mass, sometimes referred to as the Misner-Sharp mass in the literature, withinthe shell of circumferential radius R by M ( r, t ) = 4 π Z R ( r,t )0 ρ ( r, t ) R dR. (2.3)We consider the evolution of a perturbed region embedded in a flat Friedmann-Lemaitre-Robertson-Walker (FLRW) Universe with metric ds = − dt + S ( t )( dr + r dθ + r sin θdφ ) , (2.4)which is a particular case of (2.1). The scale factor in this background evolves as S ( t ) = (cid:18) tt i (cid:19) α , α ≡ γ ) , (2.5)where t i is some reference time.We denote the background solution with a suffix 0. In terms of the metric variables defined in(2.1), we find a = 1 , b = S ( t ) , R = rS ( t ) . (2.6)The background Hubble parameter is H ( t ) = ˙ R a R = ˙ SS = αt , (2.7)and the energy density is calculated from the Friedmann equation, ρ ( t ) = 3 α πGt . (2.8)The energy density perturbation is defined as δ ( t, r ) ≡ ρ ( t, r ) − ρ ( t ) ρ ( t ) . (2.9)We introduce a variable H defined by H ( t, r ) ≡ ˙ RaR = UR . (2.10) he curvature profile K ( t, r ) is defined by rewriting b as b ( t, r ) = R ′ ( t, r ) p − K ( t, r ) r . (2.11)This quantity K ( t, r ) vanishes outside the perturbed region so that the solution asymptoticallyapproaches the background FLRW solution at spatial infinity.We denote the comoving radius of a perturbed region by r i , the precise definition of which willbe given later (see eq. (2.15)), and define a dimensionless parameter ǫ in terms of the square ratioof the Hubble radius H − to the physical length scale of the configuration, ǫ ≡ (cid:18) H − S ( t ) r i (cid:19) = ( ˙ Sr i ) − = t α i t β α r , β ≡ − α ) . (2.12)When we set the initial conditions for PBH formation, the size of the perturbed region is muchlarger than the Hubble horizon. This means ǫ ≪ K (0 , r ) ≡ K i ( r ) , (2.13)where K i ( r ) is an arbitrary function of r which vanishes outside the perturbed region. Note that,from (2.11), K i ( r ) has to satisfy the condition K i ( r ) < r . (2.14)We normalize radial Lagrangian coordinate r in such a way that K i (0) = 1.In order to represent the comoving length scale of the perturbed region, we use the co-movingradius, r i , of the overdense region. We can calculate r i by solving the following equation for theenergy density perturbation defined by (2.9): δ ( t, r i ) = 0 . (2.15)Since the initial condition is taken at the superhorizon regime, when ǫ is extremely small, thefollowing lowest-order solution[46] δ ( t, r ) = 2 r r ( r K i ( r )) ′ ǫ ( t ) (2.16)suffices to calculate r i , which is obtained by solving3 K i ( r i ) + r i K ′ i ( r i ) = 0 . (2.17)Note that the physical length scale in the asymptotic Friedmann region is obtained by multiplyingby the scale factor S ( t ), the normalization of which we have not specified. We can therefore set upinitial conditions for the PBH formation with arbitrary mass scales by adjusting the normalizationof S ( t ) which appears in the expansion parameter.We also introduce the averaged over-density, denoted by ¯ δ and defined as the energy densityperturbation averaged over the over-dense region as follows:¯ δ ( t ) ≡ (cid:18) πR ( t, r od ( t )) (cid:19) − Z R ( t,r od ( t ))0 πδR dR. (2.18)Here r od ( t ) represents the comoving radius of the overdense region and is the solution of δ ( t, r od ( t )) = 0. It turns out that r od ( t ) is very close to r i calculated from (2.17), i.e. lowest-order expansion.Hereafter, we will always use the initial conditions obtained in this section. II. BASIC EQUATIONS USED IN THE NUMERICAL COMPUTATIONS
The following equations were used in [53] to analyze the gravitational collapse of sphericallysymmetric masses: ˙ U = − a (cid:18) πR Γ w P ′ + M GR + 4 πGP R (cid:19) , (3.1)˙ R = aU, (3.2)( νR ) · νR = − a U ′ R ′ , (3.3)˙ E = − P (cid:18) ν (cid:19) · , (3.4)( aw ) ′ aw = E ′ + P (1 /ν ) ′ w , (3.5) M = 4 π Z r ρR R ′ dr, (3.6)Γ = 4 πνR R ′ , (3.7) P = γρ, (3.8) w = E + P/ν, (3.9)where E ≡ ρ/ν and ν ≡ πbR . (3.10)The constraint equation reads (cid:18) R ′ b (cid:19) = Γ = 1 + U − MR . (3.11)We rewrite these equations in terms of the barred variables, which are defined by factoring out thescale factor S and the background energy density ρ as shown below:¯ R ≡ R/S, (3.12)¯ a ≡ a, (3.13)¯ U ≡ U/ ˙ S, (3.14)¯ ρ ≡ ρ/ρ , (3.15)¯ M ≡ M/ ( ρ S ) , (3.16)¯ b ≡ b/S, (3.17)¯ ν ≡ S ν, (3.18)¯Γ ≡ Γ , (3.19)¯ P ≡ P/ρ , (3.20)¯ w ≡ w/ ( ρ S ) . (3.21) irst, using these definitions, (3.1) leads to¨ S ¯ U + ˙ S ˙¯ U = − ¯ a (cid:18) πS ¯ R ¯Γ S ¯ w ¯ P ′ + Sρ ¯ M G ¯ R + 4 πGρ ¯ P S ¯ R (cid:19) . (3.22)Using the Friedmann equation ˙ SS ! = 8 πG ρ (3.23)on the right side of (3.22) we find¨ S ¯ U + ˙ S ˙¯ U = − ¯ a πS ¯ R ¯Γ S ¯ w ¯ P ′ + 38 π ˙ S S ¯ M ¯ R + 32 ˙ S S ¯ P ¯ R ! . (3.24)Dividing the both sides by ˙ S /S and noting S = (cid:18) tt i (cid:19) α , S ¨ S ˙ S = α − α , (3.25)we obtain tα ˙¯ U = − ¯ a (cid:18) π ˙ S − ¯ R ¯Γ¯ w ¯ P ′ + 38 π ¯ M ¯ R + 32 ¯ P ¯ R (cid:19) + 1 − αα ¯ U . (3.26)Defining a new time coordinate by τ ≡ αβ log ǫ (3.27)and noting the following relationship tα ∂∂t = ∂∂τ , (3.28)which is obtained from the definition of ǫ , (2.12), one finds( ¯ U ) τ = − ¯ a (cid:18) πr i exp (cid:18) βα τ (cid:19) ¯ R ¯Γ¯ w ¯ P ′ + 38 π ¯ M ¯ R + 32 ¯ P ¯ R (cid:19) + 1 − αα ¯ U , (3.29)where the subscript τ denotes a time derivative with respect to τ . Similarly, (3.2)-(3.4) are rewrittenin terms of the barred variables as follows:¯ R τ = ¯ a ¯ U − ¯ R, (3.30)(¯ ν ¯ R ) τ ¯ ν ¯ R = − ¯ a ¯ U ′ ¯ R ′ + 1 , (3.31)¯ E τ = − ¯ P (cid:18) ν (cid:19) τ . (3.32)Since the barred variables are defined by factoring out S and ρ , depending only on time, theequations (3.5)-(3.10), which do not contain time derivatives, do not change their appearance evenafter being rewritten in terms of the barred variables:(¯ a ¯ w ) ′ ¯ a ¯ w = ¯ E ′ + ¯ P (1 / ¯ ν ) ′ ¯ w , (3.33)¯ M = 4 π Z r ¯ ρ ¯ R ¯ R ′ dr, (3.34)¯Γ = 4 π ¯ ν ¯ R ¯ R ′ , (3.35)¯ P = γ ¯ ρ, (3.36) w = ¯ E + ¯ P / ¯ ν, (3.37)¯ ν = 14 π ¯ b ¯ R . (3.38)The following equation is equivalent to (3.33) and is also useful:¯ a = ¯ ρ − γ γ . (3.39)Lastly, the constraint equation (3.11) can be rewritten as¯Γ = 1 + ˙ S ( ¯ U − M π ¯ R ) . (3.40)This equation is not used for time evolution but rather is used to estimate numerical errors. Now,all the relevant equations in terms of the barred variables, (3.29), (3.30), (3.31), (3.32), (3.34),(3.35), (3.36), (3.37), (3.38), (3.39) and (3.40) have been obtained. The boundary conditions areimposed such that ¯ U = ¯ R = ¯ M = 0 and Γ = 1 at the centre, and ¯ a = ¯ ρ = 1 on the outer boudaryso that the numerical solution is smoothly connected to the FLRW solution.In the next section we discuss the numerical solutions of the basic equations presented here. IV. TYPICAL TIME EVOLUTION OF PERTURBED REGIONS IN THE COS-MIC TIME SLICING: A REVIEW
In order to test our numerical code which solves the evolution equations described in the previoussection, we first adopt the following two-parameter family of curvature profile K i ( r ) = (cid:20) B (cid:16) rσ (cid:17) (cid:21) exp (cid:20) − (cid:16) rσ (cid:17) (cid:21) , (4.1)which was also investigated in [46](hereafter PM), in order to examine the consistency with theresults obtained in PM. As mentioned previously, the amplitude of the profile is set to unity atthe origin, meaning that we use the same normalization as a spatially closed Friedmann universe.Here the parameters B and σ control the shapes of initial perturbations. The range of B is limitedto 0 ≤ B ≤ K i ( r ) is a monotonic function. Two examples of (4.1) are shown in Figure 1.In order to relate the initial curvature perturbation to the amplitude of the density perturbation,following PM, let us approximately evaluate the energy density perturbation averaged over the over-dense region, denoted by ¯ δ and defined by (2.18), at the time of the horizon crossing. Using (2.16)¯ δ ( t ) becomes ¯ δ ( t ) = (cid:18) πr (cid:19) − Z r i πr r K i ( r )) ′ ǫ ( t ) dr = 23 K i ( r i ) r ǫ ( t ) . (4.2)By setting ǫ ( t ) = 1, we define ¯ δ hc ≡ K i ( r i ) r . (4.3)Note that (4.2) is valid only when ǫ ( t ) ≪
1, so this formula gives just an approximate value of ¯ δ ( t )at the time of the horizon crossing. Still, (4.3) gives a good indicator of how strong gravity is.A profile for which ¯ δ hc is small corresponds to a small amplitude perturbation and PBH is notformed from this kind of initial configurations. The narrow profile in Figure 1 with ( B, σ ) = (0 , . δ hc = 0 .
04, while the wide profile with (
B, σ ) = (1 , .
7) corresponds to ¯ δ hc = 0 . K i ( r ) approximately tells how large ¯ δ hc is, as can be seen from (4.3). In the IG. 1: A wide and narrow initial curvature profiles K i ( r ) represented by (4.1). Note that K i ( r )has to satisfy the condition K i ( r ) < /r .following the time evolution of the initial perturbed region for these two cases is presented asillustration.First, let us look at the case of the narrow profile (¯ δ hc = 0 . ρ , normalised by the background energy density, is shown in Figure 2. In this casethe growth of the perturbation stops at some point in time after the horizon crossing and theenergy density starts to decrease in the central region, with a sound wave propagating outward,the amplitude of which gradually decreases. That is, the initial perturbation dies away and theeventual state at t = ∞ is the flat FLRW universe.Next, let us consider a case where the amplitude of the initial perturbation is sufficiently large(( B, σ ) = (1 , . , ¯ δ hc = 0 .
51) and a PBH is eventually formed. The time evolution of ¯ ρ , U and2 M/R in this case is shown in Figure 3.Near the centre the energy density increases drastically and the central perturbed region grad-ually expands outward. The central perturbed region is always surrounded by the under-denseregion.From (2.10) U is written as U = HR and it corresponds to the recession velocity in theFLRW universe. At an early stage, when the amplitude of the perturbation is small, U is positiveeverywhere-reflecting the expansion of the universe. In the central region where gravity becomesincreasingly stronger, the expansion decelerates rapidly and therefore U decreases rapidly. Thenat some point in time, the central region stops expanding and U becomes negative there, startinggravitational contraction.The mass M is defined by (2.3) and represents the total energy contained in a sphere of radius R .When the amplitude of perturbation is small, M ∼ ρ ( t ) R ( t, r ) hence 2 M/R ∝ R . As mentionedearlier, at a later stage the energy density in the central region increases dramatically as a blackhole is formed, which is contrasted with the decrease in the energy density in the surroundingregion due to the expansion of the universe. In such a situation the mass profile in the central IG. 2: An example of the time evolution of the energy density perturbation ¯ ρ , normalised bythe background energy density, in a case where ( B, σ ) = (0 , .
3) and no PBH is formed. The plotsare numbered in order of time evolution.region is steep but is relatively gentle in the surrounding region. This feature of the mass profilecan be easily understood by an analogy with a situation where a star resides in the vacuum, inwhich case the mass profile is a monotonically increasing function inside the radius of the star andis flat outside that radius. Due to this behaviour of M at a later stage, in the region away from thecentre the mass only weakly depends on R and so approximately 2 M/R ∝ R − there. This meansthat a peak appears in the profile of 2 M/R at some moment in cases where perturbation growssufficiently. Specifically, in cases where a black hole is formed, the height of this peak exceeds unityand this implies the formation of the apparent horizon from the following arguments[54].Suppose the trajectory of a photon moving outward, along which adt = bdr, (4.4)so along the geodesic of this photon, dRdt = ∂R∂t dt + ∂R∂r dr = a ( U + Γ) dt, (4.5)where the definition of U (2.2) and Γ = R ′ /b , followed from (3.7) and (3.10), have been used.From this we find dR/dt = 0 where U = − Γ holds, meaning the photon reaching that point cannotescape further away from the centre. Since we find U = − Γ when 2
M/R = 1 from the constraint(3.11), the peak of 2
M/R with the height exceeding unity means that there exist photons whichare trapped by the gravitational potential of the central black hole and cannot escape to infinity.Namely, the apparent horizon has been formed.We have calculated the time evolution with various values of B and σ and obtained results fullyconsistent with PM. In the next section, we discuss the PBH formation condition for a much widerclass of profiles than that defined by (4.1) . IG. 3: Examples of time evolution of ¯ ρ (top-left), U (top-right), M (down-left) and 2 M/R (down-right), calculated using the cosmic time slicing. These are obtained for the case (
B, σ ) = (1 , . V. TWO MASTER PARAMETERS CRUICIAL FOR PBHS FORMATION
We now proceed to our full analysis introducing the following function K i ( r ) = A " B (cid:18) rσ (cid:19) n exp " − (cid:18) rσ (cid:19) n + (1 − A ) exp " − (cid:18) rσ (cid:19) , (5.1)which can represent various shapes of profiles using the five parameters as is shown in Figure 4.Profiles with A = 1 include the profiles studied in PM but can also incorporate those with a steepertransition to the homogeneous region, whereas the inclusion of the second term makes it possibleto further realize profiles with tails (see Figure 4). It also turned out that this function with n = 2 can fit the profiles investigated in [45] sufficiently well, after translating their profiles intoour K i ( r ) since a different coordinate system was used their. Thus our function not only includesthose investigated in previous work but also enables us to investigate new shapes of profiles.For each set of five parameters, the time evolution of the perturbation is calculated to revealwhether or not a PBH is formed at the end. First, let us try to describe PBH formation conditionfor these various kinds of profiles using the quantity ¯ δ hc andΩ ≡ max r | K ′ i ( r ) | , (5.2)which corresponds to ∆ in PM and may provide a measure of the density gradient, or the pressuregradient for the case of the radiation-dominated universe. As can be seen in Figure 5, these twoquantities fail to distinguish between the profiles which lead to PBH formation and those which IG. 4: Dependence of the shapes of profiles represented by (5.1) on parameters. In each panel,the dependence of the shape on one of the parameters is shown with the rest of the parametersfixed.do not. One conceivable reason why this combination fails is because ¯ δ hc is too sensitive to theedge r = r i of configurations. To see this let us consider two profiles with almost the same shapenear the centre but one of these has slightly wider tail with a larger r i than the other. For thiswider profile, ¯ δ hc is smaller since this quantity is the energy density perturbation averaged inside r < r i . However, what is important in physically determining whether a PBH is formed or notis the shape of profile near the centre, while the information around the outer boundary is less IG. 5: The blue points correspond to profiles which lead to PBH formation and the red points toprofiles which do not. Blue and red points are not well separated by a single line in this diagram.important. Therefore, proper quantities which can clearly distinguish between the PBH-formingprofiles and the non-forming ones should take the same value for those two profiles with the sameshape in the centre. But ¯ δ hc takes different values for these two and as a result this fails, especiallywhen profiles with tails are included, as depicted in Figure 5. This also shows that investigatinglimited types of shapes only, as has been done in the previous works, is insufficient to reveal thePBH formation condition which is generally accurately applicable to various profiles.In order to obtain reliable measures to judge if PBHs are formed or not, we must introduceparameters which can separate the effects of the central part collapsing into holes and those of thesurrounding tail part. Let us define r p (0 < p <
1) by K i ( r p ) = p. (5.3)After a number of trials of combinations ∆ p,q ≡ r p − r q , (5.4) I s,j ≡ Z r s r j K i ( r ) dr, (5.5)where p < q and j = 0 , ,
2, it turned out that a relatively clear separation between configurationswhich collapse to PBHs and those which do not is obtained by the following combination:∆ ≡ r / , / = r / − r / (5.6)and I ≡ I / , = Z r / r K i ( r ) dr. (5.7)The quantity I is similar to ¯ δ hc in the sense that it is a volume integral of curvature profile closelyrelated to the density fluctuations (see eq. (4.2)), but the important difference is that it is not ffected by the tail part of the overdense region, since the upper limit of the integral is cut offat r / . On the other hand, ∆ represents the width of the transition region between K i ∼ K i ∼ ≤ A ≤
1, 0 ≤ B ≤ σ min ≤ σ ≤ σ ≤ min { σ , . / √ − A } , where σ min is chosento search only the profiles relevant to revealing the PBH formation condition and n = 1 , , , , σ is because if σ is too large compared to σ , the secondterm of the function (5.1) dominates and therefore this function reduces to a single gaussian. Ineach numerical simulation, the values of the five parameters are chosen randomly from the rangeof each parameter. The initial conditions are set up by the second order asymptotic expansion.The eventual numerical errors, estimated by the constraint equation (3.40) have been confirmedto be less than a few .As is seen there the condition for PBH formation can be quite well described by the followingfitting formula:( S (∆ − ∆ b ) + I b )Θ( − (∆ − ∆ b )) + ( S (∆ − ∆ b ) + I b )Θ(∆ − ∆ b ) < I, (5.8)where Θ denotes the unit step function and ( S , S , ∆ b , I b ) = ( − . , − . , . , . I is smaller. Thisis because when ∆ is larger, the pressure gradients are smaller and in addition gravity is relativelystronger even away from the centre, in which case gravity near the centre, measured by I , needsnot be so large compared to cases with a smaller ∆. Put differently, for I . . ≡ I cr , profileswith a smaller ∆ do not result in PBH formation because the pressure gradient is so large that thegravitational collapse is hindered.On the other hand, for I & I cr the criterion is not so sensitive to ∆ for the following reason.Since we are using a normalization condition K i (0) = 1, the configurations with a large I must alsohave a large r / . Thus those with I & I cr have such a large r / that the difference in ∆ does notaffect the formation criterion much. Note that (∆ , I ) can distinguish between the PBH-formingprofiles and the non-forming cases much more decisively than (Ω , ¯ δ hc ). This is partly because (∆ , I )does not reflect the information near the outer boundary around r = r i so much as ¯ δ hc .Some of the readers may be worried about the existence of some red points corresponding to non -PBH forming profiles in the BH formation region (the blue region), namely, false positives,and some blue points in the red region (false negatives). Note that it is in principle impossible todistinguish between two possibilities with 100 accuracy using only two quantities, while consideringthe profiles controlled by as many as five parameters. However, we argue that this choice ofquantities, (∆ , I ), separates two regions fairly well with high accuracy. The false positives and thefalse negatives constitute only around 2.1 of all the points in Figure 6. That is, eq.(5.8) can tellwhether a PBH is formed or not with around 2.1 accuracy. Note that this percentage is calculatedonly from the profiles appearing in Figure 6, which are randomly generated relatively near thethreshold line to reveal PBH formation condition. The fitting formula was obtained by minimizingthis percentage. On the other hand, the points in the gmixed regionh in Figure 5 constitute around10 of all the points. Therefore, the accuracy would be around 10 for the case of (Ω , ¯ δ hc ), if the hreshold line was drawn somehow in Figure 5, as is done in Figure 6. Note that Figures 6 and 5are depicted with exactly the same set of 5-parameter choices.The dashed line in Figure 6 corresponds to the Carr’s condition eq.(1.1). Note that the value’1/3’ was obtained in the uniform hubble slicing and this value corresponds to ¯ δ hc ≃ .
22 in thecomoving slicing [60]. It turned out that the profiles for which ¯ δ hc ≃ .
22 is distributed around thedashed line in Figure 6, so the original Carr’s condition roughly corresponds to 0 . . I .The horizontal line in Figure 6, originating from K i ( r ) < /r , is obtained as follows. First,note that the ideal profile which gives the maximum value of I for some ∆ is K i ( r ) = 1(0 < r < , /r (1 < r < p / , p / < r ), when ∆ ≤ √ − p / ≃ .
35. In addition, when p / − p / ≃ . ≤ ∆, the value of I for this ideal profile, namely, the maximum value of I forcorresponding ∆, is 1 / p / − ≃ .
62, hence the horizontal line.Let us summarize this section emphasizing the key points of what has been done. What isimportant and new is that we have investigated various shapes of profiles altogether, including thoseprofiles which were not investigated in the previous work as well as those investigated previously.We formed ”master variables” which decide if PBHs are formed or not. This means that eventhough our profiles now include 5 parameters, only two master variables are important in order todetermine if a PBH forms. Such a detailed analysis has never been done before. Indeed note thatin the previous publications the variety of shapes was restricted, where the profiles were usuallyparameterized by at most two parameters. Figure 6 tells that the PBH formation condition canbe simply described by using only two quantities for the whole variety of profiles, represented by5 parameters. To reach our conclusion depicted in Figure 6, we made a lot of trials and errors tofind the two master variables which we can use to clearly separate the domains of PBH formationand non-formation. So our result in Figure 6 is not merely a re-parameterization of the previouslyknown results but goes far beyond them. We stress once again that it is important to use variableswhich can separate the effects of the central part and the boundary region.In this section the time evolution of the perturbed region is calculated in the cosmic timeslicing, but the problem with this slicing is that computation stops because a singularity showsup at the centre soon after the apparent horizon is formed. Thus, although this slicing serveswell to determine the criteria of PBH formation, we cannot use it to investigate the final mass,which is essential in order to determine the mass spectrum. In the next section, we discuss thedetermination of the mass with a singularity avoided using what is called the null slicing. IG. 6: The PBH formation condition for the profiles represented by (5.1). The blue and redpoints correspond respectively to the profiles which lead to the black hole formation and those whichdo not. The shaded region labelled ”unphysical” corresponds to the profiles which do not satisfy K i ( r ) < /r and therefore are unphysical. The profile used as an example of the PBH-formingcases in this paper corresponds to the yellow star in this figure. The dashed line corresponds tothe Carr’s condition. A H L B H L r K i H r L C H L D H L r K i H r L E H L F H L r K i H r L FIG. 7: Profiles corresponding to the points marked as A-F in Figure 6 with values of (∆ , I ) foreach profile. I. ACCRETION ONTO PBHS AND NULL SLICING
The determination of the mass without facing a singularity by using the null slicing [48, 55–58]is discussed in this section. In this slicing, the space-time is sliced along the null geodesics ofhypothetical photons emitted from the centre and reaching a distant observer. In other words, thespace-time is sliced with the hyper-surfaces, each defined by a constant null coordinate u , the so-called observer time defined shorterly. By this construction of the null slicing, only the informationoutside the horizon is calculated, without looking into what happens inside the apparent horizon.The initial conditions are given on some hypersurface defined by u =const., which is depicted bya blue dotted line in Figure 8, and are obtained using the cosmic time slicing by calculating thenull geodesic of a hypothetical photon which reaches a distant observer after being emitted fromthe centre at some moment in time, while at the same time recording the physical quantities ofthis null geodesic [57]. As is seen in Figure 8, in this slicing, the information can be obtainedFIG. 8: Illustration of how a singularity is avoided in the null slicing.without facing a singularity until a sufficiently later time when the eventual mass of a PBH can bedetermined.Let us define the time variable u by first noting adt = bdr (6.1)along an outgoing photon. Then, u is defined by f du = adt − bdr, (6.2)where f is the lapse function and is necessary to make du a perfect differential. From this defini-tion, (6.1) holds along the hyper-surfaces each defined by u =const., meaning that these surfacescorrespond to the null geodesics of outgoing photons. Using u as the time variable then means thatthe space-time is sliced with the null slices. A boundary condition on the lapse function is imposedby setting a ( u, r = ∞ ) = f ( u, r = ∞ ) = 1, hence u = t at the surface defined by r = ∞ . The hysical meaning of this boundary condition is that u is chosen to coincide with the proper timemeasured by a distant observer residing at spatial infinity in the background FLRW universe. Forthis reason, the null slicing is also sometimes referred to as observer time slicing in the literature.Let us look at the relationship between the derivatives in the two slicings. For the spatialderivatives, one finds from (6.2) 1 b ∂∂r (cid:12)(cid:12)(cid:12) u = 1 b ∂∂r (cid:12)(cid:12)(cid:12) t + 1 a ∂∂t (cid:12)(cid:12)(cid:12) r . (6.3)Specifically for a function f ( t ), depending only on t ,1 b ∂∂r (cid:12)(cid:12)(cid:12) u f ( t ) = 1 a ∂∂t (cid:12)(cid:12)(cid:12) r f ( t ) . (6.4)For the time derivatives one immediately finds from (6.1), ∂∂u (cid:12)(cid:12)(cid:12) r = fa ∂∂t (cid:12)(cid:12)(cid:12) r . (6.5)The Einstein equations in the null slicing were obtained in [55] and were later used to simulatethe gravitational collapse followed by the formation of a black hole [56, 57] and recently to simulatethe PBH formation as well [48, 58]. We used numerical techniques similar to those used in [57, 58].The fundamental equations are as follows: U = 1 f R u , (6.6)1 f M u = − πR P U, (6.7) E u = − P (cid:18) ν (cid:19) u , (6.8) b = 14 πνR , (6.9)1 f U u = − (cid:18) π Γ R w P ′ + M + 4 πR PR (cid:19) − (cid:18) πνR U ′ + 2 U Γ R (cid:19) , (6.10)1 f (cid:18) ν (cid:19) u = 1 ν Γ (cid:18) U Γ R + 4 πνR U ′ − f U u (cid:19) , (6.11)1 b (cid:18) Γ + Uf (cid:19) ′ = − πR ρ + Pf , (6.12)where the subscript u denotes differentiation with respect to u . Let us rewrite these equationsin terms of the barred variables, defined previously. To do this, it is better to rewrite the u derivatives by the t derivatives since the barred variables were defined by factoring out S ( t ) and ρ ( t ), depending only on t . For example, using (6.5) and the definitions of the barred variables,(6.7) can be rewritten as ( S ρ ) · + S ρ ˙¯ M = − πS ˙ Sρ ¯ a ¯ R ¯ P ¯ U , (6.13)so ˙¯ M = − π αt ¯ a ¯ R ¯ P ¯ U − ( 3 αt − t ) ¯ M , (6.14)which then yields ¯ M τ = − π ¯ a ¯ R ¯ P ¯ U + 3 γ ¯ M , (6.15) ith the help of (2.5) and (3.28). Similarly, (6.6) leads to¯ R τ = ¯ a ¯ U − ¯ R, (6.16)and also, noting (3.28) and (6.5), (6.8) gives¯ E τ = − ¯ P (cid:18) ν (cid:19) τ . (6.17)Equation (6.9) leads to ¯ b = 14 π ¯ ν ¯ R . (6.18)Using (6.4) and the definitions of the barred variables, (6.11) can be transformed to give1¯ a (cid:18) S ¯ ν (cid:19) · = S ¯ ν ¯Γ S ¯ U ¯Γ S ¯ R + 1¯ a ¨ S ¯ U + ˙ SS ¯ b ¯ U ′ − a ( ˙ S ¯ U ) · ! , (6.19)which can be further rewritten using (6.5) as(log ¯ ν ) τ = 3 − a ¯ U ¯ R + 1¯Γ (cid:18) − ¯ a ¯ b ¯ U ′ + ˙ S ∂ ¯ U∂τ (cid:19) . (6.20)Using (2.5), (3.23), (6.4) and (6.5), (6.10) eventually gives¯ U τ = − ¯ U α − α −
32 ¯ a (cid:20) π ¯Γ ¯ R ¯ w (cid:18) S ¯ P ′ − α S ¯ b ¯ a ¯ P (cid:19) + 38 π ¯ M + 4 π ¯ R ¯ P ¯ R + 13 (cid:18) π ¯ νR (cid:18) ¯ b ¯ a α − α ¯ U + 1˙ S ¯ U ′ (cid:19) + 2 1˙ S ¯ U ¯Γ¯ R (cid:19)(cid:21) . (6.21)From the definitions of the barred variables and the Friedmann equation (3.23), (6.12) can berewritten as ¯Γ + ˙ S ¯ Uf ! ′ = −
32 ˙ S ¯ b ¯ R ¯ ρ + ¯ Pf . (6.22)The equations in terms of the barred variables shown above include ˙ S , which can be calculatedby ˙ S = 1 / q ǫr , the relation obtained from (2.12). Observe that t , the time variable in the cosmictime slicing, depends on both u and r ( t = t ( u, r )) in the null slicing, so do ǫ and τ . Then it isnecessary to determine ǫ and τ at a point ( u, r ). To this end, note that on the surface r =const.,we find dt ( u, r ) = f ( u, r ) a ( u, r ) du = f ( u, r ) a ( u, r ) dt ( u, r = ∞ ) , (6.23)which can be rewritten as dτ ( u, r ) = f ( u, r ) a ( u, r ) (cid:18) ǫ ( u, r = ∞ ) ǫ ( u, r ) (cid:19) /β dτ ( u, r = ∞ ) . (6.24)Here the definition of ǫ (2.12) and that of τ (3.27) have been used. The time step dτ ( u, r = ∞ )needs to be chosen to satisfy the CFL condition (see the Appendix), and the general dτ ( u, r ) canbe calculated from (6.24).The boundary conditions are imposed the same as in the cosmic time slicing by setting ¯ u = ¯ R =¯ M = 0 , Γ = 1 at the centre and ¯ a = ¯ ρ = 1 on the outer boundary so that the numerical solutioncoincides with the background FLRW solution there.We now present results of numerical computations using the null slicing. The hypersurfaces of u =const., corresponding to null geodesics, are shown in the top left panel of Figure 9. Observe hat the intervals between the null geodesics are tiny in the central region, meaning that time doesnot pass there in effect. Therefore, the formation of a singularity can be avoided in this slicingas expected. The upper lines in this figure correspond to the null geodesics of the hypotheticalphotons which are emitted from the centre at later times and feel the effects of stronger gravity,so that they need more time to reach a distant observer. In this figure there is an envelope curveof the null geodesics, which shows approximately the location of the apparent horizon. In this waythe time evolution is computed only outside the apparent horizon, so the eventual mass of a PBHcan be determined without facing a singularity. From the same figure, the apparent horizon radiuscan be confirmed to asymptote to a constant value after its formation. This means that the blackhole mass asymptotes to a constant value because R = 2 M on the apparent horizon, and thisbehavior of the mass can be confirmed by the converging curves of the mass profile in the bottomleft panel of Figure 9.The flatness of the mass profile in later time can be understood by noting that the energy densityin a region away from the centre decreases due to the expansion of the universe and also due to theexistence of an underdense region surrounding the central overdense region. As mentioned earlier,this behaviour of the mass profile is similar to that of the vacuum inhabited by a star at the centre,in which case the mass profile is a monotonically increasing function in r inside the star and is flatoutside.The time evolution of 2 M/R is shown in the bottom right panel of Figure 9. The fact thatthe height of the peak of 2
M/R saturates to unity in a late time in the null slicing is clearly seen,which is contrasted with the cosmic time slicing in which the peak is confirmed to exceed unity(see Figure 3). This feature of the time evolution being frozen near the centre in the null slicingcan also be confirmed by the behaviour of the lapse function, which is shown in the top right panelof Figure 9. Note also that the spatial derivatives of f become large near the apparent horizon,which shows the necessity of Adoptive-Mesh-Refinement(AMR), described in the Appendix. ThusFigure 9 as a whole shows that using the null slicing the Einstein equations can be solved until asufficiently late time when the eventual mass of a PBH can be determined.In this paper the mass of a PBH for one particular profile was shown just to illustrate how themass can be obtained by null slicing. In our future work, as was mentioned in the Introduction,we will investigate the mass evolution and the eventual mass of PBHs for a wide variety of profilesconsidered in this paper. VII. CONCLUSION
In this paper we have presented the results of numerical computations of the time evolution ofa perturbed region after the horizon re-entry. The initial conditions for these numerical compu-tations were given using an analytical asymptotic expansion technique developed in our previouspaper. By calculating the time evolution of various initial perturbations, the condition for PBHformation has been investigated. We have extended preceding analyses by performing many morenumerical computations of PBH formation based on the initial curvature profiles characterized byfive parameters which not only reproduce the variety of profiles near the centre but also incorporatethe possible extended features in the tail region (see eq.(5.1)).We have shown that the criterion of PBHs formation can still be expressed in terms of twocrucial (master) parameters which correspond to the averaged amplitude of over density in thecentral region and the width of transition region at outer boundary. As is shown in Figure 6, this IG. 9: Examples of time evolutions of the relevant quantities in the null slicing for the casewhere (
B, σ ) = (1 , . R , normalized bythe Hubble radius at the horizon crossing. The u =const. surfaces are shown in the top left panel.Shown in the top right panel is the lapse function, which goes to zero near the centre. The mass M normalized by the horizon mass at the time of the horizon crossing is shown in the down leftpanel. The profile of 2 M/R , the height of which asymptotes to unity, is shown in the down rightpanel. Arrows in each figure indicate the evolution in time sequence.is the case even though our profiles are characterized by as many as five parameters. We have alsoprovided a reliable physical interpretation of the two-parametric criterion.Using a null slicing we have calculated the variation of PBH mass as a result of accretion.In this paper we present only one example corresponding to K i ( r ) given in Figure 9. For thisexample, the eventual mass of the PBH is approximately equal to the horizon mass at the time ofhorizon crossing. The initial configuration in this example lies at the point marked by the yellowstar in Figure 6, which is fairly far from the critical line. However, it has been shown that theeventual mass of a PBH can be much smaller when ¯ δ hc is extremely close to δ min , the minimumvalue required to create a black hole [48, 59], which would fall on the critical line in Figure 6. Inour future work we will explore this issue for the various types of K i ( r ) investigated in this paper.Note that PBHs are formed only from extremely high peaks of perturbation, corresponding tothe tail of the probability distribution of primordial perturbation. The profiles corresponding tothese peaks have been calculated [50, 51] and turned out to be approximately spherically symmetricand monotonic near the peak. So in this paper we have considered only the spherically symmetricprofiles which are monotonic near the centre. However, deviation from sphericity to some extentis expected so we will explore effects of non-sphericity in our future work as well. CKNOWLEDGMENTS
This work was partially supported by JSPS Grant-in-Aid for Scientific Research 23340058 (J.Y.),Grant-in-Aid for Scientific Research on Innovative Areas No. 21111006 (J.Y.), Grant-in-Aid forExploratory Research No. 23654082(T.H.), and Grant-in-Aid for JSPS Fellow No. 25.8199 (T.N.).TN thanks School of Physics and Astronomy, Queen Mary College, University of London for hospi-tality received during this work. We thank B. J. Carr for useful communications. TN acknowledgesH. Kodama, K. Kohri, K. Ioka and H. Takami for helpful comments.
Appendix
In this Appendix, a few details of the techniques used in the numerical computation of thispaper are discussed.
A. Determination of time steps satisfying Courant-Friedrichs-Lewy condition
In order to guarantee numerical stability, the following Courant-Friedrichs-Lewy condition (CFLconditionjhas to be maintained during the entire computation: c s < ∆ r ∆ t , (7.1)where c s is the sound velocity, which is 1 / √ r and ∆ t represent intervals of gridpoints and timesteps respectively. In general relativity, this CFL conditionis modified as follows: c s < b ( t, r )∆ ra ( t, r )∆ t , (7.2)which is rewritten in terms of the quantities introduced previously as∆ τ < √ − α )¯ b ( t, r ) βr i √ ǫ ¯ a ( t, r ) ∆ r. (7.3)In the null slicing (7.2) is rewritten to be c s < b ( u, r )∆ rf ( u, r )∆ u + b ( u, r )∆ r . (7.4)This equation gives a condition for the time interval∆ u < − c s c s b ( u, r ) f ( u, r ) ∆ r, (7.5)which can be reexpressed noting ∆ u = ∆ t ( u, r = ∞ ) and (3.28) as a condition for ∆ τ ( u, r = ∞ ):∆ τ ( u, r = ∞ ) < √ − − α ) βr i ¯ b ( u, r ) f ( u, r ) p ǫ ( u, r = ∞ ) ∆ r. (7.6)The time step ∆ τ ( u, r = ∞ ) has to be chosen to satisfy this condition and general ∆ τ ( u, r ) canbe calculated from (6.24). . Numerical techniques employed in the calculation in the null slicing The following three numerical techniques are employed to minimize the amount of calculationswhile maintaining accuracy.The first one is to stop computing the time evolution in the central region where time is frozen.As we have seen, f → dτ → u the outer boundaryis placed at r = r b ( u ), where the following inequality is satisfied: | ¯ ρ ( u , r ) − | < c , (7.7)with c being a constant which is sufficiently small (say, 10 − ). Here, it is important to notethat the central perturbed region gradually expands outward, as is shown in Figure 3 for example.As a result, at a later time u ( > u ), r b ( u ) < r b ( u ) so (7.7) would not hold at u if we took r b = r b ( u ). If on the other hand r b is fixed to be r b ( u f ), where u f is the final time of thenumerical computation, (7.7) always holds during the computation. However, this choice leads toan unnecessarily large amount of calculations. This situation can be ameliorated by changing r b optimally by taking the following procedures. Suppose (7.7)is satisfied at u so r b is chosen to be r b = r b ( u ) and later (7.7) breaks at u . Then, at this moment r b is redefined to be r b = r b ( u ) andthe numerical solution in the region r b ( u ) < r < r b ( u ) is approximated by the background FLRWsolution and appended to the solution in the region r < r b ( u ), already calculated numerically.After u , the differential equations are solved in the region r < r b ( u ). And if (7.7) breaks againat u = u , r b is redefined to be r b = r b ( u ) and the solution is approximated by the backgroundsolution in the region r b ( u ) < r < r b ( u ), which is then appended to the numerical solution inthe region r < r b ( u ). This procedure is repeated if necessary, which guarantees (7.7) during theentire computation and as a result enables a computation with a minimal number of grids whilemaintaining accuracy in the region far from the centre .The third technique is the reduction of grid intervals in the region where the spatial derivativesbecome large, a technique known as Adoptive-Mesh-Refinement(AMR). As the perturbation grows,spatial derivatives of physical quantities become large in the central region, which can be confirmedfrom the top right panel of Figure 9 for instance. In order to maintain accuracy in such a situation,the grid intervals need to be reduced but doing so globally would drastically increase the numberof grid points. Therefore, it is optimal to reduce the grid intervals in the central region only wherethe spatial derivatives become large, while the grid intervals in the region away from the centre arekept relatively large, since the spatial derivatives are modest there.[1] Y. B. Zel’dovich and I. D. Novikov, Sov. Astron. , 602 (1967).
2] S. Hawking, Mon. Not. Roy. Astron. Soc. , 75 (1971).[3] S. Hawking, Nature , 30 (1974).[4] Y. B. Zel’dovich, A. A. Starobinskii, M. Y. Khlopov, and V. M. Chechetkin, Sov. Astron.Lett. , 110 (1977).[5] I. D. Novikov, A. G. Polnarev, A. A. Starobinskii, and Y. B. Zel’dovich, Astron. Astrophys. , 104 (1979).[6] B. V. Vainer and P. D. Naselskii, Astron. Zh. , 231 (1978), [Sov. Astron. , 138 (1978).].[7] B. V. Vainer, O. V. Dryzhakova, and P. D. Naselskii, Pis ma Astronomicheskii Zhurnal , 344(1978), [Sov. Astron. Lett. , 185 (1978).].[8] S. Miyama and K. Sato, Prog. Theor. Phys. , 1012 (1978).[9] K. Kohri and J. Yokoyama, Phys. Rev. D , 023501 (2000), astro-ph/9908160.[10] D. N. Page and S. Hawking, Astrophys. J. , 1 (1976).[11] J. H. MacGibbon, Nature , 308 (1987).[12] J. H. MacGibbon and B. J. Carr, Astrophys. J. , 447 (1991).[13] B. J. Carr, Astrophys. J. , 8 (1976).[14] B. Paczynski, Astrophys. J. , 1 (1986).[15] R. Saito and J. Yokoyama, Phys. Rev. Lett. , 161101 (2009), [ , 069901(E) (2011).],0812.4339.[16] R. Saito and J. Yokoyama, Prog. Theor. Phys. , 867 (2010), [ , 351(E) (2011).],0912.5317.[17] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D , 104019 (2010), 0912.5297.[18] K. Griest, A. M. Cieplak, and M. J. Lehner (2013), 1307.5798.[19] A. G. Polnarev, in Morphological Cosmology , edited by P. Flin and H. W. Duerbeck (1989),vol. 332 of
Lecture Notes in Physics, Berlin Springer Verlag , pp. 369–376.[20] K. Sato, Mon. Not. Roy. Astron. Soc. , 467 (1981).[21] A. H. Guth, Phys. Rev. D , 347 (1981).[22] A. A. Starobinsky, Phys. Lett. B , 99 (1980).[23] V. F. Mukhanov and G. Chibisov, Sov. Phys. JETP , 258 (1982).[24] A. H. Guth and S. Pi, Phys. Rev. Lett. , 1110 (1982).[25] S. Hawking, Phys. Lett. B , 295 (1982), revised version.[26] A. A. Starobinsky, Phys. Lett. B , 175 (1982).[27] G. Hinshaw et al. (WMAP Collaboration) (2012), 1212.5226.[28] P. Ade et al. (Planck Collaboration) (2013), 1303.5082.[29] J. Garcia-Bellido, A. Linde, and D. Wands, Phys. Rev. D , 6040 (1996).[30] H. M. Hodges and G. R. Blumenthal, Phys. Rev. D , 3329 (1990).[31] P. Ivanov, P. Naselsky, and I. Novikov, Phys. Rev. D , 7173 (1994).[32] J. Yokoyama, Astron. Astrophys. (1997).[33] J. Yokoyama, Phys. Rev. D , 083510 (1998).[34] J. Yokoyama, Phys. Rep. , 133 (1998).[35] M. Kawasaki and T. Yanagida, Phys. Rev. D , 043512 (1999).[36] J. Yokoyama, Prog. Theor. Phys. Suppl. , 338 (1999).[37] R. Saito, J. Yokoyama, and R. Nagata, J. Cosmol. Astropart. Phys. , 024 (2008).[38] A. Taruya, Phys. Rev. D , 103505 (1999).[39] B. A. Bassett and S. Tsujikawa, Phys. Rev. D , 123503 (2001).[40] A. M. Green and K. A. Malik, Phys. Rev. D , 021301 (2001).[41] M. Kawasaki, T. Takayama, M. Yamaguchi, and J. Yokoyama, Mod. Phys. Lett. A22 , 1911 , 1426 (2008), 0711.3886.[43] B. J. Carr and S. Hawking, Mon. Not. Roy. Astron. Soc. , 399 (1974).[44] B. J. Carr, Astrophys. J. , 1 (1975).[45] M. Shibata and M. Sasaki, Phys. Rev. D , 084002 (1999), gr-qc/9905064.[46] A. G. Polnarev and I. Musco, Class. Quant. Grav. , 1405 (2007), gr-qc/0605122.[47] D. K. Nadezhin, I. D. Novikov, and A. G. Polnarev, NASA STI/Recon Technical Report N , 10983 (1979).[48] J. C. Niemeyer and K. Jedamzik, Phys. Rev. D , 124013 (1999), URL http://link.aps.org/doi/10.1103/PhysRevD.59.124013 .[49] A. Polnarev, T. Nakama, and J. Yokoyama, J. Cosmol. Astropart. Phys. , 027 (2012),URL http://stacks.iop.org/1475-7516/2012/i=09/a=027 .[50] A. Doroshkevich, Astrophysics , 320 (1970), ISSN 0571-7256, URL http://dx.doi.org/10.1007/BF01001625 .[51] J. M. Bardeen, J. Bond, N. Kaiser, and A. Szalay, Astrophys.J. , 15 (1986).[52] C. W. Misner and D. H. Sharp, Phys. Rev. , B571 (1964).[53] M. M. May and R. H. White, Academic Press, New York and London (1967).[54] I. Musco and J. C. Miller, arXiv:1201.2379 [gr-qc] (2012).[55] W. C. Hernandez, Jr. and C. W. Misner, Astrophys. J. , 452 (1966).[56] J. C. Miller and S. Motta, Class. Quant. Grav. , 185 (1989), URL http://stacks.iop.org/0264-9381/6/i=2/a=012 .[57] T. W. Baumgarte, S. L. Shapiro, and S. A. Teukolsky, Astrophys. J. , 717 (1995).[58] I. Musco, J. C. Miller, and L. Rezzolla, Class. Quant. Grav. , 1405 (2005), URL http://stacks.iop.org/0264-9381/22/i=7/a=013 .[59] I. Musco, J. C. Miller, and A. G. Polnarev, Class. Quant. Grav. , 235001 (2009), 0811.1452.[60] T. Harada, C.-M. Yoo, and K. Kohri (2013), 1309.4201., 235001 (2009), 0811.1452.[60] T. Harada, C.-M. Yoo, and K. Kohri (2013), 1309.4201.