Accretion onto a small black hole at the center of a neutron star
AAccretion onto a small black hole at the center of a neutron star
Chloe B. Richards, Thomas W. Baumgarte, and Stuart L. Shapiro
2, 3 Department of Physics and Astronomy, Bowdoin College, Brunswick, ME 04011 Department of Physics, University of Illinois at Urbana-Champaign, IL 61801 Department of Astronomy and NCSA, University of Illinois at Urbana-Champaign, IL 61801
We revisit the system consisting of a neutron star that harbors a small, possibly primordial,black hole at its center, focusing on a nonspinning black hole embedded in a nonrotating neutronstar. Extending earlier treatments, we provide an analytical treatment describing the rate of secularaccretion of the neutron star matter onto the black hole, adopting the relativistic Bondi accretionformalism for stiff equations of state that we presented elsewhere. We use these accretion ratesto sketch the evolution of the system analytically until the neutron star is completely consumed.We also perform numerical simulations in full general relativity for black holes with masses up tonine orders of magnitude smaller than the neutron star mass, including a simulation of the entireevolution through collapse for the largest black hole mass. We construct relativistic initial datafor these simulations by generalizing the black hole puncture method to allow for the presence ofmatter, and evolve these data with a code that is optimally designed to resolve the vastly differentlength scales present in this problem. We compare our analytic and numerical results, and provideexpressions for the lifetime of neutron stars harboring such endoparasitic black holes.
I. INTRODUCTION AND SUMMARY
A number of authors have discussed the prospect ofusing neutron stars as dark matter detectors (see, e.g.,[1–8] and references therein). The idea is that the darkmatter may be in the form of small black holes that maybe captured by neutron stars, or may be in the form ofother particles that may coalesce to form black holes inthe interior of neutron stars. In either case, neutron starsmay end up harboring a small “endoparasitic” black holeat their center. Such a black hole will then accrete thesurrounding neutron star matter until the entire star hasbeen consumed. The current existence of neutron starpopulations therefore can be used to constrain the natureof black holes as candidates for dark matter.In principle, the above scenario could be realized in atleast two different ways. It is possible that primordialblack holes (PBHs) formed in the early Universe (see,e.g., [9, 10]), and that they contribute to, or even accountfor, the dark matter. Such PBHs can be captured bystars. Following capture, the black holes can settle closeto the stellar center, and subsequently accrete the entirestar (see, e.g., [9, 11]). The capture, settling and accre-tion process are particularly efficient for neutron stars(see, e.g., [5, 8], but see also [12]), so that the currentexistence of certain neutron star populations constrainsthe density of PBHs in certain mass ranges. In partic-ular, these arguments have been used to establish limitsin the mass range of 10 − M (cid:12) (cid:46) M BH (cid:46) − M (cid:12) [5],which is otherwise only poorly constrained (cf. [13–15];see also [16, 17] and references therein for constraintsarising from gravitational wave observations).Alternatively, it is possible that dark matter particlesare captured inside a neutron star. Under sufficientlyfavorable conditions, these particles may form a high-density object that may then collapse to form a smallblack hole (see, e.g., [1–3, 6]). The authors of [6], forexample, estimated the black hole mass to be approxi- mately M BH (cid:46) − M (cid:12) , similar to the mass range ofinterest for PBHs.Once a black hole has either been captured by the neu-tron star, or has formed inside the neutron star, andafter this black hole has settled to the center, it willaccrete the neutron star matter, ultimately consumingthe entire star. Observational signatures of this processhave been discussed by a number of authors, including[8, 18, 19], while numerical simulations have been per-formed recently by [7]. The authors of [7] consideredthree different equations of state (EOSs), as well as bothrotating and non-rotating neutron stars, but limited theirsimulations to relatively large black hole masses with M BH /M ≥ − , where M is the neutron star mass.They found that the accretion rate follows the relation˙ M BH ∝ M as suggested by Bondi accretion ([20]; seealso [21], hereafter ST, for a textbook treatment, includ-ing its relativistic generalization), and that the accretionrate is largely independent of the neutron star spin, inagreement with [22].In this paper we extend earlier work on the accretionof neutron stars by endoparasitic black holes in severalways. In Section II we provide an analytical overview de-scribing the accretion process in the core of a nonrotatingneutron star. Building on our earlier results [23, 24], weadopt a relativistic Bondi accretion model to describea nonspinning black hole accreting gas obeying a stiff,polytropic EOS to track the evolution of the system. Asdiscussed in [23], Bondi accretion for stiff EOSs with adi-abatic indicies Γ > / M BH ∝ M .We integrate the equations, crudely accounting for thequasistatic evolution of the neutron star during the ac-cretion process, to calculate the neutron star’s survivaltime.In Section III we perform numerical simulations of theTypeset by REVTEX a r X i v : . [ a s t r o - ph . H E ] F e b M BH / M M * B H Black hole growthRest-mass fluxBondi accretionMinimum accretion 10 M * B H [ M / y r ] FIG. 1. The analytical and numerical values of the accretionrate ˙ M (cid:63) BH as a function of M BH /M , computed for our fidu-cial Γ = 2 neutron star model (see Table II for details). Onthe left we express the accretion rate in geometrized units (interms of which ˙ M (cid:63) BH is dimensionless), while on the right weexpress this rate in units of solar masses M (cid:12) per year. Thesolid (blue) line represents the analytical solution (44) basedon Bondi accretion for stiff EOSs, while (red) triangles and(green) circles represent numerical measurements of the ac-cretion rate based on the growth of the black hole horizonand on the rest-mass flux, respectively (see Sections III C 2and III C 3). The dashed (red) line represents the minimumaccretion rate (45) for a Γ = 2 polytrope. (See also Table IIIfor a detailed listing of these results.) accretion process in full general relativity. Supportedby the results of [7] we focus on nonrotating, spheri-cal neutron stars and nonspinning black holes only, andadopt, for simplicity, a stiff, polytropic EOS with Γ = 2.Such a value is often used in numerical simulations to ap-proximate a stiff nuclear EOS. Using a numerical relativ-ity code implemented in spherical polar coordinates (see[25, 26]) with a logarithmic radial coordinate, we extendthe results of [7] by simulating the accretion onto blackholes with mass ratios as small as M BH /M (cid:39) − . Thisallows our simulations to extend into the mass range ofinterest from the perspective of PBHs and viable blackhole dark matter candidates, as discussed above. Suchblack hole masses thus extend down to the range of dwarfplanets.We compare our analytical and numerical results inSection IV. In particular, we compile accretion rates ob-tained both analytically and numerically for a large rangeof black hole masses in Table III. We also plot these re-sults in Fig. 1, which serves as a summary of our results.In this figure, the solid line represents the analytical ac-cretion rate quoted in Eq. (44). The dashed line showsthe minimum accretion rate for Γ = 2 as identified by[23] (see their Eq. 51). Such a minimum only arises forstiff EOSs, but is surprisingly insensitive to the EOS forΓ >
2. The open circles and filled triangles mark two dif-ferent methods for computing the accretion rate in our numerical simulation data. The former result from eval-uating the accretion rate directly from the flux of mattercrossing the black hole horizon, while the latter are basedon measuring the rate at which the area of the horizonincreases. Since the black hole grows at a rate that isproportional to the square of the black hole mass, andhence is exceedingly slow for small black hole masses, wecan accurately measure the horizon’s growth rate numer-ically only for the largest black hole masses. We note thatthe proportionality ˙ M BH ∝ M holds over many ordersof magnitude, that our numerically determined accretionrates agree well with those computed analytically fromthe relativistic Bondi accretion relation (indicating thatwe have correctly identified the constant of proportional-ity), and finally that the actual accretion rates are onlyslightly larger than the minimum accretion rate identifiedin [23]. All of the above confirms that the consumptionof a neutron star by an endoparasitic black hole is welldescribed by Bondi accretion for stiff EOSs, so that theseanalytical expressions can be used to establish the life-time of such neutron stars (see Section II C).Throughout this paper we use geometrized units with G = 1 = c unless noted otherwise. II. ANALYTICAL OVERVIEWA. Accretion rates
We assume that a nonspinning black hole with mass M BH is located at the center of a spherical neutron starwith mass M (cid:29) M BH and radius R (cid:29) M BH that wouldbe in strict hydrostatic equilibrium in the absence of theblack hole. We take the matter to be at rest far fromthe black hole, which resides at the center of a nearlyhomogeneous core. To estimate the rate at which theblack hole accretes matter from the neutron star we firstdefine the capture radius r a ≡ M BH a c . (1)Here a is the speed of sound, given by a = dPdρ (cid:12)(cid:12)(cid:12)(cid:12) s = dPdρ (cid:12)(cid:12)(cid:12)(cid:12) s ρ ρ + P , (2)and a c is its value evaluated at the center of the un-perturbed star. In (2), P is the pressure, ρ the totalmass-energy density, ρ the rest-mass density, and thederivative is taken at constant entropy. For the New-tonian estimates that we will provide in this section wewill not need to distinguish between ρ and ρ , but weintroduce them here for later reference. For a polytropicequation of state with P = Kρ Γ0 , (3)where Γ is the adiabatic index and K a constant, we maythen approximate a c (cid:39) (cid:18) Γ P c ρ c (cid:19) / (cid:39) (cid:18) Γ MR (cid:19) / . (4)Here P c and ρ c are the values of the unperturbed star’spressure and density at the center, with P c (cid:28) ρ c andwhere we have assumed hydrostatic equilibrium in writ-ing the last (rough) equality.The accretion flow onto the black hole can be approx-imated in two opposite limiting regimes, depending onwhether the neutron star mass m ( r a ) contained withinthe capture radius r a is greater or smaller than the blackhole mass M BH . In the latter case we may ignore theself-gravity of the neutron star matter, so that the ac-cretion process is described by secular Bondi accretion,i.e. adiabatic flow with asymptotically constant matterdensity and pressure and zero flow velocity ([20, 27]; seealso ST for a textbook treatment). In the former case wecannot ignore the self-gravity, and the evolution modeledas an accretion process becomes catastrophic dynamicalcollapse.Defining m ( r ) as the neutron star mass within the ra-dius r , so that M = m ( R ), we have m ( r ) (cid:39) π ρ c r (5)for sufficiently small r . We then compute the crucialmass ratio, m ( r a ) M BH (cid:39) π M ρ c a c (cid:39) π M ρ c R Γ M . (6)We now write ρ c = δ ¯ ρ = δ M πR , (7)where ¯ ρ is the unperturbed star’s mean density, and thefactor δ measures its central concentration, ρ c / ¯ ρ . Wetabulate values of δ for Newtonian polytropes of index n and adiabatic index Γ = 1 + 1 /n in Table I. Inserting (7)into (6) we now have m ( r a ) M BH (cid:39) δ Γ (cid:18) M BH M (cid:19) . (8)We see that m ( r a ) ∼ M BH only for large black holemasses M BH ∼ M and soft EOSs that result in largevalues of δ/ Γ , also listed in Table I. For black holes thatstart out with masses much smaller than the neutronstar mass, M BH (cid:28) M , almost the entire accretion pro-cess (i.e. the longest duration) will occur in the regime m ( r a ) (cid:28) M BH , and will therefore be described by qua-sistatic Bondi accretion, while only the short final epochproceeds dynamically in the regime m ( r a ) ∼ M BH . n Γ δ δ/ Γ λ s δ as well as thecombination Γ /δ for Newtonian polytropes with polytropicindex n and adiabatic index Γ = 1 + 1 /n . For Γ ≤ / λ s ; see Eq. (10). ForΓ > /
3, the relativistic accretion eivenvalues λ GR can beconstructed as in [23]; see also Section II D below.
1. Case I: m ( r a ) (cid:28) M BH – Bondi accretion In this case we can neglect the self-gravity of the neu-tron star fluid inside the capture radius r a , since the grav-itational forces are dominated by those exerted by theblack hole. In this case, the accretion is well describedby adiabatic Bondi accretion [20], the rate for which isgiven by ˙ M (cid:63) BH = − ˙ M (cid:63) = 4 πλ GR (cid:18) M BH a (cid:63) (cid:19) ρ (cid:63) a (cid:63) . (9)Here λ GR is a dimensionless “accretion eigenvalue”, typ-ically of order unity, and the (cid:63) symbols denote values asobserved by a “local asymptotic” static observer who isfar from the black hole, but still well inside the neutronstar, i.e. M BH (cid:28) r (cid:28) R . The dot in the accretion ratesdenotes a derivative with respect to time as measuredby such an observer (see also Appendix A for a detaileddiscussion). It is assumed that the density ρ (cid:63) and soundspeed a (cid:63) approach constants, and the flow speed u (cid:63) ap-proaches zero, as r > r a becomes large in this asymptoticregion, which typically resides inside the nearly homoge-neous core of the neutron star. In Appendix A, as well asSection III C 2 below, we discuss how this “local” accre-tion rate is related to the rate of mass accretion as seenby an observer far from the star. Here we simply pointout that the accretion of bound matter is not influencedby the spherical matter distribution beyond the radiusat which it becomes bound (Birkhoff’s theorem): whenself-gravity of the matter bound to the black hole can beneglected, the only gravitational source influencing theaccretion is the central black hole.The accreted mass measured by (9) is fundamentallya (baryon) rest mass. It enhances the black hole’s totalgravitational mass by a similar amount for strictly adi-abatic flow as long as the asymptotic internal energy ofthe gas is small in comparison to the rest-mass energy.A relativistic treatment of accretion for nonspinningblack holes shows that the requirement that the soundspeed be less than the speed of light demands that theflow become transonic and pass through a critical point,yielding a unique value for λ GR , as shown by ST. For softequations of state with Γ ≤ / a (cid:63) (cid:28) λ GR → λ s = 14 (cid:18) − (cid:19) (3Γ − / − (10)(see, e.g., Chapter 14, Eq. (14.3.17), and Table 14.1 inST). Unlike λ GR , the values of λ s are independent of a (cid:63) and a s , the sound speed at the transonic radius r s (see [23] for details). Stiffer equations of state require arelativistic treatment as discussed in [23].We now assume that r a is sufficiently small, r a (cid:28) R ,that we can approximate the fluid variables as seen by thelocal observer discussed above, i.e. ρ (cid:63) and a (cid:63) , by thoseat the center of the unperturbed star, i.e. ρ c and a c . Theaccretion rate then becomes˙ M (cid:63) BH = − ˙ M (cid:63) = 4 πλ GR M a c ρ c , (11)from which we can crudely estimate the accretiontimescale τ acc τ acc ≡ M BH ˙ M (cid:63) BH (cid:39) a c πλ GR M BH ρ c = Γ / δλ GR M / R / M BH . (12)Here we have used (4) and (7) in the last step. Theabove estimates the timescale for the black hole to dou-ble its mass, which is the bottleneck in the neutron starconsumption process: this epoch takes the longest be-cause both the black hole mass and hence the accretionrate (11) are the smallest they will be during the process.Dividing by the neutron star mass, we may rewrite thisresult as τ acc M (cid:39) Γ / δλ GR (cid:18) RM (cid:19) / (cid:18) MM BH (cid:19) . (13)Note that τ acc /M → ∞ as M BH /M →
0. Alternatively,we may also express the accretion timescale in terms ofthe neutron star’s dynamical (collapse) timescale τ dyn (cid:39) γ (4 πρ c / / = γδ / (cid:18) RM (cid:19) / M, (14)where γ is a factor of order unity and where we have used(7). Combining (12) and (14) we obtain τ acc τ dyn (cid:39) Γ / δ / λ GR γ (cid:18) MM BH (cid:19) . (15)We again have τ acc /τ dyn → ∞ as M BH /M →
0. We willcalculate τ acc more carefully in Section II C below.
2. Case II: m ( r a ) ∼ M BH – Dynamical accretion In this case we cannot neglect the self-gravity of theneutron star fluid inside the capture radius r a . We now generalize the definition (1) of this capture radius, anddefine a critical radius, r crit = m ( r crit ) + M BH a c , (16)inside which the gas is marginally bound by the combinedmass of the black hole and the gas itself. Using (5), wecan rewrite (16) as4 π ρ c r (cid:18) M BH m ( r crit ) (cid:19) = a c . (17)The rate at which the black hole accretes mass can nowbe expressed as the area of the sphere with the critical ra-dius, 4 πr , times the mass flux across this sphere, ρ c u c .Assuming that, at the critical radius, the fluid speed u c is comparable to the sound speed a c , as in typical Bondiflows, we then have˙ M (cid:63) BH (cid:39) πr ρ c a c = 3 a c (cid:18) M BH m ( r crit ) (cid:19) − , (18)where we have used (17) in the last equality. The corre-sponding accretion timescale is then given by τ acc = M BH ˙ M BH (cid:39) M BH / (cid:18) RM (cid:19) / (cid:18) M BH m ( r crit ) (cid:19) , (19)or τ acc τ dyn (cid:39) δ / γ Γ / M BH M (cid:18) M BH m ( r crit ) (cid:19) , (20)where we have approximated the dynamical timescale τ dyn as in (14).We now evaluate Eq. (20) in two limits. In the limit M BH ∼ m ( r crit ), in which case M BH ∼ M by (8), we no-tice that the accretion timescale τ acc becomes comparableto the dynamical timescale τ dyn , as one would expect. Inthe opposite limit, M BH (cid:29) m ( r crit ), the critical radius r crit defined in (16) reduces to r a defined in (1) and wemay approximate M BH m ( r crit ) (cid:39) M BH πρ c r a = 3 a c πρ c M = Γ δ (cid:18) MM BH (cid:19) , (21)where we have used (4) and (7) in the last equality. In-serting (21) into (20) we recover, up to factors of orderunity, the Bondi accretion timescale (15), as expected. B. Effects of stellar evolution
Our simple estimates (15) and (20) for the accretiontimescales ignore the fact that the accretion rates changeas the black hole mass M BH increases, and also ignore thefact that the neutron star structure changes as the accre-tion proceeds. We can approximate the effects of thissecular “stellar evolution” by assuming that, while thestar loses mass to the black hole, it will adjust quasistat-ically to a new equilibrium configuration while keepingits total Newtonian energy E constant. We now writethis energy as the simple functional E = − α M M BH R − − − M R , (22)where the first term accounts for the interaction betweenthe stellar gas and the black hole, with α being a constantthat depends on Γ, α = α (Γ), while the second termdescribes the neutron star’s self-energy (see Eq. 3.3.10 inST).Evaluating (22) at the initial time, denoted by (0), wehave E = − α M (0) M BH (0) R (0) − − − M (0) R (0) . (23)Since, by our assumption, expressions (22) and (23) mustbe identical, we can equate them and solve for R to find R = MM (0) αM BH + (3Γ − M/ (5Γ − αM BH (0) + (3Γ − M (0) / (5Γ − R (0) . (24)We now approximate M BH (cid:28) M , in which case (24)reduces to R (cid:39) (cid:18) MM (0) (cid:19) R (0) . (25)Using (25) in (4) then yields a c (cid:39) (cid:18) Γ M (0) R (0) (cid:19) / (cid:18) M (0) M (cid:19) / , (26)while (7) gives ρ c (cid:39) δ M (0)4 πR (0) (cid:18) M (0) M (cid:19) . (27)Inserting (26) and (27) into the Bondi accretion rate (9)then results in˙ M = − λ GR δ Γ / ( M BH (0) + M (0) − M ) M (0) / R (0) / (cid:18) MM (0) (cid:19) − / (28)where we have expressed the black-hole mass M BH interms of the evolving neutron star mass M as M BH = M BH (0) + M (0) − M. (29)Note that the last factor in (28) accounts for stellar evo-lution. C. Accretion times
We can now compute the neutron star lifetime (i.e. theaccretion time) by integrating Eq. (28). Towards thatend, it is useful to introduce the dimensionless quantities y ≡ M BH (0) + M (0) M (0) , y ≡ MM (0) , (30) and T ≡ λ GR δ Γ / (cid:18) M (0) R (0) (cid:19) / t, (31)in terms of which we may rewrite (28) as dydT = − ( y − y ) y − / . (32)As in (28), the last factor accounts for stellar evolution.
1. Without stellar evolution
We first ignore the effects of stellar evolution, so that(32) reduces to dydT = − ( y − y ) , (33)which can be integrated readily to yield[ T ] fi = (cid:20) y − y (cid:21) fi . (34)Here the square brackets serve as a reminder to insertlimits of integration. At the initial time, which we chooseto be T i = 0, we have y i = 1, so that T f = 1 y f − y − − y = M (0) M BH (0) − M (0) M BH (35)where M BH is the black hole mass at the time T f , andwhere we have used (29). We can now find the total ac-cretion time by setting the final neutron star mass equalto zero, i.e. by choosing y f = 0. Further assuming that M BH (0) (cid:28) M (0) and recalling (31) we find τ acc M (0) = Γ / δλ GR (cid:18) R (0) M (0) (cid:19) / (cid:18) M (0) M BH (0) (cid:19) , (36)which is identical to (13), as expected.
2. With stellar evolution
We now repeat the exercise, but include the last factorin (28) in order to account for stellar evolution. In thiscase the integral can be carried out as described in Ap-pendix B. Choosing, as before, y i = 1 at T i = 0, as wellas y f = 0 in order to obtain the accretion time T f = T acc ,we obtain T acc = 52 + 43 y +6 y + y y − − y / ln (cid:32) y / + 1 y / − (cid:33) (37)Alternatively, we may introduce y h = y − M BH (0) M (0) (38)and rewrite (37) as T acc = 7 y h + 493 y h + 16115 + 1 y h −
72 ( y h + 1) / ln √ y h + 1 √ y h − . (39)Taking the limit y h →
0, we see that T acc will be dom-inated by the term 1 /y h , so that we recover the sameaccretion time t acc as in (36). This is not entirely sur-prising, since most of the accretion time is spent duringearly times when neither the neutron star mass nor ra-dius change appreciably, so that stellar evolution is notimportant. At late times, however, the response of thestar to the accretion process, and the corresponding ad-justments in the stellar structure, will affect the accretiontime, as expressed by (39).As a concrete example, consider a star with M = 1 M (cid:12) .Since, in geometrized units, 1 M (cid:12) (cid:39) . (cid:39) µ s, wethen have τ acc ∼ (cid:18) R (0) M (0) (cid:19) / (cid:18) M (0) M BH (0) (cid:19) − s , (40)where we ignored factors of order unity. For a main-sequence star, with R (cid:39) M , we see that the accretiontime will exceed a Hubble time if M BH (0) (cid:46) − M (cid:12) .For a neutron star, however, R (cid:39) M , resulting in sig-nificantly smaller accretion timescales. Therefore, blackholes with masses as small as M BH (0) (cid:38) − M (cid:12) wouldbe able to consume a neutron star well within a Hubbletime. D. Fiducial neutron star model
In our comparisons with the numerical results of Sec-tion III below we consider, as a fiducial neutron starmodel, a dynamically stable equilibrium star with a cen-tral rest-mass density of ρ c = 0 . K − n governed bya polytropic equation of state (3) with Γ = 2 and n = 1, which we constructed by solving the Tolman-Oppenheimer-Volkoff (TOV) [28, 29] equations. Detailedproperties of this stellar model are listed in Table II.Since K n/ has units of length in geometrized units, wemay introduce non-dimensional quantities by rescalingany dimensional quantity with a suitable power of K ; inparticular, we define˜ ρ ≡ K n ρ, ˜ R ≡ K − n/ R, ˜ ρ ≡ K n ρ , ˜ M ≡ K − n/ M, (41)and similar for other quantities. We list “tilde” variablesthat have been rescaled with respect to K in the sec-ond column in Table II. In the third column we rescaleeach variable with respect to the neutron star’s gravita-tional mass M , while in the fourth column we rescale withrespect to the corresponding maximum-mass configura-tion. In particular we note that, for our adopted model, M/M max = 0 . M max is the maximum gravita-tional mass of a spherical star with our adopted EOS. Fi-nally, for the fifth column in Table II we assume that ourstar has a gravitational mass of M = 1 . M (cid:12) , in whichcase K takes the value K = (1 . M (cid:12) / ˜ M ) (cid:39)
156 km .From the parameters given in Table II we compute thecentral sound speed to be a c = 0 . . (42)We identify, as before, the neutron star’s central densityand sound speed with the corresponding asymptotic val-ues for the Bondi accretion onto the black hole and follow[23] to compute the accretion eigenvalue λ GR = 1 . . (43)Inserting the above values into the Bondi accretion rate(9) we obtain ˙ M (cid:63) BH = 21 .
24 ˜ M , (44)which is not significantly larger than the minimum steadystate accretion rate for a stiff polytrope with Γ = 2,˙ M (cid:63) BH , min = 9 .
29 ˜ M (45)(see Eq. 51 in [23]). Adopting the above value of K , andrecalling that, in geometrized units, M (cid:12) (cid:39) × − s,we can evaluate (45) to yield˙ M (cid:63) BH , min = 7 . × − M (cid:12) yr (cid:18) M BH − M (cid:12) (cid:19) . (46) III. NUMERICAL TREATMENTA. Initial data
We construct relativistic initial data describing a non-spinning black hole embedded at the center of a nonro-tating, spherical neutron star. Our task is to solve theHamiltonian and momentum constraints of general rela-tivity (see, e.g. [30] for discussion and references) whichwe do by generalizing the puncture method (see [31]) toallow for the presence of matter. Our approach differsfrom that adopted by [7], who constructed initial databy matching an interior black hole solution to an exte-rior neutron star solution.We assume that the initial slice is conformally flat andstatic, thereby setting K ij = 0 = S i , which trivially sat-isfies the momentum constraints. Here K ij is the extrin-sic curvature and S i is the 3-momentum density mea-sured by a normal observer n a , i.e. an observer whosefour-velocity is the normal vector n a . The Hamiltonianconstraint then becomes¯ D ψ = − πψ ρ. (47) Quantity Rescaled wrt a K Rescaled wrt M Rescaled wrt max. mass model Physical units ρ c b ˜ ρ c = 0 . M ρ c = 0 . ρ c /ρ max0 c = 0 . ρ c = 3 . × g/cm ρ c c ˜ ρ c = 0 . M ρ c = 0 . ρ c /ρ max c = 0 . ρ c = 4 . × g/cm R d ˜ R = 0 . R/M = 5 . R/R max = 1 . R = 10 . r isoe ˜ r iso = 0 . r iso /M = 4 . r iso /r maxiso = 1 . r iso = 8 .
73 km M f ˜ M = 0 . M/M = 1
M/M max = 0 . M = 2 . × g M ˜ M = 0 . M /M = 1 . M /M max0 = 0 . M = 3 . × g ψ c h ψ c = 1 . ψ c = 1 . ψ c /ψ max c = 0 . ψ c = 1 . α c i α c = 0 . α c = 0 . α c /α max c = 1 . α c = 0 . a With respect to b Central rest-mass density c Central mass-energy density d Areal radius e Isotropic radius f Gravitational mass g Rest mass h Central conformal factor i Central lapse function
TABLE II. Parameters for our fiducial Γ = 2, n = 1 polytropic neutron star model in the absence of a black hole (see textfor details). The conformal factor ψ is well defined by our assumption that, in Cartesian coordinates, the determinant of theconformally related metric is ¯ γ = 1. The lapse function α listed here is the value obtained from integrating the TOV equations,and is different from the 1+log lapse adopted in our numerical evolution calculations (see Eq. 57). Here ψ is the conformal factor, ¯ D the flat Laplace op-erator, and ρ = n a n b T ab the mass-energy density as ob-served by a normal observer. We allow for a conformalrescaling of the density, ρ = ψ m ¯ ρ, (48)where m is a yet-to-be-determined exponent. An attrac-tive choice might be m = −
6, since it leaves the properintegral over the density ρ invariant if we keep ¯ ρ fixed, (cid:90) ψ ρ d x = (cid:90) ¯ ρ d x. (49)Note, however, that this integral represents neither thegravitational nor the rest mass.Now assume that we have constructed a solution to theTOV equations [28, 29] in isotropic coordinates, so thatwe obtain radial profiles of the conformal factor ψ NS andthe mass-energy density ρ NS for the equilibrium neutronstar by itself. In particular, these functions satisfy theHamiltonian constraint (47) with¯ D ψ NS = − πψ ρ NS , (50)and we identify ¯ ρ = ψ − m NS ρ NS . (51)We now want to modify this solution so that the newsolution accounts for a black hole embedded at the cen-ter of the neutron star. Towards that end, we write theconformal factor as ψ = ψ NS + ψ BH + u, (52) where ψ BH = M r (53)is an isolated black hole’s contribution to the conformalfactor in our isotropic coordinates. We refer to M as the“puncture mass”; it has no immediate physical signifi-cance, and serves as a mass parameter only (see SectionIII C 1 and Eqs. 58 and 66 below for the black hole’sisolated horizon, or irreducible, mass). Inserting (52)into the Hamiltonian constraint (47) and observing that¯ D ψ BH = 0 we obtain¯ D ψ NS + ¯ D u = − π ( ψ NS + ψ BH + u ) m ¯ ρ, (54)or, using (50),¯ D u = − π (cid:110) ( ψ NS + ψ BH + u ) m − ψ m NS (cid:111) ¯ ρ. (55)Since ψ BH diverges as r →
0, it may be desirable tochoose m < −
5, which makes the right-hand side of (55)regular. Therefore, unless noted otherwise, we will use m = − u (see alsoFig. 9).We have now reduced the problem to finding regularsolutions u to the elliptic equation (55). Since the equa-tion is non-linear, we adopt an iterative approach. Oncewe have obtained this solution we can compute the new(physical) energy density ρ from ρ = ψ m ¯ ρ = (cid:18) ψ NS + ψ BH + uψ NS (cid:19) m ρ NS . (56) r / M M = 1e-03= 1e-04= 1e-05= 1e-06= 1e-07= 1e-08= 1e-09= 1e-10 FIG. 2. Profiles of the initial rest-mass density ρ as a func-tion of isotropic radius r , for our fiducial neutron star model(see Table II) with different black hole puncture masses M .Here and throughout we choose m = − Note that we will have ρ → r → m < M (cid:28) M . We show examples of density pro-files for black holes with different masses embedded inour fiducial neutron star model in Fig. 2.While our initial data depend on our choice of m inEq. (48), they quickly settle into a quasiequilibrium con-figuration soon after matter marginally bound to theblack hole begins to flow inward at t ∼ r a /a c ∼ M BH .The system thus relaxes to a state of quasistatic accre-tion onto the black hole that is independent of m (seeFigs. 3 as well as the discussion in III B below). B. Numerical evolution
We evolve our initial data with a code that solvesthe Baumgarte-Shapiro-Shibata-Nakamura (BSSN) for-mulation of Einstein’s equations [32–34]. We adopt areference-metric formulation (see, e.g., [35–38]) in orderto implement the equations in spherical polar coordinates(see [25, 26] for details and tests; see also [39–41]). Thelatest version of our code uses fourth-order finite differ-encing for all spatial derivatives in Einstein’s equations,together with a fourth-order Runge-Kutta time integra-tor.We impose coordinates using the “1+log” slicing con- r / M M m = 1 m = 4 m = 6 m = 910 FIG. 3. The rest-mass density ρ as a function of isotropicradius r for different values of the conformal exponent m in (48) for ˜ M = 10 − , so ˜ M BH (0) = 1 . × − , and M BH (0) /M (0) = 8 . × − . Although the initial densityprofiles, shown in the inset, clearly depend on m , they allevolve to the same density profile, shown in the large plot,once a quasi-equilibrium has been reached. Here and in sev-eral of the following graphs we suppress the innermost fewgrid points well inside the horizon, since they are affectedby numerical noise caused by the puncture singularity at theorigin. dition ( ∂ t − β i ∂ i ) α = − αK (57)(see [42]) for the lapse function α , and a “Gamma-driver”condition for the shift vector β i (see [43, 44]). On our ini-tial slice we choose a “pre-collapsed” lapse with α = ψ − and zero shift.For all simulations reported in this paper we used anumerical grid of N r = 512 radial grid points with theouter boundary at ˜ r out = 4, corresponding to about 5.7times the isotropic radius of our fiducial neutron star. Weallocate the radial grid points using a sinh function, re-sulting in a grid that becomes logarithmic asymptoticallyand allows us to resolve the vastly different length scalesassociated with the black hole and the neutron star.We similarly implement the equations of relativistichydrodynamics in spherical polar coordinates adopting areference-metric formulation [45]. We solve the resultingequations using an HLLE approximate Riemann solver[46, 47], together with a simple monotonized central-difference limiter reconstruction scheme [48].As an example of our evolution calculations, we showin Fig. 3 profiles of the density ρ for our fiducial neu-tron star model hosting a black hole with puncture mass˜ M = 1 × − . We show results for different values ofthe conformal exponent m in (48). Evidently, the ini-tial density profiles, shown in the inset of Fig. 3, showlarge differences, as one might expect. However, once theevolution reaches quasi-equilibrium, shown in the maingraph in Fig. 3, the profiles are all very similar. Thisgives us confidence that, except for a small initial tran-sient, our evolution calculations are largely independentof our choice of m , and quickly settle down into a solutiondescribing a steady-state accretion onto the endoparasiticblack hole. C. Diagnostics
We invoke a number of different diagnostics in orderto evaluate our numerical simulations.
1. Black hole Mass
A black hole’s isolated horizon, or irreducible, mass isgiven by M BH = (cid:18) A π (cid:19) / , (58)where A is the proper area of the black hole’s event hori-zon at a given instant of time.In practice, we locate in our numerical evolution cal-culations the black hole’s apparent horizon rather thanthe event horizon, since the former requires data at oneinstant of time only, rather than the entire spacetime(see, e.g., [30] for a textbook discussion). The two hori-zons should coincide for quasistatic evolution. We thencompute the apparent horizon’s proper area, and use thisvalue in (58). For many situations, in particular for sta-tionary or nearly stationary spacetimes, this yields anexcellent approximation.We can also compute the approximate initial black-hole mass as follows. Since our initial data are confor-mally flat and describe a moment of time-symmetry withzero shift, we may write the spacetime metric at that in-stant as ds = − α dt + ψ ( dr + r d Ω ) . (59)The expansion of a bundle of outgoing null geodesics or-thogonal to a spherical surface of radius r is then givenby Θ = √ rψ ddr (cid:0) rψ (cid:1) (60)(see, e.g., Eq. (7.22) in [30] with A = B = ψ ). Wewill assume that u in (52) is small compared to the otherterms (see Appendix C) so that we may approximate ψ (cid:39) ψ NS + ψ BH = ψ NS + M r . (61)For r much smaller than the neutron star radius we mayapproximate ψ NS as constant. We now find the black hole’s apparent horizon by set-ting the expansion (60) to zero, which yields ddr (cid:32) r (cid:18) ψ NS + M r (cid:19) (cid:33) (cid:39) ψ − M r = 0 (62)or r AH (cid:39) M ψ NS . (63)At the apparent horizon, we then have ψ AH ≡ ψ ( r AH ) (cid:39) ψ NS . (64)We now compute the apparent horizon’s proper area from A = 4 πψ r = 16 πψ M , (65)where have used both (63) and (64). Inserting (65) into(58) we obtain our result M BH (cid:39) ψ NS M , (66)where ψ NS may be estimated by the central value of theunperturbed neutron star’s conformal factor. We havefound that, for M (cid:28) M , Eq. (66) provides an excel-lent approximation to our numerical values for the initialblack hole masses (see Table III).
2. Accretion rates: black hole growth
We measure the rates at which the black hole accretesneutron star matter in two different ways (see also Fig. 1and Table III).In one approach, we directly measure the black holemass M BH from (58) as a function of coordinate time,and then determine the accretion rate from the slope ofthis function. We show an example in Fig. 4. We alsoincluded in Fig. 4 results for a black hole of the sameinitial mass M BH (0), but evolved in vacuum, i.e. withouta neutron star, shown as the dashed (green) line. Thelatter appears nearly horizontal, demonstrating that thegrowth observed for a black hole embedded in a neutronstar indeed results from accretion of neutron star ma-terial, rather than numerical noise. For small black holemasses, however, the accretion becomes so slow (cf. Eq. 9)that it is no longer possible to accurately determine theslope of the function M BH ( t ). Therefore we have usedthis direct measure of the black hole growth only for blackholes with ˜ M BH (0) (cid:38) − (see Fig. 1 and Table III).Measuring the slope of curves M BH ( t ) yields the ac-cretion rate ˙ M BH , where the derivative is taken with re-spect to the coordinate time t . Since this coordinatetime agrees with proper time of a static observer at in-finity, i.e. one at large distances from the neutron starwith r (cid:29) R , this measure determines the accretion rateas observed by a static observer at infinity. In SectionII A 1, we introduced the accretion rate ˙ M (cid:63) BH as observed0 t / M BH (0)0.9951.0001.0051.0101.0151.0201.0251.030 M B H / M B H ( ) M BH in NS best fit M BH in vacuum0 50 100 . . FIG. 4. Growth of the black hole’s irreducible mass (58) as afunction of time, for ˜ M BH (0) = 1 . × − (the solid orangeline). The dashed-dotted (green) line shows a linear fit, whoseslope we identify with the accretion rate ˙ M BH . For compar-ison, we also include as the dashed (red) line the irreduciblemass of a black hole with the same initial mass M BH (0) invacuum, i.e. without a black hole. by a static “local asymptotic” observer far from the blackhole, but well inside the star, with M BH (cid:28) r (cid:28) R . In˙ M (cid:63) BH , the derivative is therefore taken with respect tothis local observer’s proper time τ (cid:63) (see also AppendixA for a detailed discussion). We can relate the two ratesby recognizing that the proper time of the static, localobserver moving along a normal vector to our spacelikehypersurfaces will advance at a rate dτ (cid:63) = α (cid:63) dt , where α (cid:63) is the lapse function of this local observer. We thenhave ˙ M BH = dM BH dt = α (cid:63) dM BH dτ (cid:63) = α (cid:63) ˙ M (cid:63) BH (67)(see also Eq. A21). We may interpret this relation asstating that the rate as observed by a distant observer,˙ M BH , is red-shifted by the lapse function α (cid:63) with respectto the rate as observed by a local observer, ˙ M (cid:63) BH , as onemight expect.In Fig. 5 we show profiles of the rest-mass density ρ ,the lapse function α and the shift vector β r in one of oursimulations. This example shows how the density andlapse function “plateau” in a region M BH (cid:28) r (cid:28) R , sothat their values for a “local asymptotic” observer can beidentified to high accuracy as long as M BH (cid:28) M . Notealso that, in this region, the shift β r is very close to zero,so that in this nearly static spacetime a static observer(i.e. one whose four-velocity is aligned with the timelikeKilling vector) indeed coincides with a normal observer,as we had assumed above.We list results for accretion rates determined numeri-cally from the growth of the black hole mass in Table III. . . M ρ α r/M BH . . β r t = 0 t = 642 M BH t = 1286 M BH FIG. 5. Profiles of the rest-mass density ρ , lapse α and shiftvector β r in an evolution of our fiducial neutron star witha black hole of mass ˜ M BH (0) = 1 . × − at different(coordinate) times. At the later times, all functions havesettled down into a (quasi) equilibrium. Note also that thedensity and lapse “plateau” in a region M BH (cid:28) r (cid:28) R , wherethey take the nearly constant values ρ (cid:63) (cid:39) ρ c = 0 . α (cid:63) (cid:39) .
623 (marked by the horizontal dotted lines). Notealso that the shift is very close to zero in this region. Theblack circles denote the location of the black hole horizon.
3. Accretion rates: rest-mass flux
An alternative approach to computing the accretionrate is to measure the rate of fluid flow across the blackhole horizon. Assuming that the accretion is sufficientlyslow so that we can approximate the black hole as nearlystatic, we can compute the flux F of rest-mass accretionthrough a sphere H of radius r from F ( r ) = − (cid:90) H √− gρ u r dθdφ, (68)where g is the determinant of the spacetime metric (see,e.g., Appendix A in [49]). Assuming spherical symmetrywe may carry out the integration to obtain F ( r ) = − πα √ γρ u r r , (69)where we have used √− g = α √ γ , and where γ is thedeterminant of the spatial metric. This expression yieldsthe flux of rest-mass through any sphere of radius r , andthe accretion rate in particular when evaluated on theblack hole horizon, ˙ M BH = F ( r hor ) . (70)For stationary flow we expect F to become independentof radius. We demonstrate this behavior in Fig. 6, wherewe show profiles of F at different instants of time, for ablack hole with initial mass M BH (0) /M = 7 . × − embedded in our fiducial neutron star. Note that F takesa nearly constant value in an inner region that grows1 r / M BH / M B H t / M BH = 1.63e+02 t / M BH = 3.27e+02 t / M BH = 4.90e+02 FIG. 6. Profiles of the flux (69) at different instants of time,for a black hole with initial mass ˜ M BH (0) = 1 . × − ,˜ M = 1 × − embedded in our fiducial neutron star model.In the outer parts of the star, the non-zero flux reflects a nu-merical adjustment of the star, resulting from the fact thatthe initial data are not in perfect equilibrium on the numer-ical grid. In the inner part, the flux approaches a value thatbecomes independent of both space and time, resulting in anequilibrium accretion flow onto the black hole. For compari-son we also included, as the faint lines, profiles for an evolu-tion of the same neutron star but without a black hole, whichshows the same behavior in the outer parts of the star, butvery different behavior in the vicinity of the black hole. with time, as the fluid flow settles down into steady-stateaccretion onto the black hole.Also included in Fig. 6 are profiles of the flux F foran evolution of our fiducial neutron star model withouta black hole. While these profiles show the same behav-ior as the neutron star with the black hole in the outerparts of the star, where the flux is dominated by a nu-merical adjustment of the near-equilibrium initial data tothe numerical grid, the profiles are very different in thevicinity of the black hole. This demonstrates that theplateau in the flux observed for the evolution with theblack hole hole indeed represents steady-state accretiononto the black hole, rather than a numerical artifact.We record our numerical results for these accretionrates in Table III. Note that, just like the accretion ratecomputed in Section III C 2, the rate (70) represents arate as measured by an observer at a large distancefrom the neutron star, i.e. at r (cid:29) R . Comparing thisrate with the rates computed in Section II, which repre-sented those measured by a “local asymptotic” observerat M BH (cid:28) r (cid:28) R , again requires this local observer’slapse function α (cid:63) (see also Appendix A).Also note that the accretion rate discussed in III C 2measures changes in the black hole’s gravitational mass,while the flux (70) measures the accretion of rest mass.In our numerical simulations we find that the black hole’sgravitational mass grows at a rate somewhat larger thanthe rate of rest-mass accretion, which is presumably be- ˜ ρ W ξ Bondinumerical r/M BH . . F ? / ˜ M B H FIG. 7. Comparison of our numerical fluid flow profiles(dashed red lines) with those obtained from integrating theBondi equations (solid blue lines), for our fiducial neutronstar model hosting a black hole of initial mass ˜ M BH (0) =1 . × − , at time t = 1 . × M BH (see text for de-tails). The open black circles denote the location of the blackhole horizon, while the solid green dots mark the location ofthe critical radius (which, in the Newtonian limit, reduces tothe location at which the fluid flow becomes supersonic). cause the former includes the accretion of other forms ofenergy (in particular, internal thermal energy) in addi-tion to rest-mass energy. IV. RESULTS
In this Section we compare our numerical results inSection III for our fiducial neutron star model hostingblack holes with a wide range of different masses to ouranalytical estimates in Section II.
A. Comparison with Bondi flow
We start with a comparison of fluid flow profiles. Inour numerical simulations, we focus on data in the vicin-ity of the black hole, at sufficiently late times so thatthe fluid has had enough time to settle down into steady-state accretion. We compare these numerical results withthose resulting from a direct integration of the “relativis-tic Bondi-equations”, i.e. the equations describing spher-ically symmetric, steady-state, adiabatic fluid flow in aSchwarzschild spacetime (see, e.g., Appendix G in STand [23]).Since the coordinates used in our code are differentfrom the Schwarzschild coordinates used in the usual con-struction of the steady-state Bondi solution, only scalarquantities – for example the rest-mass density ρ – canbe easily compared directly. In order to compare an in-variant measure of the fluid four-velocity u a , we com-pute the “gamma-factor” between an observer comov-2ing with the fluid and a “Killing observer”, i.e. a staticobserver whose four-velocity is aligned with a timelikeKilling-vector ξ a = ∂/∂t , W ξ = − ξ a u a ( − ξ a ξ a ) / = − u t ( − g tt ) / . (71)In Schwarzschild coordinates, this can be expressed as W ξ = α S u t S = (cid:18) − M BH r S (cid:19) / u t S = (cid:40) (cid:18) − M BH r S (cid:19) − ( u r S ) (cid:41) / , (72)while, in our code, we evaluate W ξ = W ( α − β r v r )( α − β r β r ) / (73)with W = αu t and v r ≡ W γ ra u a = u r W + β r α , (74)where v i is the spatial projection of the four-velocity u a ,divided by W . Here we assume not only that the space-time is indeed approximately static, but also that thelapse α and the shift β r render this spacetime in a co-ordinate system that leaves the metric quantities nearlytime-independent. Note also that static observers existonly outside the black hole, so that we can evaluate W ξ only for r (cid:38) M BH .Finally, we can compare the flux (69) computed in ourcode with the accretion rate (9) as predicted from theBondi solution.We show such a comparison for ˜ M = 10 − at (coordi-nate) time t = 1 . × M BH in Fig. 7. For the rest-massdensity ρ and the gamma-factor W ξ , the curves agree sowell that they can hardly be distinguished in the figure.The graphs for the accretion appear to differ more, but,at least in part, that is because the scale of the y -axisspans a much smaller range for the (almost) constantfunctions displayed in this panel. In fact, even the accre-tion rates agree to within a small fraction of a percent inthe vicinity of the black hole, demonstrating that the ac-cretion onto an endoparasitic black hole inside a neutronstar is well-described by relativistic Bondi accretion. B. Final Fate: Total Consumption
For small initial black hole masses, the total accretiontimes are far too long for us to simulate the consumptionof the entire neutron star numerically. We have thereforeperformed simulations of such a complete consumptiononly for sufficiently large initial black hole masses.As an example, we show in Fig. 8 results from a simula-tion with an initial black hole mass of ˜ M BH (0) = 0 . t/M ADM . . . . . M B H / M A D M numericalana. est. with evol.ana. est. no evol. FIG. 8. The black hole mass M BH as a function of coordinatetime t in a simulation leading to complete collapse. The solidline represents the numerical values of the black hole mass M BH , starting with the initial black hole mass ˜ M BH (0) =0 . M BH (0) /M ADM = 0 . M ADM . The dashed line shows the analytical estimate(B16) that takes into account effects of stellar evolution, whilethe dotted line shows the estimate (35) that ignores theseeffects. In this comparison we have adopted the value α (cid:63) (cid:39) . These initial data have a total gravitational Arnowitt-Deser-Misner (ADM) mass of ˜ M ADM = 0 . M BH (0) /M ADM = 0 . t (cid:39) M ADM , however, the consumption becomes dy-namical, as described in Section II A 2. Once the blackhole has consumed the entire neutron star, its mass set-tles down to a value that agrees to within less than 0.1%with the initial ADM mass of the system, which confirmsthe accuracy of our numerical simulations.We also included in Fig. 8 the analytical estimates re-sulting from the integrations of Section II C. Specifically,the dashed line represents an estimate that includes theeffects of stellar evolution (see Section II C 2), while thedotted line ignores stellar evolution (Section II C 1).We caution that these comparisons should be consid-ered qualitative for several reasons. Perhaps most im-portantly, the initial ratio of the black hole to neutronstar mass is too large for the asymptotic region in thecore just beyond where the gas becomes bound to theblack hole to be homogeneous, as required by the Bondimodel. In addition, the final phase of accretion departsfrom secular to dynamical growth, and the Bondi rates weadopted throughout break down. Also, all the estimatesof Section II are based on rates as observed by a “localasymptotic”, static observer. Relating these “local” rates3to the “global” rates computed in the numerical simula-tions requires the lapse function of a local observer, α (cid:63) .While this can be done rather well if M BH (cid:28) M , asshown in Fig. 5, such a lapse function can no longer beidentified unambiguously if M BH (cid:46) M . In fact, the en-tire notion of a local, static observer in a region with M BH (cid:28) r (cid:28) R no longer applies during the late stages.Finally, we used simple Newtonian arguments to modeleffects of stellar evolution (see Section II B), which clearlydo not apply in the late-time, dynamical collapse. Werewe able to follow the entire evolution for an initial blackhole with M BH (cid:28) M the number of decades in bothtime and increasing mass ratio during which the analyticcurves closely match the numerical tracks in Fig 8 wouldbe considerably larger.Despite all these disclaimers, we observe in Fig. 8 thatthe qualitative agreement between the numerical resultsand analytical estimates is reasonable. In particular, wesee that including the effects of stellar evolution does im-prove this agreement. In these comparisons, we adoptedthe value α (cid:63) = 0 .
623 to convert from our analytic localto our numerical global observers. This was the valuethat we had identified for smaller black hole masses forthis fiducial neutron star. But as we discussed above, inreality this value is no longer well defined at late times.
C. Accretion rates
We summarize in Table III results from our simulationsfor a large range of initial black hole masses, spanningseven orders of magnitude in M BH /M . For each blackhole mass, we compute in our code the accretion rates˙ M BH , directly from the growth of the black hole mass (forsufficiently large initial black holes) as discussed in Sec-tion III C 2, and/or from the fluid flux, ˙ M BH = F ( r AH ),as discussed in Section III C 3. We also identify the lapse α (cid:63) of a “local asymptotic” observer far inside the neu-tron star (see Fig. 5 for an example), and then dividethe above “global” rates by this lapse in order to com-pute the rates that such a local observer would measure.The latter can also be estimated from the Bondi expres-sion (44), the result of which we list in the last columnof Table III. The entries in this Table are also shown inFig. 1.We clarify that, for the small black hole masses listedin Table III, we do not track the evolution to completion,until the entire star has been consumed, since doing sowould require following the star for many more dynam-ical timescales than is computationally feasible. Insteadwe evolve the system for a coordinate time of approxi-mately 10 M BH , by which time the accretion has settleddown into equilibrium and we can accurately measure theaccretion rate.Most importantly we observe that the different mea-sures of the accretion rates agree well with each other.In our numerical simulations of Section III, the flux of(baryon) rest mass across the horizon agrees very well with the analytical accretion rates computed from rel-ativistic Bondi accretion in Section II. Determining theaccretion rate from the growth of the black hole horizon,which provides a measure of the increase in total grav-itational mass, results in a somewhat larger value, pre-sumably because this includes internal thermal energy inaddition to rest-mass energy.Finally, we note that the initial black hole mass˜ M BH (0) agrees very well with the analytical estimate (66)with ψ c = 1 .
27 for our fiducial neutron star model (seeTable II).
V. DISCUSSION
We study in detail the process by which a small “en-doparasitic” black hole, residing at the center of a neu-tron star, consumes its host. While a number of as-pects of this problem have been studied before (see, e.g.,[5, 7, 8]), we expand on these treatments in a number ofways.Building on our previous study of Bondi accretion forstiff EOSs [23, 24] we develop a quantitative analyticaldescription of this accretion process. In particular, thisallows us to determine the constant of proportionality inthe relation ˙ M BH ∝ M that some previous authors hadadopted. We use these results to construct an approxi-mate analytic model that tracks the secular evolutionof the system as the black hole, assumed initially small,grows by accretion and ultimately consumes the entireneutron star.We also perform numerical simulations of this accre-tion process, extending previous simulations (see [7]) tosignificantly smaller ratios M BH /M and simulating longenough for the systems to achieve quasistationary accre-tion. Our numerical code adopts spherical polar coordi-nates (see [25, 26]) with a logarithmic radial coordinate,which allows us to adequately resolve the vastly differentlength scales of the black hole and neutron star.As shown in Fig. 1, our numerical results for the ac-cretion rates agree remarkably well with those computedfrom relativistic Bondi accretion over many orders ofmagnitude in M BH /M . This establishes that the ac-cretion onto small black holes at the center of neutronstars is indeed governed by secular Bondi accretion, andthat the lifetimes of such stars are determined by theseaccretion rates. In particular, this supports our findingreported in [24] that this lifetime is close to a nearly uni-versal maximum lifetime that is roughly independent ofthe properties of the neutron star and its EOS, and de-pends on the initial black hole mass M BH only.As an important application, our results corroboratearguments that use the current existence of neutron starpopulations to constrain either the contribution of pri-mordial black holes to the dark matter content of theUniverse, or that of dark matter particles that may formblack holes at the center of neutron stars after they havebeen captured (see, e.g., [1–8]). These constraints are4 ˜ M a ˜ M BH (0) b M BH (0) /M α (cid:63) d ˙ M BH /α (cid:63) e F ( r AH ) /α (cid:63) f ˙ M (cid:63) BHg − . × − . × − .
616 4 . × − . × − . × − − . × − . × − .
619 4 . × − . × − . × − − . × − . × − .
622 – 3 . × − . × − − . × − . × − .
623 – 3 . × − . × − − . × − . × − .
623 – 3 . × − . × − − . × − . × − .
623 – 3 . × − . × − − . × − . × − .
623 – 3 . × − . × − − . × − . × − .
623 – 3 . × − . × − Black hole puncture mass ˜ M = K − / M in our initial data; see Section III A. b Initial irreducible mass ˜ M BH (0) = K − / M BH (0) of the black hole; see Section III C 1. c Ratio between M BH (0) and the neutron star rest mass M . d Lapse of a “local asymptotic” static observer; see, e.g., Fig. 5 for an example. e Mass-energy accretion rate from measurements of M BH ; see Section III C 2. f Rest-mass accretion rate from flux across horizon; see Section III C 3. g Rest-mass accretion rate from Bondi expressions; see Eq. (44).
TABLE III. Accretion rates for different black hole masses embedded in our fiducial neutron star model (see Table II). Theaccretion rates ˙ M BH and F ( r AH ) represent rates as measured by a static observer at infinity. To compare these rates with thosemeasured by a “local asymptotic” static observer in the neutron star core, ˙ M (cid:63) BH , we divide the former by the lapse α (cid:63) of thelocal observer. The numerically measured rest-mass flux F ( r AH ) /α (cid:63) (see Section III C 3) agrees very well with the analyticalvalue ˙ M (cid:63) BH given by Bondi accretion (see Eq. 44). Measuring the growth of the black hole horizon ˙ M BH /α (cid:63) (see Section III C 2)results in somewhat larger values, presumably because it includes internal thermal energy in addition to rest-mass energy.Results in this table are also shown in Fig. 1. based on the notion that, given certain cosmological den-sities of these dark matter constituents and their masses,neutron stars would capture these objects and wouldthen be consumed by the black holes after times thatare in conflict with the ages of old neutron star popu-lations. In particular, these arguments have been usedto constrain the contribution of PBHs in the mass range10 − M (cid:12) (cid:46) M BH (cid:46) − M (cid:12) (see, e.g., [5, 13]). ACKNOWLEDGMENTS
We would like to thank Gordon Baym, Milton Ruiz,Maria Perez Mendoza, and Antonios Tsokaros for help-ful conversations. CBR acknowledges support throughan undergraduate research fellowship at Bowdoin Col-lege. This work was supported in parts by NationalScience Foundation (NSF) grants PHY-2010394 to Bow-doin College, and NSF grants PHY-1662211 and PHY-2006066 and National Aeronautics and Space Adminis-tration (NASA) grant 80NSSC17K0070 to the Universityof Illinois at Urbana-Champaign.
Appendix A: Black hole metric and Bondi accretioninside a star
In this Appendix we generalize the Schwarzschild met-ric so that it allows for an asymptotic region that is differ-ent from the Minkowski metric. We then adopt this formof the metric to describe approximately a small black holeresiding at the center of the neutron star, and to rederivethe equations governing Bondi flow in this context.
1. A general form of the Schwarzschild metric
Following Section 5.1 in [50], we write the sphericallysymmetric, time-independent metric in vacuum as ds = − e A dt + e B dR + R d Ω , (A1)where, in this Appendix only, R is the areal radius (see5.11 in [50], hereafter C5.11, except that we use A and B instead of α and β in order to avoid confusion with thelapse function and the shift vector). We first evaluate thecombination e B − A ) R tt + R RR , where R ab is the Riccitensor, in Einstein’s vacuum field equations. This yields ∂ R A + ∂ R B = 0 (see C5.16) and hence A = − B + C, (A2)where C is a constant of integration. We usually set thisconstant to zero (see C5.17), which results in t being theproper time of a static observer at R → ∞ , but here wewill allow C to remain non-zero, and will determine itsvalue in Section A 2 below.We next evaluate R θθ in Einstein’s equations, whichnow takes the form e A (2 R ∂ R A + 1) = e C (A3)(compare C5.18) and is solved by e A = e C − κR , (A4)where κ is another constant of integration. We usually identify this constant with 2 M , where M is the gravita-tional mass of the black hole, but we will again postponedetermining this constant until Section A 2.5Using (A4) and (A2) in (A1) we now find ds = − (cid:16) e C − κR (cid:17) dt + (cid:18) − e − C κR (cid:19) − dR + R d Ω (A5)for the general form of the Schwarzschild metric. Evi-dently, we recover the asymptotic-Minkowski form of themetric for C = 0 and κ = 2 M .
2. Identification of constants
We next identify the constants C and κ in the metric(A5) by matching to our initial data describing a blackhole embedded inside a neutron star. We assume thatthe black hole’s mass is small compared to the mass ofthe star, M BH (cid:28) M , and will restrict our analysis toa region that is large compared to M BH , but small com-pared to M (and the radius of the neutron star). We maythen approximate the neutron star’s conformal factor asconstant, ψ NS (cid:39) ψ c (see Table II and the discussion inSection III A), and may neglect u in (52) (see AppendixC). We may also neglect the neutron star’s contributionto the total mass-energy in this region, so that the metric(A5) still provides an approximate solution to Einstein’sequations, even though it was derived in vacuum.Recall that, in Section III A, we construct the initialspatial line element from dl = ψ ( dr + r d Ω ) , (A6)where r is the isotropic radius, and, by our assumptionabove, the conformal factor (52) reduces to ψ = ψ c + M r . (A7)Identifying (A6) with the spatial part of (A5), we obtain rψ = R (A8)and (cid:18) − e − C κR (cid:19) − / dR = ψ dr = R drr (A9)or dRR / ( R − e − C κ ) / = drr . (A10)Integration of both sides yields Dr = 2 (cid:18) R + (cid:112) R − Re − C κ − e − C κ (cid:19) (A11)where D is a constant of integration, which we can solvefor R to find R = Dr (cid:18) e − C κDr (cid:19) . (A12) We next use (A8) to write ψ = (cid:18) Rr (cid:19) / = D / (cid:18) e − C κDr (cid:19) (A13)and compare this with (A7) to identify D = 4 ψ c (A14)and e − C κ = 2 ψ c M = 2 M irr , (A15)where we have used (66), derived under the same assump-tions as our treatment here, in the last equality.We now evaluate the time part of the metric (A5). Inparticular, consider a “static asymptotic” observer as in-troduced in Section II A 1, i.e. an observer at areal radius R (cid:29) M irr , but far inside the neutron star so that our as-sumptions above apply. According to (A5), the propertime of such an observer advances at a rate dτ (cid:63) = e C dt ,meaning that we may identify e C with the lapse functionof such an observer, α (cid:63) .Finally, we can insert our results into (A5) to find themetric describing the spacetime in the vicinity of a blackhole embedded in a neutron star, ds = − α (cid:63) (cid:18) − M BH R (cid:19) dt + (cid:18) − M BH R (cid:19) − dR + R d Ω , (A16)where we have identified the black hole mass with itsirreducible mass, M BH = M irr .
3. Bondi accretion
The equations governing Bondi accretion, i.e. the con-tinuity equation and the Euler equation of relativistichydrodynamics for stationary, adiabatic, and sphericallysymmetric fluid flow, are usually derived assuming aSchwarzschild metric that asymptotes to a Minkowskimetric. While our metric (A16) takes a slightly differ-ent form, both the continuity equation and the Eulerequation take the exact same form. That means that theresults for relativistic Bondi accretion (as presented, forexample, in Appendix G of ST or in [23]) can be adoptedwithout change, provided that the black hole mass M BH is identified with its irreducible mass M irr . In particular,the accretion rate is given by Eq. (9).There remains one ambiguity, however, namely themeaning of the time derivative in ˙ M BH in the accre-tion rate, i.e. whose time we refer to in this derivative.In order to clarify this, we compute the radial compo-nent of the fluid flux as observed by a static observerat an arbitrary radius R outside the black hole horizon(i.e. not necessarily at R (cid:29) M BH , but well inside theneutron star). Specifically, we take the dot product be-tween the fluid flux J = ρ u (where we use bold-face6to denote a vector) with the orthonormal basis one-form˜ ω ˆ R = (1 − M BH /R ) − / ˜ ω R to find J ˆ R = ˜ ω ˆ R · J = (cid:18) − M BH R (cid:19) − / J R . (A17)We now multiply with − πR and define u ≡ − u R (which is positive for inflowing matter) to obtain the ac-cretion rate dM/dτ as measured by this observer, dM BH dτ = − πR J ˆ R = 4 πR (cid:18) − M BH R (cid:19) − / ρ u. (A18)Since, for this observer, dτ = α (cid:63) (1 − M BH /R ) / dt , wemay rewrite the above expression as dM BH dt = 4 πα (cid:63) R ρ u, (A19)which implies dM BH dτ (cid:63) = 4 πR ρ u, (A20)where τ (cid:63) is again the proper time of a static “local asymp-totic” observer. We see that the “usual” expression forthe Bondi accretion rate (e.g. 14.3.18 or G.21 in ST)refers to a time derivative with respect to the time asobserved by a static asymptotic observer, which, in ourcase, means a “local asymptotic” observer – far from theblack hole, but well inside the neutron star. In (9) andelsewhere we emphasize this by denoting the time deriva-tive with a star, ˙ M (cid:63) BH = dM BH /dτ (cid:63) . We also see that theaccretion rates observed by a local asymptotic observerand an observer at infinity are related by˙ M (cid:63) BH = dM BH dτ (cid:63) = 1 α (cid:63) dM BH dt = 1 α (cid:63) ˙ M BH , (A21)as we have observed already in (67). Appendix B: Integration of Eq. (32)
In this appendix we outline how the differential equa-tion (32) can be solved analytically.We first separate variables to obtain dT = − y / dy ( y − y ) (B1)and then integrate to find T = − I, (B2)where we have assumed that the initial time is chosen tovanish, T i = 0, and where I is given by I = (cid:90) y / dy ( y − y ) . (B3) This integral can now be integrated as follows.Using partial fractions, we rewrite (B3) as I = (cid:90) y / dy ( y − y ) = (cid:90) y / y + ( y − y ) − ( y − y ) ( y − y ) dy = (cid:90) y / ( y − y ) + 2 yy − y ( y − y ) dy = (cid:90) y / dy + 2 y (cid:90) y / y − y / y − y ) dy. (B4)Repeating the process twice more, we obtain I = (cid:90) y / dy + 2 y (cid:90) y / dy + 3 y (cid:90) y − / dy +4 y (cid:90) y − / y − y / y − y ) dy. (B5)We now split the last integral into two terms, y (cid:90) y − / y − y ( y − y ) dy (B6)= 4 y (cid:90) y / ( y − y ) dy − y (cid:90) y − / ( y − y ) dy and use a hyperbolic trig substitution y = y tanh x (B7)in both integrals, resulting in4 y (cid:90) y / ( y − y ) dy = 8 y / (cid:90) sinh x dx (B8)and 3 y (cid:90) y − / ( y − y ) dy = 6 y / (cid:90) cosh x dx. (B9)Since (cid:90) sinh x dx = 12 sinh x cosh x − x (cid:90) cosh x dx = 12 sinh x cosh x + x , (B11)as can be seen using integration by parts, we can combineresults to find y (cid:90) y − / y − y ( y − y ) dy = y / (sinh x cosh x − x ) . (B12)We now rewrite sinh x cosh x in terms of tanh x and insertthe substitution (B7) to obtain y (cid:90) y − / y − y ( y − y ) dy (B13)= y y / y − y − y / tanh − ( y/u ) / . I = 25 y / + 43 y y / + 6 y y / (B14)+ y y / y − y − y / tanh − ( y/y ) / . Combining the first four terms and usingtanh − ( x ) = 12 ln (cid:18) x − x (cid:19) (B15)we can also write this result as I = y / y + 14 y y + 70 y y − y y − y ) − y / ln (cid:32) y / + y / y / − y / (cid:33) . (B16)Recall that 1 − y = ( M BH − M BH (0)) /M (0) measures thefractional increase in the black hole mass (see Eqs. 29and 30), and that T = − I is proportional to the timeas measured by a local asymptotic static observer (seeEq. 31). Appendix C: An Approximate Analytical Solutionto the Hamiltonian constraint
In this Appendix we present an approximate but ana-lytical solution to the Hamiltonian constraint (55),¯ D u = − π (cid:110) ( ψ NS + ψ BH + u ) n − ψ n NS (cid:111) ¯ ρ. (C1)We have solved this equation to high precision in SectionIII A, and our goal here is not to reproduce that solutionquantitatively; instead, we adopt some simple argumentsthat allow us to understand the qualitative behavior ofthe solution u . In particular, we will find that u ap-proaches zero everywhere as M BH /M →
0, which justifiesneglecting u , for example, in (61) and (A7).In fact, we start by assuming that u (cid:28) ψ NS , so thatwe may neglect this term on the right-hand side of (C1),ultimately verifying the validity of this step by showingthat it is well-satisfied by our final solution. We furtherapproximate the neutron-star conformal factor ψ NS byusing its central value throughout the entire interior ofthe star, i.e. for all r < R , and we replace the density¯ ρ by its average value ¯ ρ ave = ¯ ρ c /δ , where we adopt theNewtonian central condensation δ = 3 .
29 for a Γ = 2polytrope (see Table I). While we will see that u doesindeed become arbitrarily small as M approaches zero,the latter approximation is rather crude and introducesdiscrepancies, but only of order unity.We next observe that, even in the interior of the star,the right-hand side of (C1) behaves differently dependingon whether ψ BH or ψ NS is greater. If ψ BH = M / (2 r ) (cid:29) ψ NS , and assuming n < −
5, we have( ψ NS + ψ BH ) n − ψ n NS (cid:39) ψ n BH − ψ n NS (cid:39) − ψ n NS . (C2) On the other hand, if ψ BH (cid:28) ψ NS , we may expand( ψ NS + ψ BH ) n − ψ n NS (cid:39) (5 + n ) ψ n NS ψ BH . (C3)Recognizing that ψ BH = ψ NS at the approximate locationof the apparent horizon r AH = M ψ NS (C4)identified in (63), we may approximate Eq. (C1) as¯ D u = s r < r AH ˜ s/r r AH ≤ r < R r ≥ R, (C5)where we have defined s = 2 πψ n NS ¯ ρ (C6)and ˜ s = − π (5 + n ) ψ n NS M ¯ ρ = − (5 + n ) s r AH . (C7)Since ¯ D is the flat Laplace operator in spherical sym-metry, ¯ D u = 1 r ddr (cid:18) r dudr (cid:19) , (C8)piece-wise solutions to (C5) are given by u = sr + C I r < r AH
12 ˜ sr + C II + D II r r AH ≤ r < RC III r r ≥ R. (C9)Here we have assumed regularity at the origin r = 0 and u → r → ∞ , and the four constants C I , C II , D II ,and C III are constants of integration. We can determinethe latter by requiring that both u and its first derivativeare continuous at both r AH and R . Evaluating these fourconditions yields C I = − ˜ s ( R − r AH ) − s r C II = − ˜ s RD II = 12 ˜ s r − s r C III = 12 ˜ s ( r − R ) − s r . (C10)In Fig. 9 we compare the numerical solution for u , forour fiducial neutron-star model and ˜ M = 10 − , with theapproximate analytical solution (C9). As expected, thetwo solutions do not agree quantitatively, but they nev-ertheless show very similar qualitative behavior. Basedon this qualitative agreement we may make the followingobservations about the properties of solutions u to theHamiltonian constraint (55):8 − − − r/M − − − − − u × num.approx. FIG. 9. The numerical solution u for our fiducial neutron-star model and ˜ M = 10 − with n = − • The largest value of u , in magnitude, occurs at thecenter, where it is dominated by the term u c (cid:39) − ˜ sR ∝ − ¯ ρ M R. (C11)We see that u is proportional to M , and vanishes everywhere in the limit M →