Constraining polymerized black holes with quasi-circular extreme mass-ratio inspirals

The Event Horizon Telescope made a significant breakthrough by capturing the first image of the supermassive black hole in the M87 galaxy in 2019 [1], followed by the first image of the supermassive black hole at the center of the Milky Way announced in 2022 [2]. These images provide direct visual evidence of supermassive objects at the centers of galaxies. In each galaxy, there are many stellar-mass objects orbiting the central supermassive black hole. Over time, these stellar-mass objects slowly inspiral toward the central supermassive black hole due to gravitational radiation, constituting extreme mass-ratio inspiral (EMRI) systems [3]. The gravitational waves emitted by EMRIs are expected to be detected by future space-based gravitational wave detectors, such as Laser Interferometer Space Antenna (LISA) [4], Taiji [5], TianQin [6], and DECi-hertz Gravitational-wave Observatory (DECIGO) [7].

For each EMRI, smaller objects undergo numerous orbits around the central supermassive black hole before ultimately plunging into it [3]. The orbital dynamics of the smaller objects are mainly determined by the properties of the central supermassive black hole and gravitational radiation. The gravitational waveforms emitted by an EMRI are closely related to the trajectories of the smaller objects, serving as precise mappings of the spacetime geometry around the supermassive black hole. Thus, gravitational waves from EMRIs are invaluable probes for studying the characteristics of supermassive black holes and testing gravitational theories [8−13].

Semi-relativistic methods (such as analytical and numerical kludge (NK) methods [14−18]) are effective for studying the evolutions of EMRIs. The NK method combines the exact relativistic orbital trajectory of a smaller object with an approximate expression for the gravitational radiation [15]. Typically, the orbit of the smaller object is elliptical. However, as gravitational radiation leads to a gradual loss of energy and orbital angular momentum, the eccentricity of the orbit gradually decreases [19]. Gravitational radiation from a particle in a circular orbit around a black hole was studied by the authors of Refs. [20−22]. To describe the orbital evolution of the smaller object, the quasi-circular orbit approximation is often employed [23]. In this approximation, the smaller object is assumed to remain in a circular orbit at each moment, with the radius of the orbit slowly decreasing over time due to gravitational radiation. The authors of Refs. [23−26] employed the quasi-circular orbit approximation to study the gravitational waves from a compact object orbiting a supermassive Kerr black hole. The quasi-circular orbit approximation was also used to study LISA's capability to detect new fundamental fields [27−29]. More recently, quasi-circular EMRIs have been used to investigate the properties of black holes in bumblebee gravity [30].

The existence of black holes is a key prediction of general relativity. However, the singularities predicted within black holes pose a significant theoretical challenge. As a potential solution, loop quantum gravity has proposed several black hole models that eliminate these singularities [31−33]. The authors of Refs. [34−37] explored the EMRIs around several black holes in loop quantum gravity. However, the authors of Refs. [34, 35] examined EMRIs over only a few orbital cycles, neglecting the effects of gravitational radiation on the orbital evolution of the smaller object. Through the semi-classical polymerization of the area within the interior of spherically symmetric black hole spacetimes, a quantum-corrected black hole spacetime has been proposed in Refs. [38, 39]. Various properties of this polymerized black hole have been studied, including Hawking radiation [40], quasinormal modes [41, 42], and strong gravitational lensing and shadow [42, 43].

In this study, we focus on gravitational waveforms generated by a stellar-mass object inspiraling into a supermassive polymerized black hole along the quasi-circular orbits. Treating the smaller object as a massive test particle, we derive its equations of motion and obtain the corresponding radial effective potential. We investigate the impact of the quantum corrections on stable circular orbits (SCOs), including the innermost SCOs (ISCOs), around the polymerized black hole. Then, we employ a numerical algorithm to model the quasi-circular orbit inspirals. We explore how the quantum corrections affect the evolution of the quasi-circular inspirals. Considering a stellar-mass object moving along the quasi-circular orbits around a supermassive polymerized black hole, we analyze the gravitational waveforms from this EMRI system by using the NK scheme [15] and discuss how the quantum corrections influence these waveforms. We further account for the Doppler shift induced by LISA's motion [14] and calculate the resulting gravitational wave response signal in LISA. To assess the observability of the quantum parameter's effects on gravitational waves, we calculate the mismatch between waveforms for different values of the quantum parameter and those of a Schwarzschild black hole. Finally, we investigate the critical value of the quantum parameter, below which the gravitational waves from quasi-circular EMRIs of the polymerized black hole become indistinguishable from those of a Schwarzschild black hole.

The remainder of this paper is organized as follows. In Sec. II, we analyze the geodesic of a massive test particle on the equatorial orbit around the polymerized black hole, derive the radial effective potential, and investigate the properties of the general SCOs and ISCOs. In Sec. III, we adopt a numerical algorithm to model the quasi-circular orbit inspirals. In Sec. IV, we study the gravitational waveforms of the test object along quasi-circular orbits around the supermassive polymerized black hole, investigate their dephasing, and explore their distinguishability. Finally, the conclusions of this study are discussed in Sec. V. Throughout the paper, we use the geometrized unit system with $ G = c = 1 $.

Applying the polymerization of the area operator in the interior of a black hole, a new black hole metric in loop quantum gravity has been proposed in Refs. [38, 39]. This polymerized black hole offers potential insights into the quantum dynamics inside black holes and could help in understanding the connection between quantum gravity and general relativity. The line element of the polymerized black hole is

where M is the Arnowitt–Deser–Misner mass, and k is a quantum-corrected parameter that determines the minimal radius of the spacetime. The radius of the event horizon for the polymerized black hole is $ \sqrt{4 M^2 + k^2} $. Equation (1) returns the line element of a Schwarzschild black hole when $ k = 0 $. For convenience, we use the dimensionless parameter $ \hat{k} = k/M $ instead of k in this study. Through the observations of M 87* and Sgr A* black hole shadows from the Event Horizon Telescope, Ref. [43] obtained constraints on k, with the strongest being $ \hat{k} \leq 0.36 $.

It is important to note that the true value of the parameter $ \hat{k} $ is expected to be extremely small. In an EMRI system, a smaller object orbits a central supermassive black hole over many years, and very long-duration gravitational wave signals are generated. Gravitational waves from EMRIs can map the background spacetime geometry with high precision. Although the value of the parameter $ \hat{k} $ is small, its effects may still become detectable over extended observation periods. If future space-based gravitational wave detectors detect gravitational waves from EMRIs, we could use these gravitational waves to derive precise constraints on the black hole metric parameters. In this study, we use a significantly larger value for $ \hat{k} $ so that we can better highlight its impact on the gravitational wave signal and ensure that its influence significantly exceeds numerical errors.

To explore the properties of the polymerized black hole, we study an EMRI system consisting of a stellar-mass object orbiting a central supermassive polymerized black hole. The gravitational waves from this EMRI system provide a highly accurate probe of the spacetime geometry of the polymerized black hole. The gravitational waveforms from the EMRI system are closely related to the trajectory of the smaller object. When the influence of the smaller object on spacetime is ignored, we can treat the smaller object as a massive test particle and study its orbit with a geodesic at every moment. We use the Lagrangian method [44] to derive the equations of motion for the test particle moving on the equatorial plane $ (\theta = \pi/2) $ around the polymerized black hole in Appendix A. Equation (29) can be rewritten as

where the radial effective potential is

We plot the radial effective potential (3) in Fig. 1 to investigate its properties. As shown in Fig 1(a), the radial effective potential will possess an extremum with increasing orbital angular momentum. The minimum of the radial effective potential corresponds to the stable circular orbits, and the maximum corresponds to the unstable ones. Figure 1(b) shows that the maximum of the radial effective potential decreases with the parameter $ \hat{k} $.

Figure 1.  (color online) (a) Radial effective potentials of the test particle with different values of orbital angular momentum around the polymerized black hole with $ \hat{k} =1 $. (b) Radial effective potentials of the test particle with fixed orbital angular momentum $ L/M = 4.2 $ around the polymerized black hole with different values of the parameter $ \hat{k} $.

In the late stage of an EMRI, the eccentricity of the smaller object's orbit is nearly zero [19]. Thus, the assumption of a quasi-circular orbit is a reasonable approximation in modeling the orbital evolution of EMRIs [24−26]. In this study, we assume that the smaller object is in an SCO around the supermassive polymerized black hole at every moment. The properties of a test particle moving along an SCO around the polymerized black hole should satisfy [45]

Evidently, the parameter $ \hat{k} $ affects the properties of the SCO in Eq. (4). We numerically solve Eq. (4) and plot the relations between the radius of the SCO and parameter $ \hat{k} $ and between the energy of the SCO and parameter $ \hat{k} $ in Fig. 2. One can see that, with a fixed orbital angular momentum, both the radius and energy of the test particle along the SCO decrease with the parameter $ \hat{k} $. The ISCO is the SCO with the minimal radius. In gravitational wave astrophysics, the ISCO marks the boundary between the inspiral and merger phases in a two-body system. The properties of a test particle moving along an ISCO around the polymerized black hole should satisfy the following conditions:

Figure 2.  (color online) (a) Radius of the SCO around the polymerized black hole as a function of parameter $ \hat{k} $. (b) Energy of the test particle along the SCO around the polymerized black hole as a function of parameter $ \hat{k} $.

We numerically solve Eq. (5) and plot the relations of the radius, orbital angular momentum, and energy of the ISCO with the parameter $ \hat{k} $ in Fig. 3. It can be seen that the radius, orbital angular momentum, and energy of the ISCO all increase with the parameter $ \hat{k}. $

Figure 3.  (color online) Radius (a), orbital angular momentum (b), and energy (c) of the ISCO around the polymerized black hole as functions of the parameter $ \hat{k} $.

In this section, we investigate the orbital evolution of a stellar-mass object inspiraling into a central supermassive polymerized black hole with the quasi-circular orbit approximation and study the corresponding gravitational waveforms. For the smaller object with mass m moving along a trajectory $ Z^i(t) $, we treat the Boyer-Lindquist coordinates as fictitious spherical polar coordinates and then project the object's trajectory onto the following Cartesian coordinates [46]:

Then, the symmetric and trace-free (STF) mass quadrupole of the smaller object can be defined as [47]

where

is the $ tt $-component of the stress-energy tensor for the object. Gravitational radiation carries away the smaller object's energy and orbital angular momentum. The fluxes of the energy and orbital angular momentum are given by [48]

where $ \epsilon ^{ikl} $ is the three-dimension Levi-Civita symbol, and the angled brackets denote the average of a physical quantity over several orbital periods. For a smaller object moving along a circular orbit, Eq. (9) can be simplified to

where R is the radius of the circular orbit, and Ω is the angular velocity of the object.

We assume that the stellar-mass object is in an SCO around the central supermassive polymerized black hole at every moment. Due to gravitational radiation, the energy and orbital angular momentum of the smaller object decrease, leading to a simultaneous reduction in the radius of the SCO in which the smaller object resides. Therefore, the complete orbit of the smaller object is quasi-circular. The algorithm we used to model the orbital evolution of a stellar-mass object around a supermassive black hole under gravitational radiation proceeds as follows: i) determine the small object's mass m, the supermassive black hole's mass M, and a time step $ \Delta t $ for evolution; ii) set the initial position $ (r_0, \phi_0) $ of the object; iii) calculate the initial energy $ E_0 $, orbital angular momentum $ L_0 $, and angular velocity $ \Omega_0 $ of the object; iv) calculate the mass quadrupole moment $ I_{ij} $, energy flux $ \mathrm{d}E/\mathrm{d}t $, and angular momentum flux $ \mathrm{d}L/\mathrm{d}t $ at the object's position; v) update the object's orbital angular momentum; vi) compute the new position and energy using the updated orbital angular momentum; and vii) repeat steps iv) to vi) until the object's orbit shrinks to or below the ISCO.

The algorithm described above systematically tracks the gradual decay of the orbit due to gravitational radiation, providing a detailed picture of the orbital evolution of a stellar-mass object around a supermassive black hole. Consider an EMRI system consisting of a test object with $ m = 10 M_\odot $ and a central supermassive polymerized black hole with $ M= 10^6 M_\odot $. By setting the initial position $ (r_0, \phi_0) $ of the smaller object as $ (10M,\; \pi/2) $ and assigning different values to the quantum correction parameter $ \hat{k} $, we use the algorithm described above to investigate the evolution of the smaller object's quasi-circular orbit around the central supermassive polymerized black hole. Taking $ \hat{k} = 0.1 $ as an example, we plot the complete orbit of the smaller object in Fig. 4. We also plot the evolutions of the smaller object's energy E, orbital angular momentum $ L/M $, and radius $ r/M $ of the quasi-circular orbit around the supermassive polymerized black hole with different values of the parameter $ \hat{k} $ in Fig. 5. The results show that the smaller object's orbital evolution takes a shorter time when it orbits around the central supermassive polymerized black hole with a larger parameter $ \hat{k} $. This is consistent with the radius of the ISCO increasing with the parameter $ \hat{k} $. We have numerically confirmed the validity of the orbital evolution algorithm we employ in Appendix B.

Figure 4.  (color online) Complete quasi-circular orbit of a stellar-mass object inspiraling into a central supermassive polymerized black hole with $ \hat{k} = 0.1 $. The initial position $ (r_0, \phi_0) $ of the object is $ (10 M, \; \pi/2) $.

Figure 5.  (color online) Evolutions of the smaller object's energy E, orbital angular momentum $ L/M $, and radius $ r/M $ of the quasi-circular orbit around the supermassive polymerized black hole are shown in (a), (c), and (e) within the complete orbital time and in (b), (d), and (f) within the late stage, respectively.

After obtaining the orbital evolution of a small object inspiraling into the supermassive polymerized black hole, we can investigate the corresponding gravitational waveforms to further explore the properties of the polymerized black hole. For a stellar-mass object moving around a supermassive black hole, to the order of the mass quadrupole moment, the metric perturbations describing the gravitational waves are defined as [47]

where $ D_ \text{L} $ is the luminosity distance from the EMRI system to the detector. To obtain the gravitational-wave polarizations, we should construct a new detector-adapted coordinate system $ (X,\; Y,\; Z) $ [46]. Under this coordinate system, the gravitational-wave polarizations from Eq. (11) are

The components $ h_{\zeta \zeta} $, $ h_{\iota \iota} $, and $ h_{\iota \zeta} $ are given by [15]

where ι is the inclination angle of the orbital plane of the smaller object to the $ X-Y $ plane, and ζ is the longitude of the pericenter measured in the orbital plane.

Throughout the evolution of the smaller object's quasi-circular orbit, we get the smaller object's mass quadrupole moment at each position. Then, we calculate the corresponding gravitational waveforms from Eqs. (12) and (13), with $ D_\text{L} = 2 \text{Gpc} $, $ \iota = \pi/4 $, and $ \zeta = \pi/4 $. We plot the gravitational waveforms in Fig. 6. Figures 6(a) and 6(c) show that the parameter $ \hat{k} $ has almost no effect on the gravitational waveforms during the initial 2000 s. One can find from Figs. 6(b) and 6(d) that the impact of the parameter $ \hat{k} $ becomes apparent after a one-year accumulation and the parameter $ \hat{k} $ causes an advance in the phase of the gravitational waveforms.

Figure 6.  (color online) Gravitational waveforms from a test object with $ m = 10 M_\odot $ along quasi-circular orbits around a supermassive polymerized black hole with $ M = 10^6 M_\odot $ and different values of the parameter $ \hat{k} $. (a) and (c): Gravitational waveforms within the initial 2000 s. (b) and (d): Gravitational waveforms within 2000 s one year later.

In the study of gravitational waves from EMRIs around black holes in modified gravity theories, the impact of modified gravity on the phase of the gravitational waves is more pronounced than that on amplitude. This is because the phase is closely linked to the orbital dynamics and spacetime geometry. Changes in the gravitational potential due to modified gravity can cause shifts in the orbital frequency, which in turn modify the phase evolution over time. While the amplitude of the gravitational waves also carries valuable information, the phase typically encodes the most detailed signatures of the underlying gravitational theory. As a result, analyzing the phase evolution offers a more sensitive probe for detecting the effects of modified gravity in EMRIs.

To further investigate the effect of the parameter $ \hat{k} $ on the gravitational waves' phase, we define the dephasing for the gravitational waves with $ \hat{k} $ as

where $ \Phi_0 $ is the phase of the gravitational waves in the supermassive Schwarzschild black hole case. We calculate the dephasing of the gravitational waves with different values of the parameter $ \hat{k} $ through Eq. (17). The results are shown in Fig. 7. One can see that the parameter $ \hat{k} $ induces a phase advance in the gravitational waveforms, with a larger parameter resulting in a more pronounced phase shift, and this phase advancement accumulates over time.

Figure 7.  (color online) Dephasing of the gravitational waves with different values of the parameter $ \hat{k} $.

To explore the observability of the quantum parameter's effect on gravitational waves, we should consider the response signal of the gravitational waves in a detector. The response signal of gravitational waves in a detector is related to the characteristics of that detector. The strain amplitude in LISA can be described as [49]

where the factor $ \sqrt{3}/2 $ is from the fact that the actual angle between LISA arms is $ 60^\circ $, and $ F^{+, \times}_{ \text{I,II}} $ are the ''antenna pattern'' functions. Here, we use the expressions of the “antenna pattern” functions from Ref. [14]. The “antenna pattern” functions depend on both the direction of the source $ (\theta_{S},\; \phi_{S}) $ and the orbital angular momentum direction $ (\theta_{L} = 0,\; \phi_{L} = 0) $. The motion of LISA introduces a Doppler shift correction to the phase of the gravitational waves. Following Ref. [14], we modify the phase of the gravitational waves as

where $ R_ \text{AU} $ is the astronomical unit, and $ T = 1 $ year.

For two response signals $ h_1(t) $ and $ h_2(t) $, the noise-weighted inner product between them is defined as [50, 51]

where $ \tilde{h}_i(f) $ is the Fourier transformation of $ \tilde{h}_i(t) $, $ \tilde{h}^\ast_i(f) $ is the complex conjugate of $ \tilde{h}_i(t) $, $ f_ \text{min} $ and $ f_ \text{max} $ are the boundaries of the frequency range of both $ \tilde{h}_1(f) $ and $ \tilde{h}_2(f) $, and $ S_n(f) $ is the power spectral density of LISA [52, 53]. The signal-to-noise ratio (SNR) of a response signal h is $ \sqrt{<h|h>} $. To evaluate the degree of similarity between gravitational wave signals $ h_1 $ and $ h_2 $, we introduce a faithfulness function as follows [51]:

The corresponding mismatch function is [51]

For the gravitational wave signals with SNR ρ and from a model with N intrinsic parameters, the critical value of the mismatch is $ N/(2 \rho^2) $ [51]. If the mismatch between two gravitational wave signals is smaller than the critical mismatch, these two signals are indistinguishable by the detector.

For the model of gravitational waveforms in this study, the dimension of the intrinsic parameters $ N=6 $. In fact, there are two strain signals for a gravitational wave detected by LISA, as given by Eq. (18). Thus, the final SNR of a gravitational wave is $ \rho_\text{F} = \sqrt{\rho^2_\text{I} + \rho^2_\text{II}} $. Following Refs. [27, 28], we consider the last one-year evolution of the EMRI system and rescale the luminosity distance $ D_\text{L} $ to make the final SNR $ \rho_\text{F} $ equal to 30. Then, the corresponding threshold value for the final mismatch is $ 0.0033 $. We calculate the mismatch between the gravitational waves for the supermassive polymerized black hole with different values of the parameter $ \hat{k} $ and those for the supermassive Schwarzschild black hole. The results are shown in Fig. 8. One can find that, under $ \rho_\text{F}=30 $, the mismatch between the gravitational waves for the supermassive polymerized black hole with $ \hat{k} < 0.003 $and those for the supermassive Schwarzschild black hole is smaller than the critical mismatch. This means that we can probe the parameter $ \hat{k} $ to $ {\cal{O}}(10^{-3}) $, which is tighter than the constraint $ \hat{k} < 0.36 $ from the Sgr A$ ^\ast $ black hole shadow [43], with a one-year observation of the EMRI system.

Figure 8.  (color online) Mismatch between the gravitational waves with different values of the parameter $ \hat{k} $ and those for the supermassive Schwarzschild black hole.

To ensure the robustness of the constraints on the quantum correction parameter $ \hat{k} $, it is essential to examine potential parameter degeneracies. In particular, because both the black hole mass M and quantum parameter $ \hat{k} $ influence the gravitational waveform, one must verify that their effects are distinguishable through mismatch analysis. If $ \hat{k} $ was strongly degenerated with M, the observed waveform could not uniquely determine the presence of quantum corrections. To address this, we compute the logarithmic mismatch $ \log_{10} {\cal{M}} $ over a grid of $ (\log_{10}M/M_\odot,\; \hat{k}) $ values and present the results in Fig. 9. The color map shows $ \log_{10} {\cal{M}} $, with darker shades corresponding to larger mismatches. The red dashed curve traces the logarithmic critical mismatch value at each black hole mass. As Fig. 9 shows, the location of the logarithmic critical mismatch $ \log_{10} {\cal{M}}_c = -2.482 $ varies non-trivially with M, and the mismatch exhibits a clear and continuous dependence on both parameters. This behavior indicates that the effect of the quantum correction parameter $ \hat{k} $ cannot be fully absorbed by rescaling the black hole mass, thereby demonstrating that $ \hat{k} $ and M are not degenerate. This supports the validity of our approach in constraining $ \hat{k} $ using gravitational wave signals from EMRIs.

Figure 9.  (color online) Contour plot of the logarithmic mismatch $ \log_{10} {\cal{M}} $ as a function of the supermassive black hole mass M and quantum correction parameter $ \hat{k} $. The red dashed line indicates the logarithmic critical mismatch trajectory $ \log_{10} {\cal{M}}_c = -2.482 $ across varying black hole masses.

In this study, we explored the quasi-circular EMRIs around a central supermassive polymerized black hole. By modeling the smaller object as a massive test particle, we derived its equations of motion and the associated radial effective potential. Our analysis revealed that the peak of the radial effective potential decreases with the dimensionless quantum-corrected parameter $ \hat{k} $. We examined the properties of SCOs around the polymerized black hole, finding that both the radius and energy of the particle on these orbits decrease with the parameter $ \hat{k} $, given a fixed orbital angular momentum. Conversely, for a fixed $ \hat{k} $, both the radius and energy of the particle on the SCOs increase with the particle’s orbital angular momentum. Moreover, we investigated the ISCOs around the polymerized black hole and found that the radius, orbital angular momentum, and energy of the ISCOs all increase with the parameter $ \hat{k} $.

Next, we adopted a numerical algorithm to evolve the quasi-circular orbits of a small object inspiraling into a supermassive polymerized black hole. Using the NK scheme, we analyzed the corresponding gravitational waveforms of this quasi-circular EMRI system. We found that the impact of quantum corrections on the gravitational waveforms is initially negligible but becomes apparent after one year. Additionally, we examined the influence of the parameter $ \hat{k} $ on the phase of the gravitational waveforms, discovering that $ \hat{k} $ induces a phase advance that accumulates over time.

Finally, we considered the Doppler shift caused by LISA’s motion and modified the phase of the gravitational waves accordingly. We then calculated the gravitational wave response signal detected by LISA. To assess the observability of the quantum parameter’s impact on gravitational waves, we computed the mismatch between waveforms generated with different values of the quantum parameter and those corresponding to a supermassive Schwarzschild black hole. Considering the evolution of the EMRI system over the last one year and rescaling the luminosity distance to make the SNR equal to 30, we obtained the critical value for mismatch. Our results indicate that the mismatch between the gravitational waves for the supermassive polymerized black hole with $ \hat{k} < 0.003 $ and those for the supermassive Schwarzschild black hole is smaller than the critical value. This suggests that LISA could probe the parameter $ \hat{k} $ to $ {\cal{O}}(10^{-3}) $ with a one-year observation of the EMRI system. This is tighter than the current limit of $ {\cal{O}}(10^{-1}) $ obtained from the Sgr A$ ^\ast $ black hole shadow [43]. To assess the robustness of our constraints, we investigated the potential degeneracy between the quantum correction parameter $ \hat{k} $ and black hole mass M. Our mismatch analysis showed that the effects of $ \hat{k} $ cannot be mimicked by adjusting M. This indicates that the effects of the two parameters are distinct and lead to non-degenerate states. This supports the reliability of using EMRI waveforms to constrain quantum corrections in the black hole spacetime. Exploring the constraints on the quantum correction parameter through various astronomical observations contributes to the study of quantum-corrected black hole models.

Note that we employed the quasi-circular inspiral approximation to model the orbital evolution of smaller objects in EMRIs. However, for a more comprehensive analysis, the evolution of elliptical orbits should also be considered [10]. Additionally, given that celestial objects generally possess spin, incorporating the effects of spin on orbital evolution is essential [54, 55]. Conversely, the presence of dark matter surrounding a supermassive black hole can influence both the orbital dynamics of smaller objects and the resulting gravitational waveforms [56−58]. A major challenge lies in distinguishing the effects of modified gravity from those of environmental factors on EMRIs. We will attempt to address these issues in our future work.

We thank Rui Niu for important discussions.

The Lagrangian of the test particle with mass m is [44]

where τ is the proper time. For simplicity, we set the mass of the test particle $ m = 1 $. The generalized momentum of the test particle is defined by

where the dot denotes the derivative with respect to proper time. Substituting Eqs. (1) and (23) into Eq. (24), the equations of motion for the test particle are derived as

where E and L represent the energy and orbital angular momentum of the test particle per unit mass, respectively. Here, we only consider the test particle moving on the equatorial plane $ (\theta = \pi/2) $ around the polymerized black hole. Then, from Eq. (23), we obtain

To test the validity of the numerical algorithm that we employed to evolve the quasi-circular orbits of a small object, we compare the small object's energy $ E(t) $ obtained using the numerical algorithm and energy $ \tilde{E}(t) $ obtained using the equations of motion, employing Eqs. (2), (3), and (4), over the entire orbital evolution time. For a small object with $ m = 10 M_\odot $ inspiraling into a supermassive polymerized black hole with $ M = 10^6 M_\odot $ and $ \hat{k} = 0.1 $, we set the initial position $ (r_0, \phi_0) $ of the object as $ (10 M, \; \pi/2) $ and use the orbital evolution algorithm to obtain its complete quasi-circular orbit. We also calculate the evolution of the object's energy and orbital angular momentum as $ (E(t), L(t)) $ through the orbital evolution algorithm and obtain another evolution of the object's energy $ \tilde{E}(t) $ from the equations of motion. Then, we calculate the relative numerical error between $ E(t) $ and $ \tilde{E}(t) $, which is defined as

The results are shown in Fig. B1, indicating that the evolution of the object's energy obtained from Eq. (10) agrees well with that from Eqs. (2), (3), and (4); this implies that the numerical algorithm we employed to evolve the quasi-circular orbits of a small object is reasonable.

Figure B1.  (color online) Relative numerical error between $ E(t) $ from gravitational radiation and $ \tilde{E} $ from equations of motion of the object.