Orbital Stability Study of the Taiji Space Gravitational Wave Detector

Space-based gravitational wave detection is extremely sensitive to disturbances. The Keplerian configuration cannot accurately reflect the variations in spacecraft configuration. Planetary gravitational disturbances are one of the main sources. Numerical simulation is an effective method to investigate the impact of perturbation on spacecraft orbits. This study shows that, in the context of the Taiji project, Earth's gravity is an essential factor in the change in heliocentric formation configuration, contributing to the relative acceleration between spacecrafts in the order of $\mathcal O(10^{-6})\,{\rm m\cdot s^{-2}}$. Considering 00:00:00 on 27 October 2032 as the initial orbiting moment, under the influence of Earth's gravitational perturbation, the maximum relative change in armlengths and variation rates of armlengths for Taiji is $1.6\times 10^{5}\,{\rm km}$, $32\,{\rm m\cdot s^{-1}}$, respectively, compared with the unperturbed Keplerian orbit. Additionally, by considering the gravitational perturbations of Venus and Jupiter, the armlength and relative velocity for Taiji are reduced by $16.01\%$ and $17.45\%$, respectively, compared with when only considering that of Earth. The maximum amplitude of the formation motion indicator changes with the orbit entry time. Results show that the relative velocity increase between the spacecrafts is minimal when the initial orbital moment occurs in July. Moreover, the numerical simulation results are inconsistent when using different ephemerides. The differences between ephemerides DE440 and DE430 are smaller than those between DE440 and DE421.


I. INTRODUCTION
On 11 February 2016, the Laser Interference Gravitational-Wave Observatory (LIGO) [1] announced the direct observation of gravitational wave signals, opening a new era of physics [2,3].Compared with ground-based gravitational wave detection, space-based detection can address the limitations of ground-based noise [4][5][6] and interferometer scale, enabling the detection of gravitational waves in the medium-and low-frequency bands [7][8][9].The success of the Laser Interferometer Space Antenna Pathfinder (LPF) validates the feasibility of gravitational wave detection in space [10][11][12][13].The subsequent detection plan, the Laser Interferometer Space Antenna (LISA) project, is currently under pre-research and simulation, with a targeted launch date of around 2035 [14].
In 2016, the Chinese Academy of Sciences proposed the Taiji Space Gravitational Wave Detection Program [15,16].The program aims to use Taiji spacecrafts to form a formation of three spacecrafts in Earth-like orbits around the Sun, describing an equilateral triangle with a side length of 3 million km [17].These three spacecrafts would be distributed approximately 20 • in front of (or behind) the Earth, and the angle between the spacecraft formation plane and the ecliptic plane would be 60 • [18,19].
Space-based gravitational wave detection detects gravitational waves through the use of laser interferometry [20].As a gravitational wave passes through, it weakly distorts space, causing a small change in the length of the interferometer arms.By measuring the interference effect of the laser beam in the interferometer, this small change in arm length can be detected, and a gravitational wave signal is obtained.In space-based gravitational wave detection, the length of the laser interference arm changes with time [21].The instability of the laser frequency is one of the main sources of noise.Time-delay interferometry (TDI) technology [22,23] is proposed to reduce this measurement noise.The implementation of TDI requires precise spacecraft absolute range measurement and real-time communication.Thus, the stability of spacecraft orbits should be investigated.Spacecrafts operate in a complicated gravitational environment.The orbit calculations must consider a variety of disturbances, such as planetary gravity, planetary tidal force, lunar gravity, solar light pressure, and post-Newtonian effects [24].Planetary gravity is one of the dominant factors, especially the gravities of Earth, Venus, and Jupiter, which has a significant impact on orbital stability [8,25].Furthermore, to meet the requirements for gravitational wave detection within targeted frequency bands, stringent constraints are imposed on the parameters governing the configuration of the Taiji formation [26].
• To mitigate significant shifts within the sensitive frequency band, it is essential to limit the arm length of the Taiji constellation within a range of 3 ± 0.035 million kilometers.
• The Taiji constellation does not need to keep a rigid geometry, but the change rate of arm-length is limited to under 10 m/s.
• Due to the design requirements of the optical system, the breathing angle of the Taiji constellation should be maintained at 60 ± 1 degrees.
To meet the aforementioned detection requirements, the configuration of the Taiji formation needs to maintain orbital stability.Gravitational perturbations from celestial bodies induce changes in spacecraft velocities, thereby altering their orbits.Over time, these variations may result in the armlength of spacecraft, breathing angle, and velocities exceeding the maximum acceptable limit.Once this threshold is surpassed, we will be unable to effectively probe the target frequency band for the gravitational wave detection of Taiji.Hence, assessing the impact of planetary gravitational perturbations on satellite formations is crucial for calculating experimental duration and facilitating subsequent orbit optimization efforts.
This study calculates the armlength, armlength variation rate, and breathing angle variation of the Taiji constellation under planetary gravitational perturbations.Employing ephemeris DE440 and based on the Keplerian orbit of Taiji, these variations under the influence of the gravities of Venus, Earth, and Jupiter were evaluated.Our analysis showed that the variations in the heliocentric formation configuration are primarily influenced by Earth's gravitational perturbations.Numerical simulations based on ephemeris DE440 show that the contribution of the Earth's gravity to the relative acceleration between spacecrafts is approximately O(10 −6 ) m • s −2 .Considering 00:00:00 on 23 October 2032 as the initial orbiting moment as an example, when the Earth's gravitational perturbation is determined, within six years the maximum amplitude of the relative variation in the armlength is 1.6×10 5 km, maximum amplitude of the relative variation in breathing angle is 3.2 • , and maximum amplitude of the relative variation in the armlength variation rates is 32 m • s −1 , compared with the unperturbed Keplerian configuration.By adding the gravitational perturbations of Venus and Jupiter, the armlength and relative velocity are reduced by 16.01% and 17.45%, respectively, compared with when only the Earth's gravitational perturbation is considered.Results showed that the maximum amplitude of Taiji constellation armlengths, armlength variation rate, and breathing angles varied with orbital entry moments.The maximum amplitude of the interspacecraft armlength change rate is the smallest when the initial orbital moment occurs in July.When the orbital entry moment is 00:00:00 on 1 July 2032, the maximum amplitude of the reduced armlengths is approximately 8 × 10 4 km, the maximum amplitude of the reduced breathing angles is approximately 1.6 • , and the maximum armlength variation rate is approximately 18.1 m • s −1 within six years.
We show that by using different ephemerides for the calculations, the constellation armlengths deviate in the order of O(1) m, constellation armlength variation rate deviates in the order of O(10 −7 ) m • s −1 , and relative accelerations deviate in the order of O(10 −13 ) m • s −1 within six years.The differences between ephemerides DE440 and DE430 are smaller than those between DE440 and DE421.
The remainder of the paper is presented as follows: In Section II, we derive the expressions for the positions and velocities of spacecrafts in the Keplerian configuration.The armlengths, armlength variation rate, inter-spacecraft relative acceleration, and breathing angles are determined.In Section III, we discuss the effect of planetary gravitational perturbations on the heliocentric formation configuration and the effect of different entering orbit moments.In Section IV, we compare the calculation results using different ephemerides.In Section V, we present the conclusions.

II. THE KEPLERIAN ORBITAL CONFIGURATION OF TAIJI A. The Keplerian Orbit
Several orbitals have been designed for the Taiji heliocentric formation configuration [8,9,27].We adopted the Keplerian orbital configuration model proposed in ref. [28] as part of the Taiji program preview, as shown in Figure 1.In the Keplerian configuration, the heliocentric reference system ECLIPJ2000 is chosen, the argument of periapsis is denoted as ω, and the longitude of the ascending node is denoted as Ω.When t = 0, spacecraft 1 (SC1) is at the highest point of the orbit and is located in the XZ plane.This can be associated with the determination of ω and Ω of SC1 as 270 • [29].As the three orbits have rotational symmetry about the Z-axis with a rotation angle of 120 • , the ω values of spacecraft 2 (SC2) and spacecraft 3 (SC3) are the same as that of SC1 with ω = 270 • .SC2 has Ω = 30 • , and SC3 has Ω = 150 • .Notably, the orbital semi-major axis of the three spacecrafts is a, eccentricity is e, and inclination is ϵ.The spacecraft position vectors and where E i is the eccentric anomaly, which can be obtained from the Kepler equation where ω = G(ms+m i ) is the average angular velocity of the spacecraft.
In the Taiji program, the half-length axis is a = 1 AU, the orbital eccentricity is e ≈ 5.789×10 −3 , and the orbital inclination ε is given by Equation (4) [28] cos ε = √ 3 3 where ϕ is the angle between the plane of spacecraft formation and ecliptic plane, and α is a small parameter for the expansion, e 2 + 2e + cos 2 ϕ − cos ϕ .The spacecraft position r i and velocity v i varied with time t, as expressed in the above equation.

B. Configuration Parameter
The armlength l ij , armlength variation rate v i , relative acceleration between the spacecrafts lij , and breathing angle β i of the constellation are as follows: where θ is the angle between vector ṙi − ṙj and vector r i − r j and i, j = 1, 2, 3.
By defining the reduced armlength as lij = l ij − 3 × 10 6 km and the reduced breathing angle as βi = β i −60 • , we obtain the variation of Taiji constellation reduced armlengths, armlength variation rate, and reduced breathing angles within six years, as shown in Figures 2 and 3.In the Keplerian orbital configuration, the reduced armlength amplitude is approximately 1 × 10 4 km, the rate of the armlength amplitude variation is 1.2 m/s, breathing angle amplitude is −0.24 • ∼0.24 • , and the relative acceleration amplitude is approximately 4.5 × 10 −7 m/s 2 .

III. TAIJI HELIOCENTRIC FORMATION CONFIGURATION A. Influence of the Gravitational Field of the Solar System Bodies
As we aim to evaluate the influence of the gravitational fields of celestial bodies on the configuration of spacecraft formation, we consider only the Newtonian gravitational forces of the Sun and other 9 bodies, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, and the Moon, in the solar system (k = 1, 2, ..., 9).The motion equation of the spacecraft in heliocentric ECLIPJ2000 can be expressed as where r i is the spacecraft position vector, µ s = GM s , µ k = GM k , s denotes the Sun, and G is the gravitational constant.
For this specific formation configuration, the initial mission time is 00:00:00 on 27 October 2032 (heliocentric ecliptic coordinate system), and the initial position and velocity are presented in Table I.The acceleration of each planet toward SC1 is We obtain the planet position from the DE440 ephemeris, substitute it into the Equation ( 6), and obtain the result as shown in Figure 4.The Figure shows that the gravitational attractions exerted by the Earth, Venus, and Jupiter on spacecrafts are greater than those of other celestial bodies in the solar system and may have a greater variation on the spacecraft formation configuration.The contribution of the gravitational disturbances of the Earth, Venus, and Jupiter to the relative acceleration between spacecrafts obtained by numerical simulation is shown in Figure 5.As shown in Figures 4 and 5, the gravitational perturbations of Venus, Earth, and Jupiter cause a variation in the spacecraft orbit.The maximum variation is of O(10 −7 ) m/s 2 .The Earth's gravitational perturbation has the greatest variation in the heliocentric formation configuration of O(10 −6 ) m/s 2 .Subsequently, we will examine variations in spacecraft formation configuration parameters perturbed solely by Earth's gravity.

B. Influence of the Gravitational Perturbation of the Earth
Considering the gravitational perturbation of the Earth, the motion equation of the spacecraft can be expressed as Substituting the initial conditions in Table I into Equation ( 8), the spacecraft position vector r i (t) can be obtained.Substituting r i (t) into Equation ( 5) and comparing it with the Keplerian orbits, the result is shown in Figure 6.Considering the Earth's gravitational perturbations, the orbits of spacecrafts in the heliocentric formation deviate from the Kepler orbit.This deviation changes the spacecraft formation configuration and increases the magnitude of the configuration parameters.Specifically, the maximum amplitude of the relative change in armlength between spacecrafts is 1.6 × 10 5 km, the maximum amplitude of the relative change in relative velocity between spacecrafts is 32 m • s −1 , and the maximum relative change in the breathing angle is 3.2 • within six years.

C. Influence of the Gravitational Perturbation of Venus and Jupiter
To obtain accurate configuration parameter variations, the influence of the gravities of Venus and Jupiter on the spacecraft formation configuration must be considered.Choosing k = 2, 3, 5 to represent Venus, Earth, and Jupiter in the equation, respectively, for a single spacecraft in the heliocentric formation configuration, the following acceleration can be obtained.
Substituting the initial conditions in Table .I into Equation ( 9), the spacecraft position vector r i (t) can be obtained.Substituting r i (t) into Equation ( 5) and comparing it with the perturbation of the spacecraft formation configuration parameters caused exclusively by the Earth's gravity, the results are shown in Figure 7. On can see that in addition to the gravity perturbation of the Earth, the configuration parameters will change when the gravity perturbation from Venus, Earth, and Jupiter is included.The maximum relative variations in inter-spacecraft armlength and relative velocity are 5×10 4 km and 9 m • s −1 , respectively, within six years.The maximum relative variation in breathing angle is approximately 0.8 • .However, the relative magnitude of the changes does not accurately reflect the influence of the perturbations of Venus and Jupiter.Figure 8 compares the changes in spacecraft formation configuration parameters in the three cases: the unperturbed Keplerian configuration, adding only the gravitational perturbation of the Earth, and adding the perturbations of Venus and Jupiter.As shown in Figure 8, when considering the gravitational disturbances of Venus and Jupiter in addition to the gravitational disturbances of the Earth, the reduced armlengths and armlength variation rate decrease by 16.01% and 17.45%, respectively.
The findings indicate that when the orbit entry time is set to 00:00:00 on 27 October 2032, the Earth's gravity amplifies the variations in the constellation armlength, change rate, and breathing angle.However, the gravitational perturbations of Venus and Jupiter reduce the magnitude of the constellation configuration parameters.

D. Orbital entry moment
The orbital entry moment is an optimization variation that affects the orbit design, and the different initial orbiting moments change the initial position and velocity of the spacecraft.In this section, different initial positions and velocities are obtained by changing the initial orbital entry moment of the formation.Further, the changes in configuration parameters at different orbital entry moments are analyzed.The specified orbiting time is set as 00:00:00 on the first day of each month in 2032.The simulation assesses the variations in the armlengths, breathing angles, and armlength variation rates over six years following the different orbiting times, which are shown in Figure 9. Notably, variations in the maximum amplitudes of the configuration parameters are observed at different orbital moments.Specifically, in July, the maximum amplitudes of the armlength and relative velocity between the spacecrafts are the smallest over one year.Consequently, selecting July as the moment of orbit entry is conducive to extending the mission duration.
Through sensitivity analysis of the initial orbital entry moment of the heliocentric formation, the orbital moment chosen was 00:00:00 on 1 July 2032.The initial positions and velocities of the spacecraft are presented in Table II.By substituting the parameters in Table II into Equation ( 9), the position vector r i (t) of the spacecraft can be determined.The vector in Equation ( 1) is used to calculate the variation in spacecraft formation configuration parameters, as shown in Figure 10.The findings indicate that selecting the orbital entry moment as 00:00:00 on 1 July 2032, while TABLE II.Initial spacecraft position and velocity (00:00:00 on 1 July 2032).accounting for the gravitational perturbation of Venus, Earth, and Jupiter, results in a maximum amplitude of approximately 8 × 10 4 km in the reduced armlengths, a maximum armlength variation rates of 18.1 m • s −1 , and a maximum amplitude of approximately 1.6 • in the reduced breathing angles over six years.

IV. SIMULATION WITH DIFFERENT EPHEMERIS
The accuracy of the numerical simulation results depends on that of the planetary position data.Differences in planetary orbital positions may occur in various ephemerides due to variations in the chosen fitting data and the integration method used [30].These deviations can affect the numerical simulation of spacecraft configuration parameters.Therefore, the biases in numerical simulations caused by different ephemerides as input should be extensively investigated.
We obtained planetary position data from the ephemerides DE440, DE430 and DE421 [31][32][33].In the heliocentric equatorial coordinate system J2000, the maximum disparities in the coordinates of Venus, Earth, and Jupiter relative to the solar center of mass are presented in Table III across various ephemerides.Comparatively, the maximum coordinates of the ephemerides of the three planets differ significantly: Venus and Earth differ at the 100 m level, while Jupiter differs at the 10 km level.
Discrepancies in the planetary coordinates between ephemerides can cause variations in the forces acting on the spacecrafts.Consequently, these variations may alter the relative acceleration between the spacecrafts, thereby affecting the overall configuration of the formation.Three ephemerides, DE421, DE430, and DE440, were used for numerical simulations with an initial entry time of 00:00:00 on 1 July 2032.Figure 11 shows the difference in spacecraft armlengths, armlength FIG. 10.The variations of the reduced armlength, reduced breathing angle, and the inter-spacecraft relative velocity when the orbital insertion time is 00:00:00 on 1 July 2032.
variation rates, and relative accelerations within six years.
Figure 11 shows differences in the spacecraft configuration parameters obtained via numerical simulation using different ephemerides.Specifically, the variation between DE440 and DE430 is relatively small, with a difference of approximately 0.5 m in armlength over six years, 1 × 10 −7 m/s in armlength variation rates, and 3 × 10 −14 m/s 2 in relative acceleration.Conversely, the variation between DE421 and DE440 is significant, with a difference of approximately 2.5 m in armlength, 5 × 10 −7 m/s in armlength variation rates, and 1 × 10 −13 m/s 2 in relative acceleration within six years.
Notably, the variable period of the difference in spacecraft structural parameters obtained via the numerical simulation of different ephemerides is approximately one year, and the frequency is relatively low.This is relatively small compared to the change in the gravitational wave frequency band targeted by the Taiji project.

V. CONCLUSIONS
Using the DE series of ephemerides, we investigated the influence of celestial gravitational perturbations on the heliocentric formation configuration.Specifically, variations in the major planets in the solar system on the orbital stability of the Taiji space gravitational wave detector were investigated.Ephemeris DE440 was used to interpolate the planetary orbital position during the mission, and the high-order Runge-Kutta numerical integration method was used to calculate planetary gravity.The trajectory of each spacecraft in the heliocentric formation configuration under the influence of perturbation was numerically simulated.In our investigation, to streamline the model, solar system bodies are regarded as point masses, with ancillary factors such as their geometric characteristics and rotational dynamics being omitted.These aspects will be meticulously scrutinized in forthcoming research.
Results showed that the gravitational acceleration affecting an individual spacecraft within a spacecraft formation due to planetary gravity is approximately 10 −7 m • s −2 .The relative acceleration between spacecrafts due to planetary gravity is approximately 10 −6 m • s −2 .The difference in the heliocentric formation configuration is mainly affected by the gravity of the Earth.Considering 00:00:00 on 27 October 2032 as the initial orbit entry moment as an example, when only the Earth's gravitational perturbation is considered, the maximum amplitude of the relative change in the armlength within six years was 1.6 × 10 5 km, the maximum amplitude of the relative change in breathing angle was 3.2 • , and the maximum amplitude of the relative change in armlength variation rates was 32 m • s −1 .Further, we discuss the spacecraft formation configuration changes with the addition of gravitational perturbations from Venus and Jupiter.The results showed that compared with the case where only the Earth's gravitational perturbation was considered, after the gravitational perturbations of Venus and Jupiter were added, within six years the maximum amplitude of the relative change in armlength was approximately 5 × 10 4 km, maximum amplitude of the breathing angle was approximately 0.8 • , and maximum amplitude of the relative change in armlength variation rates was 9 m/s.After adding the gravitational perturbations of Venus and Jupiter, the armlength and relative velocity were reduced by 16.01% and 17.45%, respectively, compared with when only the Earth's gravitational perturbation was considered.
Variations in the maximum amplitudes of the formation configuration parameters were observed when entering orbit at different times.The smallest increase in armlength variation rates between the spacecrafts occurred when the initial orbiting time was in July.Consequently, selecting July as the initial orbital insertion time extends the experiment duration.By considering 00:00:00 on 1 July 2032 as the initial orbit entry moment, under the gravitational perturbations of Venus, Earth, and Jupiter, using the DE440 ephemeris, the maximum amplitude of the reduced armlength was 8 × 10 4 km, the maximum amplitude of the reduced breathing angle was 1.6 • , and the maximum armlength variation rates was 18.1 m • s −1 within six years.
Ephemeris data were used to perform numerical simulations.By comparison, the maximum disparity in the orbital positions of Venus and Jupiter was 300 m and 60 km, respectively.The armlength deviation over the six-year simulation period was approximately 1 m.The deviation in armlength variation rates was in the order of O(10 −7 ) m • s −1 , and the deviation in relative accelerations was approximately 1×10 −13 m • s −2 .The differences between ephemerides DE440 and DE430 are smaller than those between DE440 and DE421.Nonetheless, the impact of discrepancies in high-frequency bands induced by factors such as planetary rotation remains unexplored and warrants consideration in future studies.

FIG. 2 .
FIG. 2. Armlength, rate of change in armlength, and relative acceleration in the Keplerian orbital configuration.

FIG. 4 .FIG. 5 .
FIG. 4. Acceleration of SC1 by the planets and the Moon in the solar system.

FIG. 6 .
FIG.6.Considering the gravitational perturbation of the Earth, the variations in configuration parameters relative to the Kepler orbit within six years: (top) variation in the armlength, (middle) variation in the breathing angle, (bottom) variation in the armlength variation rate, when the orbital insertion time is 00:00:00 on 27 October 2032.

FIG. 7 .
FIG. 7. Considering the gravitational perturbations of Venus, Earth and Jupiter, the configuration parameters variations in 6 years relative to the Earth-only perturbations: (top) change in armlength; (middle) change in breathing angle; and (bottom) change in armlength variation rate, when the orbital insertion time is 00:00:00 on 27 October 2032.

FIG. 9 .
FIG.9.From top to bottom, the three rows represent the maximum amplitude of the reduced armlength, the maximum amplitude of the reduced breathing angle, and the maximum amplitude of the inter-spacecraft relative velocity.