Gravitational-wave Emission from a Primordial Black Hole Inspiraling inside a Compact Star: a Novel Probe for Dense Matter Equation of State

Primordial black holes of planetary masses captured by compact stars are widely studied to constrain their composition fraction of dark matter. Such a capture may lead to an inspiral process and be detected through gravitational wave signals. In this Letter, we study the post-capture inspiral process by considering two different kinds of compact stars, i.e., strange stars and neutron stars. The dynamical equations are numerically solved and the gravitational wave emission is calculated. It is found that the Advanced LIGO can detect the inspiraling of a $10^{-5}$ solar mass primordial black hole at a distance of 10 kpc, while a Jovian-mass case can even be detected at megaparsecs. Promisingly, the next generation gravitational wave detectors can detect the cases of $10^{-5}$ solar mass primordial black holes up to ${\sim}1$ Mpc, and can detect Jovian-mass cases at several hundred megaparsecs. Moreover, the kilohertz gravitational wave signal shows significant differences for strange stars and neutron stars, potentially making it a novel probe to the dense matter equation of state.


INTRODUCTION
Primordial black holes (PBHs; Zel'dovich & Novikov 1967) can be formed through a vast range of mechanisms. For example, they may be generated due to the density inhomogeneity in the early universe (for a review, see Khlopov 2010). Although having not been directly detected yet, PBHs are considered to be a candidate for dark matter (for a review, see Carr & Kühnel 2020). To constrain PBHs' composition fraction of dark matter, the event rate of collisions between planetary-mass PBHs and compact stars has been widely discussed and the corresponding electromagnetic emissions have been extensively studied (e.g., Capela et al. 2013;Pani & Loeb 2014;Fuller et al. 2017;Abramowicz et al. 2018;Génolini et al. 2020). Interestingly, during such a collision, the PBH-compact star system will also emit gravitational waves (GWs) as the PBH accretes matter from the compact star. Kurita & Nakano (2016) and East & Lehner (2019) studied the GW signals from an accreting PBH, which is assumed to have plunged into its compact companion and stays at the exact center of the neutron star (NS). Horowitz & Reddy (2019) studied the GW emission during the inspiraling process of a PBH inside an NS, but their treatment of the dynamics and the compact star structure is still preliminary. Génolini et al. (2020) analyzed the process of a PBH plunging into an NS by considering dynamical friction, accretion, and GW emission. For a trapped PBH, they found that both the frequency and amplitude of the GWs are constant during the inspiral. However, their dynamical equations are mainly appropriate for PBHs in deeply subsonic motions. Also, they assumed that the NS has a homogeneous structure when calculating the motion of the trapped PBH, causing additional deviation in the derived GW waveforms.
On the other hand, the equation of state (EoS) of dense matter determines the structure of compact stars. According to the strange-quark matter hypothesis (Itoh 1970;Bodmer 1971;Farhi & Jaffe 1984;Witten 1984), pulsars may actually be strange stars (SSs) consisting of strange-quark matter (Alcock et al. 1986). The strange-quark matter is self-bound, thus strange dwarfs and even strange planets can stably exist. However, a 1.4 M SS has a radius very similar to that of a normal NS with comparable mass, thus it is hard to distinguish between these two types of compact stars via observations (Geng et al. 2015). An interesting method to identify strange quark objects is to search for very close-in binary systems containing a strange planet and a compact star (Geng et al. 2015;Kuerban et al. 2019Kuerban et al. , 2020. When the orbital radius of a planet is less than ∼ 5.6 × 10 10 cm or the orbital period is less than ∼ 6100 s, then it cannot be a normal matter planet, but should be a strange planet, because the tidal force is too strong to allow any kind of normal matter planets to stably exist there (Huang & Yu 2017;Kuerban et al. 2019Kuerban et al. , 2020. However, in these close-in "planetary" systems, there is still a possibility that the planetary-mass object is actually a PBH. We thus need to further scrutinize its nature. GW astronomy may shed new light on the study of compact star structure. Tidal deformability and maximal mass of compact stars, which can be hinted at through GW signals from double compact star mergers (for a review, see Guerra Chaves & Hinderer 2019), may reflect the internal composition and structure of compact stars (e.g., . In fact, the observations of GW170817 have already put useful constraints on the tidal deformability of "neutron stars" (Annala et al. 2018;Abbott et al. 2018Abbott et al. , 2019. However, these constraints are still too weak to pin down the EoS (Lai et al. 2019;Shibata et al. 2019;Mondal & Gulminelli 2021). In this Letter, we will study the GW emission produced by a PBH inspiraling inside a compact star. Different from a strange planet, a PBH can inspiral and tunnel inside the compact star, producing special GW signals. Two kinds of EoSs will be assumed for the compact stars in our calculations, i.e., a strange-quark matter EoS and a hadronic matter EoS. The results will help us judge whether the close-in planetary-mass object is a strange planet or a PBH.
The structure of this Letter is as follows. In Section 2, we describe the model by setting up the equations of motion and compact star structure. Numerical results on the dynamics are then presented in Section 3. In Section 4, the GW emissions are calculated and analyzed. The effects of orbital eccentricity are considered in Section 5, and the event rate is calculated in Section 6. Finally, we summarize and discuss our results in Section 7. Throughout this Letter, a natural unit system of c = G = = 1 is adopted.
2. MODEL SETUP As a PBH inspirals inside a compact star, it gradually loses its orbital kinetic energy through interactions with compact star matter. We focus on three main channels -dynamical friction, accretion, and GW emission (Génolini et al. 2020). Moreover, we assume that the trajectory of the PBH within the compact star is quasi-circular and do not consider the effects of ellipticity in most of our calculations. This condition can be guaranteed by the tidal dissipation (Ogilvie 2014) and various circularization effects (Macedo et al. 2013). Anyway, a brief discussion on the effects of eccentricity will be presented in Section 5.

Dynamical Friction
A wake is produced when the PBH (of mass m PBH ) moves through the surrounding medium (with a density of ρ and the sound speed c s ). The wake will then exert a gravitational drag force on the PBH, known as the dynamical friction force (Ostriker 1999). For a circular orbit, the force can be written on the basis of the radial (r) and azimuthal (φ) components as (Kim & Kim 2007) where v is the relative speed of the PBH with respect to the medium. The coefficients I r and I ϕ are functions of the Mach number (M = v/c s ) and the distance (r) between the PBH and the center of the compact star (Kim & Kim 2007): ( and where r m = √ em PBH /(2v 2 ) (Cantó et al. 2011). However, the above expression of I ϕ fails for M → 0 because it gives an unphysically positive azimuthal drag force. In fact, in this deep subsonic phase, the circular-orbit dynamical friction is expected to be similar to the case of a straight trajectory (Ostriker 1999), i.e., I ϕ = ln [(1 + M)/(1 − M)] /2 − M → M 3 /3. To account for this asymptotic behavior, we use the following polynomial expansion for M 1, This expression reduces to the solution of Ostriker (1999) for M → 0. When M increases, it smoothly transitions to Equation 3 at M = 0.08588. So when M < 0.08588, Equation 4 is applied instead of Equation 3 in our calculations.

Accretion
As the PBH accretes matter from the compact star, it also accumulates a negative momentum, resulting in a drag force of (Edgar 2004;Génolini et al. 2020) For a PBH with a finite velocity, the accretion rate is (Bondi 1952;Shima et al. 1985;Edgar 2004) where λ is the accretion eigenvalue depending on the EoS of the accreted matter. Equation 6 is applicable for m PBH M CS (Richards et al. 2021), where M CS is the mass of the compact star.

Structure of Compact Stars
The dynamical friction and accretion explicitly depend on ρ and c s . Therefore, we shall obtain the compact star structure before further studying the equations of motion of the PBH. The structure of a non-rotating compact star is available by solving the Tolman-Oppenheimer-Volkoff equation (Tolman 1939;Oppenheimer & Volkoff 1939) where P is the pressure, and M r is the mass within the radius r so that dM r /dr = 4πr 2 ρ. To solve Equation 7, one needs the EoS of compact star matter. Two typical kinds of compact stars, SSs and normal NSs, are considered in our study. For SSs, we adopt the simple bag model with massless quarks (Farhi & Jaffe 1984), whose EoS reads P = (ρ − 4B)/3. The bag constant B is taken as 57 MeV fm −3 . For NSs, we use the hadronic BSk 24 EoS 1 (Potekhin et al. 2013;Pearson et al. 2018). The sound speed is calculated from c s = dP/dρ. The density and sound speed profiles of a 1.4 M SS and NS are shown in Figure 1. The SS has a smaller radius (R SS = 11.0 km) than the NS (R NS = 12.6 km). Moreover, the SS has a rather uniform density 0.01 0.1 1 10 100 10 -6 10 -5 10 -4 10 -3 10 -2 Figure 2. Evolution of the PBH mass. For the dash-dotted, dotted, dashed, and solid curves, the initial mass of the PBH is taken as 10 −3 , 10 −4 , 10 −5 , and 10 −6 M , respectively. The mass of the compact star is 1.4 M , which is assumed to be either an SS (thick curves) or an NS (thin curves). and sound speed profile with a sharp edge, while the NS's density and sound speed vary significantly from its center to surface. These differences indicate that the drag force exerted on the inspiraling PBH may be different in these two cases.
For the parameter λ, although there are analytical expressions for EoSs when the adiabatic index is Γ 1 (Richards et al. 2021;Aguayo-Ortiz et al. 2021), no simple expression is available when Γ < 1 (at the phase transition region of a NS; Potekhin et al. 2013) and Γ → ∞ (at the surface of a SS; Xia et al. 2021). So, we use a commonly-used value of Γ = 4/3 in our calculations, which naturally gives λ = 1/ √ 2 (Génolini et al. 2020;Aguayo-Ortiz et al. 2021). Such a value is appropriate for most of the density range in NSs and SSs (Potekhin et al. 2013;Xia et al. 2021).

Equation of Motion
During the inspiral, we assume that the compact star structure inside r stays unchanged and the whole compact star remains spherically symmetric to its center. Thus the equation of motion can be written in the relative-motion frame as where M = M CS + m PBH is the total mass of the system. The last post-Newtonian term (Blanchet 2014) is introduced to account for the GW energy-loss. The PBH is initially assumed to be located at the surface of the compact star, with a Keplerian velocity ( M/R X , X stands for SS or NS) in a circular orbit. We terminate our calculation when m PBH /r = 1/12, which means r reduces to a value comparable to the innermost stable circular orbit of the PBH and the post-Newtonian method loses its accuracy (Blanchet et al. 2011).

BINARY EVOLUTION
We have solved Equation 8 numerically. For the initial mass of the PBH, we take four typical values, i.e., 10 −3 , 10 −4 , 10 −5 , and 10 −6 M . As for the compact star, we assume that it is either an SS or an NS, both with an initial mass of M 0 = 1.4 M . Figure 2 shows the evolution of the PBH mass. As expected, a PBH with a larger initial mass accretes faster. Moreover, the mass of the PBH interacting with an SS increases faster than that of the corresponding NS case. This is due to the large surface density of the SS. Note that in all our calculations, the condition of m PBH M CS is satisfied so that Equation 6 is applicable. Figure 3 shows the evolution of the separation between the PBH and the center of the compact star. We see that the PBHs quickly inspiral inward as its mass increases. At the end of our calculation, the PBH almost settles at the center of the compact star, which means the remnant matter of the compact star will soon be completely swallowed by the black hole (Blanchet et al. 2011). In the NS cases, the separation decreases sharply in the range of 9 km r 12 km. This is because the sound speed and density of the NS increase quickly in the region. As a result, the Mach number drops to an intermediate value, leading to an enhancement in the dynamical friction (Kim & Kim 2007). Figure 4 illustrates the evolution of the quantity m PBH r 2 . It is interesting to note that Génolini et al. (2020) suggested that m PBH r 2 is an adiabatic invariant during the inspiral. Our Figure 4 shows that m PBH r 2 is roughly constant only at the very early stage of the inspiral, while it generally decreases in the later process. The later decreasing of m PBH r 2 may be caused by two reasons. First, the constant-m PBH r 2 approximation is only valid in the deeply subsonic regime, but the PBH moves at a speed comparable to or larger than the sound speed. Second, when the PBH decelerates to a subsonic speed, it becomes very massive and it is also near the center of the compact star. As a result, the post-Newtonian term dominates the dynamical evolution.
It would be helpful to compare the individual contributions of different energy-losing channels. The energy dissipation rates of the damping forces areĖ DF = F DF · v andĖ acc = F acc · v, while the quadrupole GW emission has a power of (Maggiore 2007) Taking the initial PBH mass as 10 −6 M , we have plot −Ė of the dynamical friction, accretion, and GW emission against the PBH orbital radius in Figure 5. We can see that at small radii where the PBH is deeply subsonic, the accretion dominates, but the dynamical friction is still non-negligible. At larger radii, the dynamical friction dominates over accretion, especially in the NS case where the PBH has an intermediate Mach number. In the NS crust, where the density is very low, only the GW emission has a significant contribution, because the dynamical friction and accretion are very ineffective. Note that Equation 9 underestimates the GW power when r approaches the innermost stable circular orbit of the PBH (Maggiore 2007;Blanchet et al. 2011). For all the three energy-losing channels, the dissipation rate generally scales asĖ ∝ m 2 PBH . Note that when M 1, an additional logarithmic term presents in the dynamical friction, but the dissipation rate still does not deviate too much from theĖ ∝ m 2 PBH pattern. As a result, the relative dissipation rates of the three energy-losing channels are almost irrelevant to m PBH .

GRAVITATIONAL WAVE SIGNAL
Under the quadrupole approximation, the leading-order GWs of a point mass moving inside a spherical object are similar to the case of two point masses (Nazin & Postnov 1995;Ginat et al. 2020). As a result, the GW waveforms of a PBH moving inside a compact star are (Creighton & Anderson 2011)  where µ = m PBH M CS /M is the reduced mass, D L is the luminosity distance, ϕ is the orbital phase. Here the system is assumed to be face-on. Taking the distance as D L = 1 kpc, the GW strain amplitude h = 4µv 2 /D L is shown in Figure 6. Generally, the strain amplitude h decays with time. But it is interesting to note that there is obvious difference between the two curves corresponding to SS and NS. Thus the GW signal may be used to probe the EoS of dense matter. To decide whether the GW signal can be detected by a particular GW detector or not, it is useful to calculate the GW strain spectral amplitude (Moore et al. 2015;Zou et al. 2020) and compare it with the sensitivity of GW detectors. Hereh(f ) is an average of the Fourier transforms of h + and h × at frequency f . Actually, h f is the square root of the GW power spectral density. In Figure 7, we illustrate h f in a frequency range of M 0 /R 3 X /π f 4ρ c /(3π), where ρ c is the central density of the compact star. Beyond this main frequency range (3 -5 kHz) of the GW signals, our calculation of h(f ) is significantly polluted by spectral leakage. Comparing h f with the sensitivity curves of various GW detectors, we can see that the GW amplitude from such a PBH system containing a 10 −5 M black hole at 1 kpc away is about one magnitude higher than the sensitivity curve of the current Advanced LIGO 2 . It means that the Advanced LIGO and the upcoming LIGO A+ upgrade can safely detect the inspiraling of a 10 −5 solar mass primordial black hole at a distance of 10 kpc. For a Jovian-mass PBH, due to the much stronger GW emission, the detection horizon can even be pushed to megaparsecs. The next generation detectors such as Einstein Telescope 3 and Cosmic Explorer 4 , with sensitivities almost two magnitudes better, can even detect the GWs from an inspiraling 10 −5 M PBH at ∼ 1 Mpc, and detect Jovian-mass cases at several hundred megaparsecs. Figure 7 clearly shows that although GWs in the NS case and the SS case have similar amplitudes, the shape of the detailed h f curves actually is quite different for these two EoSs. For example, h f in the SS cases generally decreases homogeneously with the increase of the frequency, while h f in the NS cases shows a concavity at intermediate frequencies.
As a result, GW detectors working in kilohertzes can hopefully help us probe the EoS of dense matter through signals from PBHs inspiraling inside compact stars.

EFFECTS OF ECCENTRICITY
The spatial speeds of PBHs may distribute in a wide range. High-speed PBHs could not be captured by a compact star. Only those PBHs with a relative speed small enough would be successfully captured by a compact star. For simplicity, the PBH is assumed to be initially in a circular orbit in our calculations. However, a PBH could be accelerated to a speed significantly larger than the circular Keplerian speed when it is gravitationally captured by the compact star, thus an orbit with a large eccentricity is essentially possible (Génolini et al. 2020). At such an early stage of the capture, the system does not produce inspiraling GW signals as described in Section 4, but generates intermittent GW emissions at each passage of the perihelion as pointed out by Génolini et al. (2020).
If the perihelion is located outside the compact star, the eccentricity will gradually decrease due to GW emission (Peters 1964;Blanchet et al. 2011;Creighton & Anderson 2011;Blanchet 2014). Moreover, the tidal dissipation can also circularize the orbit (Ogilvie 2014). In most cases, the circularization will be sufficient and the PBH enters the compact star in a quasi-circular orbit.
On the other hand, if the PBH enters the compact star with a substantial eccentricity, or even its perihelion is initially located inside the compact star, it will lose part of its kinetic energy due to the interaction with the compact star matter near the perihelion, and will be finally trapped inside the compact star after several passages (Abramowicz et al. 2018). When the PBH moves completely inside the compact star, the accretion process can efficiently circularize the orbit (Macedo et al. 2013), especially for the extreme mass ratio cases considered in this work (m PBH /M CS ∼ 10 −6 − 10 −3 ). Therefore, our previous calculations of the dynamics and GW emissions are still valid for the later stage of a PBH captured by a compact star. For the earlier stage when the trapped PBH still has a non-negligible eccentricity, Equation 8 can be used to calculate the motion of the PBH, but the dynamical friction term and the resulting GW waveform will be very different.

EVENT RATE
If PBHs account for a fraction X of the dark matter, the detectable GW event rate of compact stars capturing PBHs can be calculated by (Capela et al. 2013;Fuller et al. 2017) where ρ DM is the density of dark matter, N is the number of compact stars within GW detectors' horizon, R S is the Schwarzschild radius of compact stars, and E loss ∝ m 2 PBH is the energy loss during the PBH-compact star interaction. Here we have assumed a Maxwellian distribution for the PBH velocities, with a dispersion ofv. For the PBH mass, we take a typical value of m PBH = 10 −5 M in our calculations below.
There are 281 pulsars within the 1 kpc range of the Sun 5 (Manchester et al. 2005). Considering that pulsars have an average lifetime of ∼ 10 7 years (compared to the lifetime of ∼ 10 10 years of the Milky Way) and beaming factor of ∼ 3.5 (Lyne & Graham-Smith 2012), the true number of compact stars within 1 kpc of the Sun can be estimated as ∼ 10 6 . Since the GW amplitude of an inspiraling 10 −5 M PBH at 1 kpc away is about one magnitude higher than the sensitivity curve of the Advanced LIGO (see Figure 7 (a)), the Advanced LIGO's horizon is therefore ∼ 10 kpc, which nearly includes the whole Milky Way. As a result, there are N ∼ 10 9 compact stars in the detection horizon. This is a logical number, which actually coincides with the total number of compact stars in the Milky Way. For PBHs with m PBH ∼ 10 −5 M , one has 3E loss m PBHv 2 (Capela et al. 2013), thus the exponential term in Equation 13 can be neglected. Adopting a radius of R X ∼ 12 km for a typical 1.4 M compact star (see Section 2.3) and taking the PBH velocity dispersion asv ∼ 105 km s −1 (Fuller et al. 2017), the detectable event rate within a 10 kpc horizon (by the Advanced LIGO) is where ρ DM = 0.4 GeV cm −3 is the local dark matter density (de Salas & Widmark 2021). This rate seems to be very small. However, note that in Equation 14 we have assumed that both the PBHs and the compact stars distribute homogeneously in a spherical volume with a radius of 10 kpc. In a realistic case, both the dark matter and the compact stars concentrate at the Galactic center, and the high densities will markedly increase the event rate. The Galactic center is ∼ 8 kpc away from us. In fact, the dark matter density at the Galactic center might be as high as ρ DM = 10 6 GeV cm −3 (Bertone & Merritt 2005), and a population of up to N ∼ 10 8 compact stars may reside there (Cordes & Lazio 1997). In such a dense environment, multi-body interactions can further increase the capture rate by a factor of ∼ 3.5 (Brayeur & Tinyakov 2012). Moreover, note that Equation 13 only takes the dynamical friction into account (Capela et al. 2013), while the energy loss through GW emission can additionally enhance the capture rate by up to a factor of 10 for PBHs with m PBH 10 −6 M (Génolini et al. 2020). Considering all the above factors and assuming PBHs account for a portion of 0.01 -0.1 of the dark matter at the Galactic center (Niikura et al. 2019), the event rate solely from the Galactic center can be estimated as F ∼ 10 −5 -10 −4 year −1 . This rate is rather low, but still nonzero. It is also consistent with the nondetection of such events by the LIGO-Virgo Collaboration up until now.
Future GW detectors like the Einstein Telescope and the Cosmic Explorer have sensitivities nearly two magnitudes better than the advanced LIGO (see Figure 7), thus may have a horizon of ∼ 1 Mpc. About 10 -100 galaxies may reside in this range, leading to a total detectable event rate of F ∼ 10 −4 -10 −2 year −1 . Note that there are still other factors that can further enhance the event rate. For example, globular clusters are expected to have a relatively high dark matter density as well as a low velocity dispersion. Thus they are also interesting places for PBH captures (Das et al. 2022). Moreover, the tidal effect (Génolini et al. 2020) and possible underestimation of ρ DM and N at the Galactic center can also increase the event rate, which are however beyond the scope of this study. Anyway, the upcoming next-generation GW detectors may hopefully succeed in detecting such events.

SUMMARY AND DISCUSSION
The process of a planetary-mass PBH inspiraling inside a compact star is investigated in detail. Such a process is of great interest for constraining PBH's composition fraction of dark matter. The effects of dynamical friction, accretion, and GW emission on the motion of the PBH are taken into account. Two kinds of compact stars are considered, i.e., SSs and normal NSs. It is found that the resulting GW signals show significant difference between the two cases. Encouragingly, the current Advanced LIGO detector can detect a 10 −5 M PBH case at 10 kpc away from us, and detect a Jovian-mass PBH case up to megaparsecs. The planned Einstein Telescope and Cosmic Explorer can detect 10 −5 M PBHs at ∼ 1 Mpc, and detect Jovian-mass PBHs at several hundred megaparsecs. Note that the GW frequencies are near the high end of the sensitivity curve for most ground-based interferometers. Future ad hoc detectors working at kilohertz frequencies like Neutron Star Extreme Matter Observatory (Ackley et al. 2020) will be powerful tools for probing the EoSs of dense matter.
Our dynamical equation (Equation 8) is invalid at the final stage when the whole compact star is to be swallowed by the PBH to form a stellar mass black hole. In our study, we did not calculate the GWs emitted at this final stage. However, since the PBH has moved to the central region of the compact star at the end of our calculation, the accretion by the black hole should be highly spherical and the GWs are expected to be much weaker according to Birkhoff's theorem. In our modeling, we have omitted the spinning of the compact star for simplicity. But note that the rotation of the compact star would not affect the results significantly, because even a millisecond pulsar rotates much slower than its Keplerian velocity (Génolini et al. 2020). Another approximation is that a constant λ is adopted during our calculations. To overcome this weakness, the GW signals should be calculated through full general relativistic hydrodynamic simulations, which is beyond the scope of the current study.
Before the PBH contacts the NS/SS surface and begins its journey inside the compact star, it orbits around its host as a very close-in object, which will also emit strong GWs. Such GWs may be recognized as originating from a "planet"-compact star system. Geng et al. (2015) and Kuerban et al. (2020) have suggested that GWs from such "planet"-compact star systems could be efficiently used to identify strange quark planets, because a normal matter planet will be tidally disrupted by the compact star when it is still further away so that no GWs are available (Geng et al. 2015;Kuerban et al. 2020). Here, we would like to further remind that the possibility that the "planet" is actually a PBH should be further repelled before finally identifying it as a strange quark planet. The GWs emitted by the PBH inspiraling inside its host compact star, as calculated by us in this study, can help us with this task. The merging of a strange quark planet with an SS will only produce some kind of ringing-down in GWs, while a PBH will tunnel through its host and produce complicated GW patterns as demonstrated in this study.
PBHs colliding with compact stars are also hypothetically associated with many other interesting phenomena, including the formation of solar mass black holes (Takhistov et al. 2021) and fast radio bursts (Abramowicz et al. 2018;Kainulainen et al. 2021). Future observations of the GWs will help constrain these hypotheses.