Revisiting Primordial Black Hole Capture by Neutron Stars

A sub-solar mass primordial black hole (PBH) passing through a neutron star, can lose enough energy through interactions with the dense stellar medium to become gravitationally bound to the star. Once captured, the PBH would sink to the core of the neutron star, and completely consume it from the inside. In this paper, we improve previous energy-loss calculations by considering a realistic solution for the neutron star interior, and refine the treatment of the interaction dynamics and collapse likelihood. We then consider the effect of a sub-solar PBH population on neutron stars near the Galactic center. We find that it is not possible to explain the lack of observed pulsars near the galactic center through dynamical capture of PBHs, as the velocity dispersion is too high. We then show that future observations of old neutron stars close to Sgr A* could set stringent constraints on the PBHs abundance. These cannot however be extended in the currently unconstrained asteroid-mass range, since PBHs of smaller mass would lose less energy in their interaction with the neutron star and end up in orbits that are too loosely bound and likely to be disrupted by other stars in the Galactic center.


I. INTRODUCTION
Primordial black holes (PBHs) have long been recognized as viable dark matter candidates [1][2][3] (for recent reviews, see e.g.Refs.[4][5][6][7]).The interest in PBHs goes beyond their cosmological relevance, as their presence has been invoked to explain a variety of astrophysical observations (c.f.Ref. [8]), such as 1) microlensing events by planetary-mass objects [9]; 2) microlensing of quasars [10]; 3) microlensing events by objects between 2 and 5 solar masses [11]; 4) correlations in the X-ray and cosmic infrared background fluctuations [12]; 5) nonobservations of certain ultrafaint dwarf galaxies [13]; 6) masses, spins and coalescence rates for black holes found by LIGO/Virgo [14]; 7) relationship between the mass of a galaxy and that of its central supermassive black hole [15]; 8) so-called G objects near the Galactic centre (GC) [16,17].Many constraints exist on the PBH abundance [4][5][6][7].Viable scenarios however exist in which PBHs can in principle be compatible with all constraints and observations [5], for instance when the change in degrees of freedom as the Universe cools from the QCD scale to the electron mass leads to a multi-modal PBH mass function [18].
In this article, we assess constraints on the abundance of sub-stellar mass PBHs, by performing a state-of-the-art analysis of dynamical PBH capture by neutron stars (NS), extending and refining previous studies [19][20][21][22][23][24].Once captured, the PBH sinks into the core of the NS, accretes * rcaiozzo@sissa.itthe stellar matter, and leaves behind a black hole of mass around that of the initial NS.
We specifically focus on the innermost parsec of our Milky Way.This region, which we refer to as the Galactic center (GC), is of particular interest for the process of PBH capture by NSs, due to both the expected high density of dark matter and the predicted large population of pulsars [25,26].It is estimated that thousands of radio pulsars should be present in the GC [25].However, a number of surveys have struggled to find any within the innermost 25 parsec [27].This was thought to be mainly due to strong temporal scatterings, altering the discrimination of pulsar signals from the inner parsecs [28].However, the 2013 discovery of the magnetar SGR J1745-29 within 0.1 pc from Sagittarius A* [29], suggests that the temporal scattering cannot be as strong as previously thought [30].Several mechanisms have been put forward to explain this missing-pulsar problem in the GC, including efficient magnetar formation, and disruption of the expected pulsar population [27].We investigate here the possibility of NS disruption by PBHs (see also Ref. [24]).
Millisecond pulsars (MSPs) would be excellent candidates for PBH capture at the GC.Their presence could hence be used to formulate abundance constraints on PBHs.Because of their energy spectrum, MSPs are a likely candidate for the observed GeV excess from the Galactic bulge [31], which suggests a large population of these objects extending into the Galactic center.Ref. [32] predicts from numerical simulations 67 long lived NS-Xray binaries in the innermost parsec of the Milky Way.Given that these binaries are the progenitor of millisecond pulsars, it is reasonable to expect a population of MSPs extending in the central parsec.Despite the fact that their detection in GC is technically challenging, the temporal scattering measured from the SGR J1745-29 suggests that such a detection might be possible [27].
In the article, we investigate the rate of dynamical capture of PBHs by NSs in the GC, revising and extending existing calculations, considering in particular more realistic star interiors.We calculate the likelihood of collapse for a NS of a given age in that region, and assess the disruption rate of pulsars in the GC.
The paper is structured as follows: in Section II we discuss the formalism for energy loss of a PBH crossing a NS.In Section III we explain the NS interior solution used for our calculations and the effect this has on the energy loss.Section IV explores the relevant properties of the GC.In Section V we explain how to calculate the PBH capture rate from NSs.We then expand on how to calculate the likelihood of collapse for a NS, accounting for the time of orbital collapse for the PBH (Section VI) as well as dense environment effects (Section VII), and how to use this to obtain constraints on the abundance of PBHs from the observation of old NS in the GC (Section VIII).Throughout this paper we use units in which c = 1.

II. PATH OF PBHS AND ENERGY LOSS
In order to calculate the energy lost by the PBH through gravitational interaction with the NS, we start by deriving an expression for the PBH trajectory inside and outside the star.Following previous work from Génolini et al. [21], we use the Newtonian approximation and assume that the path of the PBH is not significantly perturbed by the NS interactions during the first crossing.This is a valid assumption since the kinetic energy of the PBH in the star is much larger then the energy lost by the PBH during one crossing.
Path of PBHs.For a NS of mass M and radius R, all possible approach orbits are fully characterized by just two parameters: the asymptotic velocity at infinity, v i , and the impact parameter b.The path of the PBH outside the star is described by a hyperbolic orbit, which is parameterized by a semi-minor axis equal to the impact parameter and a semi-major axis a and eccentricity e given by  where b = b/R is the reduced impact parameter and v ⋆ = (GM/R)1/2 is the typical velocity of the star.If b is smaller then the critical impact parameter b , the PBH will cross the NS and a new expression becomes necessary to describe the path of the PBH inside the star.Approximating the NS as a sphere of constant density, the path of the PBH, while inside the star, will be described by an ellipse centered at the NS core, being parameterized by a semi-major and semi-minor axis given by Patching the two regimes together, we get [21] where the angles ϕ i and ψ i are given in Table I and ε = 1 − α − /α + is the eccentricity of the ellipse parametrizing the internal path.This is the same parametrization as used in Ref. [21], but with corrected typographical errors in the equations for the total deflection angle, the semi-major and semiminor axis inside the star.
Energy losses.Black holes, if not electrically charged, can only interact gravitationally with surrounding matter, and PBHs are no exception.Assuming that charge and angular momentum are negligible, as recent analyses [33,34] seem to suggest for a large class of PBH formation scenarios 1 , there are only four relevant mechanisms for a PBH to lose energy while orbiting close or through a NS: • gravitational waves (GW) emission; • dynamical friction (DF); • accretion of NS matter onto the PBH; and • production of surface waves (SW).
If the total energy lost, ∆E tot , through the initial hyperbolic orbit is larger than the initial energy E i = 1 2 mv 2 i , the PBH will be stuck in a bound orbit and, excluding external forces, will eventually sink to the NS's center and consume it [40][41][42].
Gravitational Waves emission.A PBH orbiting close or passing through a NS, will lose energy due to GW emission.Given a PBH of mass m, the GW energy loss per unit angle outside the NS is given by [43] and inside the NS it will be [21] The energy loss per unit angle for the whole orbit can then be expressed as a piece-wise function in the form It is then straightforward to compute the energy loss from gravitational waves for a given orbit as In the case where the PBH does not cross the star, this can be solved analytically as [43] ∆E gw = 8 15 + Surface waves.Tidal effects are mostly negligible in PBH-NS interactions where the PBH does not cross into the star [44].However, if the PBH passes through the NS, ripples are induced on the surface, leading to an associated energy loss from the PBH.These ripples are referred to in the literature as surface waves, there is no exact solution for the associated energy loss from a realistic star, but Ref. [45] has shown that the result should be of order being a suitable approximation given that SW will be subdominant to accretion and dynamical friction for most of the impact parameters smaller than b c .
Accretion.While the PBH is inside the NS, some of the star's mass will be accreted by the black hole.This mass will have to be accelerated to the local PBH velocity, which will cause an accretion force in the direction opposite to the PBH's motion.The internal temperature of the NS at formation is around 10 11 K.The star then rapidly cools down through neutrino emissions, and during the lifetime of a pulsar, the internal temperature is T (10 Myr) ∼ 10 7 K [46], with an associated thermal velocity v th ≈ 10 −3 .Given that the PBH velocity within the NS is much higher than the local thermal velocity of the neutrons, we treat them as static in the NS reference frame.In the limit v ≪ c s (PBH velocity much smaller then the sound speed), including relativistic effects, the rate of mass capture for a PBH with speed v (with associated Lorentz factor γ) in a perfect fluid at equilibrium with density ρ and pressure P is [47] where d c is the critical impact parameter for a Schwarzschild black hole at the local PBH velocity, which is given by solving the set of equations [40] where is the Schwarzschild radius of the PBH.The collision with the captured mass can then be treated as perfectly inelastic, the consequent drag force experienced by the PBH will be Approximating the path of the PBH through the NS as described before, the total loss of kinetic energy is where ϕ 0 and ϕ 1 are the entrance and exit angle, respectively.
Dynamical friction.The nucleons approaching at a distance larger than d c from the center will not be captured but will still be deflected and hence take momentum from the PBH.This will also cause a force on the PBH, known as dynamical friction.In the NS reference frame, the change in momentum of a deflected neutron reads [40] where θ is the total deflection experienced by the neutron in the PBH reference frame.Expressing the integral in terms of x = R s /r, the deflection angle is given by θ All components of ∆p not parallel to the direction of motion of the PBH will cancel out by symmetry if one assumes that the density is approximately constant on the scales of interest.Adding a small correction by considering the relativistic momentum density for a perfect fluid to the result from Ref. [40], the force experienced by the PBH is given by The upper bound of the integral d max is defined as the distance at which the scattered nucleon no longer gains enough energy to be ejected from the Fermi sea [21], this being equivalent to the local binding energy µ(r).Assuming that all the NS matter is composed by neutrons, the local value of d max will be the impact parameter d which solves the equation where we have updated the equation given in Ref. [19] with the appropriate relativistic effects.The total energy loss from dynamical friction is then For both DF and accretion the calculations have been made under the assumption that the gravitational effects of the PBH dominate the neutron's behaviour locally.
Given that d max < O(10R s ) in the core of the NS -where most of the energy loss from DF occurs -we expect this approximation to hold.
The dynamical friction and accretion calculations are done under the assumption that the neutrons in the PBH's wake are approximately collisionless.This is a valid assumption in the limit c s ≪ v, which however starts to break down in the core of the NS, and in the case of the 2.12 M ⊙ NS does not hold at all since the sound speed is larger than the PBH velocity for most of the interior (although this is dependent on the chosen equation of state).
Implementing a fully collisional treatment might significantly enhance the predicted energy loss from both channels.Previous analyses have shown that there is a strong enhancement of dynamical friction when v ≈ c s , since a gravitational "sonic boom" creates a divergence in the force when the speed of the perturber is equal to the medium's speed of sound [48].This process has been extended to the case of relativistic dynamical friction in a collisional fluid [49] under the assumption of a weak gravitational field coming from the medium, which is not valid in the case of a NS, since the gravitational wake would be disrupted by the star's gravitational effect.We also stress that current collisional dynamical friction models are not suited to take into account the suppression coming from the binding energy of neutrons, which reduces the DF effects up to a factor of 10 in our calculations.

III. ENERGY LOSS FOR A REALISTIC NEUTRON STAR INTERIOR
The PBH energy loss in the NS interior is highly dependent on the local density, PBH velocity and binding energy per nucleon.Therefore, it is important to have an accurate description of the star's interior.
We use the BSk-20 model [50] to solve the TOV equations to find a numerical solution for the star interior; some of the relevant values are shown in Table II.The Brussels-Montreal (BSk) models for energy-density functionals are a series of analytical models based on zero range effective interactions between nucleons [51].These allow to obtain analytical approximations for the equation of state of a neutron star.We chose to work with theBSk-20 model since it has the best fit to the massradius data collected from the LIGO/VIRGO data from the GW170817 NS merger [52,53].
We use our interior solutions to numerically evaluate the Newtonian paths inside a NS.To these we superimpose the local GR velocities and interior values to compute the energy loss from dynamical friction and accretion on any given path.The GW energy loss in the interior is computed using the constant-density approximation as discussed in Section II, since ∆E gw is subdominant to the other energy losses in the star's interior, this should not have a significant effect on the final results.
The energy loss from each component and the total energy loss are computed for two NS interiors: one which represents a NS of mass 1.52 M ⊙ (which is close to the expected mean mass for MSP) and one which represents an extremal 2.12 M ⊙ NS (which is close to the maximum allowed mass for MSP) [54].The result for the 1.52 M ⊙ NS are shown in Fig. 1.
The inclusion of a complete description enhances the dynamical friction and accretion for PBHs at low values of the impact parameter b when passing through the star's core, but decreases the energy loss for PBHs moving only through the less dense crust around b ≈ b c .For a radial orbit and for a realistic interior, the result is a more then doubled in energy loss.Overall this results in a higher capture rate (by approximately 30%) for a NS with realistic interior in environments with high velocity dispersion (like the GC).
The higher density and mass of the heavier star also have a strong impact on the total energy lost by the PBH, allowing the extremal star to capture more PBHs.The difference in energy loss reaches a factor of 6 for radial orbit, this together with a slightly higher critical impact parameter causes a significantly higher capture rate for the more massive star.
We also run analyses for the BsK-19 and BsK-21 models [50], which represent the two extrema of our uncertainty on the mass-radius function for NS, representing denser and less dense NS models.Using these equations of state, we compute interior solutions for NS of the same mass (1.52 M ⊙ ), to keep the main parameter for the PBH energy loss constant among the different models.As expected, denser stars have a higher maximal energy loss but a smaller critical impact parameter.While the velocity dispersion is small, the capture rate is dominated by GWs, and are hence the same for the different equations of state (since the NS masses are the same), as the velocity dispersion increases, the effects from the changes in the critical impact parameter and maximal energy loss compete against each other, with the denser star leading to a higher capture rate when the velocity dispersion is high, similarly to the previously discussed cases.Overall, the associated changes in capture rate are at most ∼ 30%.

IV. CENTRAL PARSEC
In order to analyze the PBH capture undergone by NSs in the GC one must have a good description of the dark matter density, stellar mass and velocity distribution within this region.In this section we collect these relevant pieces of information: stellar distribution, dark matter  distribution and the expected velocity distribution of dark matter in the reference frame of local NSs.
Stellar density.The stellar density near the Galactic center can be modelled from high-resolution adaptive optics imaging [55], which gives as the best-fit description of the central stellar cusp where r is the radial distance from Sagittarius A*.The function is split in two regimes with α = 1.4 for r < 0.39 pc and α = 2.0 for r ≥ 0.39 pc.From this one can derive an approximate stellar number density as n ⋆ (r) = ρ ⋆ (r)/⟨M ⋆ ⟩, where ⟨M ⋆ ⟩ = 1.3 M ⊙ is the average mass for the Salpeter mass function.
Pulsars age.An important parameter for dynamical capture in pulsars is the age of the NS, given that likelihood of PBH capture is directly proportional to age of the NS.Regular isolated pulsars have an expected lifetime of 10 -100 Myr [56], meaning their expected age should be of order ∼ 10 Myr (although younger pulsars are easier to detect [56]).
It is harder to directly estimate the age of MSPs given that the standard mechanisms.The characteristic age, does not work for MSPs, since they have received angular momentum from an external source [57].The expected age of a MSP in the GC can still be estimated by looking at other old stars in the GC, this process gives an expected age of ⟨T GC ⟩ = 10.4Gyr [58].
PBH density.We assume that the radial PBH density and velocity profiles matches those of the dark matter.In the case where they compose only a fraction of the total dark matter, we simply rescale the density profile normalization.The dark matter distribution in the central parsec is unfortunately hard to infer from observations, since the gravitational potential is completely dominated by baryons.The best strategy to estimate the density profile is to extract it from numerical simulations which include baryonic effects, and which successfully reproduce the rotation curve and morphological properties of the Milky Way [59][60][61].Here, we consider two models for the dark matter density profile.The first is the standard Navarro--Frenk-White (NFW) model, which provides an excellent fit in dark matter simulations without baryons.
The second is the "adiabatic-growth" model, labelled "A" below.The latter assumes that the baryons contracted adiabatically and symmetrically within the pre-existing dark matter halo, compressing the dark matter and increasing its density [62].When applied to an halo with an initial NFW profile, the result is a steeper profile inside the solar circle, with power-law index 1.5 [63].For both the NFW and the adiabatic profiles, we considered the case with (subscript 's') and without (subscript 'n') a central 'spike' induced by the adiabatic growth of the central black hole [64,65], as described in Ref. [66] and shown in Fig. 2. The general form of the Galactic dark matter density inward of the solar system reads where ρ ⊙ represents the local density in the solar system, γ c is the power-law index for the cusp, r s = 0.4 pc is the distance at which the DM spike behaviour begins, r in = 8 × 10 −7 pc is the smallest distance at which a particle can orbit the black hole, and γ s is the power-law index for the spike.For the non-spiked case γ s = γ c .The relevant values for all models are shown in Table III.
PBH velocity distribution.The velocity distribution for a spherically-symmetric isothermal system, as NFW,n NFW,s A,n A,s ρ ⊙ [GeV/cm 3 ] 0.3 0.3 0.5 0.5 γ c 1 1 1.5 1.5 γ s 1 2.33 1.5 2.4 is the dark matter halo, should approximately follow a Maxwellian distribution [67].From the virial theorem, one has where M c (r) is the total mass contained in the sphere of radius r centered at Sagittarius A*, and v rms is the root-mean-square velocity.M c is obtained by integrating the stellar-DM density up to a sphere of radius r, and then adding the mass of Sagittarius A* mass, M Sag = 4.1 × 10 6 M ⊙ [68] as the central-mass value.The velocity distribution can then be written as The approximate velocity distribution obtained in this way gives results close to the observed velocity dispersion near the galactic center [69], although the data has a slight suppression of radial velocity, likely due to interaction with the dense environment and with the supermassive black hole [70].This bias might be larger at closer distances from Sagittarius A* but in the distances considered in this paper this should at most add O(10%) corrections to the velocity distribution, and more likely negligible corrections if v rms is not substantially changed.There is little difference in expected velocity between different dark matter distribution models, especially in the central parsec, which is of interest here.For r < 0.1 pc, the velocity dispersion is dominated by the supermassive black hole mass such that v rms ∝ r −1/2 .Equation ( 24) describes the velocity distribution of dark matter for a static observer and is the one usually used in most of the literature on dynamical capture.However, any neutron star in the central parsec will also be orbiting around Sagittarius A*.We calculate the velocity distribution by moving to the frame of reference of a NS moving with velocity v NS = v rms through the dark matter halo.This is the case for a NS moving on a circular orbit in the GC, which gives the velocity distribution This velocity distribution reduces the capture rate by a factor 4.5 with respect to the standard Maxwell distri- bution in the small-v/v rms limit, which is the relevant regime for low-mass PBH capture.

V. CAPTURE RATE
Assuming a locally-uniform PBH distribution, it is possible to calculate the NS capture rate G, which is the rate at which PBHs are captured into a bound orbit at any point in the Galaxy [21] where b G (v) is the maximum radius at which capture can occur for a given velocity v (see Fig. 3 for numerical evaluation), its distribution f (v) (in the NS reference frame), and the local number density of the PBHs n PBH = ρ PBH /m.The capture rate might be enhanced by a factor of 3.5 for tight binaries according to numerical simulations [71].This is relevant for MSPs which are expected to be in tight binaries for a large portion of their existence, thereby accelerate their angular velocities up to the observed high values [72].Including this effect, it is then possible, for a MSP moving through a circular orbit in the Galactic center, and using the velocity distribution discussed before, to plot the capture rate/velocity-dispersion graph for different masses as shown in Fig. 4. The PBH capture rate is systematically enhanced for the heavier neutron star (2.12 M ⊙ ), where the difference is larger (up to a factor of five) for higher velocity dispersions, since the larger maximum possible capture velocity v max = 2∆E max /m of the denser star becomes more relevant to the capture statistics in such environments.The capture rates for different masses tend to converge at high velocities.This can be explained by analysing the behaviour of the capture radius-velocity relation as given in Fig. 3.For all masses, b G converges to the critical impact parameter before rapidly collapsing to zero as the maximum possible energy loss of a PBH passing through a radial trajectory, ∆E max = ∆E(b = 0), becomes insufficient to disperse the original kinetic energy of the PBH.This first convergence is clearly caused by the sudden jump in energy loss as the PBH crosses the NS's surface, which means that, for velocities close to v max , the capture radius is approximately constant and similar to b c .
In high-velocity-dispersion environments like the GC, most of v max for different masses are concentrated at the beginning of the velocity distribution where the velocity is much smaller than the velocity dispersion.This property extends to more masses the higher the value of v rms .Expanding around v ≪ v rms to first order, one can approximate the velocity distribution of PBHs as The same proportionality holds also for the regular Maxwell distribution discussed in the literature.The largest contribution to G comes from velocities close to v max where b G ≈ b c ≈ R v ⋆ /v, meaning the integrand is roughly proportional to v. The capture rate will then be proportional to where in the last step, one needs to use the property ∆E ∝ m 2 discussed in Appendix B. Therefore, the convergence of the capture rate is a general feature of the capture dynamics for PBHs with Maxwell-like velocity distributions, and it is more pronounced in environments with larger velocity dispersions.As a consequence of this convergence, the final PBH abundance constraints exhibit a flat behaviour as a function of mass.

VI. TIME TO COLLAPSE
Once the PBH is captured in a bound orbit with a NS, it will take some time for the PBH to lose enough energy to fall to the center of the star and consume it.If the time required for the PBH to destroy the star is longer than the age of the star, then the capture might not yield any observable consequence and such captures must be ignored in the calculations.
Immediately after capture, the PBH will be in an elliptical orbit of energy where a 0 is the semi-major axis of the initial captured orbit, and ∆E tot is the total energy lost by the PBH by summing over all the energy loss channels from Section III.
Assuming that the orbit is such that the distance of closest approach r p is the same [21], and approximating the approaching orbit as nearly parabolic, yields b 0 ≈ b a 0 /a .The collapse of the orbit can then be split in two cases: the GW-induced collapse for b > b c .and the collapse dominated by the energy lost inside the NS for b ≤ b c .
In the first case (b > b c ), the orbital decay will be completely dominated by the GW energy loss.The decay time for a highly eccentric binary is given by [73] which can be approximated as If the energy lost by the PBH through gravitational waves is much larger than its initial energy, one can utilise Eq. ( 8) such that E orb ≈ −∆E GW .Under the assumption which is similar to the result presented in Ref. [21] but with slightly different exponents.
In the second case (b ≤ b c ), the energy loss will be dominated by the extremely fast intervals when the PBH crosses through the NS.Under the assumption that the first few orbits after capture follow a similar path through the neutron star as the initial hyperbolic orbit (which should be the case given the small effect on the kinetic energy of the PBH while inside the NS), the energy loss for each orbit around the NS will be approximately equivalent.One can then approximate the total decay time as the sum of the orbital periods T ell = 2π a 3 n /GM that the PBH goes through as the orbit decays.This gives where ζ is the Hurwitz zeta function.T decay is roughly proportional to m −3/2 , and will therefore be considerably higher for lower masses.This puts constraints on the minimum PBH mass which can effectively disrupt a NS within the total age of the star.For an old MSP of age 10 Gyr and mass 1.52 M ⊙ , the minimum mass able to cause a collapse is m min ≈ 10 −16 M ⊙ .Once the orbit has completely decayed, the mass accretion rate of the black hole can be modeled by the spherical Bondi accretion, which gives as the time for the NS to collapse into a black hole [21] T In the last stages of stellar collapse, the fast rotation of the NS will slow down the accretion and the new regime will be Eddington-like.However, this will be relevant only in the last moments of the star's collapse [74].This does not impact our analysis on the scale of pulsar lifetimes.T collapse is then an order of magnitudes smaller than T decay for all PBH masses of interest and can therefore be ignored.

VII. EXPECTED COLLAPSE NUMBER AND SURVIVAL LIKELIHOODS
Given a NS of age of T ⋆ , it is possible to find the expected number of captures leading to collapse within its lifetime.We call this the expected collapse number C. Of course, no star can undergo more than one collapse, but this number just represents the average number of PBH encounters capable of leading to collapse that the star should have undergone in its lifetime.We define the critical radius b C (m, v) for collapse, as the maximum distance such that a PBH of mass m and velocity v is captured and falls into the NS within a time T < T ⋆ .PBH masses for which any capture orbit takes more than the age of the MSPs to collapse will have b C = 0 for any velocity.As an example in Figs. 5 and 6, PBHs with mass m min ≈ 10 −16 M ⊙ or smaller cannot cause a collapse within the age of the star.Hence, they have b C = C = 0 and do not appear in the log-log graph.
The orbital decay time T decay (b, v) can be comparable to the star's life.In that case, there is a shorter time window for capture and collapse T eff = T ⋆ − T decay (b, v).In order to account for this, one can start by defining an effective area of capture for a PBH with mass m and velocity v as This gives the expected collapse number, In the following calculations, we approximate the decay time as the radial orbit decay time T d , such that T d = T decay (b = 0, v), which is a good approximation in highvelocity-dispersion environments for masses much larger than m min .The time dampening is strongest for lower mass PBHs with larger orbits.An important question is if, the case of MSPs, the companion star might negatively impact the orbital collapse.It is possible that the PBH might be slingshot out of the system through the gravitational interaction with the other two bodies.Numerical simulations are necessary to check the actual effect on the expected collapse number.
The expected number of collapses C for a given NS can be interpreted as the mean number of PBHs, with the correct characteristics leading to the star's collapse, which should have crossed the NS in its lifetime.Given that PBH captures are independent of each other, the likelihood for a NS to have encountered a given number of PBHs with the correct mass, velocity and impact parameter to lead to a collapse will be given by a Poisson distribution with mean C. The probability of the NS's survival will then be given by P (Collapse) = 1 − e −C .In general the expected suppression rate of an initial population of N 0 neutron stars at a distance r from the GC after a time t will be in the form Since C is roughly proportional to t, the suppression will be much higher for older populations of NSs.We find that even assuming an adiabatic and spiked dark matter halo, a NS mass of 1.52 M ⊙ (which is larger than the average mass of regular pulsars [75]), and ignoring dense environment issues which will be discussed later in this section, the survival rate after the average pulsar lifetime is N (10 7 yr) ≈ N 0 .This holds even at distances of 10 −3 pc, meaning it is impossible for the PBHs to disrupt a significant fraction of pulsars at galactocentric distances larger than that.This clearly cannot explain the lack of pulsars detection in the inner 1 -10 parsecs [27].
This means that dynamical capture of PBHs cannot be the dominant channel to explain the missing pulsar problem, unless other processes significantly increase the PBH capture rate (e.g.capture from adiabatic contraction see Ref. [41] and more recently in Ref. [76]).

VIII. CONSTRAINTS AND DENSE ENVIRONMENT EFFECTS
As mentioned in Ref. [77], the assumption that the PBH-neutron star system is gravitationally isolated from the surrounding distribution of stars breaks down in dense stellar environments.The apoastron r a for a PBH of mass m, initial velocity v and impact parameter b is given by where a 0 is the semi-major axis of the initial bound orbit given by Eq. (29).By approximating the local stars as uniformly distributed, we get an expected distance ) as the cutoff of the gravitational influence of the NS, where M r = ⟨M ⋆ ⟩/M .In the case of the central few parsecs, it is important to consider the pull of the central mass when compared to the NS at the PBH's apoastron.The distance at which the pull of the NS and the central mass are equivalent is r c = r M/M c .We can then impose the condition r a ≤ min(r ⋆ , r c ) to all captures to make sure our two-body-orbital-collapse formalism is a valid approximation.In practice, at the Galactic centre we always find r c < r ⋆ .Redefining b C to include this property for the bound orbit, allows us to re-evaluate it numerically for different radii.This turns out to exclude collapses for a wide range of low masses.The closer a NS gets to the GC, the more masses are excluded from capture due to the gravitational effects of the central mass.
Assuming a Dirac delta for the mass function of PBHs, we set for a given PBH mass an upper limit on the PBH abundance, such that the likelihood of such an observation to occur for the given PBH population is under 5%.We consider the possibility of setting such constraints from future observations of one or more MSPs in the central parsec.As mentioned before, we take the age of these objects to be T MSP = 10.4Gyr.
In the case of a single observation it is straightforward to use C ∝ ρ PBH to find f PBH = ρ PBH /ρ DM , which would give a likelihood of survival for the MSP under 5%, where C DM is the expected collapse number for ρ PBH = ρ DM .Using this, we can plot the constraints on different PBH masses from the detection of a MSP at different distances from the GC, for different DM distributions and for MSPs of different masses.The possible scenarios for constraints from a single MSP observation are shown in Fig. 8. Constraints above the upper limit f PBH (m) = 1 are for unphysical scenarios in which the PBH abundance exceeds the DM abundance in the Milky Way.No physical constraint can be put from the observation of a single MSP in the NFW case, with or without a spike.However, as discussed before it seems that the adiabatic case should be the more likely physical scenario, in which case some constraints might be possible, especially in the presence of a spike around Sagittarius A*.There is a substantial enhancement of constraints from the observation of an extremal 2.12 M ⊙ MSP.In the adiabatic-spik case both an observation of an average mass or an extremal MSPs are enough to put physical constraints, while in the adiabatic non-spiked case only the observation of an extremal MSP would be enough.The most optimistic case for constraints from a single observation are shown together with previous constraints in Fig. 7.
The possibility to set constraints with future observations of pulsars located very close to the Galactic center depends on the detectability of these objects despite the strong inter-stellar scattering expected in this region [78].Under our approximations, the above constraints effectively vanish in the currently unconstrained "asteroid mass" window (10 −16 -10 −10 M ⊙ ), due to the fact that Sagittarius A* and the dense stellar environment disrupt the orbits of captured PBHs at lower masses.It might however be possible to constrain this mass range by observing an old MSP in an ultra-faint dwarf with similar dark matter density and velocity dispersion but with lower baryon density.However, this does not seem to be a promising avenue both because we don't expect to see many pulsars in these regions, due to the low baryonic fraction.Furthermore, even if we could observe one, it would be hard to reconstruct the distance from the center of the dwarf with high precision.
In the future, it may become possible to set interesting constraints in this mass region, by estimating the global evolution of the PBH phase space due to interactions with the global compact-star population in the central parsec, and three-body gravitational interactions with other stars, in the gravitational potential of Sagittarius A*, and how this may impact PBH capture in the region.

IX. CONCLUSIONS
In this work we have studied PBH-NS interactions and orbital dynamics.We have improved previous energy-loss calculations by treating the interaction of PBHs with NS using a realistic solution for its interior, and refining the treatment of the interaction dynamics and collapse likelihood.FIG.7: Constraints on the PBH abundance arising from the most optimistic case of an observation of a single 2.12 M ⊙ MSP at 10 −3 pc in the Galactic centre, assuming an adiabatically compressed dark matter profile with a spike around Sgr A* (green).Existing constraints on the PBHs abundance are also shown for reference.Data and code for this figure are taken from Ref. [6].
We found that dynamical capture of PBHs cannot explain the missing pulsar problem.The velocity dispersion is simply too high for the PBH to efficiently disrupt the star within the pulsar's lifetime (∼ 10 Myr), even under the most optimistic assumptions.For the lower masses of the PBH's unconstrained region there are still an average of 10 3 PBH-NS crossing in the pulsar lifetime at 10 pc from the GC, with more crossings closer to the GC.We have also shown that it is possible to set stringent constraints on sub-solar mass PBHs from MSP observations in the Galactic center.The main uncertainty in the calculation arises from the many-body dynamics of the central parsec.This affects the dynamics and capture of PBHs, including hyperbolic encounters with other stars and the effect of binary companions.These are expected to be present, allowing for the formation of MSPs, on the PBH's orbital collapse.It is important to investigate how all these processes will impact the orbital collapse of the captured PBHs.More constraints might become possible with a better understanding of the interactions between PBHs and the global population of pulsars at the Galactic center.
These constraints would be significantly enhanced if multiple MSPs were observed near the Galactic center.This would be indicative of the presence of a large population of objects in a region with large dark matter density.
We discuss for completeness the case of multiple pulsars detection in Appendix A.
In our analysis, both the dynamical friction and accretion calculations are performed under the assumption that the neutrons in the PBH's wake are approximately collisionless.We have argued that a fully collisional treatment might significantly enhance the energy loss from both channels.A detailed investigation of the collisional effects would likely strengthen the constraints derived in this paper and extend them to lower PBH masses.

FIG. 1 :
FIG.1: Total energy loss and energy loss for each mechanism versus impact parameter for a PBH of mass 10 −9 M ⊙ and initial velocity 10 −3 , and a NS of mass 1.52 M ⊙ .

3 ]FIG. 2 :
FIG. 2: Dark matter density versus distance from Sagittarius A*, the black line represents the ISCO.Spiked (full lines), spikeless (dotted lines), NFW (blue) and Adiabatic (red) profiles are shown.See main text for further details.

FIG. 3 :
FIG. 3: Capture radius b G as a function of PBH's initial speed for a NS of mass 1.52 M ⊙ .

FIG. 5 :
FIG.5: Critical radius for collapse as function of PBH initial velocity for a NS of mass 1.52 M ⊙ and age ⟨T GC ⟩.

FIG. 6 :
FIG.6: Expected number of collapses as a function of the local v rms , with local PBH density of 1 GeV/cm 3 , for a MSP of mass 1.52 M ⊙ and age ⟨T GC ⟩.

TABLE I :
Angles for the orbital formalism.

TABLE II :
Relevant values for the two NS interior solutions.

TABLE III :
Relevant values for the four dark matter models.
FIG.4: Expected number of captures as a function of the local v rms , with local PBH density of 1 GeV/cm 3 , for a MSP of mass 1.52M ⊙ and age < T GC >.