Rational orbits and gravitational waves in static spherical spacetimes: an open-source numerical framework

From the detection of gravitational wave signals emitted by binary black hole mergers by the Laser Interferometer Gravitational-Wave Observatory (LIGO) [1], to the successive imaging of the supermassive black holes in the M87 galaxy and at the Galactic Center (Sgr ${\rm{A}} ^{*} $) by the Event Horizon Telescope (EHT) [2, 3], and further to the observational testing of Hawking's area law [4], the continuously accumulating evidence has made the existence of black holes practically indisputable. Consequently, serving both as natural sources of extreme gravitational fields and as crucial building blocks for refining the theory of gravity, black holes have attracted extensive attention from the scientific community.

Owing to the intense gravitational fields of black holes, small celestial bodies can be gravitationally bound to orbit around them. These orbital motions exhibit a rich variety of patterns, including standard Keplerian orbits, precessing orbits, chaotic motions [5−13], off-equatorial motions [14], and extremely accelerated motions in polar regions [15−17]. Regardless of the specific dynamical behavior, the characteristics of these orbits inherently encode the fundamental properties of the background spacetime, making them a premier probe for investigating black holes. Interestingly, by studying the geodesic motion of timelike particles in black hole spacetimes, Levin and Perez-Giz established a numerical framework for calculating strictly closed timelike orbits [18]. Within a complete orbital cycle, the azimuthal and radial periods of such orbits exhibit an integer ratio, and hence they are referred to as rational orbits. Specifically, the configuration of a rational orbit can be strictly characterized by three integers $ (z,w,v) $, which makes it possible to extract richly structured orbital information within curved spacetimes. Subsequently, this methodology has been widely applied by researchers to various modified theories of gravity [19−64], establishing a direct connection between spacetime parameters and orbital configurations. This ultimately provides an extensive theoretical catalog for observing timelike orbits, refining theories of gravity, and constraining spacetime parameters.

It is worth noting that rational orbits are not merely mathematical constructs; they not only encode spacetime information but also play a significant role in orbital resonance and nonlinear theories. Therefore, further exploring the potential observational signatures of rational orbits holds substantial practical importance. Interestingly, a small celestial body moving on a rational orbit and the central black hole can constitute an extreme-mass-ratio inspiral system, whose mass quadrupole moment evolves with time, thereby generating gravitational radiation. Due to the breakdown of the adiabatic evolution condition at resonances, these radiation signals carry unique fingerprint features, providing a crucial means of extracting the intrinsic properties of spacetime. More importantly, the gravitational wave signals of such systems typically fall within the mHz frequency window, making them core scientific targets for future space gravitational wave observatories, such as LISA [65, 66], TianQin [67−69], and Taiji [70, 71]. This further underscores the need to establish a library of gravitational wave signals from rational orbits and to utilize them as a unique probe for testing gravity in the strong field regime.

Despite the extensive research devoted to rational orbits in curved spacetimes, a universal and comprehensive open-source numerical tool capable of efficiently and conveniently simulating these orbits and their associated gravitational wave signals has not yet been available. This deficiency implicitly impedes further detailed exploration within the field. To address this gap, we have developed a computational platform based on Mathematica and the OpenMP computing framework to realize the aforementioned functionalities, making it seamlessly applicable to any static spherically symmetric spacetime. Specifically, with only the covariant metric defined by the user, the code directly returns circular orbits, rational orbits, irrational orbits, and the specific energy range corresponding to the stable orbital regions, alongside the gravitational wave signals and characteristic strains for the designated parameters.

The remainder of this paper is organized as follows. In the next section, we detail the programming logic and fundamental formulas of the open-source code, mathematically demonstrating how to derive rational orbits, the corresponding gravitational wave signals, and characteristic strains starting from a static spherically symmetric spacetime line element. In Section III, we adopt the Schwarzschild spacetime as a test case to showcase the versatile functionalities of the code, focusing on the observational prospects of gravitational wave signals emitted by an extreme-mass-ratio inspiral system composed of the Galactic Center black hole and a small celestial body. In the final section, we draw our conclusions and provide a discussion. Throughout the mathematical derivations in this paper, we adopt geometrized units, where the black hole mass M, the gravitational constant G, and the speed of light c are all set to unity ($ M=G=c=1 $). When considering practical applications, we restore the relevant quantities to the International System of Units.

We consider a static spherically symmetric spacetime, whose line element in the coordinates $ x^{\mu}=(t,r,\theta,\varphi) $ is written as

where $ g_{\mu\nu} $ the covariant metric, with $ f(r) = -g_{tt} $ representing the metric potential. We introduce the four-velocity $ \dot{x}^{\mu}=(\dot{t},\dot{r},\dot{\theta},\dot{\varphi}) $. Then, the Lagrangian $ {\cal{L}} $ governing the particle motion in the spacetime is written as

Since the spacetime is static and spherically symmetric, the line element (1) does not explicitly depend on time t or the azimuthal angle φ, which leads to two conserved quantities for the particle motion, namely the specific energy E and the specific angular momentum L. According to the Euler-Lagrange equations, the mapping between these conserved quantities and the velocities is given by

Here, $ p_{\mu} $ represents the conjugate momentum corresponding to the generalized coordinates of the particle. Subsequently, we have

Based on the Legendre transformation, we can obtain the Hamiltonian governing the particle motion in the spacetime as

where $ g^{\mu\nu} $ denotes the contravariant metric, and in the spherically symmetric case, $ g^{\mu\nu}=1/g_{\mu\nu} $. Given the initial conditions of the timelike particle $ (t,r,\theta,\varphi) $ and $ (p_{t},p_{r},p_{\theta},p_{\varphi}) $, the particle orbit can be obtained from the canonical equations

We consider a static spherically symmetric spacetime, whose line element in the coordinates $ x^{\mu}=(t,r,\theta,\varphi) $ is written as

where $ g_{\mu\nu} $ the covariant metric, with $ f(r) = -g_{tt} $ representing the metric potential. We introduce the four-velocity $ \dot{x}^{\mu}=(\dot{t},\dot{r},\dot{\theta},\dot{\varphi}) $. Then, the Lagrangian $ {\cal{L}} $ governing the particle motion in the spacetime is written as

Since the spacetime is static and spherically symmetric, the line element (1) does not explicitly depend on time t or the azimuthal angle φ, which leads to two conserved quantities for the particle motion, namely the specific energy E and the specific angular momentum L. According to the Euler-Lagrange equations, the mapping between these conserved quantities and the velocities is given by

Here, $ p_{\mu} $ represents the conjugate momentum corresponding to the generalized coordinates of the particle. Subsequently, we have

Based on the Legendre transformation, we can obtain the Hamiltonian governing the particle motion in the spacetime as

where $ g^{\mu\nu} $ denotes the contravariant metric, and in the spherically symmetric case, $ g^{\mu\nu}=1/g_{\mu\nu} $. Given the initial conditions of the timelike particle $ (t,r,\theta,\varphi) $ and $ (p_{t},p_{r},p_{\theta},p_{\varphi}) $, the particle orbit can be obtained from the canonical equations

We discuss the orbital solutions of timelike particles in the equatorial plane, so we can safely set $ \theta=\pi/2 $ and $ \dot{\theta}=p_{\theta}=0 $, along with the additional conserved quantity of the particle $ {\cal{L}}=-1/2 $. Substituting these conditions into the Lagrangian, we obtain

where $ V_{\rm{eff}} $ represents the particle effective potential ¹, written as

For standard circular orbits in the spacetime, the condition $ \dot{r}=0 $ must be satisfied. Specifically, there exist two types of special circular orbits in the spacetime, namely the innermost stable circular orbit (ISCO) and the marginally bound circular orbit (MBO). The former corresponds to the minimum point of the effective potential and can be determined by solving

to obtain the orbital radius $ r_{\rm{isco}} $ and the corresponding orbital angular momentum $ L_{\rm{isco}} $. Furthermore, by substituting the roots of Eqs. (11) and (12) into Eq. (9), we can obtain the specific energy $ E_{\rm{isco}} $ corresponding to the ISCO. It is worth pointing out that for metrics of the type $ g_{tt}=-1/g_{rr} $, the $ V_{\rm{eff}} $ corresponding to $ r_{\rm{isco}} $ is exactly the specific energy, which is clear from the radial evolution Eq. (9).

As for the MBO, it corresponds to the maximum point of the effective potential and is therefore unstable; any perturbation can cause a timelike particle on the MBO to plunge into the black hole or escape to infinity. Based on condition (11) and (9), with $ E_{\rm{mbo}}=1 $, we can obtain the orbital radius $ r_{\rm{mbo}} $ and $ L_{\rm{mbo}} $ of the MBO.

For any arbitrary angular momentum $ L_{i} $ belonging to the interval $ (L_{\rm{isco}},L_{\rm{mbo}}) $, the effective potential (10) possesses both a maximum point and a minimum point, corresponding to the unstable circular orbit (UnSCO) and the stable circular orbit (SCO) for $ L_{i} $. Their respective radii can be obtained by individually solving Eq. (11) for a given $ L_{i} $. Subsequently, by substituting the two radii and $ L_{i} $ into Eq. (9), we can acquire the specific energies corresponding to the UnSCO and the SCO, denoted as $ E_{\rm{max}} $ and $ E_{\rm{min}} $, respectively. More importantly, by arbitrarily selecting $ E_{i} $ within the interval $ (E_{\rm{min}},E_{\rm{max}}) $, the combination of $ E_{i} $ and $ L_{i} $ corresponds to a stable precessing orbit. In other words, substituting $ E_{i} $ and $ L_{i} $ into the equation $ \dot{r}^{2}=0 $ yields three real roots. Arranged in descending order as $ r_{2} $, $ r_{1} $, and $ r_{0} $, the orbit is confined to oscillate within the range $ [r_{1},r_{2}] $, where $ r_{1} $ and $ r_{2} $ correspond to the periastron and apastron of the orbit, respectively, while $ r_{0} $ lacks physical significance and is merely mathematically valid.

We discuss the orbital solutions of timelike particles in the equatorial plane, so we can safely set $ \theta=\pi/2 $ and $ \dot{\theta}=p_{\theta}=0 $, along with the additional conserved quantity of the particle $ {\cal{L}}=-1/2 $. Substituting these conditions into the Lagrangian, we obtain

where $ V_{\rm{eff}} $ represents the particle effective potential ¹, written as

For standard circular orbits in the spacetime, the condition $ \dot{r}=0 $ must be satisfied. Specifically, there exist two types of special circular orbits in the spacetime, namely the innermost stable circular orbit (ISCO) and the marginally bound circular orbit (MBO). The former corresponds to the minimum point of the effective potential and can be determined by solving

to obtain the orbital radius $ r_{\rm{isco}} $ and the corresponding orbital angular momentum $ L_{\rm{isco}} $. Furthermore, by substituting the roots of Eqs. (11) and (12) into Eq. (9), we can obtain the specific energy $ E_{\rm{isco}} $ corresponding to the ISCO. It is worth pointing out that for metrics of the type $ g_{tt}=-1/g_{rr} $, the $ V_{\rm{eff}} $ corresponding to $ r_{\rm{isco}} $ is exactly the specific energy, which is clear from the radial evolution Eq. (9).

As for the MBO, it corresponds to the maximum point of the effective potential and is therefore unstable; any perturbation can cause a timelike particle on the MBO to plunge into the black hole or escape to infinity. Based on condition (11) and (9), with $ E_{\rm{mbo}}=1 $, we can obtain the orbital radius $ r_{\rm{mbo}} $ and $ L_{\rm{mbo}} $ of the MBO.

For any arbitrary angular momentum $ L_{i} $ belonging to the interval $ (L_{\rm{isco}},L_{\rm{mbo}}) $, the effective potential (10) possesses both a maximum point and a minimum point, corresponding to the unstable circular orbit (UnSCO) and the stable circular orbit (SCO) for $ L_{i} $. Their respective radii can be obtained by individually solving Eq. (11) for a given $ L_{i} $. Subsequently, by substituting the two radii and $ L_{i} $ into Eq. (9), we can acquire the specific energies corresponding to the UnSCO and the SCO, denoted as $ E_{\rm{max}} $ and $ E_{\rm{min}} $, respectively. More importantly, by arbitrarily selecting $ E_{i} $ within the interval $ (E_{\rm{min}},E_{\rm{max}}) $, the combination of $ E_{i} $ and $ L_{i} $ corresponds to a stable precessing orbit. In other words, substituting $ E_{i} $ and $ L_{i} $ into the equation $ \dot{r}^{2}=0 $ yields three real roots. Arranged in descending order as $ r_{2} $, $ r_{1} $, and $ r_{0} $, the orbit is confined to oscillate within the range $ [r_{1},r_{2}] $, where $ r_{1} $ and $ r_{2} $ correspond to the periastron and apastron of the orbit, respectively, while $ r_{0} $ lacks physical significance and is merely mathematically valid.

It is not difficult to find that the combinations of $ E_{i} $ and $ L_{i} $ are infinite, and in most cases, the corresponding precessing orbits are not closed. That is to say, after the timelike particle completes one period in the azimuthal direction, its radial coordinate cannot return to the initial value. As the orbital evolution time increases, the envelope of the precessing orbit eventually forms a circle. Interestingly, within these infinite combinations, there exist certain special cases where closed precessing orbits emerge, which are referred to as rational orbits. Specifically, the period $ T_{\varphi} $ of the azimuthal angle φ of the timelike particle is exactly proportional to the period $ T_{r} $ of the radial coordinate, such as $ T_{\varphi}/T_{r}=q $, where q is an arbitrary rational number.

The authors in [18] point out that the orbital rational number q can reflect the orbital characteristics from two aspects. On one hand, it is defined as

Here, z represents the number of leaves the orbit can exhibit within a complete orbital period, which can also be understood as the number of times the timelike particle passes through the apastron. w counts the nearly circular loops the orbit completes around the black hole near the periastron. In the equatorial plane described by the Cartesian coordinate system $ (x,y) $, a larger w implies that there are more intersection points between the orbit and the x-axis within a complete orbital period. v is used to indicate the direction of the next apastron toward which the timelike particle heads after being released from the current apastron. Usually, along the direction of increasing φ, the apastrons are sequentially labeled as $ v=0,1,2,3\cdots z-1 $. Specifically, for a single leaf orbit ($ z=1 $), v can take the value of 0 or 1. Overall, once $ (z,w,v) $ is given, the orbital configuration is determined and corresponds to the rational number q.

On the other hand, during the process in which the orbit evolves from the apastron to the periastron and then reaches the next adjacent apastron, the accumulation of the azimuthal angle $ \Delta\varphi $ is related to the rational number q:

According to the geodesic integration, we have

where $ \varphi_{1} $ and $ \varphi_{2} $ represent the azimuthal angles of the periastron and apastron, respectively, and the derivatives are written as

By combining Eqs. (14) and (15), we establish the connection between the rational orbital configuration $ (z,w,v) $ and the constants of motion $ L_{i} $ and $ E_{i} $. In addition, we can obtain the proper time elapsed along the rational orbit between adjacent periastron and apastron

Then, the period corresponding to a complete orbit is written as $ T_{q} = 2z\Delta\tau $.

With the theoretical framework described above, the procedure for simulating rational timelike orbits in spherically symmetric spacetimes is clear. First, based on the conditions for the ISCO and the MBO, we calculate $ L_{\rm{isco}} $ and $ L_{\rm{mbo}} $. Second, we select $ L_{i} $ within the open interval $ (L_{\rm{isco}}, L_{\rm{mbo}}) $, use Eq. (11) to calculate the radii of the stable and unstable circular orbits corresponding to this angular momentum, and obtain the specific energies corresponding to these two radii according to Eq. (9). Thus, we acquire the specific energy range $ (E_{\rm{min}},E_{\rm{max}}) $ for stable precessing orbits for a given $ L_{i} $. Third, by scanning $ E_{i} $ within $ (E_{\rm{min}},E_{\rm{max}}) $, we calculate the rational number q corresponding to $ (L_{i},E_{i}) $ according to Eq. (14), and obtain the mapping relationship between q and $ E_{i} $. Fourth, given the configuration parameters $ (z,w,v) $ of the target orbit, we utilize Eq. (13) to calculate the corresponding rational number q, and employ interpolation methods to obtain $ E_{i} $. By substituting $ E_{i} $ and $ L_{i} $ into Eq. (9), we obtain the periastron and apastron radii $ r_{1} $ and $ r_{2} $. Consequently, the initial conditions $ (t,r,\theta,\varphi) $ and $ (p_{t},p_{r},p_{\theta},p_{\varphi}) $ for the rational orbit corresponding to $ (z,w,v) $ are read as $ (0,r_{2},\pi/2,0) $ and $ (-E_{i},0,0,L_{i}) $. Finally, by utilizing a numerical integration algorithm to integrate the canonical Eqs. (8) over the proper time $ T_{q} $, we can obtain the evolution of the target rational orbit within a complete period.

It is not difficult to deduce that the key to obtaining a rational orbit lies in confirming the precise values of the orbital specific energy and specific angular momentum corresponding to $ (z,w,v) $. If these values deviate, the closed nature of the orbit is broken, and the rational orbit transitions into an irrational orbit. To obtain the irrational orbit, one can add a minute perturbation to q. For instance, the $ q=1+1/2 $ corresponding to $ (2,1,1) $ can be modified to $ q=1+1/2+1/200=1+101/200 $. Then, the modified q is subjected to the aforementioned interpolation procedure to obtain the initial conditions of the irrational orbit. The irrational orbit simulated in this manner can still exhibit the characteristics of $ (z,w,v) $ over a short evolution time. However, as the evolution time increases, the envelope of the orbit eventually evolves into a circle.

It is not difficult to find that the combinations of $ E_{i} $ and $ L_{i} $ are infinite, and in most cases, the corresponding precessing orbits are not closed. That is to say, after the timelike particle completes one period in the azimuthal direction, its radial coordinate cannot return to the initial value. As the orbital evolution time increases, the envelope of the precessing orbit eventually forms a circle. Interestingly, within these infinite combinations, there exist certain special cases where closed precessing orbits emerge, which are referred to as rational orbits. Specifically, the period $ T_{\varphi} $ of the azimuthal angle φ of the timelike particle is exactly proportional to the period $ T_{r} $ of the radial coordinate, such as $ T_{\varphi}/T_{r}=q $, where q is an arbitrary rational number.

The authors in [18] point out that the orbital rational number q can reflect the orbital characteristics from two aspects. On one hand, it is defined as

Here, z represents the number of leaves the orbit can exhibit within a complete orbital period, which can also be understood as the number of times the timelike particle passes through the apastron. w counts the nearly circular loops the orbit completes around the black hole near the periastron. In the equatorial plane described by the Cartesian coordinate system $ (x,y) $, a larger w implies that there are more intersection points between the orbit and the x-axis within a complete orbital period. v is used to indicate the direction of the next apastron toward which the timelike particle heads after being released from the current apastron. Usually, along the direction of increasing φ, the apastrons are sequentially labeled as $ v=0,1,2,3\cdots z-1 $. Specifically, for a single leaf orbit ($ z=1 $), v can take the value of 0 or 1. Overall, once $ (z,w,v) $ is given, the orbital configuration is determined and corresponds to the rational number q.

On the other hand, during the process in which the orbit evolves from the apastron to the periastron and then reaches the next adjacent apastron, the accumulation of the azimuthal angle $ \Delta\varphi $ is related to the rational number q:

According to the geodesic integration, we have

where $ \varphi_{1} $ and $ \varphi_{2} $ represent the azimuthal angles of the periastron and apastron, respectively, and the derivatives are written as

By combining Eqs. (14) and (15), we establish the connection between the rational orbital configuration $ (z,w,v) $ and the constants of motion $ L_{i} $ and $ E_{i} $. In addition, we can obtain the proper time elapsed along the rational orbit between adjacent periastron and apastron

Then, the period corresponding to a complete orbit is written as $ T_{q} = 2z\Delta\tau $.

With the theoretical framework described above, the procedure for simulating rational timelike orbits in spherically symmetric spacetimes is clear. First, based on the conditions for the ISCO and the MBO, we calculate $ L_{\rm{isco}} $ and $ L_{\rm{mbo}} $. Second, we select $ L_{i} $ within the open interval $ (L_{\rm{isco}}, L_{\rm{mbo}}) $, use Eq. (11) to calculate the radii of the stable and unstable circular orbits corresponding to this angular momentum, and obtain the specific energies corresponding to these two radii according to Eq. (9). Thus, we acquire the specific energy range $ (E_{\rm{min}},E_{\rm{max}}) $ for stable precessing orbits for a given $ L_{i} $. Third, by scanning $ E_{i} $ within $ (E_{\rm{min}},E_{\rm{max}}) $, we calculate the rational number q corresponding to $ (L_{i},E_{i}) $ according to Eq. (14), and obtain the mapping relationship between q and $ E_{i} $. Fourth, given the configuration parameters $ (z,w,v) $ of the target orbit, we utilize Eq. (13) to calculate the corresponding rational number q, and employ interpolation methods to obtain $ E_{i} $. By substituting $ E_{i} $ and $ L_{i} $ into Eq. (9), we obtain the periastron and apastron radii $ r_{1} $ and $ r_{2} $. Consequently, the initial conditions $ (t,r,\theta,\varphi) $ and $ (p_{t},p_{r},p_{\theta},p_{\varphi}) $ for the rational orbit corresponding to $ (z,w,v) $ are read as $ (0,r_{2},\pi/2,0) $ and $ (-E_{i},0,0,L_{i}) $. Finally, by utilizing a numerical integration algorithm to integrate the canonical Eqs. (8) over the proper time $ T_{q} $, we can obtain the evolution of the target rational orbit within a complete period.

It is not difficult to deduce that the key to obtaining a rational orbit lies in confirming the precise values of the orbital specific energy and specific angular momentum corresponding to $ (z,w,v) $. If these values deviate, the closed nature of the orbit is broken, and the rational orbit transitions into an irrational orbit. To obtain the irrational orbit, one can add a minute perturbation to q. For instance, the $ q=1+1/2 $ corresponding to $ (2,1,1) $ can be modified to $ q=1+1/2+1/200=1+101/200 $. Then, the modified q is subjected to the aforementioned interpolation procedure to obtain the initial conditions of the irrational orbit. The irrational orbit simulated in this manner can still exhibit the characteristics of $ (z,w,v) $ over a short evolution time. However, as the evolution time increases, the envelope of the orbit eventually evolves into a circle.

Under general circumstances, small celestial bodies orbit around supermassive black holes at the centers of galaxies, and the mass ratio between them can exceed $ 10^{4} $. The orbits of these small celestial bodies can be characterized by geodesics and inherently encode spacetime information. More importantly, the extreme-mass-ratio inspiral systems composed of small celestial bodies and central black holes serve as observational targets for space gravitational wave projects. To this end, our code further integrates the computation of gravitational wave waveforms and characteristic strains radiated by precessing orbits, including both rational and irrational orbits.

In [72], the authors provided approximate analytical formulas for calculating the gravitational wave polarization states of extreme-mass-ratio inspiral systems, which are sufficient to capture the essential gravitational wave features, written as

Here, η is the symmetric mass ratio parameter of the system, given by $ \eta=mM/(m+M)^{2} $, where m is the mass of the small celestial body. $ D_{\rm{L}} $ represents the luminosity distance from the extreme-mass-ratio inspiral system to the Earth. ι and ζ denote the inclination angle and the longitude of the periastron of the orbit, respectively. After obtaining r and φ at any given moment of the orbit, substituting them into Eqs. (19) and (20) allows the time-domain information of the two gravitational wave polarization states to be obtained directly.

Based on the fast Fourier transform, we can convert the gravitational wave signals calculated from Eqs. (19) and (20) into the frequency-domain, yielding $ \hat{h}_{+}\left(f\right) $ and $ \hat{h}_{\times}\left(f\right) $, where f represents the frequency. Furthermore, by combining the frequency-domain signals of the two polarization states according to [73]

we can obtain the characteristic strain $ h_{\rm{c}} $ of the gravitational wave. By comparing the $ h_{\rm{c}} $ of different orbits with the sensitivity curves of space gravitational wave detectors, the observational prospects of the orbital gravitational radiation can be thoroughly discussed.

Under general circumstances, small celestial bodies orbit around supermassive black holes at the centers of galaxies, and the mass ratio between them can exceed $ 10^{4} $. The orbits of these small celestial bodies can be characterized by geodesics and inherently encode spacetime information. More importantly, the extreme-mass-ratio inspiral systems composed of small celestial bodies and central black holes serve as observational targets for space gravitational wave projects. To this end, our code further integrates the computation of gravitational wave waveforms and characteristic strains radiated by precessing orbits, including both rational and irrational orbits.

In [72], the authors provided approximate analytical formulas for calculating the gravitational wave polarization states of extreme-mass-ratio inspiral systems, which are sufficient to capture the essential gravitational wave features, written as

Here, η is the symmetric mass ratio parameter of the system, given by $ \eta=mM/(m+M)^{2} $, where m is the mass of the small celestial body. $ D_{\rm{L}} $ represents the luminosity distance from the extreme-mass-ratio inspiral system to the Earth. ι and ζ denote the inclination angle and the longitude of the periastron of the orbit, respectively. After obtaining r and φ at any given moment of the orbit, substituting them into Eqs. (19) and (20) allows the time-domain information of the two gravitational wave polarization states to be obtained directly.

Based on the fast Fourier transform, we can convert the gravitational wave signals calculated from Eqs. (19) and (20) into the frequency-domain, yielding $ \hat{h}_{+}\left(f\right) $ and $ \hat{h}_{\times}\left(f\right) $, where f represents the frequency. Furthermore, by combining the frequency-domain signals of the two polarization states according to [73]

we can obtain the characteristic strain $ h_{\rm{c}} $ of the gravitational wave. By comparing the $ h_{\rm{c}} $ of different orbits with the sensitivity curves of space gravitational wave detectors, the observational prospects of the orbital gravitational radiation can be thoroughly discussed.

In geometrized units, the covariant metric of the Schwarzschild spacetime is written as

Our code requires only the definition of the covariant metric to operate with full functionality.

According to the conditions for the ISCO and the MBO, we obtain $ r_{\rm{isco}}=6 $ and $ r_{\rm{mbo}}=4 $ for the Schwarzschild spacetime, with the corresponding orbital angular momenta being $ L_{\rm{isco}}=3.4641 $ and $ L_{\rm{mbo}}=4 $. To quantify the selection of $ L_{i} $, we introduce the angular momentum parameter $ \varepsilon \in (0,1) $, such that

We take $ \varepsilon=0.5 $, which corresponds to $ L_{i}=3.7321 $, as an example to illustrate the results. Figure 1 displays the evolution of the effective potential $ V_{\rm{eff}} $ with respect to the radial coordinate r. It can be observed that the effective potential curve exhibits two extremum points, where the maximum point corresponds to the unstable circular orbit and the minimum point corresponds to the stable circular orbit. Their radii can be obtained by solving Eq. (11), yielding $ 4.37294 $ and $ 9.55526 $ respectively. According to Eq. (9), we further determine the specific energies corresponding to the maximum and minimum values to be $ E_{\rm{max}}=0.968443 $ and $ E_{\rm{min}}=0.954626 $, respectively. In other words, for $ L_{i}=3.7321 $, there exists a specific energy range $ (E_{\rm{min}},E_{\rm{max}}) $ that allows stable precessing orbits to exist.

Figure 1.  (color online) Effective potential $ V_{\rm{eff}} $ as a function of r in the case of $ L_{i}=3.7321 $. The red and blue star markers correspond to the maximum and minimum points, respectively, representing the locations of the unstable circular orbit and the stable circular orbit.

We select $ E_{1}=0.959 $, $ E_{2}=0.963 $, and $ E_{3}=0.967 $ within the specific energy range and plot the variation of $ \dot{r}^{2} $ with respect to r, as illustrated in Fig. 2, where the black dash-dotted lines represent $ E_{\rm{min}} $ and $ E_{\rm{max}} $. We find that the curves corresponding to the three different specific energies $ E_{i} $ are entirely located between the two critical dash-dotted lines, and they intersect the line $ \dot{r}^{2}=0 $ three times, indicating the existence of precessing orbits. The timelike particles can perform oscillatory motion within the lobes of their respective curves where $ \dot{r}^{2}>0 $ (between the two larger roots of $ \dot{r}^{2}=0 $). This result not only demonstrates the high precision of the code but also verifies the feasibility of the previously mentioned method for finding stable precessing orbits.

Figure 2.  (color online) $ \dot{r}^{2} $ as a function of r. The red, green, and blue solid lines correspond to $ E_{1}=0.959 $, $ E_{2}=0.963 $, and $ E_{3}=0.967 $ respectively. The upper and lower dash-dotted lines correspond to $ E_{\rm{max}} $ and $ E_{\rm{min}} $, respectively. It can be observed that the red, blue, and green lines intersect the line $ \dot{r}^{2}=0 $ three times, indicating the existence of stable orbits, and the two larger roots correspond to the apastron and periastron of the orbits. In contrast, the two dash-dotted lines intersect the line $ \dot{r}^{2}=0 $ only twice, corresponding to circular orbits.

Our code also supports searching for the conditions for stable precessing orbits in a broader parameter space. Figure 3 illustrates the evolution of $ E_{\rm{min}} $ (red solid line) and $ E_{\rm{max}} $ (blue solid line) with respect to the specific angular momentum $ L_{i} $. It can be found that both critical values of E increase as $ L_{i} $ increases, and the region between them resembles a triangle. Arbitrarily selecting a combination of specific energy and specific angular momentum within the shaded region yields a stable bound orbit.

Figure 3.  (color online) Evolution of the maximum specific energy $ E_{\rm{max}} $ and the minimum specific energy $ E_{\rm{min}} $ corresponding to stable orbits with respect to the orbital angular momentum $ L_{i} $. The shaded region between the two curves represents the parameter space where stable orbits exist.

Next, we return to the case of $ \varepsilon=0.5 $. Figure 4 illustrates the evolution of the rational number q with respect to $ E_{i} $, obtained by combining Eqs. (14) and (15), where $ E_{i} \in (E_{\rm{min}},E_{\rm{max}}) $. It can be observed that the variation of q initially exhibits a slow increase as $ E_{i} $ grows, and subsequently experiences a sharp increase when the specific energy reaches $ E_{\rm{max}} $. This occurs because when $ E_{i} $ approaches its maximum value, the orbit tends to circularize, which naturally leads to a rapid accumulation of $ \Delta\varphi $. Furthermore, we find that in this scenario, q possesses a lower bound of approximately $ q=0.637 $ (the left endpoint of the curve). Combined with the connection between q and the orbital configuration shown in Eq. (13), this phenomenon indicates that the dynamics naturally filter out certain structures; for instance, the configuration $ q=w+v/z=0+1/2 $ is strictly prohibited.

Figure 4.  Rational number q as a function of the specific energy $ E_{i} $. It is clear that q exhibits a lower limit since $ E_{i} $ has a minimum value.

Fixing $ w=1 $, we simulate the rational orbits for $ z=1, 2, 3, 4, 5 $ for one period each, where $ v=1 $ for $ z=1 $, and $ v=z-1 $ for $ z \gt 1 $. As illustrated in the first row of Fig. 5, the starting point of each orbit is located at the rightmost apastron. It can be found that the number of orbital leaves, namely the number of apastrons, is perfectly consistent with z. The process of the orbit evolving from one apastron to the adjacent apastron produces $ w+1 $ loops around the central black hole; therefore, as z increases, the structure in the central region of the orbit becomes increasingly complex. Subsequently, we add a minute quantity to the rational number q to observe the behavior of irrational orbits. As shown in the second row of Fig. 5, each orbit is simulated for three times the period of the original rational orbit. It can be observed that the irrational orbits lack strict periodicity and cannot perfectly overlap.

Figure 5.  Rational orbits (first row) and irrational orbits (second row) corresponding to different orbital configurations. From left to right, the numbers of orbital leaves are $ z=1 $, 2, 3, 4, and 5, and we have $ v=z-1 $ for $ z \gt 1 $, while $ v=1 $ for $ z=1 $. The irrational orbits are obtained by adding $ 1/100z $ to q. Here, we fix $ w=1 $. The simulated proper time for rational orbits is exactly one orbital period $ T_{q} $, whereas the simulation time for irrational orbits is three times $ T_{q} $.

Upon obtaining the orbits, we fix ι and ζ to $ \pi/4 $. By utilizing Eqs. (19) and (20), the code returns the two gravitational wave polarization states $ h_{+} $ and $ h_{\times} $ corresponding to the orbits. Meanwhile, since the gravitational waves are potentially observable, in this step we restore the geometrized units to the International System of Units, which requires specifying the specific parameters of the source. We consider the supermassive black hole at the Galactic Center, with a black hole mass $ M \approx 4 \times 10^{6} $ $ M_{\odot} $ and a distance $ D_{\rm{L}} \approx 8000 $ pc, and assume the mass of the timelike small celestial body to be $ m = 100 $ $ M_{\odot} $. Then, the intensity scaling factor of the gravitational waves is $ \eta GM/(D_{\rm{L}}c^{2}) \approx 5.9887 \times 10^{-16} $, and the time scaling factor is $ GM/c^{3} \approx 19.7088 $. Figure 6 presents the gravitational wave polarization states corresponding to the orbits of various configurations in Fig. 5, where the simulation time for rational orbits (left column) is exactly one orbital period, and that for irrational orbits is three times the period of the rational orbits. It is not difficult to find that the oscillatory behaviors of the two gravitational wave polarization states are similar, and the evolution of their amplitudes is also close. The gravitational wave waveforms are all composed of two components: one is the high-frequency oscillatory component, such as the segment where $ \tau= 4000 \sim 5000 $ in panel (a); the other is the relatively flat region containing broad peaks that connects the high-frequency oscillations, such as the segment where $ \tau= 5000 \sim 12000 $ in panel (c). Some of the authors of this paper pointed out in [74] that the high-frequency oscillatory component originates from the orbital motion near the periastron, while the flat region arises from the contribution of the motion of the small celestial body near the apastron. Therefore, we can observe that the number of high-frequency oscillations in the waveform perfectly matches the number of orbital leaves z. On the other hand, we find that the gravitational wave waveforms of rational orbits and irrational orbits are almost identical within a relatively short period of time.

Figure 6.  (color online) Time-domain gravitational wave signals of rational orbits (left column) and irrational orbits (right column). From top to bottom, z sequentially increases from 1 to 5. In each panel, red and blue correspond to the $ h_{+} $ and $ h_{\times} $ polarization states, respectively. The left column shows that within one orbital period, the number of high-frequency oscillations in the gravitational wave signal is consistent with the number of orbital leaves.

Our code employs the fast Fourier transform to process the time-domain signals of the gravitational waves to obtain frequency-domain information. Figure 7 displays the frequency-domain signals of the two gravitational wave polarization states for the rational orbits in the first row of Fig. 5, where red corresponds to the “plus” polarization state and blue corresponds to the ''cross'' polarization state. It is found that the signal windows of both polarization states range from $ 10^{-4} \sim 10^{-1.5} $ Hz, within which the signal strength at the mHz level and below remains robust. We also find that the amplitude range of $ h_{\times} $ is larger than that of $ h_{+} $, and the former can even drop as low as $ 10^{-25} $ when the frequency f is higher than the order of $ 10^{-2} $ Hz. Furthermore, it is difficult to identify the orbital characteristics using solely the frequency-domain signals of the gravitational waves.

Figure 7.  (color online) Frequency-domain information obtained via the fast Fourier transform from the time-domain gravitational wave signals of the rational orbits in Fig. 6. From left to right, the number of orbital leaves z sequentially increases from 1 to 5. It can be observed that the frequency-domain amplitudes of the gravitational waves are relatively strong below 0.01 Hz. Meanwhile, it is almost impossible to establish a connection between the frequency-domain gravitational wave signals and the orbital configurations.

Utilizing Eq. (21), our numerical tool further obtains the characteristic strain $ h_{\rm{c}} $ of the gravitational waves for each orbit, as shown in Fig. 8, where the black curve represents the sensitivity curve of LISA [65]. It can be observed that, although these portions are relatively small, the characteristic strain curves of the gravitational waves for the five rational orbits exhibit portions that are higher than the LISA sensitivity curve. This indicates that the gravitational radiation from rational orbits has the potential to be observed. On the other hand, our model describes a dynamical system composed of an intermediate mass black hole and the Galactic Center supermassive black hole. The fact that the gravitational wave signals of this system can be detected by LISA implies that there is a theoretical foundation for utilizing LISA to search for intermediate mass black holes near the Galactic Center.

Figure 8.  (color online) Comparison between the characteristic strains of gravitational waves for different orbital configurations and the LISA sensitivity curve (black). It can be found that the characteristic strains of all the displayed orbital configurations have portions exceeding the LISA sensitivity curve, and these correspond to the mHz frequency window. This demonstrates that the gravitational radiation from the rational orbits investigated in this paper has the potential to be observed in the future.

In geometrized units, the covariant metric of the Schwarzschild spacetime is written as

Our code requires only the definition of the covariant metric to operate with full functionality.

According to the conditions for the ISCO and the MBO, we obtain $ r_{\rm{isco}}=6 $ and $ r_{\rm{mbo}}=4 $ for the Schwarzschild spacetime, with the corresponding orbital angular momenta being $ L_{\rm{isco}}=3.4641 $ and $ L_{\rm{mbo}}=4 $. To quantify the selection of $ L_{i} $, we introduce the angular momentum parameter $ \varepsilon \in (0,1) $, such that

We take $ \varepsilon=0.5 $, which corresponds to $ L_{i}=3.7321 $, as an example to illustrate the results. Figure 1 displays the evolution of the effective potential $ V_{\rm{eff}} $ with respect to the radial coordinate r. It can be observed that the effective potential curve exhibits two extremum points, where the maximum point corresponds to the unstable circular orbit and the minimum point corresponds to the stable circular orbit. Their radii can be obtained by solving Eq. (11), yielding $ 4.37294 $ and $ 9.55526 $ respectively. According to Eq. (9), we further determine the specific energies corresponding to the maximum and minimum values to be $ E_{\rm{max}}=0.968443 $ and $ E_{\rm{min}}=0.954626 $, respectively. In other words, for $ L_{i}=3.7321 $, there exists a specific energy range $ (E_{\rm{min}},E_{\rm{max}}) $ that allows stable precessing orbits to exist.

Figure 1.  (color online) Effective potential $ V_{\rm{eff}} $ as a function of r in the case of $ L_{i}=3.7321 $. The red and blue star markers correspond to the maximum and minimum points, respectively, representing the locations of the unstable circular orbit and the stable circular orbit.

We select $ E_{1}=0.959 $, $ E_{2}=0.963 $, and $ E_{3}=0.967 $ within the specific energy range and plot the variation of $ \dot{r}^{2} $ with respect to r, as illustrated in Fig. 2, where the black dash-dotted lines represent $ E_{\rm{min}} $ and $ E_{\rm{max}} $. We find that the curves corresponding to the three different specific energies $ E_{i} $ are entirely located between the two critical dash-dotted lines, and they intersect the line $ \dot{r}^{2}=0 $ three times, indicating the existence of precessing orbits. The timelike particles can perform oscillatory motion within the lobes of their respective curves where $ \dot{r}^{2}>0 $ (between the two larger roots of $ \dot{r}^{2}=0 $). This result not only demonstrates the high precision of the code but also verifies the feasibility of the previously mentioned method for finding stable precessing orbits.

Figure 2.  (color online) $ \dot{r}^{2} $ as a function of r. The red, green, and blue solid lines correspond to $ E_{1}=0.959 $, $ E_{2}=0.963 $, and $ E_{3}=0.967 $ respectively. The upper and lower dash-dotted lines correspond to $ E_{\rm{max}} $ and $ E_{\rm{min}} $, respectively. It can be observed that the red, blue, and green lines intersect the line $ \dot{r}^{2}=0 $ three times, indicating the existence of stable orbits, and the two larger roots correspond to the apastron and periastron of the orbits. In contrast, the two dash-dotted lines intersect the line $ \dot{r}^{2}=0 $ only twice, corresponding to circular orbits.

Our code also supports searching for the conditions for stable precessing orbits in a broader parameter space. Figure 3 illustrates the evolution of $ E_{\rm{min}} $ (red solid line) and $ E_{\rm{max}} $ (blue solid line) with respect to the specific angular momentum $ L_{i} $. It can be found that both critical values of E increase as $ L_{i} $ increases, and the region between them resembles a triangle. Arbitrarily selecting a combination of specific energy and specific angular momentum within the shaded region yields a stable bound orbit.

Figure 3.  (color online) Evolution of the maximum specific energy $ E_{\rm{max}} $ and the minimum specific energy $ E_{\rm{min}} $ corresponding to stable orbits with respect to the orbital angular momentum $ L_{i} $. The shaded region between the two curves represents the parameter space where stable orbits exist.

Next, we return to the case of $ \varepsilon=0.5 $. Figure 4 illustrates the evolution of the rational number q with respect to $ E_{i} $, obtained by combining Eqs. (14) and (15), where $ E_{i} \in (E_{\rm{min}},E_{\rm{max}}) $. It can be observed that the variation of q initially exhibits a slow increase as $ E_{i} $ grows, and subsequently experiences a sharp increase when the specific energy reaches $ E_{\rm{max}} $. This occurs because when $ E_{i} $ approaches its maximum value, the orbit tends to circularize, which naturally leads to a rapid accumulation of $ \Delta\varphi $. Furthermore, we find that in this scenario, q possesses a lower bound of approximately $ q=0.637 $ (the left endpoint of the curve). Combined with the connection between q and the orbital configuration shown in Eq. (13), this phenomenon indicates that the dynamics naturally filter out certain structures; for instance, the configuration $ q=w+v/z=0+1/2 $ is strictly prohibited.

Figure 4.  Rational number q as a function of the specific energy $ E_{i} $. It is clear that q exhibits a lower limit since $ E_{i} $ has a minimum value.

Fixing $ w=1 $, we simulate the rational orbits for $ z=1, 2, 3, 4, 5 $ for one period each, where $ v=1 $ for $ z=1 $, and $ v=z-1 $ for $ z \gt 1 $. As illustrated in the first row of Fig. 5, the starting point of each orbit is located at the rightmost apastron. It can be found that the number of orbital leaves, namely the number of apastrons, is perfectly consistent with z. The process of the orbit evolving from one apastron to the adjacent apastron produces $ w+1 $ loops around the central black hole; therefore, as z increases, the structure in the central region of the orbit becomes increasingly complex. Subsequently, we add a minute quantity to the rational number q to observe the behavior of irrational orbits. As shown in the second row of Fig. 5, each orbit is simulated for three times the period of the original rational orbit. It can be observed that the irrational orbits lack strict periodicity and cannot perfectly overlap.

Figure 5.  Rational orbits (first row) and irrational orbits (second row) corresponding to different orbital configurations. From left to right, the numbers of orbital leaves are $ z=1 $, 2, 3, 4, and 5, and we have $ v=z-1 $ for $ z \gt 1 $, while $ v=1 $ for $ z=1 $. The irrational orbits are obtained by adding $ 1/100z $ to q. Here, we fix $ w=1 $. The simulated proper time for rational orbits is exactly one orbital period $ T_{q} $, whereas the simulation time for irrational orbits is three times $ T_{q} $.

Upon obtaining the orbits, we fix ι and ζ to $ \pi/4 $. By utilizing Eqs. (19) and (20), the code returns the two gravitational wave polarization states $ h_{+} $ and $ h_{\times} $ corresponding to the orbits. Meanwhile, since the gravitational waves are potentially observable, in this step we restore the geometrized units to the International System of Units, which requires specifying the specific parameters of the source. We consider the supermassive black hole at the Galactic Center, with a black hole mass $ M \approx 4 \times 10^{6} $ $ M_{\odot} $ and a distance $ D_{\rm{L}} \approx 8000 $ pc, and assume the mass of the timelike small celestial body to be $ m = 100 $ $ M_{\odot} $. Then, the intensity scaling factor of the gravitational waves is $ \eta GM/(D_{\rm{L}}c^{2}) \approx 5.9887 \times 10^{-16} $, and the time scaling factor is $ GM/c^{3} \approx 19.7088 $. Figure 6 presents the gravitational wave polarization states corresponding to the orbits of various configurations in Fig. 5, where the simulation time for rational orbits (left column) is exactly one orbital period, and that for irrational orbits is three times the period of the rational orbits. It is not difficult to find that the oscillatory behaviors of the two gravitational wave polarization states are similar, and the evolution of their amplitudes is also close. The gravitational wave waveforms are all composed of two components: one is the high-frequency oscillatory component, such as the segment where $ \tau= 4000 \sim 5000 $ in panel (a); the other is the relatively flat region containing broad peaks that connects the high-frequency oscillations, such as the segment where $ \tau= 5000 \sim 12000 $ in panel (c). Some of the authors of this paper pointed out in [74] that the high-frequency oscillatory component originates from the orbital motion near the periastron, while the flat region arises from the contribution of the motion of the small celestial body near the apastron. Therefore, we can observe that the number of high-frequency oscillations in the waveform perfectly matches the number of orbital leaves z. On the other hand, we find that the gravitational wave waveforms of rational orbits and irrational orbits are almost identical within a relatively short period of time.

Figure 6.  (color online) Time-domain gravitational wave signals of rational orbits (left column) and irrational orbits (right column). From top to bottom, z sequentially increases from 1 to 5. In each panel, red and blue correspond to the $ h_{+} $ and $ h_{\times} $ polarization states, respectively. The left column shows that within one orbital period, the number of high-frequency oscillations in the gravitational wave signal is consistent with the number of orbital leaves.

Our code employs the fast Fourier transform to process the time-domain signals of the gravitational waves to obtain frequency-domain information. Figure 7 displays the frequency-domain signals of the two gravitational wave polarization states for the rational orbits in the first row of Fig. 5, where red corresponds to the “plus” polarization state and blue corresponds to the ''cross'' polarization state. It is found that the signal windows of both polarization states range from $ 10^{-4} \sim 10^{-1.5} $ Hz, within which the signal strength at the mHz level and below remains robust. We also find that the amplitude range of $ h_{\times} $ is larger than that of $ h_{+} $, and the former can even drop as low as $ 10^{-25} $ when the frequency f is higher than the order of $ 10^{-2} $ Hz. Furthermore, it is difficult to identify the orbital characteristics using solely the frequency-domain signals of the gravitational waves.

Figure 7.  (color online) Frequency-domain information obtained via the fast Fourier transform from the time-domain gravitational wave signals of the rational orbits in Fig. 6. From left to right, the number of orbital leaves z sequentially increases from 1 to 5. It can be observed that the frequency-domain amplitudes of the gravitational waves are relatively strong below 0.01 Hz. Meanwhile, it is almost impossible to establish a connection between the frequency-domain gravitational wave signals and the orbital configurations.

Utilizing Eq. (21), our numerical tool further obtains the characteristic strain $ h_{\rm{c}} $ of the gravitational waves for each orbit, as shown in Fig. 8, where the black curve represents the sensitivity curve of LISA [65]. It can be observed that, although these portions are relatively small, the characteristic strain curves of the gravitational waves for the five rational orbits exhibit portions that are higher than the LISA sensitivity curve. This indicates that the gravitational radiation from rational orbits has the potential to be observed. On the other hand, our model describes a dynamical system composed of an intermediate mass black hole and the Galactic Center supermassive black hole. The fact that the gravitational wave signals of this system can be detected by LISA implies that there is a theoretical foundation for utilizing LISA to search for intermediate mass black holes near the Galactic Center.

Figure 8.  (color online) Comparison between the characteristic strains of gravitational waves for different orbital configurations and the LISA sensitivity curve (black). It can be found that the characteristic strains of all the displayed orbital configurations have portions exceeding the LISA sensitivity curve, and these correspond to the mHz frequency window. This demonstrates that the gravitational radiation from the rational orbits investigated in this paper has the potential to be observed in the future.

In this paper, based on Mathematica and the OpenMP framework, we have developed a numerical code for solving rational orbits and irrational orbits applicable to any static spherically symmetric spacetime, which also supports the further calculation of gravitational wave signals in the time-domain and frequency-domain, as well as the characteristic strains ². Users only need to define the covariant metric of the target spacetime and set the parameters of the source to operate it smoothly. Meanwhile, our code exhibits high execution efficiency. On a platform with an Intel i7–14700 processor and 8 GB RAM, simulating and plotting the data from Figs. 3$ -$8 takes only 65 seconds. This numerical tool possesses considerable guiding value in the fields of orbital dynamics in curved spacetimes and gravitational wave astronomy.

Utilizing this code, we have demonstrated the rational orbits and irrational orbits of various orbital configurations in the Schwarzschild spacetime and discovered that the number of high-frequency oscillations in the time-domain gravitational wave waveforms perfectly matches the number of orbital leaves. More importantly, by comparing the characteristic strains of the gravitational radiation from each orbit with the LISA sensitivity curve, we find that the gravitational waves radiated by extreme-mass-ratio inspiral systems composed of an intermediate mass black hole and the Galactic Center supermassive black hole are promising targets for future observation. This further demonstrates that utilizing LISA to search for intermediate mass black holes near the Galactic Center is feasible.

Furthermore, it is worth mentioning that astrophysical black holes often inherit the angular momentum of their progenitor stars and consequently possess spin. Their spacetimes are therefore predominantly characterized by axisymmetric geometries. Thus, it is necessary to extend the applicable scope of our code to rotating spacetimes. In this procedure, the crucial step is to adapt the calculations for the timelike particle radial Eq. (9), the effective potential (10), and the azimuthal angle accumulation $ \Delta\varphi $ to their counterparts for rotating black holes. In future work, we plan to expand the capability of this code to automatically search for rational orbits within rotating spacetimes.

In this paper, based on Mathematica and the OpenMP framework, we have developed a numerical code for solving rational orbits and irrational orbits applicable to any static spherically symmetric spacetime, which also supports the further calculation of gravitational wave signals in the time-domain and frequency-domain, as well as the characteristic strains ². Users only need to define the covariant metric of the target spacetime and set the parameters of the source to operate it smoothly. Meanwhile, our code exhibits high execution efficiency. On a platform with an Intel i7–14700 processor and 8 GB RAM, simulating and plotting the data from Figs. 3$ -$8 takes only 65 seconds. This numerical tool possesses considerable guiding value in the fields of orbital dynamics in curved spacetimes and gravitational wave astronomy.

Utilizing this code, we have demonstrated the rational orbits and irrational orbits of various orbital configurations in the Schwarzschild spacetime and discovered that the number of high-frequency oscillations in the time-domain gravitational wave waveforms perfectly matches the number of orbital leaves. More importantly, by comparing the characteristic strains of the gravitational radiation from each orbit with the LISA sensitivity curve, we find that the gravitational waves radiated by extreme-mass-ratio inspiral systems composed of an intermediate mass black hole and the Galactic Center supermassive black hole are promising targets for future observation. This further demonstrates that utilizing LISA to search for intermediate mass black holes near the Galactic Center is feasible.

Furthermore, it is worth mentioning that astrophysical black holes often inherit the angular momentum of their progenitor stars and consequently possess spin. Their spacetimes are therefore predominantly characterized by axisymmetric geometries. Thus, it is necessary to extend the applicable scope of our code to rotating spacetimes. In this procedure, the crucial step is to adapt the calculations for the timelike particle radial Eq. (9), the effective potential (10), and the azimuthal angle accumulation $ \Delta\varphi $ to their counterparts for rotating black holes. In future work, we plan to expand the capability of this code to automatically search for rational orbits within rotating spacetimes.