Primordial black holes (PBHs), first proposed by [1,2], may have formed in the early Universe, and may contribute to or even make up most of its dark-matter content (see also [3]). While observational constraints on PBHs limit their possible contribution to the dark matter in some mass ranges, they remain viable candidates in other mass windows, including between about 10 -16 M ⊙ and 10 -10 M ⊙ as well as around 10 -6 M ⊙ (see, e.g., [4][5][6] for reviews and details).

If PBHs exist, some of them are likely to interact with stars and other celestial objects. Such interactions have been invoked as possible origins of several astrophysical phenomena, including the 1908 Tunguska event in Siberia [7] (but see [8]), neutron star (NS) implosions and "quiet supernovae" [9,10], fast radio bursts [9,11,12], the formation of low-mass stellar black holes [13][14][15][16], microquasars [17], and the origin of supermassive black holes (e.g., [18]), possibly via the formation of PBH clusters [19,20]. Gravitational-wave (GW) signatures of PBHs have been surveyed recently in [21], and the prospect of detecting PBHs using solar-system ephemerides has been discussed in [22,23] and references therein.

A collision with a star results in the PBH being gravitationally bound if it loses a sufficient amount of energy in the encounter, which is most likely to happen in collisions with NSs (see, e.g., [11,[24][25][26]). The PBH may still emerge from the star, but can no longer escape to infinity. Losing more energy in subsequent passages, the PBH at some point remains completely inside the star, settles down toward its center, accretes stellar material, and ultimately induces the dynamical collapse of the host star (see [27][28][29] for numerical simulations). While the expected event rates are small (see, e.g., [24,25,[30][31][32] as well as Sec. I in the Supplemental Material [33] for estimates) they depend strongly on a number of assumptions and may be more favorable in special environments, e.g., globular clusters and galactic centers. Small black holes may also form inside neutron stars from the collapse of other dark-matter particles (e.g., [10,27,34,35]), or be captured by neutron stars by other processes (e.g., [18,26,31,36]).

While the PBH spirals toward the center of the NS it emits gravitational radiation that-at least in principlemay be observable by next-generation GW detectors, and that would reveal information about the stellar structure [31]. The authors of [32,37], for example, examined this scenario assuming circular orbits. Since the PBH typically enters the host star on a noncircular orbit, and since the retarding forces inside the star may not circularize the orbit (see, e.g., [38,39]), the PBH's orbit is likely to remain eccentric (see also [40] for a numerical demonstration). In this Letter we discuss a qualitatively new feature of such noncircular orbits, namely continuous, quasiperiodic GW beats. These beats are caused by a precession of the PBH's orbit inside the star, the GW frequency for which is superimposed on the higher frequency arising from a single orbit. The resulting GW envelopes for the two GW polarizations are exactly out of phase, so that the GW signal alternates between being dominated by one or the other polarization. In the stellar interior, both Newtonian and relativistic effects contribute to this precession, but we find that the latter dominate in NSs. As we demonstrate both numerically and analytically, the rate of the precession, and hence the beat frequency, depends rather strongly on the NS structure, so that a future observation of such a GW beat could provide strong constraints on the nuclear equation of state (EOS), let alone confirmation of the capture by a NS of a smaller and lower-mass intruder.

As a numerical demonstration we show in Fig. 1 characteristic orbits of PBHs inside NSs governed by three different EOSs varying in stiffness, together with their associated GW signals. We adopt a simple particle testmass approximation to describe the PBH moving on geodesics in the gravitational field of a relativistic star. As dynamical friction and accretion drag forces are small perturbations that operate on secular timescales much longer than orbital times (e.g., [26,30]), we can probe the precession by neglecting these forces and examining a few orbits via geodesics. We assume the stars to be governed by a polytropic EOS

Here P is the pressure, ρ 0 the rest-mass density, K a constant, and the adiabatic exponent Γ ¼ 1 þ 1=n may be expressed in terms of the polytropic index n. We construct the stellar models by solving the Oppenheimer-Volkoff (OV) equations, adjusting the central density so that the stellar compaction is always given by

where M Ã is the star's total gravitational mass and R Ã is areal radius (see Sec. II in the Supplemental Material [33] for details). In Fig. 1 we show examples for a star with constant total mass-energy density ρ (corresponding to Γ ¼ ∞), as well as Γ ¼ 3 and Γ ¼ 2 polytropes, which serve as examples of both highly and moderately stiff candidates for the NS EOS.

We then solve the relativistic geodesic equations in order to track the PBH's orbit. As discussed in Sec. III of the Supplemental Material [33], we choose to solve these equations in terms of the isotropic metric and radius r rather than the interior Schwarzschild metric and areal radius R, although the conversion is straightforward. We always start orbits with vanishing radial speed u r ¼ 0 at an initial (areal) radius Rð0Þ ¼ R frac R Ã and with angular momentum l ¼ l frac l circ , where l circ is the angular momentum corresponding to a circular orbit at radius Rð0Þ. For the examples shown in Fig. 1 we used R frac ¼ 0.6 and l frac ¼ 0.1. All orbits are confined to a plane, which we arbitrarily take to be the xy plane. We also evaluate the leading-order GW signals h þ and h × along the z axis using the quadrupole formalism [42][43][44].

As can be seen in the left column of Fig. 1, the rate at which the PBH's orbit precesses depends strongly on the structure of the host star, and hence its EOS. For all three examples we show the orbits for a time span Δt ¼ 2000M Ã , which corresponds to about 14 ms for M Ã ¼ 1.4M ⊙ . All orbits start out on the positive x axis. For the constantdensity star in the top row, the orbit has rotated by just over 45°during this time, while for Γ ¼ 3 it has rotated by a little less than 180°, and a little over 360°for Γ ¼ 2.

If the GW signal is dominated by one polarization initially (in our case h þ ), then it will be dominated by the other polarization after the orbit has rotated by 45°, and will return to the original polarization after a rotation through 90°. The resulting GW beats, and their dependence on the NS structure, can be seen in the right panels of Fig. 1. Specifically, we observe that the signal has shifted from being dominated by h þ to being dominated by h × for the constant density star, while for Γ ¼ 3 it has gone from h þ to h × and then back to h þ almost twice, and just over four times for Γ ¼ 2.

The above behavior can be understood in part in the context of Bertrand's theorem ( [45], see also [46]), which states that for potentials VðRÞ ¼ V 0 þ kR m , where V 0 and k are constants, only the exponents m ¼ -1 and m ¼ 2 will always result in closed orbits.

The former exponent, m ¼ -1, corresponds to a Newtonian point-mass potential, which is the leading order term in the potential in the exterior of a star. In general, deviations from this exterior potential result both from relativistic corrections as well as several Newtonian effects, including tidal and rotational deformations of the star (which are not present for our static spherical stars). These deviations result in precession of the orbit, including the well-known relativistic perihelion advance of Mercury (see Sec. IV.A of the Supplemental Material [33]).

The latter exponent, a harmonic-oscillator potential with m ¼ 2, is realized in the interior of a Newtonian constant-density star, for which MðRÞ ¼ 4πρR 3 =3 and hence VðRÞ ¼ V 0 þ 2πGρR 2 =3. In this case, deviations result from both density nonuniformities and relativistic corrections. While we do not expect any precession, either in the interior or exterior, for a homogeneous spherical star in Newtonian gravitation, relativistic effects will cause precession for such a star both in the exterior and interior.

Even for an inhomogeneous star, the stellar core becomes increasingly homogeneous as R → 0. In Newtonian gravitation, the PBH's orbit therefore starts with a nearly closed orbit with m ¼ -1 far outside the NS, and ends with a tighter, nearly closed orbit with m ¼ 2 well inside the NS-quite remarkably realizing both cases of Bertrand's theorem as extreme limits. Moreover, since the approximately homogeneous core extends to larger radii for stiffer EOSs than for softer EOSs, we expect that, for a given orbital radius, the precession will be slower for a stiffer EOS than for a softer EOS in Newtonian theory. This precession rate variation is also found in general relativity, as revealed in the numerical examples of Fig. 1.

We can gain analytical insight into the above effects by considering small perturbations of circular orbits. Geodesics in static and spherically symmetric spacetimes possess two conserved constants of motion, namely the energy per unit mass e ¼ -u t and the angular momentum per unit mass l ¼ u φ . Using the normalization of the fourvelocity u a , g ab u a u b ¼ -1, a first integral of the equations of motion can be written in the form

where the constant E ≡ ðe 2 -1Þ=2 plays the role of the kinetic energy at infinity for velocities v ∞ ≪ 1, and where we split the effective potential V eff ðRÞ into the two terms

In many cases (e.g., a Newtonian point-mass) VðRÞ is independent of e and l, but we now allow this term to depend on these two constants, which is the case for the orbits in general relativity considered here. In the following we assume the orbit to be in the equatorial plane, so that θ ¼ π=2, and we provide details of how VðRÞ can be determined in Sec. IV of the Supplemental Material [33]. For a stable circular orbit at (areal) radius R 0 the effective potential V eff ðRÞ must take a minimum there, so that we have V 0 ðR 0 Þ≡ðdV=dRÞ R 0 ¼l 2 =R 3 0 . Since l¼g φφ u φ ¼R 2 dφ=dτ, the proper time τ φ needed to complete one orbit, i.e., to advance from an angle φ 0 to φ 0 þ 2π, is given by

In the vicinity of a stable orbit we may approximate the effective potential as a parabola

, where k ≡ V 00 eff ðR 0 Þ > 0 and η ≡ R -R 0 . Inserting these into (2) and taking a derivative with respect to proper time τ results in the harmonic-oscillator equation η þ kη ¼ 0 for η, where the double dot denotes a second derivative with respect to τ. Accordingly, the proper time τ R needed to travel from the orbit's pericenter to the apcenter and back to the pericenter is given by

and the ratio between the two times τ R and τ φ is

For a Newtonian point-mass potential VðRÞ ¼ pR -1 , where p is a constant, we have τ R ¼ τ φ , as expected for Kepler orbits. For a harmonic-oscillator potential VðRÞ ¼ pR 2 we have τ R ¼ τ φ =2, so that the orbit, which is centered at and symmetric about the origin, features two pericenters in each revolution. According to Bertrand's theorem, these two cases are the only potentials that lead to closed orbits. While perturbations of the point-mass potential in the stellar exterior are familiar-yielding, for example, the relativistic perihelion advance of Mercury-we now focus on the stellar interior. Specifically, we show in Sec. IV.B of the Supplemental Material [33] that, in the vicinity of the center, the potential VðRÞ can be written in the form

where V 0 , p, and q are constants. Inserting ( 7) into ( 6) and assuming qR 2 0 ≪ p we find, to leading order,

As the black hole advances from one pericenter to the next, its positional angle φ therefore advances by

beyond the angle π that would result in a closed orbit. We refer the reader to Sec. IVof the Supplemental Material [33] where the constants p and q are evaluated in general relativity for nearly circular orbits. Here we present a more transparent Newtonian treatment in order to illustrate the key ingredients. In the vicinity of the stellar center we may approximate the density as ρðRÞ ≃ ρ c þ ρ ð2Þ R 2 =2, where ρ ðnÞ c ≡ ðd n ρ=dR n Þ R¼0 . Integrating once we find the enclosed mass MðRÞ, and integrating again we obtain the potential

Comparing with (7) we identify both p and q and compute

Evidently, the Newtonian pericenter advance is related to the degree of inhomogeneity, consistent with our discussion above.

We may evaluate the term ρ ð2Þ c in (11) using the Newtonian equations of hydrostatic equilibrium. For a polytropic EOS (1) we find

where a ¼ ðΓP=ρÞ 1=2 is the (Newtonian) speed of sound. Finally we may use the central condensation δ ≡ ρ c =ρ, where ρ ¼ 3M Ã =ð4πR 3 Ã Þ is the average density, to rewrite the Newtonian pericenter advance as

For Newtonian polytropes δ depends on Γ only (see Table I for specific values). For smaller Γ, δ is larger and a 2 c smaller [for a given compaction GM Ã =ðc 2 R Ã Þ]; we therefore see that the pericenter advance is larger for a softer EOSs (for a given value of R 0 =R Ã ).

While the above analysis captures the leading-order Newtonian terms, we have found that the pericenter advance in the NSs considered here is dominated by relativistic terms. However, the pericenter advance's dependence on the EOS's stiffness is similar to that observed from the above Newtonian analysis even in the context of general relativity-namely, a softer EOS will lead to a more rapid precession of the orbit, and therefore to higher-frequency GW beats. This can be observed in Fig. 1 as well as in the Table I, where we list pericenter adnvances Δφ for nearly circular orbits (l frac ¼ 0.99) close to the center (R frac ¼ 0.05) for a range of polytropic indices. We compare numerical results from the integration of the geodesic equation with analytical results from the perturbation of nearly circular orbits and find excellent agreement.

The range of polytropic exponents Γ listed in Table I roughly covers values adopted in piecewise-polytropic approximations for candidate nuclear EOSs (see Table III in [47]; note in particular the larger range of values for Γ 3 , which governs the high-density core). The resulting values of Δφ show significant variation, suggesting that a potential of the resulting GW beats would provide a sensitive probe of the EOS.

the pericenter advance Δϕ we may also compute the precession frequency. Since we defined Δφ as the (excess) advance from one pericenter to the next, and since, to leading order, orbits in the stellar interior feature two pericenters per orbit, the angular precession frequency measured locally is given by Ω prec ¼ 2Δφ=τ φ , where τ φ is the orbital (proper) period. Related to the precession frequency is the (proper) precession period τ prec that it takes either of the two GW polarization amplitudes to go through one complete cycle, i.e., for the pericenter to advance by an angle π,foot_0

The number of orbits completed during a precession period τ prec is therefore N orbit ¼ τ prec =τ φ ¼ π=ð2ΔφÞ. Since, during one revolution, the GW signal completes two cycles, the number of such GW cycles completed as either GW polarization goes through a full beat cycle associated with the GW envelope is given by

For the orbits in Fig. 1, for example, we found Δϕ ¼ 0.0255, 0.0597, and 0.127 for Γ ¼ ∞, 3, and 2, respectively, resulting in N GW ¼ 123, 52.6, and 24.8. In Table I we also provide data for N GW , τ φ , and τ prec for the nearly circular orbits considered there.

The above periods are proper times as measured by an observer comoving with the PBH. For the nearly circular orbits in Table I we may simply multiply these periods with u t in order to obtain the corresponding coordinate time periods. In Table I we list the resulting GW frequency associated with a single orbit, f GW orb ≃ 2=ðu t τ φ Þ, and the GW beat frequency f GW beat ≃ 1=ðu t τ prec Þ, both as measured by a distant observer.

To summarize, we discuss quasiperiodic GW beats as a qualitatively new feature of continuous GW signals emitted by PBHs captured inside NSs. The beats are due to orbital precession, which is caused both by relativistic effects and density nonuniformity. Adopting a polytropic EOS we demonstrate both numerically and analytically that the beat frequency depends quite strongly on the structure of the NS and hence the stiffness of the EOS. For the NSs considered here the precession rate and beat frequency are largely due to relativistic gravitation, so that a Newtonian treatment would significantly underestimate the effect. If such beats were to be observed by next-generation GW detectors, e.g., the Einstein Telescope [48], the Cosmic Explorer [49], or the Neutron Star Extreme Matter Observatory (NEMO) [50], they would therefore provide valuable constraints on the nuclear EOS.

Clearly, the beat frequency also depends on the radius and eccentricity of the PBH's orbit, which would have to be found independently. The latter is related to the relative maximum and minimum amplitudes in each one of the GW polarizations, and it may be possible to determine the radius from the prior inspiral signal. Knowing these orbital parameters, as well as the host star's compaction, an observed beat frequency could then be compared with those found for orbits inside general relativistic stellar models constructed for TABLE I. Numerical and analytical data for polytropic stellar models and pericenter advances for nearly circular orbits close to the center of a stellar host with compaction GM Ã =ðc 2 R Ã Þ ¼ 1=6. We list, for different values of Γ ¼ 1 þ 1=n, Newtonian values of the central condensation δ ¼ ρ c =ρ (which depends on the Γ alone) and the Newtonian estimate Δφ=ðR 0 =R Ã Þ 2 , adopting numerical solutions to the Lane-Emden equation, together with Δφ for R 0 =R Ã ¼ 0.05. For the relativistic data we computed δ from solutions to the OV equations. Adopting R 0 =R Ã ¼ 0.05 again we computed Δφ num from numerical solutions to the geodesic equations, using l frac ¼ 0.99 for nearly circular orbits, and Δφ ana analytically as presented in the Supplemental Material [33]. We also list the number of GW wave cycles N GW completed during a beat cycle [see Eq. ( 15)], and, assuming a host star with mass M Ã ¼ 1.4M ⊙ , the orbital time τ φ as well as the precession time τ prec . In the last two columns we provide the corresponding GW frequencies f GW orb ≃ 2ðu t τ φ Þ -1 and f GW beat ≃ ðu t τ prec Þ -1 as measured by a distant observer.

Newton GR

2 Δφ δ Δφ num Δφ ana N GW τ φ [ms] τ prec [ms] f GW orb [kHz] f GW beat [Hz] 1.75 1.33 4.89 1.26 0.00315 7.46 0.0161 0.0162 194 0.194 18.9 5.56 28.7 2.0 1.0 3.29 0.775 0.00193 3.98 0.00841 0.00844 372 0.273 51.0 4.48 12.1 2.25 0.8 2.60 0.556 0.00139 2.94 0.00608 0.00612 513 0.321 82.4 4.01 7.81 2.5 0.67 2.23 0.432 0.00108 2.43 0.00494 0.00498 630 0.354 111 3.73 5.91 2.75 0.57 2.00 0.352 0.00088 2.14 0.00427 0.00430 731 0.379 138 3.54 4.85 3.0 0.5 1.84 0.298 0.00075 1.94 0.00382 0.00385 816 0.398 163 3.41 4.18 3.25 0.44 1.72 0.257 0.00064 1.80 0.00351 0.00353 890 0.414 184 3.31 3.72 ∞ 0 1 0 0 1 0.00164 0.00165 1904 0.553 527 2.62 1.38