{'fraction_non_alphanumeric': 0.040743476120137866, 'fraction_numerical': 0.019366486131626455, 'mean_word_length': 4.418630502445532, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accretion onto black holes inside neutron stars with piecewise-polytropic equations of state: analytic and numerical treatments

We consider spherically symmetric accretion onto a small, possibly primordial, black hole residing at the center of a neutron star governed by a cold nuclear equation of state (EOS). We generalize the relativistic Bondi solution for such EOSs, approximated by piecewise polytropes, and thereby obtain analytical expressions for the steady-state matter profiles and accretion rates. We compare these rates with those found by time-dependent, general relativistic hydrodynamical simulations upon relaxation and find excellent agreement. We consider several different candidate EOSs, neutron star masses and central densities and find that the accretion rates vary only little, resulting in an accretion rate that depends primarily on the black hole mass, and only weakly on the properties of the neutron star.


I. INTRODUCTION
Multiple authors have suggested that neutron stars may act as dark-matter detectors (e.g. [1][2][3][4][5][6][7][8]). In one scenario, primordial black holes (PBHs), which may either contribute to or even account for the dark matter in the Universe, may be captured by neutron stars and subsequently accrete the entire star [9,10]. The observed existence of neutron star populations has been invoked to constrain primordial black holes in a mass window of about 10 −15 M M BH 10 −9 M (see [5]) that is poorly constrained by other arguments and observations (see, e.g. [11][12][13][14][15]). Even this constraint assumes that both the capture and accretion processes are sufficiently fast. While different authors have arrived at different estimates for the rates of the former [5,8,16,17], we will derive analytical rates for the latter in this paper, even for neutron stars governed by realistic nuclear equations of state (EOSs).
In an alternative scenario, other candidate dark-matter particles, possibly including axions, also can be captured by neutron stars. Under sufficiently favorable conditions, these particles may coalesce to form a high-density object that then collapses to a small black hole (e.g. [1][2][3]6]), thereby resulting in the same accretion process as the scenario above. Independently of the precise scenario, it is of interest to explore the rate at which a central "endoparasitic" black hole disrupts its host neutron star.
Spherically symmetric, steady-state accretion onto a point mass of a fluid that is homogeneous and at rest far from the mass is described by the Bondi solution in Newtonian physics [18]. Most treatments of Bondi accretion, or its relativistic counterpart describing accretion onto a Schwarzschild black hole [19] (see also Appendix G in [20], hereafter ST), focus on soft EOSs with adiabatic indices 1 ≤ Γ ≤ 5/3, which is suitable for most astrophysical plasmas. While the EOS governing the cores of neutron stars is not known, most realistic candidates for the EOS at nuclear densities are stiff. As we discussed in [21] (hereafter RBS), accretion for stiff EOSs with Γ > 5/3 shows some qualitative differences from that of soft EOSs. In particular, there exists a minimum steadystate accretion rate for stiff EOSs. As shown in Appendix A of [22], these results also hold for a black hole inside a neutron star, as long as the black hole mass M BH is much smaller than that of the neutron star, M BH M NS . In this case the relativistic Bondi formalism yields the accretion rate as measured by a "local asymptotic observer", i.e. one who is far from the black hole, but deep inside the neutron star. The above minimum accretion rate therefore results in a maximum survival time for neutron stars harboring a black hole [23].
A number of authors have also performed numerical simulations of accretion onto endoparasitic black holes inside neutron stars. East and Lehner [7] considered black holes with masses M BH ≥ 10 −2 M NS and adopted nuclear EOSs (modeled with the "piecewise-polytrope" approximation described below). In particular, they observed that the accretion rate is proportional to M 2 BH , as expected from the relativistic Bondi formalism (see eq. 11 below), and that the effects of rotation are small (see also [24]). Focusing on nonrotating configurations and stiff Γ-law EOSs, we performed simulations for much smaller black holes with masses M BH 10 −9 M NS , and found that the accretion rates agree very well with those predicted by the Bondi formalism, even quantitatively [22].
In this paper we generalize the results of RBS and [22] to allow for realistic, nuclear EOSs, approximated by piecewise polytropes (PWPs). We review the PWP treatment of EOSs in Section II, and then derive analytical expressions for the stationary accretion rates onto small black holes inside neutron stars governed by such EOSs in Section III. We perform time-dependent, numerical simulations in full general relativity for these EOSs as described in Section IV. Once these numerical solutions have relaxed into a quasi-stationary solution they agree very well with the analytical solutions, as shown in Section V. Briefly summarizing in Section VI, we find that, for these realistic nuclear EOSs, the accretion rates depend only weakly on the EOS and the neutron star's density, and hence mainly on the black hole mass. These rates are just slightly larger than the minimum accretion rates that RBS and [23] computed under the assumption of Γ-law EOSs. Unless noted otherwise we use geometrized units with G = 1 = c in this paper, and adopt the solar mass 1M = 1.99 × 10 33 g = 1.47 × 10 5 cm = 4.9 × 10 −6 s as the fundamental unit.

II. NUCLEAR EOSS: APPROXIMATION BY PIECEWISE POLYTROPES
As demonstrated by [25], a large class of candidates for realistic, cold nuclear EOSs can be approximated remarkably well with piecewise polytropes (PWPs). Specifically, we write the pressure P as a function of the rest-mass density ρ 0 as where the constants K i and Γ i take different values in four different density "regions" labeled by the index 0 ≤ i ≤ 3. The different regions are separated by three "boundary densities" that can be chosen to be the same for all EOSs (see Table I). We follow [26] in our implementation of the PWP EOS; in particular, we model the low-density crust with only one piece rather than the four pieces adopted by [25] (see their Table II), and we also choose the values for the boundary densities as in [26]. The different values Γ i , as well as one of the constants K i , can be found from fits to each EOS; the remaining values of K i are then related to each other by imposing continuity of P across the boundaries between different regions. The specific internal energy then takes the form (see eq. 5 in [25]) where the b i are constants of integration that are chosen to make continuous across the boundaries between different regions. In the lowest-density region the constant b 0 vanishes, but in general b i = 0 for i > 0. 1 The total mass-energy density ρ is then 1 Note that [25] used symbols a i for these constants; we choose b i here in order to avoid confusion with the sound speed a below.
While P , and ρ are continuous across density boundaries, the speed of sound a, computed from where the subscript s denotes constant entropy, in general is not continuous -in fact, a may not even grow monotonically with ρ 0 . In each region i we may invert (4) to find but, since a may neither be a continuous nor a monotonic function of ρ 0 , we may not be able to invert (4) globally. A list of the Γ i for a large number of EOSs is given in Table III of [25]. In this paper we consider representatives of four different families of EOSs, namely the SLy [27], AP3 and AP4 [28], MS1 [29], and H4 [30] EOSs, which represent different theoretical approaches to constructing realistic nuclear EOSs. For these EOSs we provide all the above PWP parameters in Table II. All these EOSs result in maximum allowed masses that are consistent with all observed neutron star masses, including the largest currently known mass of 2.08 +0.07 −0.07 M (where the errors represent a 68.3% credibility interval), reported for the millisecond pulsar J0740+6620 from measurements of the relativistic Shapiro effects (see [31,32]; see also [33] for NICER and XMM analysis of the same pulsar, resulting in a consistent value for its mass, as well as [34] for constraints on the EOS resulting from NICER radius measurements of J0740+6620). Other high-mass neutron stars include PSR J1614-2230 with a mass of M = 1.93M (see [35]), and PSR J0348+0432 with a mass of M = 2.01M (see [36]). With the exception of MS1 and H4, the above EOSs are also consistent with the neutron star masses and tidal distortions inferred from the gravitational wave signal GW170817 and electromagnetic follow-up observations [37]; while both MS1 and H4 appear to be ruled out based on tidal distortions, we include them regardless as examples of stiff EOSs.

III. BONDI ACCRETION FOR PIECEWISE POLYTROPES: ANALYTICAL SOLUTION
We now generalize the Bondi solution [18][19][20], describing stationary, adiabatic, spherically symmetric fluid flow onto a Schwarzschild black hole, for PWPs. We follow the derivation in Appendix G of ST up to their eq. (G.22) (hereafter ST.G.22) unchanged, since it does not yet make any assumptions about the EOS. In particular (ST.G.21), the integrated continuity equation  , which we express in units of solar mass. We use Γ0 = 1.35692 in the lowest density region for all EOSs, adopting the value for the highest-density crust piece in [25] (see their Table II). We also provide the maximum gravitational mass Mmax of non-rotating neutron stars, and the corresponding central density ρ max 0c .
as well as (ST.G.22), the integrated Euler equation remain valid, as does (ST.G.17) for the relations at the critical radius, (provided this critical radius exists; see our discussion below). In the above equations r is the areal radius and u ≡ |u r | the inward radial component of the fluid's four-velocity; note that the above expressions employ Schwarzschild coordinates. Rather than adopting a single polytrope, as in ST and RBS, we now adopt the PWPs described above. In particular, the first factor on the left-hand side of (7) then takes the form ρ + P ρ 0 Note also that, using (5), we may relate the density at the critical point to that in the "local asymptotic region", denoted by a subscript * , by 10) (cf. eq. 10 in RBS, hereafter RBS.10). Inserting this into (6), and using (8), we may now write the accretion rate as observed by a local asymptotic observer (hence the superscript * ) aṡ which assumes the same form as, e.g., (ST.G.33) or (RBS.11), except that the dimensionless accretion eigenvalue λ is now given by . (12) Note that this reduces to eq. (RBS.12), as expected, when the critical point is in the same density region as the asymptotic observer, so that Γ s = Γ * , K s = K * and b s = b * . Note also that M BH on the right-hand side of (11) is the black hole's gravitational mass, while the lefthand side is the rate at which rest mass crosses the black hole's horizon.
In order to evaluate λ for given asymptotic values we need to relate a s to a * , which we will do using eq. (7). We start by inserting (5) into (9) to obtain Evaluating the left-hand side of (7) both at r s and in the local asymptotic region, where r * M and u * 1, then yields (cf. ST.G.29), or, taking the inverse of both sides and using (8), As for a single polytropic EOS, the relation (15) forms a cubic equation for x = a 2 s that we may write as  Table I); note the small discontinuities in the sound speed a at the corresponding locations.
(cf. RBS. 19). Unlike in RBS, however, we now need to evaluate the constants b i and Γ i in the regions corresponding to ρ 0 * and ρ 0s (or, equivalently, a * and a s ). For a given value of ρ 0 * , we know how to choose b * and Γ * , but, unless the local asymptotic values are in the highest-density region i = 3 already, we do not know a priori in which density region the critical point will be. Stated differently, solving (15) for a s requires values b s and Γ s , but choosing those depends on what region a s ends up in. We can solve this problem as follows. Say ρ 0 * is in density region j. Assuming that ρ 0s ≥ ρ 0 * we then consider all regions i ≥ j, and solve eq. (15), using Cardano's formula as described in RBS, to obtain candidate solutions a si for each one, disregarding unphysical solutions for which a 2 si < 0. We then evaluate (5) for remaining solutions a si in region i, and keep only those candidate solutions for which the corresponding rest-mass density ρ 0i is indeed in region i. In some cases, for low values of a * , we still find viable solutions in multiple regions from this procedure. For each one of these remaining solutions we can then construct fluid profiles by integrating eqs. (ST.G.10) both inwards and outwards away from the critical radius r si . For the examples that we considered, at most one solution a si resulted in global solutions with exactly one critical point. In the following we always adopt this solution as the analytical Bondi profile. We show an example of such a profile for the SLy EOS, extending over all four density regions, in Fig. 1. The above approach reduces to the simpler single Gamma-law case treated in [22], of course, if ρ 0 * is in the highest-density region already.
Also note that the above procedure is not guaranteed to yield solutions. As demonstrated in Appendix G in ST, the existence of a critical radius is guaranteed if all functions are continuous. Specifically, ST argue that, since their function (see ST.G.13) is negative for large r but positive for small r approaching the black hole horizon, it must have a root, which then provides the condition (8). In our treatment here, D must still change sign, but since, for piecewise polytropes, it can no longer be assumed to be continuous everywhere, this does not imply that it necessarily has a root.
For most examples that we considered we were able to find critical points without any problems. For small densities for the MS1 EOS, however, the above procedure did not yield any solutions, which we believe is related to the large discontinuity in the sound speed resulting from the large difference between Γ 2 and Γ 3 . In Fig. 2 below we show numerical results demonstrating that, in this case, a discontinuity in the sound speed leads to a "jump" across the critical point at which the last two terms in eq. (8) are equal, so that a strict root of D does not exist. Clearly, this behavior is an artifact of the piecewise-polytropic treatment of the EOS, which results in these discontinuities.

IV. NUMERICAL TREATMENT
We complement our analytical results by performing numerical simulations of the accretion onto black holes at the center of neutron stars as follows.
We construct initial data from a solution to the Tolman-Oppenheimer-Volkoff equations [38,39] for a given EOS and a given central density. Following [22] we then adopt a generalized puncture method to place a black hole with puncture mass M at the center of the neutron star, and solve the Hamiltonian constraint for the conformal factor ψ, assuming a moment of time symmetry, to obtain solutions to Einstein's constraint equations. As demonstrated in Section III.C.1 of [22], for M M NS the black hole's gravitational mass is well approximated by M BH ψ NS M, where ψ NS is the conformal factor at the center of the unperturbed neutron star.
Even though we start with initial data describing cold fluids, we allow for heating (e.g, by shocks) by adding to the cold pressure (1) thermal contributions. Specifically, we compute thermal contributions to the internal energy density U ≡ ρ 0 from where ρ 0 and are computed from the dynamically evolved quantities, and cold is given by (2). We then write U th = U nucl + U rad , where U nucl accounts for finitetemperature corrections to an ideal, nonrelativistic, nucleon Fermi gas, and U rad for contributions from radiation, (compare, e.g., [56]). In the above equations m B is the baryon rest mass (which we take to be equal to the neutron rest mass), n = ρ 0 /m B the baryon number density, k B the Boltzmann constant, Planck's constant, T the temperature, a rad the radiation constant, and the non-dimensional constant η depends on which particles contribute to the radiation. Allowing for photons (η ph = 1), three flavors of neutrinos (η ν = 3 × 7/8), as well as electron-positron pairs (η pairs = 7/4) we have η = η ph + η ν + η pairs = 43/8. We insert (20) and (21) into (19) and use a root-finding method to find the temperature T . Knowing T , we can finally compute the pressure P = P cold + P th using where Γ th = 5/3 for nonrelativistic nucleons. In all our simulations we observe a transition from our (astrophysically artificial) initial data to a nearly time-independent equilibrium solution describing accretion onto the black hole. Especially for initial data with smaller initial densities we see that this transition launches an outgoing shock wave; the accretion solution is then attained inside this shock wave. While this shock wave does lead to some heating, we find that this heating is small, especially for large initial densities, and presumably transient, so that we still find good agreement between our numerical and analytical solutions, even though the latter has been constructed for a cold gas.
In order to measure the accretion rate we monitor the flux F of rest mass through spheres S of radius r, where g is the determinant of the spacetime metric (see, e.g., Appendix A in [57]). The accretion rate is then given by the flux F evaluated on the horizoṅ The discontinuity in the sound speed a prevents an equality of the two terms in this case, so that the analytical procedure of Section III does not yield a critical point (compare the discussion at the end of Section III).
As in (11), this rate measures the accretion of rest mass rather than the change of the black hole's gravitational mass. Note also that (24) measures the accretion rate as seen by an observer at infinity, i.e. at large distances r R NS from the neutron star, while (11) measures that as seen by a local asymptotic observer with M BH r R NS . As discussed in Section III.C.2 of [22], we may compute the former from the latter usinġ where α * is the lapse function in the local asymptotic region.

V. RESULTS
We perform numerical simulations for the EOSs listed in Table II. For our analytical treatment, we adopt as the central density of the unperturbed neutron star a value just below that of the maximum mass configuration, as well as some smaller densities, for each EOS (see Table  III). We then compute the analytical accretion rates from (11). We compare these rates to our numerical simulation values from (24), adopting neutron star models with the central densities in the absence of the black hole identical to those chosen above. We provide details of all our results in Table III; in particular we list the analytical and numerical values for the accretion rates as measured by an observer at infinity.
As we discussed in Section III, the analytical approach described there does not yield analytical solutions for the   (11), the resulting analytical accretion rateṀ * 0 as measured by a local asymptotic observer, this observer's value of the lapse α * and the corresponding analytical accretion rate as measured by an observer at infinity α * Ṁ * 0 (see 25), as well as the numerical accretion rateṀ0 as computed from (24). For the MS1 EOS and for the smaller central densities the analytical approach of Section III did not yield solutions; see text for details.
MS1 EOS for our smaller central densities. Recall that this approach relies on identifying a critical point defined by equality between the three terms in eq. (8). The first of these three points, u s , is a gauge-dependent quantity, but the last two terms are gauge-invariant. In Fig. 2 we therefore show numerical profiles of these two terms for the MS1 EOS and ρ 0 = 7.5 × 10 −4 M −2 = 4.65 × 10 14 g cm −3 at a sufficiently late time for the evolution to have settled into a stationary equilibrium solution close to the black hole. We see that the difference between the two terms, M BH /(2r) − a 2 /(1 − 3a 2 ), does indeed change sign; however, this difference does not have a root because of the discontinuity of a, as we had discussed in Section III. This discontinuity, and hence the absence of a point at which equality in (8) holds, is an artifact of the PWP representation of the EOS. We note, however, that the numerical solution with its finite grid effectively interpolates between the discontinuity, thereby passing through a critical point and achieving a smooth flow and well-defined accretion rate, as indicated in Table III. In the last column of Table III we list the ratioṡ M 0 /M 2 BH and note that, for each EOS, these values depend only weakly on the central density or, equivalently, the mass of the neutron star host. Even between different EOSs these values do not vary significantly. We may therefore approximate the accretion rate, for any of the EOSs and central densities considered here, aṡ where, within about 30% or so, χ 0.06. Using M = 4.9 × 10 −6 s = 1.6 × 10 −13 yr we may write (26) aṡ which is just slightly larger than the minimum accretion rate reported in [23] (where it was computed under the assumption of single Gamma-law EOSs). 2 We also show examples of accretion profiles, computed both numerically and analytically, in Fig. 3. We note that for many of the examples that we considered these profiles feature superluminal sound speeds a > c in regions close to the black hole. This behavior is not unexpected; as discussed in RBS, it is unavoidable when Γ ≥ 3, and will also occur for softer EOS with 2 < Γ < 3 if the asymptotic densities are sufficiently high. Ref. [7] also noted the appearance of superluminal sound speeds for some EOSs, and artificially adjusted those EOSs in the corresponding high-density regimes to ensure that a < c. We instead allow the sound speed to exceed the speed of light, both in our analytical and numerical treatments. 3 The possibility of the sound speed becoming superluminal in ultradense matter has certainly been discussed in the past (see, e.g., [58]). While the equations admit such solutions without breaking down, there are strong arguments for rejecting such behavior on causality grounds, as it violates a basic principle of relativity, as emphasized by [59].
The appearance of superluminal sound speeds may be a consequence of either the underlying EOSs or their PWP representation, of course. To justify some EOSs and their PWP fits for treatments involving stable neutron stars, it is sometimes argued that the sound speed at the stellar core remains subluminal even for the maximum-mass configuration. However, the solutions that we present here provide examples of stationary equilibrium solutions in which the densities significantly exceed superluminal values for such configurations containing small black holes, even outside black-hole horizons. These solutions provide motivation for constructing nu-where M 0 is the initial black hole mass, and where we have adopted χ 0.06 in the last equality. The survival time (28) is slightly smaller, but close to the maximum survival time reported by [23], t max 1.0 × 10 6 s (10 −10 M /M 0 ). This timescale has been invoked to set constraints on the masses of primordial black holes [5,8] (but see also [16,17]).
Finally, we caution that many of these solutions feature superluminal sound speeds, and emphasize the requirement that any viable EOS and its PWP representation must avoid this unphysical behavior at the supranuclear densities encountered here.

ACKNOWLEDGMENTS
It is a pleasure to thank Charles Gammie, Chloe Richards, Lunan Sun, and Antonios Tsokaros for numerous helpful conversations. SCS acknowledges support through an undergraduate research fellowship at Bowdoin College. This work was supported in part by National Science Foundation (NSF) grants PHY-1707526 and PHY-2010394 to Bowdoin College, and NSF grants PHY-1662211 and PHY-2006066 and National Aeronautics and Space Administration (NASA) grant 80NSSC17K0070 to the University of Illinois at Urbana-Champaign.