Time-Delay Interferometry for ASTROD-GW

{CJK}

UTF8gbsn

\englishtitle

Time-Delay Interferometry for ASTROD-GW \englishauthorGang Wang \englishadvisorProf. Wei-Tou Ni \englishinstitutePurple Mountain Observatory
Chinese Academy of Sciences \englishdegreeMaster \englishmajorAstrophysics-Gravitational Waves Detection

\makeenglishtitle

{englishabstract}

In the detection of gravitational waves in space, the arm lengths between spacecraft are not equal due to their orbital motion. Consequently, the equal arm length Michelson interferometer used in Earth laboratories is not suitable for space. To achieve the necessary sensitivity for space gravitational wave detectors, laser frequency noise must be suppressed below secondary noise sources such as optical path noise and acceleration noise. To suppress laser frequency noise, time-delay interferometry (TDI) is employed to match the two optical paths and retain gravitational wave signals. Since planets and other solar system bodies perturb the orbits of spacecraft and affect TDI performance, we simulate the time delay numerically using the CGC2.7 ephemeris framework. To examine the feasibility of TDI for the ASTROD-GW mission, we devised a set of 10-year and a set of 20-year optimized mission orbits for the three spacecraft starting on June 21, 2028, and calculated the path mismatches in the first- and second-generation TDI channels. The results demonstrate that all second-generation TDI channels meet the ASTROD-GW requirements. A geometric approach is used in the analysis and synthesis of both first-generation and second-generation TDI to clearly illustrate the construction process.

\englishkeywords

ASTROD, Time-Delay Interferometry, Gravitational Waves Detection, Numerical Calculation, Geometric Construction

Chapter 1 Introduction

¹¹1This chapter is derived from review [1, 2].

In 1905, Poincaré [3] and Einstein [4] proposed the theory of special relativity. Poincaré [3] attempted to establish a relativistic theory of gravity, mentioned gravitational waves (GWs), and inferred their propagation speed to be the same as that of light based on Lorentz invariance. Subsequently, physicists (including Einstein himself) attempted to establish a relativistic theory of gravity [5, 6]. It was not until 1915 that Einstein proposed the general theory of relativity, which successfully explained the anomalous precession of Mercury’s perihelion [7]. After proposing the general theory of relativity in 1915, Einstein predicted the existence of GWs and estimated their strength [8]. To this day, Einstein’s general theory of relativity has become the standard theory of gravity, widely applied in the global positioning system, planetary and lunar ephemeris calculations, solar system space navigation and exploration, astrophysics, and cosmology.

Maxwell’s electromagnetic theory predicts the existence of electromagnetic waves. Einstein’s general theory of relativity and other relativistic theories of gravity predict the existence of GWs. GWs propagate through spacetime, forming ripples. The status of GWs in gravitational physics is analogous to that of electromagnetic waves in electromagnetic physics. The existence of GWs is a direct consequence of general relativity and is an inevitable result of all gravitational theories with a finite propagation speed.

The importance of GW detection is twofold:

: (i) To explore fundamental physics and cosmology, especially black hole physics and the early universe;
: (ii) To serve as a tool for astronomical and astrophysical research, studying compact celestial bodies and calculating their density distribution, complementing electromagnetic wave astronomy and cosmic ray (including neutrino) observations.

The evolution of the binary pulsar orbit demonstrates the existence of GW radiation [9]. In general relativity, a binary system emits energy in the form of GWs. The energy loss leads to a shrinking orbit and a shorter orbital period. Thirty years of observations of the relativistic binary B1913+16 show a cumulative advance in periastron time by 35 s. In relativity, the orbital decay rate of a binary pulsar can be calculated from the pulsar system parameters determined by pulsar timing observations. After correcting for the relative acceleration between the solar system and the pulsar binary system, Weisberg and Taylor [9] found that the measured orbital decay rate matches the predicted rate due to GW radiation from general relativity to within $(1.3\pm 2.1)\times 10^{-3}$ . Hulse and Taylor, who discovered this binary system, were awarded the 1993 Nobel Prize in Physics.

Similar to how electromagnetic waves are divided into radio waves, microwaves, infrared rays, light waves, ultraviolet rays, X-rays, and gamma rays based on frequency, GWs can also be categorized into different frequency bands [1, 2, 10, 11, 12, 13]:

: (i) Ultra-high frequency band ( $>$ 1 THz): Detection methods include terahertz resonant cavities, optical resonant cavities, and innovative methods yet to be invented.
: (ii) Very high frequency band (100 kHz - 1 THz): This is the most sensitive frequency band for laboratory detection of GWs using microwave resonant systems and short-arm length laser interferometers.
: (iii) High frequency band (10 Hz - 100 kHz): This is the most sensitive frequency band for ground-based detection of GWs using cryogenic resonators and laser interferometers.
: (iv) Mid frequency band (0.1 Hz - 10 Hz): This is the most sensitive frequency band for space-based laser interferometers with short arm lengths ( $10^{3}-10^{5}$ km).
: (v) Low frequency band (100 nHz - 0.1 Hz): This is the most sensitive frequency band for space-based laser interferometers with long arm lengths ( $10^{6}-10^{9}$ km).
: (vi) Very low frequency band (300 pHz - 100 nHz): This is the most sensitive frequency band for pulsar timing experiments.
: (vii) Ultra-low (sub-Hubble) frequency band (10 fHz - 300 pHz): This is the frequency band between the very low frequency and extremely low frequency bands, and is most sensitive for precise measurements of quasar and radio source proper motion experiments.
: (viii) Hubble (extremely low) frequency band (1 aHz - 10 fHz): This is the most sensitive frequency band for cosmic background radiation anisotropy and polarization experiments.
: (ix) (Infrared) Hubble frequency band ( $<$ 1 aHz): Inflationary cosmological models predict GWs in this band. Confirmation of these inflationary models can indirectly provide evidence for GWs in this band.

The main activities for high-frequency GW detection are in ground-based long-arm laser interferometers. The TAMA 300 m arm length interferometer[14], GEO 600 m arm length interferometer[15], and kilometer-scale laser interferometer GW detectors LIGO[16] (two 4 km arm lengths, one 2 km arm length) and Virgo[17] have essentially achieved their original design sensitivity goals. Around 100 Hz, the sensitivity of LIGO and Virgo has reached $10^{-23}/\sqrt{\text{Hz}}$ . Currently, both LIGO and Virgo are undergoing upgrade plans for the next generation detectors — AdLIGO[18] and AdVirgo[19], which will improve sensitivity by ten times, increasing the number of detectable GW sources by about 1000 times. The 3 km cryogenic laser interferometer GW detector LCGT has begun construction, with the first phase being a room temperature detector, which will be converted directly to a cryogenic third-generation long-arm interferometer after completion, commissioning, and observation[20], with sensitivity comparable to AdLIGO and AdVirgo. The European third-generation long-arm interferometer ET has begun planning[21]. It is expected that within about five years, humans will be able to directly detect GWs for the first time.

Space-based GW detection laser interferometers (LISA[22], ASTROD[23, 24], ASTROD-GW[1, 11, 13], Super-ASTROD[25], DECIGO[26], and Big Bang Observer (BBO)[27, 28]) offer the high signal-to-noise ratio, and are crucial for studying astrophysics, cosmology, and fundamental physics. We review space-based GW detectors in the next chapter. Achieving the required target sensitivity for laser interferometry space GW detectors necessitates reducing laser frequency noise. If the time delay matching of the two beams is achieved, their interference signal laser frequency noise can be subtracted, achieving the goal. The better the matching (smaller time delay difference), the better the noise reduction. In the third chapter, we discuss the principles and methods of time-delay interferometry (TDI). Designing appropriate mission orbits and establishing a suitable ephemeris is necessary. In the fourth chapter, we establish the CGC2.7 ephemeris framework. In the fifth chapter, we further discuss and select the mission orbits for ASTROD-GW based on previous work. In the sixth chapter, we introduce numerical methods for TDI and calculate the path mismatches in TDI for dynamic orbit. In the following seventh chapter, we analyze and construct various interferometry paths when the three interferometer arms work simultaneously and calculate the time delays for some of these paths for ASTROD-GW. Finally, in the eighth chapter, we provide conclusions and a brief discussion.

Chapter 2 Space-borne Gravitational Wave Detectors

¹¹1This chapter is derived from references [1, 2, 13].

The gravitational field of the solar system is determined by three factors: the dynamic distribution of matter within the solar system, the dynamic distribution of matter outside the solar system (Milky Way, extragalactic systems, universe, etc.), and GWs passing through the solar system. Different relativistic gravitational theories predict different gravitational fields for the solar system; thus, precise measurements of the solar system’s gravitational field can test these relativistic gravitational theories. Additionally, these measurements can detect GWs, determine the distribution of matter in the solar system, and observe the measurable (testable) effects of the Milky Way and the universe. We can determine the gravitational field of the solar system by measuring/monitoring various natural and artificial celestial bodies. In the solar system, the motion of celestial bodies or spacecraft follows the following astronomical dynamics equation:

\mathbf{a}=\mathbf{a}_{\text{N}}+\mathbf{a}_{\text{1PN}}+\mathbf{a}_{\text{2PN% }}+\mathbf{a}_{\text{Gal-Cosm}}+\mathbf{a}_{\text{GW}}+\mathbf{a}_{\text{non-% grav}}

(2.1)

where $\mathbf{a}$ is the acceleration of the celestial body or spacecraft, $\mathbf{a}_{\text{N}}$ is the acceleration due to the Newtonian gravitational theory from the mass distribution in the solar system, $\mathbf{a}_{\text{1PN}}$ is the first post-Newtonian correction acceleration, $\mathbf{a}_{\text{2PN}}$ is the second post-Newtonian correction acceleration, $\mathbf{a}_{\text{Gal-Cosm}}$ is the acceleration due to the mass distribution in the Milky Way and the universe, $\mathbf{a}_{\text{GW}}$ is the acceleration due to GWs, and $\mathbf{a}_{\text{non-grav}}$ is the acceleration from all non-gravitational sources. The distance between spacecraft/celestial bodies is determined by the gravitational field of the solar system (including the gravitational effects of solar oscillations), the fundamental gravitational theories, and the GWs passing through the solar system. Accurately measuring these time-varying distances can identify the causes of these variations. Some orbital combinations are better for testing relativistic gravity; some are easier for measuring solar system parameters; some are easier for detecting GWs. These are all integral parts of mission design.

Space-borne GW detectors mostly use drag-free spacecraft, so in using Equ. (2.1), $\mathbf{a}_{\text{non-grav}}$ is considered noise, and generally, local gravitational variations are also considered part of this noise.

2.1 LISA

LISA[22] (Laser Interferometer Space Antenna) has an interferometer arm length of approximately 5 million kilometers. Its goal is to detect GWs in the $10^{-4}$ to 1 Hz frequency band, primarily in the low-frequency region, with some coverage in the mid-frequency region. Its strain detection sensitivity at 1 mHz is $4\times 10^{-21}/\sqrt{\text{Hz}}$ . LISA, ASTROD, and ASTROD-GW have a wealth of GW sources: galactic compact binaries (neutron stars, white dwarfs, etc.) and extragalactic sources. Extragalactic targets include supermassive black hole binaries, the formation of supermassive black holes, and the cosmic GW background. The LISA space mission is hoped to be launched by 2020.

Refer to caption — Figure 2.1: Orbit of the LISA mission [Credit: ESA]

2.2 ASTROD

The general concept of ASTROD (Astrodynamical Space Test of Relativity Using Optical Devices) is to use drag-free spacecraft flying in formation within the solar system, using optical ranging between them, to map the gravitational field of the solar system, measure related solar system parameters, test relativistic gravity, observe solar g-mode oscillations, and detect GWs. The baseline plan for ASTROD was proposed in 1993, with conceptual and laboratory studies initiated simultaneously.

2.3 ASTROD-GW Space Gravitational Wave Detection Program

Considering the need to optimize ASTROD missions for GW detection, the spacecraft forming the space interferometer implement nearly equal arm lengths. The three spacecraft are located near the Sun-Earth Lagrangian points L3, L4, and L5, forming an almost equilateral triangular array, as shown in Figure 2.2, with an arm length of about 260 million kilometers (1.732 astronomical units). The three spacecraft conduct laser interferometric ranging between each other.

For ASTROD-GW, which focuses on GW detection, the orbit design scheme can be optimized. For the Earth-Sun Lagrangian points L3, L4, and L5, L4 and L5 are stable, while L3 is unstable, but with an instability timescale of about 50 years, making it effectively a stable point for missions lasting 10-20 years. With this choice, the spacecraft are about 260 million kilometers apart, 52 times the arm length of LISA, allowing detection of GWs at frequencies 52 times lower than LISA. Considering ASTROD-GW will be after LISA, the requirement for acceleration noise is assumed to be as same as LISA’s requirements, and for the Doppler shift between the three spacecraft, the requirements are less stringent than for LISA, allowing the use of LISA’s Doppler frequency synthesis and related technologies.

2.4 Super-ASTROD, DECIGO, BBO Space Gravitational Wave Detection Programs

The Super-ASTROD [25] mission concept involves using 3-5 spacecraft in 5 AU orbits and one spacecraft at the Sun-Earth Lagrangian point L1/L2. These spacecraft will use optical ranging to detect primordial GWs in the 100 nHz - 1 mHz frequency band, test fundamental of space-time, and map the mass distribution and dynamics of the outer solar system.

DECIGO [26] (DECi-hertz Interferometer Gravitational Wave Observatory) is a future Japanese space GW detector. Its goal is to detect various GWs in the 1 mHz - 100 Hz range, opening a new observation window for GW astronomy, primarily studying primordial GWs from the Big Bang. Its concept involves using 3 drag-free spacecraft separated by 1000 km, using Fabry-Perot Michelson interferometers to measure the relative displacement between the spacecraft.

The BBO [27, 28] (Big Bang Observer) has an orbit similar to LISA, but BBO’s arm length is 2-5 $\times 10^{4}$ km. As a follow-up to LISA, BBO aims to detect GWs in the 0.01 - 10 Hz band, filling the gap between ground and space GW detectors. Its primary goal is to study primordial GWs from the Big Bang.

Chapter 3 Principle of Time-Delay Interferometry

To achieve the required sensitivity for laser interferometer space GW detectors, it is essential to reduce laser frequency noise. The closer the time delays of the two interfering laser beams match, the better the laser frequency noise can be reduced, and the closer the sensitivity can approach the target. In this section, we discuss the principles and methods of TDI.

ASTROD employs TDI in its interferometric scheme. In space interferometers, due to long distances, the laser reaching another spacecraft must be amplified before it can be transmitted back and to other spacecraft. The method of amplification involves phase-locking the local laser with the incoming weak light before transmitting it. In the development of phase-locking with weak light, National Tsing Hua University (Hsinchu) first achieved phase-locking with 2 pW of incoming light in 2000 [29]. The Jet Propulsion Laboratory (JPL) at Caltech further improved this to 40 fW of incoming light in 2008 [30]. In 1996, Wei-Tou Ni et al. [31, 32] proposed the following two TDI schemes:

: (i) Path 1: S/C 3 $\rightarrow$ S/C 1 $\rightarrow$ S/C 3 $\rightarrow$ S/C 2 $\rightarrow$ S/C 3
Path 2: S/C 3 $\rightarrow$ S/C 2 $\rightarrow$ S/C 3 $\rightarrow$ S/C 1 $\rightarrow$ S/C 3
After phase-locking and amplification, the laser beams following Path 1 and Path 2 interfere at S/C 3. If the optical paths of Path 1 and Path 2 are equal (i.e., the distances between spacecraft are constant), then the phase of the two beams will be the same at the start and after traveling equal paths, thus canceling out the laser noise. If GWs, Lense-Thirring effects, or other factors cause the optical paths to differ, the difference is relatively small. Due to the variation in spacecraft distances over time, the corresponding noise from these changes must be considered when stabilizing the laser frequency.
: (ii) Path 1: S/C 3 $\rightarrow$ S/C 1 $\rightarrow$ S/C 2 $\rightarrow$ S/C 3
Path 2: S/C 3 $\rightarrow$ S/C 2 $\rightarrow$ S/C 1 $\rightarrow$ S/C 3
After phase-locking and amplification, the laser beams following Path 1 and Path 2 interfere at S/C 3.

These two schemes were later referred to as first-generation TDI for LISA mission, with second-generation schemes providing better cancellation [33]. In the following, we provide a detailed explanation of first- and second-generation TDI as discussed in [33].

In a traditional Michelson interferometer, the two interferometric arms can be made precisely equal, ensuring that laser passing through both arms have the same time delay, effectively canceling out the laser phase noise. However, in space-based GW missions, due to the orbital dynamics of the spacecraft, the distances between the three spacecraft are continuously changing, preventing the formation of an equal-arm space interferometer. Consequently, laser frequency noise cannot be canceled out at the traditional receiver, necessitating the development of methods to handle unequal-arm space interferometers and remove laser noise from the signal.

For an unequal-arm Michelson interferometer, as shown in Fig. 3.1, the two interferometric arms have different lengths, $L_{1}$ and $L_{2}$ . The light beams passing through these arms do not interfere directly at a single photodetector but rather interfere with their respective outgoing beams upon return. The Doppler measurement results at their respective photodetectors are denoted as $y_{1}(t)$ and $y_{2}(t)$ , with $C(t)$ representing the laser frequency noise. The GW signals entering each Doppler measurement are denoted as $h_{1}(t)$ and $h_{2}(t)$ , and the remaining noise components as $n_{1}(t)$ and $n_{2}(t)$ . The measurements $y_{1}(t)$ and $y_{2}(t)$ can be expressed as:

y_{1}(t)=C(t-2L_{1})-C(t)+h_{1}(t)+n_{1}(t),

(3.1)

y_{2}(t)=C(t-2L_{2})-C(t)+h_{2}(t)+n_{2}(t).

(3.2)

From equations (3.1) and (3.2), it is important to note the time-varying stochastic process $C(t)$ in the Doppler measurements $y_{1}$ and $y_{2}$ . The laser signals entering the photodetector at arm 1 at time $t$ include noise $C(t-2L_{1})$ from $2L_{1}$ earlier (assuming $c=1$ ), from which the current noise $C(t)$ must be subtracted, as shown in Eq. (3.1). A similar result can be obtained for arm 2.

By comparing the difference between the measurements in Eqs. (3.1) and (3.2), we obtain:

y_{1}(t)-y_{2}(t)=C(t-2L_{1})-C(t-2L_{2})+h_{1}(t)-h_{2}(t)+n_{1}(t)-n_{2}(t).

(3.3)

Considering how laser noise enters Eq. (3.3) rather than Eqs. (3.1) and (3.2), a time-shifting $y_{1}(t)$ could be implemented before the laser traverses arm 2, i.e., $y_{1}(t-2L_{2})$ , and shifting $y_{2}(t)$ to before the laser traverses arm 1, i.e., $y_{2}(t-2L_{1})$ , we can derive:

y_{1}(t-2L_{2})=C(t-2L_{1}-2L_{2})-C(t-2L_{2})+h_{1}(t-2L_{2})+n_{1}(t-2L_{2}),

(3.4)

y_{2}(t-2L_{1})=C(t-2L_{2}-2L_{1})-C(t-2L_{1})+h_{2}(t-2L_{1})+n_{2}(t-2L_{1}).

(3.5)

Comparing the difference between Eqs. (3.4) and (3.5), we obtain:

\begin{split}y_{1}(t-2L_{2})-y_{2}(t-2L_{1})=&C(t-2L_{1})-C(t-2L_{2})+h_{1}(t-% 2L_{2})-h_{2}(t-2L_{1})\\ &+n_{1}(t-2L_{2})-n_{2}(t-2L_{1}).\end{split}

(3.6)

Comparing Eqs. (3.3) and (3.6), we can see that both contain the same laser frequency noise. This implies that by subtracting Eq. (3.6) from Eq. (3.3), we obtain new data $X$ that is free from laser frequency noise $C(t)$ :

\begin{split}X&\equiv[y_{1}(t)-y_{2}(t)]-[y_{1}(t-2L_{2})-y_{2}(t-2L_{1})]\\ &=h_{1}(t)-h_{2}(t)+n_{1}(t)-n_{2}(t)\\ &-[h_{1}(t-2L_{2})-h_{2}(t-2L_{1})+n_{1}(t-2L_{2})-n_{2}(t-2L_{1})].\end{split}

(3.7)

From the expression of $X$ , we can see that by performing appropriate time delay operations in time domain and combining different Doppler measurements, it becomes feasible to eliminate the laser frequency noise. This is the essence of TDI.

Chapter 4 CGC2.7 ephemeris framework

For the purposes of orbital design and numerical calculations for TDI, having a planetary ephemeris with sufficient accuracy is crucial. Therefore, it is necessary to introduce the ephemeris framework we used, CGC2.7 (CGC: Center for Gravitation and Cosmology). The early ephemeris framework CGC1.0 was established by Dah-Wei Chiou and Wei-Tou Ni [36, 37], and the subsequent improved CGC2.0 was developed by Chien-Jen Tang and Wei-Tou Ni [38, 39]. Our ephemeris framework, based on CGC2.0, is numbered as CGC2.7 ephemeris framework. The main interactions considered in the CGC2.7 ephemeris framework include:

•

Newtonian and post-Newtonian interactions among major celestial bodies (the Sun, the nine planets, the Moon, Ceres, Pallas, and Vesta);
•

The effects of the Sun’s second zonal harmonic and Earth’s second to fourth zonal harmonics on other celestial bodies and on themselves;
•

Newtonian perturbations of 349 asteroids on major celestial bodies.

Additionally, we provide a brief introduction to the numerical integration algorithm of the ephemeris framework and compare the accuracy of the ephemeris of JPL’s DE405 in this section.

4.1 Newtonian and Post-Newtonian Interactions Between Celestial Bodies

Based on the relativistic gravitational theory with two PPN (Parameterized Post-Newtonian) parameters, $\gamma$ and $\beta$ , Brumberg[34] derived the equations for Newtonian and post-Newtonian acceleration corrections that need to be considered for a celestial body $i$ under the influence of other celestial bodies in the adopted barycentric coordinate system, as shown in Eqs. (4.1-4.3).

\ddot{\mathbf{r}}_{i}=-\displaystyle\sum_{j\neq{i}}\frac{GM_{j}}{r_{ij}^{3}}% \mathbf{r}_{ij}+\displaystyle\sum_{j\neq{i}}m_{j}(A_{ij}\mathbf{r}_{ij}+B_{ij}% \dot{\mathbf{r}}_{ij})

(4.1)

\begin{split}A_{ij}=&\frac{{\dot{\mathbf{r}}_{i}}^{2}}{r^{3}_{ij}}-(\gamma+1)% \frac{\dot{\mathbf{r}}^{2}_{ij}}{r_{ij}^{3}}+\frac{3}{2r_{ij}^{5}}(\mathbf{r}_% {ij}\dot{\mathbf{r}}_{j})^{2}+G[(2\gamma+2\beta+1)M_{i}+(2\gamma+2\beta)M_{j}]% \frac{1}{r_{ij}^{4}}\\ &+\displaystyle\sum_{k\neq{i,j}}GM_{k}[(2\gamma+2\beta)\frac{1}{r_{ij}^{3}r_{% ik}}+(2\beta-1)\frac{1}{r_{ij}^{3}r_{jk}}+\frac{2(\gamma+1)}{r_{ij}r_{jk}^{3}}% \\ &-(2\gamma+\frac{3}{2})\frac{1}{r_{ik}r_{jk}^{3}}-\frac{1}{2r_{jk}^{3}}\frac{% \mathbf{r}_{ij}\mathbf{r}_{ik}}{r_{ij}^{3}}]\end{split}

(4.2)

B_{ij}=\frac{1}{r_{ij}^{3}}[(2\gamma+2)(\mathbf{r}_{ij}\dot{\mathbf{r}}_{ij})+% \mathbf{r}_{ij}\dot{\mathbf{r}}_{j}]

(4.3)

In Eq. (4.1), the first term on the right-hand side represents the Newtonian interaction, while the second term accounts for the post-Newtonian interaction. Here, $\mathbf{r}_{i},\dot{\mathbf{r}}_{i},\ddot{\mathbf{r}}_{i}$ denote the position, velocity, and acceleration vectors of celestial body $i$ in the barycentric coordinate system of the solar system, respectively. $G$ is the gravitational constant, $M_{i}$ is the mass of celestial body $i$ , and $m_{i}=GM_{i}/c^{2}$ ; the vector $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ and $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ represent the position vector and distance between the centers of mass of celestial bodies $i$ and $j$ , respectively; $c$ is the speed of light. For Einstein’s General Relativity, the PPN parameters $\gamma=\beta=1$ . With $\alpha=0$ in Eqs. (4.2) and (4.3), we obtain the relativistic post-Newtonian acceleration correction.

4.2 Effects from Extended Bodies

For non-spherical extended celestial bodies, their interaction with other bodies, which are treated as point masses, is handled in the $O-\xi\eta\zeta$ coordinate system. The origin is placed at the center of mass of the non-spherical perturbing body, with the body’s rotation axis as the polar axis and the equatorial plane as the coordinate plane. The $\xi$ axis points from the origin to the point mass, and the perturbing body’s rotation axis lies in the $\xi\zeta$ plane. The $\eta$ axis is determined by the right-hand rule.

For any point K outside the perturbing body, its coordinates are expressed in spherical coordinates as $(r,\lambda,\theta)$ . Let $\phi=\frac{\pi}{2}-\theta$ , and the unit coordinate vector at point K is given by:

\mathbf{u}=\frac{\mathbf{r}}{r}=\begin{pmatrix}\cos\phi\cos\lambda\\ \cos\phi\sin\lambda\\ \sin\phi\end{pmatrix}.

The effect of the shape of a celestial body on the acceleration of a point mass is given by (see, e.g., Moyer, 1971),

\begin{split}\begin{pmatrix}\ddot{\xi}_{1}\\ \ddot{\eta}_{1}\\ \ddot{\zeta}_{1}\end{pmatrix}=&\frac{GM}{r^{2}}\left\{\displaystyle\sum_{i=2}^% {\infty}J_{i}(\frac{a}{r})^{i}\begin{pmatrix}(i+1)P_{i}(\sin\phi)\\ 0\\ -\cos\phi P^{{}^{\prime}}_{i}(\sin\phi)\end{pmatrix}\right\}\\ &+\frac{GM}{r^{2}}\left\{\sum_{i=2}^{\infty}(\frac{a}{r})^{i}\sum_{j=1}^{i}% \begin{pmatrix}-(i+1)P^{j}_{i}(\sin\phi)(C_{ij}\cos{j\lambda}+S_{ij}\sin{j% \lambda})\\ j\sec\phi P^{j}_{i}(\sin\phi)(-C_{ij}\sin{j\lambda}+S_{ij}\cos{j\lambda})\\ \cos\phi P^{\prime j}_{i}(\sin\phi)(C_{ij}\cos{j\lambda}+S_{ij}\sin{j\lambda})% \end{pmatrix}\right\}\end{split}

where $G$ is the gravitational constant, $M$ is the mass of the non-spherical perturbing body, $a$ is the radius of the perturbing body, $r$ is the distance between the center of mass of the non-spherical perturbing body and the point mass, $J_{i}$ are the zonal harmonic coefficients of the body, $P_{i}(\sin\phi)$ are the Legendre polynomials of degree $i$ , $P^{j}_{i}(\sin\phi)$ are the $j$ th power of $P_{i}(\sin\phi)$ , and $C_{ij}$ and $S_{ij}$ are the spherical harmonic coefficients of the perturbing body[35].

Let $\mathbf{p}=\cos\phi\mathbf{\Phi}^{0}+\sin\phi\mathbf{u}$ be the unit vector in the direction of the polar axis, where $\mathbf{\Phi}^{0}$ is the unit gradient vector of $\phi$ . Thus, the expression for the effect of the 2-4 order zonal harmonics of the non-spherical perturbing body on the acceleration of the point mass can be respectively obtained:

\ddot{\mathbf{r}}_{pm}=\frac{GMa^{2}J_{2}}{r^{4}}[(7.5\sin^{2}\phi-1.5)\mathbf% {u}-3\sin\phi\mathbf{p}],

(4.4)

\ddot{\mathbf{r}}_{pm}=\frac{GMa^{3}J_{3}}{r^{5}}[(17.5\sin^{3}\phi-7.5\sin% \phi)\mathbf{u}-(7.5\sin^{2}\phi-1.5)\mathbf{p}],

(4.5)

\ddot{\mathbf{r}}_{pm}=\frac{GMa^{4}J_{4}}{r^{6}}[(39.375\sin^{4}\phi-26.25% \sin^{2}\phi+3)\mathbf{u}-(17.5\sin^{3}\phi-7.5\sin\phi)\mathbf{p}],

(4.6)

where $J_{i}$ represents the $i$ th zonal harmonic coefficient of the celestial body.

In this ephemeris framework, zonal harmonic terms of celestial bodies consider interactions between extended body and a point mass, including the following cases:

•

Solar second-order zonal harmonics interacting with point masses of other celestial bodies. According to the formula for second-order zonal harmonics, the solar second-order zonal harmonic coefficient $J_{2}=2\times 10^{-7}$ , and the solar radius is $a=696000$ km.
•

Earth’s 2nd to 4th order zonal harmonics interacting with point masses of other celestial bodies. The Earth’s zonal harmonic coefficients are $J_{2}=0.1082626\times 10^{-2}$ , $J_{3}=-0.2533\times 10^{-5}$ , and $J_{4}=-0.1616\times 10^{-5}$ , with the Earth’s equatorial radius $a=6378.137$ km. To compute Earth’s polar direction, precession and nutation of the Earth need to be considered.

The unit rotation matrices are defined as:

\begin{split}\mathbf{R}_{x}(\theta)=&\begin{pmatrix}1&0&0\\ 0&\cos\theta&\sin\theta\\ 0&-\sin\theta&\cos\theta\end{pmatrix},\\ \mathbf{R}_{y}(\theta)=&\begin{pmatrix}\cos\theta&0&-\sin\theta\\ 0&1&0\\ \sin\theta&0&\cos\theta\end{pmatrix},\\ \mathbf{R}_{z}(\theta)=&\begin{pmatrix}\cos\theta&\sin\theta&0\\ -\sin\theta&\cos\theta&0\\ 0&0&1\end{pmatrix}.\end{split}

4.2.1 the Sun’s quadrupole moment harmonics

When calculating the solar quadrupole moment, the formula for the solar axis vector $\mathbf{p}$ in Eq. (4.4) is given by:

\mathbf{p}=\mathbf{R}_{x}(-\varepsilon)\begin{pmatrix}\cos\Omega\sin I\\ \sin\Omega\sin I\\ \cos I\end{pmatrix}=\begin{pmatrix}1&0&0\\ 0&\cos\varepsilon&-\sin\varepsilon\\ 0&\sin\varepsilon&\cos\varepsilon\end{pmatrix}\begin{pmatrix}\cos\Omega\sin I% \\ \sin\Omega\sin I\\ \cos I\end{pmatrix}

where $\Omega$ and $I$ are respectively the longitude of the ascending node and inclination of the solar equatorial plane relative to the ecliptic plane; $I=7^{\circ}15^{\prime}$ , $\Omega=75^{\circ}46^{\prime}+84^{\prime\prime}\cdot T$ , and $\varepsilon=23^{\circ}26^{\prime}21^{\prime\prime}.448$ is the obliquity of the ecliptic for the J2000 epoch.

T=\frac{JD(t)-JD(2000)}{36525.0}=\frac{JD(t)-2451545.0}{36525.0},

(4.7)

where $t$ is the dynamical time, JD represents the Julian Date, JD(t) denotes the Julian Date corresponding to the dynamical time, and JD(2000) represents the Julian Date corresponding to the J2000 epoch, which is 2451545.0.

4.2.2 Earth’s 2nd to 4th degree zonal harmonics

To obtain the Earth’s pole vector $\mathbf{p}$ , we need to compute the Earth’s precession and nutation. In CGC2.7, instead of numerical integration, we calculate the precession and nutation using the theoretical framework outlined by J. Wahr and H. Kinoshita in the IAU 1980 theory [41, 42].

Earth’s Precession

Precession refers to the transformation between the epoch mean equatorial geocentric coordinates and the true equatorial geocentric coordinates, representing the difference between these two coordinate systems. The transformation method is outlined as follows [43]:

\mathbf{r}_{m}=(\mathbf{PR})\mathbf{R},

where $\mathbf{PR}$ is the precession matrix, composed of three rotation matrices:

\mathbf{PR}=\mathbf{R}_{z}(-z_{A})\mathbf{R}_{y}(\theta_{A})\mathbf{R}_{z}(-% \xi_{A}),

where $\xi_{A}$ , $z_{A}$ , and $\theta_{A}$ are the equatorial precession angles, computed as:

\left\{\begin{aligned} &\xi_{A}=2306^{\prime\prime}.2181T+0^{\prime\prime}.301% 88T^{2}+0^{\prime\prime}.017998T^{3}\\ &z_{A}=2306^{\prime\prime}.2181T+1^{\prime\prime}.09468T^{2}+0^{\prime\prime}.% 018203T^{3}\\ &\theta_{A}=2004^{\prime\prime}.3109T-0^{\prime\prime}.42665T^{2}-0^{\prime% \prime}.041833T^{3}\end{aligned}\right.

The corresponding right ascension precession $m_{A}$ and declination precession $n_{A}$ are:

\left\{\begin{aligned} &m_{A}=\xi_{A}+z_{A}=4612^{\prime\prime}.4362T+1^{% \prime\prime}.39656T^{2}+0^{\prime\prime}.036201T^{3}\\ &n_{A}=\theta_{A}\end{aligned}\right.

Nutation of the Earth

Instantaneous mean equatorial geocentric coordinate system and instantaneous true equatorial geocentric coordinate system conversion. The difference between these two coordinate systems is known as nutation. The conversion method is according to [43]:

\mathbf{r}_{t}=(\mathbf{NR})\mathbf{r}_{m}

where $\mathbf{NR}$ is the nutation matrix, composed of three rotation matrices:

\mathbf{NR}=\mathbf{R}_{x}(-\Delta\varepsilon)\mathbf{R}_{y}(\Delta\theta)% \mathbf{R}_{z}(-\Delta\mu),

where $\Delta\mu$ , $\Delta\theta$ , and $\Delta\varepsilon$ represent nutation in right ascension, nutation in declination, and nutation in obliquity, respectively. The nutation series adopted is from IAU(1980), and for meter precision, the first 20 terms of this series are used. The computation formulas are as follows:

\left\{\begin{aligned} &\Delta\psi=\sum_{j=1}^{20}(A_{0j}+A_{1j}t)\sin(\sum_{i% =1}^{5}k_{ji}\alpha(t))\\ &\Delta\varepsilon=\sum_{j=1}^{20}(B_{0j}+B_{1j}t)\cos(\sum_{i=1}^{5}k_{ji}% \alpha(t))\end{aligned}\right.\Rightarrow\left\{\begin{aligned} &\Delta\mu=% \Delta\psi\cos\varepsilon\\ &\Delta\theta=\Delta\psi\sin\varepsilon\end{aligned}\right.

where $\psi$ is nutation in longitude, and $\Delta\mu$ and $\Delta\theta$ are nutation in right ascension and nutation in declination, respectively. The formula for the obliquity of the ecliptic $\varepsilon$ is:

\varepsilon=23^{\circ}26^{\prime}21^{\prime\prime}.448-46^{\prime\prime}.8150T.

The formulas for calculating the five angular terms in the nutation series are:

\left\{\begin{aligned} &\alpha_{1}=134^{\circ}57^{\prime}46^{\prime\prime}.733% +(1325^{r}+198^{\circ}52^{\prime}02^{\prime\prime}.633)T+31^{\prime\prime}.310% T^{2}\\ &\alpha_{2}=357^{\circ}31^{\prime}39^{\prime\prime}.804+(99^{r}+359^{\circ}03^% {\prime}01^{\prime\prime}.224)T-0^{\prime\prime}.577T^{2}\\ &\alpha_{3}=93^{\circ}16^{\prime}18^{\prime\prime}.877+(1342^{r}+82^{\circ}01^% {\prime}03^{\prime\prime}.137)T-13^{\prime\prime}.257T^{2}\\ &\alpha_{4}=297^{\circ}51^{\prime}01^{\prime\prime}.307+(1236^{r}+307^{\circ}0% 6^{\prime}41^{\prime\prime}.328)T-6^{\prime\prime}.891T^{2}\\ &\alpha_{5}=125^{\circ}02^{\prime}40^{\prime\prime}.280-(5^{r}+134^{\circ}08^{% \prime}10^{\prime\prime}.539)T+7^{\prime\prime}.455T^{2}\end{aligned}\right.

where $1^{r}=360^{\circ}$ , and the relevant coefficients for the first 20 terms of the nutation series are shown in Table 4.1.

Table 4.1: IAU1980 coefficients for the first 20 terms of the nutation

$j$	period (Days)	$k_{j1}$	$k_{j2}$	$k_{j3}$	$k_{j4}$	$k_{j5}$	$A_{0j}$ ¹¹1In, Table 4.1, $A_{0j},A_{1j},B_{0j},B_{1j}$ are in unit of $0^{\prime\prime}.0001$	$A_{1j}$	$B_{0j}$	$B_{1j}$
1	6798.4	0	0	0	0	1	-171996	-174.2	92025	8.9
2	182.6	0	0	2	-2	2	-13187	-1.6	5736	-3.1
3	13.7	0	0	2	0	2	-2274	-0.2	977	-0.5
4	3399.2	0	0	0	0	2	2062	0.2	-895	0.5
5	365.2	0	1	0	0	0	1426	-3.4	54	-0.1
6	27.6	1	0	0	0	0	712	0.1	-7	0.0
7	121.7	0	1	2	-2	2	-517	1.2	224	-0.6
8	13.6	0	0	2	0	1	-386	-0.4	200	0.0
9	9.1	1	0	2	0	2	-301	0.0	129	-0.1
10	365.3	0	-1	2	-2	2	217	-0.5	-95	0.3
11	31.8	1	0	0	-2	0	-158	0.0	-1	0.0
12	177.8	0	0	2	-2	1	129	0.1	70	0.0
13	27.1	-1	0	2	0	2	123	0.0	-53	0.0
14	27.7	1	0	0	0	1	63	0.1	-33	0.0
15	14.8	0	0	0	2	0	63	0.0	-2	0.0
16	9.6	-1	0	2	2	2	-59	0.0	26	0.0
17	27.4	-1	0	0	0	1	-58	-0.1	32	0.0
18	9.1	1	0	2	0	1	-51	0.0	27	0.0
19	205.9	2	0	0	-2	0	48	0.0	1	0.0
20	1305.5	-2	0	2	0	1	46	0.0	-24	0.0

Based on the above discussion, the transformation from the epoch mean equatorial geocentric coordinates to the instantaneous true equatorial geocentric coordinates is given by:

\mathbf{r}_{t}=(\mathbf{NR})(\mathbf{PR})\mathbf{r}.

(4.8)

4.3 Perturbations of Asteroids

Asteroid data numbered up to April 14, 2010, were downloaded from the Lowell Asteroid Database [44]. Asteroids with diameter data and classified density were selected to calculate asteroid masses. Due to the large number of asteroids, the handling in the ephemeris is as follows:

•

Select asteroids with mass $\mu>10^{-15}$ for (1) Ceres, (2) Pallas, and (4) Vesta, along with the Sun, the eight major planets, Pluto, and the Moon, for numerical integration to obtain state data.
•

For the remaining selected 349 asteroids, numerical integration is not used to compute their state at a given dynamical time. Instead, their states at specific times are calculated based on their orbital elements and Keplerian motion, using the method of undisturbed two-body motion, followed by calculation of the asteroids’ perturbation on the other 11 major integrated bodies.

Table 4.2 presents the classification statistics of the selected 349 asteroids. Density values for C, S, and M classifications are taken from JPL DE405 (Standish, 1998). The density for E classification is based on the average value from Wasson (1974). The G classification is considered a subtype of C, therefore using the same density value as C. For U-type asteroids, the density is calculated as a weighted average of densities from the other types [39].

Table 4.2: Classification and Statistics of Asteroids(Bowell, 1999)

Type	Symbol	Density (g/ $\text{cm}^{2}$ )	Number	Mass $(\text{AU}^{3}/\text{day}^{2})$
Carbonaceous chondrite	C	1.8	141	$6.36\times 10^{-14}$
Silicate (S-type)	S	2.4	93	$2.37\times 10^{-14}$
Metallic (M-type)	M	5.0	30	$1.39\times 10^{-14}$
Enstatite achondrite	E	3.65	3	$1.88\times 10^{-16}$
Extremely ultraviolet	G	1.8	5	$3.07\times 10^{-15}$
Others	U	2.16	77	$2.11\times 10^{-14}$

Newtonian Perturbations of Selected Asteroids[43]
The Kepler orbit for asteroids is described as

E-e\sin E=n(t-\tau),

(4.9)

where $\tau$ is the time when the asteroid passes through perihelion, and $n=\frac{2\pi}{T}$ . By using Eq. (4.9), we can determine $E$ , and subsequently, with known orbital elements of the asteroid, derive its position and velocity vectors. The transformation of asteroid orbital elements to position and velocity vectors in a Cartesian coordinate system is given by:

\mathbf{r}=a(\cos E-e)\widehat{P}+a\sqrt{1-e^{2}}\sin E\widehat{Q},

(4.10)

\dot{\mathbf{r}}=\frac{\sqrt{\mu a}}{r}[-\sin E\widehat{P}+\sqrt{1-e^{2}}\cos E% \widehat{Q}],

(4.11)

in which

\widehat{P}=\begin{pmatrix}\cos\Omega\cos\omega-\sin\Omega\sin\omega\cos i\\ \sin\Omega\cos\omega+\cos\Omega\sin\omega\cos i\\ \sin\omega\sin i\end{pmatrix},

\widehat{Q}=\begin{pmatrix}-\cos\Omega\sin\omega-\sin\Omega\cos\omega\cos i\\ -\sin\Omega\sin\omega+\cos\Omega\cos\omega\cos i\\ \cos\omega\sin i\end{pmatrix},

where the orbital elements are defined as follows: $a$ is the semi-major axis of the orbit; $e$ is the eccentricity of the orbit; $\omega$ is the argument of periapsis; $\Omega$ is the longitude of the ascending node; $i$ is the inclination of the orbital plane; $M$ is the mean anomaly.

4.4 Numerical Integration Methods

For a second-order differential equation

\frac{d^{2}r}{dt^{2}}=f(t,r,\dot{r}),

(4.12)

the classical fourth-order Runge-Kutta numerical integration method is employed [45]:

\left\{\begin{aligned} &r_{n+1}=r_{n}+h\dot{r}_{n}+\frac{h}{6}(M_{1}+M_{2}+M_{% 3})\\ &\dot{r}_{n+1}=\dot{r}_{n}+\frac{1}{6}(M_{1}+2M_{2}+2M_{3}+M_{4})\\ &M_{1}=hf(t_{n},r_{n},\dot{r}_{n})\\ &M_{2}=hf(t_{n}+\frac{h}{2},r_{n}+\frac{h}{2}\dot{r}_{n},\dot{r}_{n}+\frac{M_{% 1}}{2})\\ &M_{3}=hf(t_{n}+\frac{h}{2},r_{n}+\frac{h}{2}\dot{r}_{n}+\frac{h}{4}M_{1},\dot% {r}_{n}+\frac{M_{2}}{2})\\ &M_{4}=hf(t_{n}+h,r_{n}+h\dot{r}_{n}+\frac{h}{2}M_{2},\dot{r}_{n}+M_{3})\end{% aligned}\right.

(4.13)

where $h$ is a step, and $\dot{r}$ is the first-order differential with respect to $t$ . CGC2.7 ephemeris utilizes the fourth-order Runge-Kutta numerical integration method, implemented in C++ programming language. The integration step size is 0.02 days. The starting time is June 21, 2028, 12:00 (JD 2461944.0). Initial positions for 14 celestial bodies are provided by JPL’s DE405 ephemeris [46], and the masses and harmonic parameters for each body are obtained from the DE405 ephemeris header files.

4.5 Accuracy of CGC2.7 comparing to DE405

To assess the accuracy of CGC2.7 ephemeris framework, we compared it with the DE405 ephemeris [46]. which released by JPL. The comparison results are shown in Figure 4.1. We compared the data of inner planets (Mercury, Venus, Earth, and Mars) over 10 years. Using the Sun as the origin and employing ecliptic coordinates from DE405, we plotted the differences between the two calculated orbits for each planet. For example, in the first row of comparisons for Mercury in Figure 4.1, the first column shows the variation over time of the difference in heliocentric distances between DE405 and CGC2.7. The second column displays the difference in Mercury’s ecliptic longitude over time, and the third column shows the difference in ecliptic latitude. Similarly, the second row depicts comparison for Venus’s orbit, the third row shows Earth’s orbit, and the fourth row shows results of Mars’s orbit.

Chapter 5 ASTROD-GW Mission Orbit

The ASTROD-GW mission employs three spacecraft located near the Sun-Earth Lagrange points L3, L4, and L5. These spacecraft orbit the Sun in nearly circular orbits, forming an approximate equilateral triangle as shown in Figure 5.1. Each side of the triangle is approximately 260 million kilometers (1.732 AU). Lagrange points L4 and L5 are stable, while L3 is an unstable point with an instability timescale of about 50 years. Given that the expected mission duration of ASTROD-GW is 20 years, the spacecraft can maintain stable positions near L3, L4, and L5 throughout the mission period. Additionally, each spacecraft is equipped with micropropulsion thrusters for orbit corrections to mitigate disturbances and maintain precise positioning.

The LISA spacecraft orbits the Sun in a triangular formation with arms of approximately 5 million kilometers each. Given that ASTROD-GW has an interferometric arm length 52 times longer than LISA, if the relative arm length difference of ASTROD-GW can be controlled to be less than 1/52 of LISA’s relative arm length difference and the relative velocity between spacecraft is not greater than that of LISA, then the technological requirements for laser interferometry will not be higher than those already developed for LISA. With its longer arm length, ASTROD-GW exhibits higher sensitivity to lower-frequency GWs. This chapter optimizes the mission orbit according to established methods [47, 48].

5.1 Initial Orbit Selection

The initial time for ASTROD-GW mission orbit is set at 12:00 on June 21, 2028 (JD2461944). Considering the Solar System as a restricted three-body problem in an elliptical orbit on the ecliptic plane, we determine the locations of the five Lagrange points. Initially, the spacecraft are positioned near points L3, L4, and L5. The positions of the three spacecraft in the heliocentric ecliptic coordinates are as follows: L3 is located near (0, 1 AU, 0), L4 is near ( $\sqrt{3}/2\ \text{AU},-1/2\ \text{AU},0$ ), and L5 is near ( $-\sqrt{3}/2\ \text{AU},-1/2\ \text{AU},0$ ). This forms an approximately equilateral triangle with sides of about $\sqrt{3}$ AU, and each spacecraft is about 1 AU from the Sun. The initial velocities of the spacecraft are calculated based on the stellar year of the Earth’s orbital period (365.26536 days). The orbital speed of the spacecraft in circular motion around the Sun is $v_{0}=0.01720209895\ \text{AU/day}$ . The initial velocity vectors of the spacecraft in the ecliptic plane are perpendicular to their initial position vectors: at L3, the spacecraft has an initial velocity of $(-v_{0},0,0)$ ; at L4, the initial velocity is $(1/2\ v_{0},\sqrt{3}/2\ v_{0},0)$ ; and at L5, the initial velocity is $(1/2\ v_{0},-\sqrt{3}/2\ v_{0},0)$ .

After determining the initial states of the spacecraft in the heliocentric ecliptic coordinates, the coordinates are transformed to the J2000 Solar System Barycentric (SSB) coordinates used for ephemeris calculations. The J2000.0 obliquity of the ecliptic is $\varepsilon=23^{\circ}26^{\prime}21^{\prime\prime}.448$ . The transformation from spacecraft positions and velocities in heliocentric ecliptic coordinates to Solar System Barycentric coordinates is given by Eq. (5.1), where $\mathbf{R}_{SC,SSB}$ , $\mathbf{V}_{SC,SSB}$ represent the spacecraft positions and velocities in J2000 Solar System Barycentric coordinates, $\mathbf{R}_{SC,Sun}$ , $\mathbf{V}_{SC,Sun}$ represent the spacecraft positions and velocities in heliocentric ecliptic coordinates, and $\mathbf{R}_{\odot,SSB}$ , $\mathbf{V}_{\odot,SSB}$ represent the Sun’s position and velocity in J2000 Solar System Barycentric coordinates.

\begin{split}&\mathbf{R}_{SC,SSB}=\left[\begin{array}[]{rrr}1&0&0\\ 0&\cos{\varepsilon}&-\sin{\varepsilon}\\ 0&\sin{\varepsilon}&\cos{\varepsilon}\end{array}\right]\cdot\mathbf{R}_{SC,Sun% }+\mathbf{R}_{\odot,SSB},\\ &\mathbf{V}_{SC,SSB}=\left[\begin{array}[]{rrr}1&0&0\\ 0&\cos{\varepsilon}&-\sin{\varepsilon}\\ 0&\sin{\varepsilon}&\cos{\varepsilon}\end{array}\right]\cdot\mathbf{V}_{SC,Sun% }+\mathbf{V}_{\odot,SSB}.\end{split}

(5.1)

5.2 Orbital Optimization

Through the initial orbit selection, we determine the initial conditions of the spacecraft. By computing the spacecraft’s orbit over a 20-year mission period and analyzing the spacecraft orbit data, we obtain graphs showing the orbital periods and heliocentric distances of the three spacecraft over time. On one hand, we adjust the average 20-year period of the spacecraft to align it as closely as possible with the Earth’s orbital period. On the other hand, we adjust the spacecraft’s eccentricity to approach a nearly circular orbit. However, the optimal period may not necessarily be exactly 1 stellar year. In necessary cases, we may slightly deviate the spacecraft’s 20-year average period from 1 stellar year to achieve better optimization results.

5.2.1 Orbital Optimization Methodology

In this optimization process [47, 48, 49], we utilize the heliocentric coordinate system. The total energy of a planet’s motion around the Sun is given by:

E=\frac{mv^{2}}{2}-\frac{G(M+m)m}{r}=-\frac{G(M+m)m}{2a},

(5.2)

where $M$ is the mass of the Sun, $m$ is the mass of the planet, $G$ is the gravitational constant, $v$ is the relative velocity between the planet and the Sun, $r$ is the heliocentric distance of the planet, and $a$ is the semi-major axis of the planet’s elliptical orbit around the Sun.
The orbital period of a planet around the Sun is given by:

T=\frac{2\pi a^{3/2}}{\sqrt{G(M+m)}}.

(5.3)

To derive Eq. (5.2), take the total differential:

mv^{2}\frac{dv}{v}+\frac{G(M+m)m}{r}\frac{dr}{r}=\frac{G(M+m)m}{2a}\frac{da}{a}.

(5.4)

Assuming a nearly circular orbit for the spacecraft, where $r\approx a$ in Eq. (5.2), we get:

mv^{2}\approx\frac{G(M+m)m}{a}\approx\frac{G(M+m)m}{r}.

(5.5)

Substituting Eq. (5.5) into Eq. (5.4), we obtain the relationship:

\frac{dv}{v}+\frac{dr}{r}=\frac{1}{2}\frac{da}{a}.

(5.6)

Taking the logarithm of both sides of Eq. (5.3) and differentiating gives:

\frac{dT}{T}=\frac{3}{2}\frac{da}{a}\quad\Rightarrow\quad da=\frac{2}{3}\frac{% dT}{T}\cdot a.

(5.7)

5.2.1.1 Optimization of Mission Orbital Period

Optimize the spacecraft’s orbital period by adjusting the velocity and modifying the heliocentric distance. Taking the logarithm and differentiating Eq. (5.5), we obtain:

2\frac{dv}{v}\approx-\frac{da}{a}\approx-\frac{dr}{r}.

(5.8)

From Eqs. (5.7) and (5.8), combining these two equations yields:

dv\approx-\frac{1}{3}\frac{dT}{T}v,\quad dr\approx\frac{2}{3}\frac{dT}{T}r.

(5.9)

In practical adjustment processes, the above equations should be written in vector form:

\delta\mathbf{V}\approx-\frac{1}{3}\frac{dT}{T}\mathbf{V},\quad\delta\mathbf{R% }\approx\frac{2}{3}\frac{dT}{T}\mathbf{R}.

(5.10)

The formulas for calculating the adjusted velocity and position after adjusting the period are:

\begin{split}&\mathbf{V}_{new}=\mathbf{V}+\delta\mathbf{V}\approx\left(1-\frac% {1}{3}\frac{dT}{T}\right)\mathbf{V},\\ &\mathbf{R}_{new}=\mathbf{R}+\delta\mathbf{R}\approx\left(1+\frac{2}{3}\frac{% dT}{T}\right)\mathbf{R}.\end{split}

(5.11)

5.2.1.2 Optimization of Mission Orbital Eccentricity

As following, we attempt to adjust the spacecraft’s orbital eccentricity to make it closer to a circular orbit. From Eq. (5.3), it is evident that as long as the semi-major axis of the orbit remains unchanged, the orbital period will not change. Therefore, in the subsequent optimization steps, we initially set $da=0$ . From Eq. (5.6), we have:

\frac{dv}{v}=-\frac{dr}{r}

(5.12)

Typically, based on the variation of heliocentric distance over time, we determine the adjustment $dr$ , and then adjust the velocity and position accordingly. The formulas used in the actual adjustment process are:

\delta\mathbf{V}=-\frac{dR}{R}\mathbf{V},\quad\delta\mathbf{R}=\frac{dR}{R}% \mathbf{R},

(5.13)

where R denotes the initial heliocentric distance. The formulas for calculating the adjusted velocity and position after adjusting the eccentricity are:

\begin{split}&\mathbf{V}_{new}=\mathbf{V}+\delta\mathbf{V}\approx\left(1-\frac% {dR}{R}\right)\mathbf{V},\\ &\mathbf{R}_{new}=\mathbf{R}+\delta\mathbf{R}\approx\left(1+\frac{dR}{R}\right% )\mathbf{R}.\end{split}

(5.14)

5.2.2 Process of Orbital Optimization

In the process of optimizing the spacecraft’s orbit, we use in the heliocentric ecliptic coordinate system. And the orbit is computed using the CGC2.7 ephemeris framework, followed by plotting and analysis. The graphs include variations over time in the arm lengths between spacecraft, variations in arm length differentials over time, changes in heliocentric distances of spacecraft over time, and changes in Doppler relative velocities between spacecraft over time.

The optimization method for spacecraft orbit periods involves adjusting the spacecraft periods based on the variation of arm length differentials over time to minimize these differences. Due to various influencing factors, the average periods of the spacecraft change continuously over time. The period at L3 tends to decrease over time, while those at L4 and L5 tend to increase. Arm lengths Arm12 initially decrease and then increase, Arm13 increase initially and then decrease, while Arm23 shows relatively smaller changes, as shown in Figure 5.2. According to Eq. (5.11), adjusting the periods involves reducing the initial velocity of S/C1 to lengthen its initial period, and increasing the initial velocities of S/C2 and S/C3 to shorten their initial periods. This approach compensates for deviations in the spacecraft periods during operation, thereby controlling the differences in periods among the three spacecraft within a certain range and appropriately managing the arm length differences. The initial conditions of the spacecraft after adjusting the mission orbit periods are summarized in Table 5.1, with corresponding time-varying graphs shown in Figure 5.3.

Table 5.1: Initial position and velocities after optimizing the orbital periods at epoch JD2461944.0.

( J2000 solar-system-barycentric equatorial, in AU and AU/day ) x(AU)/vx(AU/day) y(AU)/vy(AU/day) z(AU)/vz(AU/day) S/C1 position 1.15400625657242E-3 9.15261701184544E-1 3.96854368692135E-1 S/C1 velocity -1.72008300146903E-2 4.88112077380618E-6 2.07014410548162E-6 S/C2 position 8.67179419751799E-1 -4.60958432557285E-1 -1.99809748470332E-1 S/C2 velocity 8.60233950755836E-3 1.36730847422802E-2 5.92795834111164E-3 S/C3 position -8.64871407921666E-1 -4.60958431471891E-1 -1.99809747999756E-1 S/C3 velocity 8.60233303252539E-3 -1.36633074545593E-2 -5.92381152958984E-3

The method of optimizing the orbital period can reduce the arm length difference to a certain extent, but usually such orbits cannot fully meet mission requirements, which requires us to adopt other methods for further optimization. According to the spacecraft’s variation in heliocentric distance, we see that the spacecraft’s orbital heliocentric distance changes within a range of approximately $(1\pm 6\times 10^{-5})$ AU. According to Eq. (5.14), we attempt to reduce the orbit eccentricity by adjusting the initial heliocentric distance with a variation of less than $6\times 10^{-5}$ AU. As for whether to increase or decrease the initial heliocentric distance, both cases are calculated and compared to find the optimized results. When the optimized initial heliocentric distance approaches the periapsis or apoapsis distance of the spacecraft’s orbit, further optimization of the orbit eccentricity is not feasible.

In most cases, the two orbit optimization methods cannot achieve satisfactory results in once optimization. Therefore, multiple iterations of orbit period and orbit eccentricity optimization are needed to obtain the optimized results that meet the requirements of the ASTROD-GW mission. The final optimized results for the ASTROD-GW mission orbit are shown in Table 5.2, the average changes in spacecraft orbit periods are shown in Table 5.3, and the variations in interference arm lengths, arm length differences, interference arm angles, and Doppler velocities over time are shown in Figure 5.4.

Table 5.2: Initial conditions after final optimization at epoch JD2461944.0

( J2000 solar-system-barycentric equatorial) x(AU)/vx(AU/day) y(AU)/vy(AU/day) z(AU)/vz(AU/day) S/C1 position 1.15400625657242E-3 9.15289225648841E-1 3.96866302001196E-1 S/C1 velocity -1.72003163872199E-2 4.88112077380618E-6 2.07014410548162E-6 S/C2 position 8.67153438989685E-1 -4.60944670325136E-1 -1.99803781815802E-1 S/C2 velocity 8.60259754371050E-3 1.36734947883888E-2 5.92813611775755E-3 S/C3 position -8.64862747667628E-1 -4.60953844061175E-1 -1.99807759114913E-1 S/C3 velocity 8.60242121981651E-3 -1.36634452887228E-2 -5.92387128797985E-3

Table 5.3: Periods of three S/C

	5-yr average	10-yr average	15-yr average	20-yr average
S/C1(day)	365.25662	365.25767	365.25636	365.25636
S/C2(day)	365.25591	365.25564	365.25620	365.25646
S/C3(day)	365.25624	365.25420	365.25656	365.25721

Chapter 6 Numerical Simulation for Time-Delay Interferometry

In Chapter 3, we briefly introduced the principles of unequal arm-length TDI. However, in a realistic space mission, due to the dynamics of spacecraft orbits, the lengths of interferometric arms vary over time. This variability makes the laser noise cannot be canceled using traditional methods. In this chapter, we will further analyze TDI and numerically calculate the optical paths of TDI to eliminate laser noise.

6.1 Algorithm of Calculation

For ASTROD-GW, the distance between spacecraft remains approximately $\sqrt{3}$ AU, and the light signal takes about 15 minutes to travel from one spacecraft to another. Therefore, the signal received by a spacecraft is emitted by the other spacecraft approximately 15 minutes earlier, necessitating time delay calculation. During laser interferometry between spacecraft, due to the continuous motion of the transmitting and receiving spacecraft in the coordinate system, and the time delay between emission and reception, it is crucial to accurately determine the position of the receiving spacecraft at the time of signal reception. The CGC2.7 ephemeris framework can calculate a series of discrete states of each spacecraft during the mission period, but it does not include all states at the moment of signal reception. To obtain the state of a spacecraft at any arbitrary time, we employ the Chebyshev polynomial interpolation method based on the spacecraft’s existing state data [51, 52].

In computing TDI, the instantaneous state of the spacecraft emitting the laser is known, but the state of the receiving spacecraft at the reception time is unknown in advance. It requires an iterative approach to approximate the state at reception time. The specific calculation method is as follows: Since the propagation of light in spacetime is not perfectly straight, the accuracy of the calculation requires consideration of post-Newtonian corrections in the solar system. Taking the post-Newtonian effect of the Sun as an example, considering the solar center as the origin of the coordinate system, point P1 stationary in the reference frame, with position vector $\mathbf{R_{1}}$ , and point P2 with position vector $\mathbf{R_{2}}$ , from which a light signal is emitted from P1 to P2. The propagation time of the light signal $T_{Travel}$ is divided into two parts: one part is the calculation under flat spacetime, and the other part is the correction due to the post-Newtonian effect.

T_{Travel}=T_{Newton}+\Delta T_{PN},

(6.1)

where $T_{Newton}$ is the propagation time in Newtonian flat spacetime, and $\Delta T_{PN}$ is the corresponding time correction due to the post-Newtonian effect. (1) Calculation of propagation time in Newtonian gravity approximation [53]:

T_{Newton}=\frac{R}{c};

(6.2)

(2) Calculation of time correction for PN (Post-Newtonian Light Propagation) additional time delay:

\begin{split}\Delta T_{PN}=&\frac{2GM}{c^{3}}\ln\left(\frac{R_{1}+R_{2}+R}{R_{% 1}+R_{2}-R}\right)\\ &+\frac{G^{2}M^{2}}{c^{5}}\frac{R}{R_{1}R_{2}}\left[\frac{15}{4}\frac{\arccos(% \mathbf{N}_{1}\cdot\mathbf{N}_{2})}{|\mathbf{N}_{1}\times\mathbf{N}_{2}|}-% \frac{4}{1+\mathbf{N}_{1}\cdot\mathbf{N}_{2}}\right],\end{split}

(6.3)

where $GM$ is the product of the mass of the celestial body for which PN effect is to be calculated and the gravitational constant, $c$ denotes the speed of light, $R_{1}=|\mathbf{R}_{1}|,R_{2}=|\mathbf{R}_{2}|$ are the distances of P1 and P2 from the origin, $\mathbf{N}_{1}=\mathbf{R}_{1}/R_{1},\mathbf{N}_{2}=\mathbf{R}_{2}/R_{2}$ are unit vectors pointing from the origin to P1 and P2 respectively, and $R=|\mathbf{R}_{2}-\mathbf{R}_{1}|$ represents the distance between P1 and P2. Here we consider only the post-Newtonian effect of the Sun.

For calculating the propagation time between spacecraft, the method is as follows. Without loss of generality, assume at time $T_{0}$ , S/C1 emits a laser signal towards S/C2, and S/C2 receives the signal at time $T^{r}_{0}$ . The propagation time in the coordinate system is $T^{r}_{0}-T_{0}$ . The iterative method to compute the time when S/C2 receives the signal is as follows:

$\displaystyle T^{r}_{0}$	$\displaystyle=T_{0}+T_{1}+T_{2}+T_{3}+...$	(6.4)
$\displaystyle(i)\ T_{1}$	$\displaystyle=\frac{\|r_{2}(T_{0})-r_{1}(T_{0})\|}{c}+\Delta T_{1,PN}$
$\displaystyle(ii)\ T_{1}+T_{2}$	$\displaystyle=\frac{\|r_{2}(T_{0}+T_{1})-r_{1}(T_{0})\|}{c}+\Delta T_{2,PN}$
$\displaystyle(iii)\ T_{1}+T_{2}+T_{3}$	$\displaystyle=\frac{\|r_{2}(T_{0}+T_{1}+T_{2})-r_{0}(T_{0})\|}{c}+\Delta T_{3,PN}$
	$\displaystyle......,$

where $\Delta T_{i,PN}$ represents the correction to the light propagation time due to PN effect computed in the $i$ -th iteration. Perform multiple iterations until $|T_{i}|<10^{-11}$ s (3 mm). The value of $r_{2}(T)$ is obtained using Chebyshev polynomial interpolation [52]. We employ a 14th-order Chebyshev polynomial interpolation over an interval of 8 days, achieving interpolation accuracy on the order of $10^{-14}$ AU.

6.2 Preliminary Analysis of Michelson-type TDI

Due to the mutual laser interferometric links between the three spacecraft, each spacecraft has two optical platforms, resulting in six links when all six optical platforms are functional. However, there may be cases where optical platforms do not function properly. In the following discussion, we consider two scenarios: (1) when interference is performed using 2 out of the 3 interferometric arms in this section; (2) when all 3 interferometric arms are simultaneously operational in next chapter.

6.2.1 Michelson-type TDI

In space missions, it is possible that one of the three interferometric arms may not function properly. In such cases, we can refer to the unequal arm Michelson interference and optimize the interference method to partially meet the requirements of interferometry and achieve the corresponding sensitivity goals.

The path of light propagation involves path $a$ , where spacecraft S/C1 emits a laser to S/C2 and then returns to S/C1; and path $b$ , where S/C1 emits a laser to S/C3 and then returns to S/C1 via S/C3. When the two laser follow their respective paths and return to S/C1, the propagation times of the two laser are examined. When two paths exhibit symmetry, the combination of paths (will described by $a$ and $b$ forms the paths) traveled by the two laser, effectively minimizing the mismatch of two paths and thus reducing the laser noise.

: Path a : $\text{SC1}\rightarrow\text{SC2}\rightarrow\text{SC1}$ ;
: Path b : $\text{SC1}\rightarrow\text{SC3}\rightarrow\text{SC1}$ .

According to theoretical calculations, for the 1st-generation TDI paths, static arm length differences can be eliminated, and for the 2nd-generation TDI paths, differences in velocities with the same relative velocities can be eliminated. There are several groups of TDI paths [54]:

	$\displaystyle\text{1st generation}:$	$\displaystyle 1\rightarrow 2\rightarrow 1\rightarrow 3\rightarrow 1(ba-ab)$
	$\displaystyle\text{2nd generation}:$	$\displaystyle 1\rightarrow 2\rightarrow 1\rightarrow 3\rightarrow 1\rightarrow 3% \rightarrow 1\rightarrow 2\rightarrow 1(baab-abba)$

	$\displaystyle(i)\qquad n=1,$	$\displaystyle[ab,ba]\equiv abba-baab$
	$\displaystyle(ii)\qquad n=2,$	$\displaystyle[a^{2}b^{2},b^{2}a^{2}],[abab,baba],[ab^{2}a,ba^{2}b]$
	$\displaystyle(iii)\qquad n=3,$	$\displaystyle[a^{3}b^{3},b^{3}a^{3}],[a^{2}bab^{2},b^{2}aba^{2}],[a^{2}b^{2}ab% ,b^{2}a^{2}ba],[a^{2}b^{3}a,b^{2}a^{3}b],$
		$\displaystyle[aba^{2}b^{2},bab^{2}a^{2}],[ababab,bababa],[abab^{2}a,baba^{2}b]% ,[ab^{2}a^{2}b,ba^{2}b^{2}a],$
		$\displaystyle[ab^{2}aba,ba^{2}bab],[ab^{3}a^{2},ba^{3}b^{2}],\text{% lexicographic (binary) order}$

6.3 Results of Numerical Simulation

According to the numerical calculation method, the results of mismatches in TDI channels in Section 6.2 are shown in Figures 6.3 and 6.4. The horizontal axis represents mission time (in days), and the vertical axis represents path mismatch (time difference in seconds).

ASTROD-GW requires that the path difference between two paths be within 150 ns (50 m). Both [ab, ba] and [abba, baab] satisfy this requirement. Experimentally, after demonstrating compliance with the path difference requirement, ASTROD-GW also meets the noise requirements. In 2010, de Vine et al [57] from the JPL demonstrated the implementation of LISA’s TDI experiment in the laboratory. Through TDI, laser frequency noise was reduced by approximately $10^{9}$ , and clock phase noise was reduced by $6\times 10^{4}$ , restoring the system’s intrinsic displacement noise on the laboratory test bench. ASTROD-GW should also undergo relevant laboratory experiments, with the feasibility of the principle being similar.

Chapter 7 Geometry Analysis and Numerical Simulation for TDI

To facilitate discussion, we label the interferometric arms of the three spacecraft as follows: spacecraft $i$ ’s interferometric arms are denoted as $L_{i},L_{i^{\prime}}$ ; clockwise, they are labeled as $L_{1^{\prime}},L_{2^{\prime}},L_{3^{\prime}}$ ; counterclockwise, they are $L_{1},L_{2},L_{3}$ . This is illustrated in Figure 7.1. Using $y_{ij}$ to denote interferometric measurement from spacecraft $j$ originating from spacecraft $i$ .

7.1 First-Generation TDI

The first-generation TDI eliminates static differences in arm lengths through specific paths. Therefore, in this analysis, we initially consider the arm lengths as static. Commas denote time delays, where $y_{ij,k}=y_{ij}(t-L_{k})$ and $y_{ij,kr}=y_{ij}(t-L_{k}-L_{r})$ [56]. Because the arm lengths are static, the time delays following the commas can be exchanged. We proceed with a detailed analysis for first-generation TDI as the first case.

7.1.1 TDI with Six Links

7.1.1.1 Sagnac

The Sagnac effect includes three channels, designated as $\alpha,\beta,\gamma$ depending on which spacecraft serves as the starting point [55]. Sagnac ( $\alpha$ ) employs six interferometric arms, illustrated in Figure 7.3. Interference involves two laser beams: $\alpha$ -Beam1 departs from spacecraft SC1, passes through SC3 and SC2, and returns to SC1; $\alpha$ -Beam2 departs from SC1, passes through SC2 and SC3, and returns to SC1. Virtual interference occurs at SC1 where these two laser beams meet. The formula expressions for the $\alpha$ is given by Eq. (7.1),

\begin{split}\alpha=y_{13,13}+y_{32,3}+y_{21}-y_{31}-y_{23,2^{\prime}}-y_{12,1% ^{\prime}2^{\prime}}.\end{split}

(7.1)

Through analysis, we can express the specific process as Eq. (7.2). We need to calculate the time difference $dt$ between when these two laser beams reach SC1 after traveling their respective paths. The results for Sagnac ( $\alpha$ ) are depicted in Figure 7.4. Due to the influence of the Sagnac effect, the computed results at this stage are not meet the requirement. In subsequent analyses, we will refine the Sagnac interference type to achieve more satisfactory results.

\begin{split}\alpha-Beam1:\hskip 20.0pt\overbrace{1\overset{L_{2}}{% \longrightarrow}3\overset{L_{1}}{\longrightarrow}2\overset{L_{3}}{% \longrightarrow}1}^{y_{13,13}+y_{32,3}+y_{21}}\ ||\ t&\\ \alpha-Beam2:\hskip 20.0pt\underbrace{1\overset{L_{3^{\prime}}}{% \longrightarrow}2\overset{L_{1^{\prime}}}{\longrightarrow}3\overset{L_{2^{% \prime}}}{\longrightarrow}1}_{y_{12,1^{\prime}2^{\prime}}+y_{23,2^{\prime}}+y_% {31}}\ ||\ t&\end{split}

(7.2)

During the calculation, the time difference is obtained by subtracting the time taken by $\alpha$ -Beam2, which travels through $L_{3^{\prime},1^{\prime}2^{\prime}}+L_{1^{\prime},2^{\prime}}+L_{2^{\prime}}$ , from the time taken by $\alpha$ -Beam1, which travels through $L_{2,12}+L_{1,3}+L_{3}$ . After the light beam along $\alpha$ -Beam1 arrives at SC1 at time $t$ , we calculate backward along the negative time direction starting from time $t$ for the process of $\alpha$ -Beam2. This negative time is then added to the time taken by $\alpha$ -Beam1, resulting in the time difference $dt$ as given by Eq. (7.3). If we denote the positive direction of time with ” $\rightarrow$ ” and the negative direction with ” $\leftarrow$ ”, we can represent the interference process of Sagnac ( $\alpha$ ) as $\overrightarrow{213}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}$ , where the numbers under the arrows denote the arm labels.

\begin{split}dt&=L_{2,12}+L_{1,3}+L_{3}-(L_{2^{\prime}}+L_{1^{\prime},2}+L_{3^% {\prime},1^{\prime}2^{\prime}})\\ &=L_{2,12}+L_{1,3}+L_{3}+(-L_{2^{\prime}})+(-L_{1^{\prime},2^{\prime}})+(-L_{3% ^{\prime},1^{\prime}2^{\prime}})\end{split}

(7.3)

7.1.2 TDI with Eight Links

7.1.2.1 Unequal-Arm Michelson

Regarding the 1st-generation Michelson TDI, detailed discussions have been previously conducted, and thus will not be further elaborated upon here.

7.1.2.2 Relay

Relay involves three observables based on different spacecraft starting points, named U, V, and W, respectively. In Relay (U), four interferometric arms $L_{1},L_{1^{\prime}},L_{2^{\prime}}$ and $L_{3^{\prime}}$ are used, with interferometry employing two laser paths. The interferometric path of Relay (U) between spacecraft is illustrated in Figure 7.6. The formula expression of Relay-U is given by Eq. (7.4) [55],

U=y_{31,3^{\prime}1^{\prime}1}+y_{12,1^{\prime}1}+y_{23,1}+y_{32}-y_{12}-y_{31% ,3^{\prime}}-y_{23,2^{\prime}3^{\prime}}-y_{32,1^{\prime}2^{\prime}3^{\prime}}.

(7.4)

Beam U-Beam1 starts from SC3, passes through SC1 and SC2 back to SC3, then through $L_{1}$ to SC2; Beam U-Beam2 starts from SC3, first reaches SC2, then sequentially passes through SC3 and SC1, finally returning to SC2. The interference of two beams occurs at SC2, and the specific process is represented as Eq. (7.5). If the relative positions of the three spacecraft are represented by worldlines, and the interference paths of TDI are reflected across spacecraft, as shown in Figure 7.6, where the blue lines indicate the interference paths.

\begin{split}U-Beam1:\hskip 20.0pt\overbrace{3\overset{L_{2^{\prime}}}{% \longrightarrow}1\overset{L_{3^{\prime}}}{\longrightarrow}2\overset{L_{1^{% \prime}}}{\longrightarrow}3\overset{L_{1}}{\longrightarrow}2}^{y_{31,3^{\prime% }1^{\prime}1}+y_{12,1^{\prime}1}+y_{23,1}+y_{32}}\ ||\ t&\\ U-Beam2:\hskip 20.0pt\underbrace{3\overset{L_{1}}{\longrightarrow}2\overset{L_% {1^{\prime}}}{\longrightarrow}3\overset{L_{2^{\prime}}}{\longrightarrow}1% \overset{L_{3^{\prime}}}{\longrightarrow}2}_{y_{32,1^{\prime}2^{\prime}3^{% \prime}}+y_{23,2^{\prime}3^{\prime}}+y_{31,3^{\prime}}+y_{12}}\ ||\ t&\end{split}

(7.5)

From the path diagram of Relay (U) on the spacecraft worldlines, it can be seen that the virtual TDI path meet at SC2 at time $t$ , indicating that the two laser paths are continuous at $t$ , thereby considering the two paths of the entire interference process as continuous, enclosing an almost closed path. In the calculation process, U-Beam1 first reaches SC2 in the positive time direction at time $t$ , then returns from SC2 along U-Beam2 in the negative time direction to SC3. The resulting time difference $dt$ for Relay (U) is calculated through this process, expressed as Eq. (7.6).

\overrightarrow{{}_{3}2^{\prime}_{1}3^{\prime}_{2}1^{\prime}_{3}1}_{2}% \overleftarrow{3^{\prime}_{1}2^{\prime}_{3}1^{\prime}_{2}1_{3}}

(7.6)

The interferometric arms used are denoted by $i,i^{\prime}(i=1,2,3)$ , with subscripts indicating the spacecraft through which they pass. We select different times within the entire mission duration to simulate TDI and obtain the variation of $dt$ with time as shown in Figure 7.7.

Further analysis reveals that if we first calculate along the path of U-Beam2, moving from SC3 in the positive time direction to SC2, and then return along the path of U-Beam1 in the negative time direction back to SC3, the result obtained will be opposite in sign to the calculation where we first move along U-Beam1 and then return along U-Beam2 in reverse. This reversal in results is easy to understand, and it precisely facilitates the subsequent self-splicing construction of second-generation TDI paths.

dt[\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}\ \overleftarrow{3^{\prime}% 2^{\prime}1^{\prime}1}]=-dt[\overrightarrow{11^{\prime}2^{\prime}3^{\prime}}\ % \overleftarrow{11^{\prime}3^{\prime}2^{\prime}}]

(7.7)

7.1.2.3 Beacon

The three observables include Beacon, labeled P, Q, and R respectively depending on the starting spacecraft, involve different interference paths. The expression for the interference is given by Eq. (7.8) [55],

P=y_{13,11^{\prime}3^{\prime}}+y_{32,1^{\prime}3^{\prime}}+y_{23,3^{\prime}}+y% _{12,2}-y_{13,3^{\prime}}-y_{32,2}-y_{23,12}-y_{12,1^{\prime}12}.

(7.8)

And the analysis based on different beams is detailed in Eq. (7.9).

\begin{split}P-Beam1-1:\hskip 20.0pt\overbrace{1\overset{L_{3^{\prime}}}{% \longrightarrow}2\overset{L_{1^{\prime}}}{\longrightarrow}3\overset{L_{1}}{% \longrightarrow}2}^{y_{12,1^{\prime}12}+y_{23,12}+y_{32,2}}\ |\ 1\overset{L_{2% }}{\dashrightarrow}3\ ||\ t&\\ P-Beam1-2:\hskip 80.0pt\underbrace{1\overset{L_{2}}{\longrightarrow}3}_{y_{13,% 3^{\prime}}}\ |\ 1\overset{L_{3^{\prime}}}{\dashrightarrow}2\ ||\ t&\\ P-Beam2-1:\hskip 20.0pt\overbrace{1\overset{L_{2}}{\longrightarrow}3\overset{L% _{1}}{\longrightarrow}2\overset{L_{1^{\prime}}}{\longrightarrow}3}^{y_{13,11^{% \prime}3^{\prime}}+y_{32,1^{\prime}3^{\prime}}+y_{23,3^{\prime}}}\ |\ 1% \overset{L_{3^{\prime}}}{\dashrightarrow}2\ ||\ t&\\ P-Beam2-2:\hskip 80.0pt\underbrace{1\overset{L_{3^{\prime}}}{\longrightarrow}2% }_{y_{12,2}}\ |\ 1\overset{L_{2}}{\dashrightarrow}3\ ||\ t&\end{split}

(7.9)

In Beacon (P), four interferometric arms are used: $L_{1},L_{1^{\prime}},L_{2},$ and $L_{3^{\prime}}$ , with four laser beams employed for interference. Beam P-Beam 1-1 starts from SC1, passes through $L_{3^{\prime}}$ to SC2, then through $L_{1^{\prime}}$ to SC3, and finally through $L_{1}$ to SC2 at time $t$ , accounting for the delay $L_{2}$ . Beam P-Beam 1-2 starts from SC1 to SC3, where it arrives with a delay $L_{3^{\prime}}$ at time $t$ . Similar analysis applies to Beams P-Beam 2-1 and 2-2 as per Eqs. (7.9). From Eq. (7.9), it’s evident that Beam P-Beam 1-1 and Beam P-Beam 2-2 reach SC2 simultaneously at time $t-L_{2}$ , allowing their virtual interference at SC2. Similarly, Beam P-Beam 1-2 and Beam P-Beam 2-1 reach SC3 simultaneously at time $t-L_{3}$ , enabling their interference at SC3. Their combined results are then simulated.

The interference paths of Beacon (P) are represented on the spacecraft’s relatively separations, as depicted in Figure 7.9. When calculating the paths of TDI, for ease of processing, all beam paths can be viewed as continuous. We use ” $\rightarrow$ ” to indicate the positive direction of time and ” $\leftarrow$ ” to indicate the negative direction of time. Thus, Beam P-Beam 1-1 and 2-2 interference can be represented as $\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{\prime}}$ , and Beam P-Beam 1-2 and 2-1 interference as $\overrightarrow{2}\ \overleftarrow{1^{\prime}12}$ , where the numbers below the arrows denote the interferometric arm labels. This approach effectively addresses simultaneous interference arrival at the spacecraft during computation. These are then combined to form Eq. (7.10). For clarity, spacecraft labels are subscripted under the numbers. The computed results are shown in Figure 7.10.

\overrightarrow{{}_{1}3^{\prime}_{2}1^{\prime}_{3}1}_{2}\overleftarrow{3^{% \prime}}_{1}\overrightarrow{2}_{3}\overleftarrow{1^{\prime}_{2}1_{3}2_{1}}

(7.10)

7.1.2.4 Monitor

In Monitor, three observables are labeled E, F, and G, respectively [55]. In first channel (E), four interferometric arms $L_{1},L_{1^{\prime}},L_{2^{\prime}}$ , and $L_{3}$ are used, with four laser beams employed for interference, and the expression is given by Eq. (7.11) [55].

E=y_{32,1^{\prime}2^{\prime}}+y_{23,2^{\prime}}+y_{31}+y_{21,11^{\prime}}-y_{3% 1,1^{\prime}1}-y_{21}-y_{32,3}-y_{23,13}.

(7.11)

According to the analysis with different beams as shown in Eq. (7.12), Beam E-Beam 1-1 departs from SC2, passes through $L_{1^{\prime}}$ to reach SC3, then passes through $L_{1}$ to return to SC2, and finally arrives at SC1 via $L_{3}$ at time $t$ . Beam E-Beam 1-2 departs from SC3, passes through $L_{2}$ to reach SC1, and upon reaching SC3, it experiences a time delay of $L_{1^{\prime}}+L_{1}$ until virtual interference time $t$ . Similar analysis applies to Beams E-Beam 2-1 and 2-2 based on the following two Eqs. from (7.12).

\begin{split}E-Beam\ 1-1:\hskip 40.0pt\overbrace{2\overset{L_{1^{\prime}}}{% \longrightarrow}3\overset{L_{1}}{\longrightarrow}2\overset{L_{3}}{% \longrightarrow}1}^{y_{23,13}+y_{32,3}+y_{21}}\ ||\ t&\\ E-Beam\ 1-2:\hskip 30.0pt\underbrace{3\overset{L_{2^{\prime}}}{\longrightarrow% }1}_{y_{31,1^{\prime}1}}\ |\ 2\overset{L_{1^{\prime}}}{\dashrightarrow}\ 3% \overset{L_{1}}{\dashrightarrow}2\ ||\ t&\\ E-Beam\ 2-1:\hskip 40.0pt\overbrace{3\overset{L_{1}}{\longrightarrow}2\overset% {L_{1^{\prime}}}{\longrightarrow}3\overset{L_{2^{\prime}}}{\longrightarrow}1}^% {y_{32,1^{\prime}2^{\prime}}+y_{23,2^{\prime}}+y_{31}}\ ||\ t&\\ E-Beam\ 2-2:\hskip 30.0pt\underbrace{2\overset{L_{3}}{\longrightarrow}1}_{y_{2% 1,11^{\prime}}}\ |\ 3\overset{L_{1}}{\dashrightarrow}\ 2\ \overset{L_{1^{% \prime}}}{\dashrightarrow}3\ ||\ t&\end{split}

(7.12)

In Eq. (7.12), we observe that Beams E-Beam 1-1 and E-Beam 2-1 arrive at SC1 simultaneously at time $t$ . Therefore, we can arrange for both to interfere at SC1 at time $t$ , and Beams E-Beam 1-2 and E-Beam 2-2 to interfere at SC1 at time $t-L_{1}-L_{1^{\prime}}$ . The results of both interferences are then combined. The TDI paths are depicted in Fig. 7.12. When calculating the TDI, we use a method similar to Beacon (P), where we consider all beams as continuous for ease of calculating, progressing forward or backward in time. We use ” $\rightarrow$ ” to denote the positive direction of time and ” $\leftarrow$ ” to denote the negative direction. Thus, we can represent the interference of E-Beam 1-1 and 2-1 as $\overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}1}$ , and the interference of E-Beam 1-2 and 2-2 as $\overrightarrow{2^{\prime}}\ \overleftarrow{3}$ , connecting them to obtain Eq. (7.13). The numerical results are shown in Fig. 7.13.

\overrightarrow{{}_{2}1^{\prime}_{3}1_{2}3}_{1}\overleftarrow{2^{\prime}_{3}1^% {\prime}_{2}1}_{3}\overrightarrow{2^{\prime}}_{1}\overleftarrow{3_{2}}

(7.13)

7.2 Second-Generation TDI

For the second-generation TDI, the goal is to eliminate spacecraft with the same relative velocity. We use the conventions and starting approach in [56] in this section. The semicolons is used to define the delay symbols. In this case, we consider the variation of arm lengths with time ( $L_{l}(t)\neq L_{l^{\prime}}(t)$ ), and after the semicolon, delay factors can no longer be exchanged. From Eqs. (7.14)-(7.16), it can be seen that each delay factor $k$ generates a first-order differential term for each delay factor to its right. For example, in Eq. (7.16), $k$ generates a first-order differential term $\dot{L}_{k}(L_{m}+L_{n})$ for the two delay factors $m$ and $n$ to its right, and $m$ generates a first-order differential term $\dot{L}_{m}L_{n}$ for the delay factor $n$ to its right. For the 2nd-generation, it is possible to find suitable paths to cancel out all these first-order differential terms.

y_{ij;k}(t)\equiv y_{ij}(t-L_{k}(t))

(7.14)

\begin{split}y_{ij;km}(t)&\equiv y_{ij}(t-L_{m}(t)-L_{k}(t-L_{m}(t)))\\ &\simeq y_{ij}(t-L_{m}(t)-L_{k}(t)+\dot{L}_{k}L_{m})\\ &\simeq y_{ij,km}+\dot{y}_{ij,km}\dot{L}_{k}L_{m}\end{split}

(7.15)

\begin{split}y_{ij;kmn}(t)&\equiv y_{ij}(t-L_{n}(t)-L_{m}(t-L_{n}(t))-L_{k}[t-% L_{n}(t)-L_{m}(t-L_{n}(t))])\\ &\simeq y_{ij}(t-L_{n}-L_{m}-L_{k}+\dot{L}_{m}L_{n}+\dot{L}_{k}(L_{n}+L_{m}))% \\ &\simeq y_{ij,kmn}+\dot{y}_{ij,kmn}[\dot{L}_{k}(L_{m}+L_{n})+\dot{L}_{m}L_{n}]% \end{split}

(7.16)

We continue to use the combination of ” $\rightarrow$ ”, ” $\leftarrow$ ”, and arm length labels to represent the paths of laser interference. For 1st-generation TDI observables, when satisfying $\textbf{N}[\overrightarrow{l}]=\textbf{N}[\overleftarrow{l}]$ and returning to the starting spacecraft after completing the entire path, we can eliminate fixed unequal arm length differences [56]. We refer to this case as $L$ closed. For 2nd-generation TDI paths, based on satisfying $L$ closed, further satisfying the conditions in Eq. (7.17) (the path should at least simultaneously include $\overrightarrow{l}\overrightarrow{\dot{m}},\overleftarrow{l}\overleftarrow{% \dot{m}},\overrightarrow{l}\overleftarrow{\dot{m}},\overleftarrow{l}% \overrightarrow{\dot{m}}$ ) can eliminate the first-generation differentials of arm lengths. We refer to this case as $\dot{L}$ closed [56]. Here, $\textbf{N}[i]$ indicates the number of occurrences of $i$ in the string. The count of $l\dot{m}$ is statistically calculated as follows: each $\overrightarrow{l}$ with itself, and all $\overrightarrow{m}$ and $\overleftarrow{m}$ to its right; each $\overleftarrow{l}$ with all $\overrightarrow{m}$ and $\overleftarrow{m}$ to its right. For example, for $\overrightarrow{3^{\prime}3}\ \overleftarrow{2^{\prime}2}$ , all counted double-letter combinations include $\overrightarrow{3^{\prime}3^{\prime}},\overrightarrow{3^{\prime}3},% \overrightarrow{3^{\prime}}\overleftarrow{2^{\prime}},\overrightarrow{3^{% \prime}}\overleftarrow{2},\ \overrightarrow{33},\ \overrightarrow{3}% \overleftarrow{2^{\prime}},\ \overrightarrow{3}\overleftarrow{2}$ , and $\overleftarrow{2^{\prime}2}$ .

\textbf{N}[\overrightarrow{l}\overrightarrow{\dot{m}},\overleftarrow{l}% \overleftarrow{\dot{m}}]=\textbf{N}[\overrightarrow{l}\overleftarrow{\dot{m}},% \overleftarrow{l}\overrightarrow{\dot{m}}]\hskip 20.0pt(l,m=1,1^{\prime},2,2^{% \prime},3,3^{\prime})

(7.17)

For $L$ closed interference paths, when we appropriately connect another closed interference path at a suitable point, the newly formed TDI path is at least still $L$ closed. The term ”appropriate point” refers to the spacecraft at the junction point, which, after following the newly connected path, returns to the original junction spacecraft. During the joining process, encountering pairs like $\overrightarrow{l}\overleftarrow{l}$ or $\overleftarrow{l}\overrightarrow{l}$ ( $l=1,1^{\prime},2,2^{\prime},3,3^{\prime}$ ) that have no physical significance can be eliminated. For example, in the case of an unequal arm length Michelson TDI path ( $L$ closed):

\overrightarrow{{}_{1}3^{\prime}_{2}3_{1}2_{3}2^{\prime}}_{\textbf{1}}% \overleftarrow{3_{2}3^{\prime}_{1}2^{\prime}_{3}2_{1}}

(7.18)

We connect an interference path at SC1 as follows: $\overrightarrow{{}_{1}2_{3}2^{\prime}_{1}3^{\prime}_{2}3}_{\textbf{1}}% \overleftarrow{2^{\prime}_{3}2_{1}3_{2}3^{\prime}_{1}}$ , resulting in:

\overrightarrow{{}_{1}3^{\prime}_{2}3_{1}2_{3}2^{\prime}}_{\textbf{1}}% \overrightarrow{[2_{3}2^{\prime}_{1}3^{\prime}_{2}3}_{\textbf{1}}% \overleftarrow{2^{\prime}_{3}2_{1}3_{2}3^{\prime}]}_{\textbf{1}}\overleftarrow% {3_{2}3^{\prime}_{1}2^{\prime}_{3}2_{1}}

(7.19)

The resulting Eq. (7.19) is the previously computed second-generation unequal arm length Michelson interference, and according to (7.17), it is also $\dot{L}$ closed.

7.2.1 Self-Splicing to Construct Second-generation TDI Observables

Further analysis reveals that for a first-generation TDI, reversing the order of all arm lengths also reverses the direction of time, for example $\overrightarrow{213}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}% \Rightarrow\overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312}$ . The TDI paths before and after transformation are essentially the same, with their computed results differing by approximately a negative sign ( $dt[\overrightarrow{213}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}]\simeq% -dt[\overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312}]$ ). By appropriately splicing these two interference modes before and after transformation, a new TDI observables is constructed that is $\dot{L}$ closed, constituting a second-generation TDI observable. For first-generation TDI, while maintaining $L$ closure, we can flexibly transform them to facilitate splicing construction.

Taking Sagnac ( $\alpha$ ) type as an example, let’s explain the splicing construction process. The original Sagnac ( $\alpha$ ) interference path is $\overrightarrow{213}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}$ . After the reversal transformation, we obtain $\overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312}$ . Then, we transform $\overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312}$ to have SC1 as both the starting and ending points, resulting in paths starting and ending at SC2: $\overrightarrow{1^{\prime}2^{\prime}}\ \overleftarrow{312}\ \overrightarrow{3^% {\prime}}$ , and starting and ending at SC3: $\overrightarrow{2^{\prime}}\ \overleftarrow{312}\ \overrightarrow{3^{\prime}1^% {\prime}}$ . These results are then spliced into the original interference path corresponding to the spacecraft, resulting in three independent second-generation TDI paths named $\alpha$ 12-1, $\alpha$ 12-2, and $\alpha$ 12-3, as shown in Eq. (7.20). Sagnac ( $\alpha$ )-type:

\begin{split}\overrightarrow{213}\ \overleftarrow{2^{\prime}1^{\prime}3^{% \prime}}\overset{\text{\tiny{inversed}}}{\Longleftrightarrow}&\ % \overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312}\\ \overrightarrow{213}_{1}\overleftarrow{2^{\prime}1^{\prime}3^{\prime}}+% \overrightarrow{{}_{1}3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{312_{1}}% \Longrightarrow&\ \overrightarrow{213[3^{\prime}1^{\prime}2^{\prime}}\ % \overleftarrow{312]2^{\prime}1^{\prime}3^{\prime}}\ (\alpha\text{12-1})\\ \overrightarrow{21_{2}3}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}+% \overrightarrow{{}_{2}1^{\prime}2^{\prime}}\ \overleftarrow{312}\ % \overrightarrow{3^{\prime}_{2}}\Longrightarrow&\ \overrightarrow{21[1^{\prime}% 2^{\prime}}\ \overleftarrow{312}\ \overrightarrow{3^{\prime}]3}\ % \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}\ (\alpha\text{12-2})\\ \overrightarrow{2_{3}13}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}}+% \overrightarrow{{}_{3}2^{\prime}}\ \overleftarrow{312}\ \overrightarrow{3^{% \prime}1^{\prime}_{3}}\Longrightarrow&\ \overrightarrow{2[2^{\prime}}\ % \overleftarrow{312}\ \overrightarrow{3^{\prime}1^{\prime}]13}\ \overleftarrow{% 2^{\prime}1^{\prime}3^{\prime}}\ (\alpha\text{12-3})\end{split}

(7.20)

Using this method, based on existing first-generation TDI paths Michelson (X), Relay (U), Beacon (P), and Monitor (E), the second-generation TDI paths are constructed by splicing the original paths with their reversed counterparts. Subsequently, numerical calculations are performed on the resulting second-generation TDI observables. Michelson (X)-type:

\begin{split}\left.\begin{split}&\overrightarrow{3^{\prime}322^{\prime}}\ % \overleftarrow{33^{\prime}2^{\prime}2}\\ &\overrightarrow{22^{\prime}3^{\prime}3}\ \overleftarrow{2^{\prime}233^{\prime% }}\end{split}\right\}\Rightarrow\left\{\begin{split}\text{X16-1: }&% \overrightarrow{3^{\prime}322^{\prime}[22^{\prime}3^{\prime}3}\ \overleftarrow% {2^{\prime}233^{\prime}]33^{\prime}2^{\prime}2}\\ \text{X16-2: }&\overrightarrow{3^{\prime}322^{\prime}[3^{\prime}3}\ % \overleftarrow{2^{\prime}233^{\prime}}\ \overrightarrow{22^{\prime}}]% \overleftarrow{33^{\prime}2^{\prime}2}\\ \text{X16-3: }&\overrightarrow{22^{\prime}3^{\prime}3[22^{\prime}}\ % \overleftarrow{33^{\prime}2^{\prime}2}\ \overrightarrow{3^{\prime}3}]% \overleftarrow{2^{\prime}233^{\prime}}\end{split}\right.\end{split}

(7.21)

Relay (U)-type:

\begin{split}\left.\begin{split}&\overrightarrow{2^{\prime}3^{\prime}1^{\prime% }1}\ \overleftarrow{3^{\prime}2^{\prime}1^{\prime}1}\\ &\overrightarrow{11^{\prime}2^{\prime}3^{\prime}}\ \overleftarrow{11^{\prime}3% ^{\prime}2^{\prime}}\end{split}\right\}\Rightarrow\left\{\begin{split}\text{U1% 6-1: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}\ \overleftarrow{3^{% \prime}2^{\prime}1^{\prime}[3^{\prime}2^{\prime}}\ \overrightarrow{11^{\prime}% 2^{\prime}3^{\prime}}\ \overleftarrow{11^{\prime}]1}\\ \text{U16-2: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}\ % \overleftarrow{3^{\prime}2^{\prime}}[\overrightarrow{11^{\prime}2^{\prime}3^{% \prime}}\ \overleftarrow{11^{\prime}3^{\prime}2^{\prime}]1^{\prime}1}\\ \text{U16-3: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1[1^{\prime}2^{% \prime}3^{\prime}}\ \overleftarrow{11^{\prime}3^{\prime}2^{\prime}}\ % \overrightarrow{1}]\overleftarrow{3^{\prime}2^{\prime}1^{\prime}1}\end{split}% \right.\end{split}

(7.22)

Beacon (P)-type:

\begin{split}\left.\begin{split}&\overrightarrow{3^{\prime}1^{\prime}1}\ % \overleftarrow{3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}12}\\ &\overrightarrow{211^{\prime}}\ \overleftarrow{2}\ \overrightarrow{3^{\prime}}% \ \overleftarrow{11^{\prime}3^{\prime}}\end{split}\right\}\Rightarrow\left\{% \begin{split}\text{P16-1: }&\overrightarrow{3^{\prime}1^{\prime}1[1^{\prime}}% \ \overleftarrow{2}\ \overrightarrow{3^{\prime}}\ \overleftarrow{11^{\prime}3^% {\prime}}\ \overrightarrow{21]}\ \overleftarrow{3^{\prime}}\ \overrightarrow{2% }\ \overleftarrow{1^{\prime}12}\\ \text{P16-2: }&\overrightarrow{211^{\prime}[1}\ \overleftarrow{3^{\prime}}\ % \overrightarrow{2}\ \overleftarrow{1^{\prime}12}\ \overrightarrow{3^{\prime}1^% {\prime}]}\ \overleftarrow{2}\ \overrightarrow{3^{\prime}}\ \overleftarrow{11^% {\prime}3^{\prime}}\\ \text{P16-3: }&\overrightarrow{211^{\prime}}\ \overleftarrow{2}\ % \overrightarrow{[3^{\prime}1^{\prime}1}\ \overleftarrow{3^{\prime}}\ % \overrightarrow{2}\ \overleftarrow{1^{\prime}12]}\ \overrightarrow{3^{\prime}}% \ \overleftarrow{11^{\prime}3^{\prime}}\end{split}\right.\end{split}

(7.23)

Monitor (E)-type:

\begin{split}\left.\begin{split}&\overrightarrow{1^{\prime}13}\ \overleftarrow% {2^{\prime}1^{\prime}1}\ \overrightarrow{2^{\prime}}\ \overleftarrow{3}\\ &\overrightarrow{3}\ \overleftarrow{2^{\prime}}\ \overrightarrow{11^{\prime}2^% {\prime}}\ \overleftarrow{311^{\prime}}\ \\ \end{split}\right\}\Rightarrow\left\{\begin{split}\text{E16-1: }&% \overrightarrow{1[1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}1}\ % \overrightarrow{2^{\prime}}\ \overleftarrow{3]}\ \overrightarrow{1^{\prime}2^{% \prime}}\ \overleftarrow{311^{\prime}}\ \overrightarrow{3}\ \overleftarrow{2^{% \prime}}\\ \text{E16-2: }&\overrightarrow{1^{\prime}[11^{\prime}2^{\prime}}\ % \overleftarrow{311^{\prime}}\ \overrightarrow{3}\ \overleftarrow{2^{\prime}]}% \ \overrightarrow{13}\ \overleftarrow{2^{\prime}1^{\prime}1}\ \overrightarrow{% 2^{\prime}}\ \overleftarrow{3}\\ \text{E16-3: }&\overrightarrow{11^{\prime}2^{\prime}}[\overleftarrow{3}\ % \overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}1}\ % \overrightarrow{2^{\prime}}]\overleftarrow{311^{\prime}}\ \overrightarrow{3}\ % \overleftarrow{2^{\prime}}\end{split}\right.\end{split}

(7.24)

7.2.2 Cross-Splicing to Construct Second-Generation TDI Observables

7.2.2.1 Joint Beacon (P) and Monitor (E)

As shown in Figure 7.20, if we mirror Beacon (P) in the time direction $\bar{\text{P}}$ (shown as dashed lines), we observe that the TDI paths of $\bar{\text{P}}$ aligns with Monitor (E). With first-generation precision, $dt[\text{Beacon (P)}]=-dt[\text{Monitor (E)}]$ , hence we proceed to construct second-generation TDI paths by splicing Beacon (P) and Monitor (E). During computation, we start from the initial spacecraft and proceed clockwise, sequentially numbering spacecraft encountered as primary indices $i(i=1,2,3)$ , and using secondary labels $a,b,c$ for each spacecraft visited. In Figure 7.20, for clarity, we denote the spacecraft sequence in Monitor (E) as $i\bar{a},i\bar{b},i\bar{c}(i=1,2,3)$ .

As depicted in Figure 7.20, similar to the self-splicing method described earlier, we splice another closed interference path at appropriate spacecraft points. We attempt to splice Monitor (E) into Beacon (P) at suitable spacecraft points. Theoretically, spacecraft $i(i=1,2,3)$ in Monitor (E) can splice into any spacecraft $i$ in Beacon (P). The number of times E passes through SC1, SC2, and SC3 are 2, 3, and 3 respectively, while P passes through SC1, SC2, and SC3 2, 3, and 3 times as well. Therefore, there are 22 (4+9+9) possible connections between them. However, not all configurations are independent. Here, we employ a geometric approach to analyze and select appropriate configurations for computation.

During splicing, we shift Monitor (E)’s path on the world line to align with Beacon (P)’s path at spacecraft $i$ with matching primary indices. This indicates that when traveling along path P to the corresponding spacecraft, we switch to path E, traverse E’s closed path, return to spacecraft $i$ , and then continue along the remaining path of P. To denote the splicing point, we use the primary index of the splicing spacecraft followed by the secondary indices from P and E’s respective splicing points. For instance, $1ab$ signifies splicing at $1a$ on path P into E, then departing from $1\bar{b}$ on E after encircling it, and returning to path P. To facilitate further description of the various splicing configurations of Beacon (P) and Monitor (E) resulting in TDI paths, we denote them as PE + actual connection number + splicing points. For example, PE16-1ab indicates the splicing of P and E at SC1a on P and SC1 $\bar{\text{b}}$ on E, resulting in 16-link paths.

When several splicing configurations yield identical path structures, although described differently, their computed results at the second-generation are similar. From Figure 7.22, it can be observed that Beacon (P) and Monitor (E) form the same relative positions on the world line at 1ab, 1ba, 2ab, and 3cc, with differences in computed results being below 1%. At the second-generation, we consider their interference paths to be inherently identical. Eq. (7.25) lists several overlapping cases in this splicing process, where we select one TDI observable as independent and present it. We identify all independent second-generation TDI paths formed by P and E as shown in (7.26), with their respective computed results illustrated in Figure 7.23.

During the splicing process, we discovered that certain paths cancel out in pairs, reducing the overall number of paths traversed. For instance, paths such as PE14-2ac interference, as depicted in Figure 7.22. When path E splices from P’s 2a point into its own $2\bar{\text{c}}$ point, the final segment of E’s path, moving counterclockwise through arm 1 to return to $2\bar{\text{c}}$ in reverse time direction, immediately proceeds clockwise through arm 1 to 3a or $3\bar{\text{b}}$ , these two segments cancel each other completely in numerical computations, and have no physical meaning. Consequently, the actual interference path results in only 14 connections after cancellation.

\begin{split}&\text{PE16-1aa;}\hskip 90.0pt\text{PE16-1ab: 1ba,2ab,3cb;}\\ &\text{PE16-1bb: 2aa,2bb,2cc,3aa,3bb,3cc;}\\ &\text{PE16-2ac: 2cb,3ab,3ca;}\hskip 25.0pt\text{PE16-2ba: 3bc;}\\ &\text{PE16-2bc: 2ca,3ac,3ba. }\end{split}

(7.25)

Beacon (P) and Monitor (E)-type:

\begin{split}\left.\begin{split}&\overrightarrow{{}_{1a}3^{\prime}_{2a}1^{% \prime}_{3a}1}_{2b}\overleftarrow{3^{\prime}}_{1b}\overrightarrow{2}_{3b}% \overleftarrow{1^{\prime}_{2c}1_{3c}2_{1a}}\\ &\overrightarrow{{}_{2a}1^{\prime}_{3a}1_{2b}3}_{1a}\overleftarrow{2^{\prime}_% {3b}1^{\prime}_{2c}1}_{3c}\overrightarrow{2^{\prime}}_{1b}\overleftarrow{3_{2a% }}\end{split}\right\}\Rightarrow\left\{\begin{split}\text{PE16-1aa: }&% \overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{\prime}}\ % \overrightarrow{2}\ \overleftarrow{1^{\prime}12[2^{\prime}1^{\prime}1}\ % \overrightarrow{2^{\prime}}\ \overleftarrow{3}\ \overrightarrow{1^{\prime}13}]% \\ \text{PE16-1ab: }&\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{% \prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}12[3}\ \overrightarrow{% 1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}1}\ \overrightarrow{2^{% \prime}}]\\ \text{PE16-1bb: }&\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{% \prime}[3}\ \overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}% 1}\ \overrightarrow{2^{\prime}]2}\ \overleftarrow{1^{\prime}12}\\ \text{PE14-2ac: }&\overrightarrow{3^{\prime}}[\overleftarrow{1}\ % \overrightarrow{2^{\prime}}\ \overleftarrow{3}\ \overrightarrow{1^{\prime}13}% \ \overleftarrow{2^{\prime}\mathbf{1^{\prime}}}]\overrightarrow{\mathbf{1^{% \prime}}1}\ \overleftarrow{3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^{% \prime}12}\\ =\ &\overrightarrow{3^{\prime}}\ \overleftarrow{1}\ \overrightarrow{2^{\prime}% }\ \overleftarrow{3}\ \overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}% }\ \overrightarrow{1}\ \overleftarrow{3^{\prime}}\ \overrightarrow{2}\ % \overleftarrow{1^{\prime}12}\\ \text{PE16-2ba: }&\overrightarrow{3^{\prime}1^{\prime}1[1^{\prime}13}\ % \overleftarrow{2^{\prime}1^{\prime}1}\ \overrightarrow{2^{\prime}}\ % \overleftarrow{3]3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}12}% \\ \text{PE14-2bc: }&\overrightarrow{3^{\prime}1^{\prime}\mathbf{1}}[% \overleftarrow{\bf{1}}\ \overrightarrow{2^{\prime}}\ \overleftarrow{3}\ % \overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}}]% \overleftarrow{3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}12}\\ =\ &\overrightarrow{3^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{3}\ % \overrightarrow{1^{\prime}13}\ \overleftarrow{2^{\prime}1^{\prime}}\ % \overleftarrow{3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}12}\\ \end{split}\right.\end{split}

(7.26)

7.2.2.2 Joint Relay (U) and $\bar{\text{U}}$

Similar to Beacon (P), Relay (U) is mirrored in time to create $\bar{\text{U}}$ , as shown in Figure 7.25. For first-generation calculations, $dt[\text{U}]=-dt[\bar{\text{U}}]$ . By utilizing both paths, we can construct corresponding second-generation TDI paths. We continue to calculate in a clockwise manner for clarity, assigning sequential identifiers to each spacecraft. Relay (U) departs from spacecraft 3a, travels clockwise around U’s path, and returns to 3a; $\bar{\text{U}}$ departs from spacecraft $2\bar{\text{a}}$ , travels clockwise around $\bar{\text{U}}$ ’s path, and returns to $2\bar{\text{a}}$ .

According to the previous analysis method, U and $\bar{\text{U}}$ can have 22 (4+9+9) connections, but these 22 connections are not completely independent. The categorization of cases with inherent correlations that we analyzed is shown in Eq. (7.27). It should be noted that in the case of U $\bar{\text{U}}$ -1ab, the path formed by U and $\bar{\text{U}}$ is depicted in Figure 7.25. It appears similar to the path discussed earlier in PE-1ab (Figure 7.22), but the actual computation process for the paths is different. However, the difference in the computation results for the second-generation TDI is less than 1%, as shown in the bottom-right graph of Figure 7.26. There are 7 independent cases in the splicing of U and $\bar{\text{U}}$ , and the paths are represented as Eq. (7.28). The numerical computation results for each are shown in Figure 7.26 (where the labels in the figure use UU instead of U $\bar{\text{U}}$ ).

\begin{split}&\text{U}\bar{\text{U}}\text{16-1aa: 2cc, 3ac, 3cb;}\hskip 25.0pt% \text{U}\bar{\text{U}}\text{16-1ab: 1ba,2ab,3cc;}\\ &\text{U}\bar{\text{U}}\text{16-1bb: 2aa,2bb,3ba;}\hskip 28.0pt\text{U}\bar{% \text{U}}\text{16-2ac: 2cb,3aa,3bb;}\\ &\text{U}\bar{\text{U}}\text{16-2ba;}\hskip 100.0pt\text{U}\bar{\text{U}}\text% {16-2bc: 2ca,3bc,3ca;}\\ &\text{U}\bar{\text{U}}\text{16-3ab.}\end{split}

(7.27)

Relay (U) and $\bar{\text{U}}$ -type:

\begin{split}\left.\begin{split}&\overrightarrow{{}_{3a}2^{\prime}_{1a}3^{% \prime}_{2a}1^{\prime}_{3b}1}_{2b}\overleftarrow{3^{\prime}_{1b}2^{\prime}_{3c% }1^{\prime}_{2c}1_{3a}}\\ &\overrightarrow{{}_{2a}1^{\prime}_{3a}1_{2b}3_{1a}2}_{3b}\overleftarrow{1^{% \prime}_{2c}1_{3c}2_{1b}3_{2a}}\end{split}\right\}\Rightarrow\left\{\begin{% split}\text{U}\bar{\text{U}}\text{16-1aa: }&\overrightarrow{2^{\prime}[2}\ % \overleftarrow{1^{\prime}123}\ \overrightarrow{1^{\prime}13]3^{\prime}1^{% \prime}1}\ \overleftarrow{3^{\prime}2^{\prime}1^{\prime}1}\\ \text{U}\bar{\text{U}}\text{16-1ab: }&\overrightarrow{2^{\prime}}[% \overleftarrow{3}\ \overrightarrow{1^{\prime}132}\ \overleftarrow{1^{\prime}12% }]\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{\prime}2^{\prime}1% ^{\prime}1}\\ \text{U}\bar{\text{U}}\text{16-1bb: }&\overrightarrow{2^{\prime}3^{\prime}1^{% \prime}1}\ \overleftarrow{3^{\prime}[3}\ \overrightarrow{1^{\prime}132}\ % \overleftarrow{1^{\prime}12]2^{\prime}1^{\prime}1}\\ \text{U}\bar{\text{U}}\text{14-2ac: }&\overrightarrow{2^{\prime}3^{\prime}}[% \overleftarrow{123}\ \overrightarrow{1^{\prime}132}\ \overleftarrow{\mathbf{1^% {\prime}}}]\overrightarrow{\mathbf{1^{\prime}}1}\ \overleftarrow{3^{\prime}2^{% \prime}1^{\prime}1}\\ =\ &\overrightarrow{2^{\prime}3^{\prime}}\ \overleftarrow{123}\ % \overrightarrow{1^{\prime}1321}\ \overleftarrow{3^{\prime}2^{\prime}1^{\prime}% 1}\\ \text{U}\bar{\text{U}}\text{16-2ba: }&\overrightarrow{2^{\prime}3^{\prime}1^{% \prime}1[1^{\prime}132}\ \overleftarrow{1^{\prime}123]3^{\prime}2^{\prime}1^{% \prime}1}\\ \text{U}\bar{\text{U}}\text{14-2bc: }&\overrightarrow{2^{\prime}3^{\prime}1^{% \prime}\mathbf{1}}[\overleftarrow{\mathbf{1}23}\ \overrightarrow{1^{\prime}132% }\ \overleftarrow{1^{\prime}]3^{\prime}2^{\prime}1^{\prime}1}\\ =\ &\overrightarrow{2^{\prime}3^{\prime}1^{\prime}}\ \overleftarrow{23}\ % \overrightarrow{1^{\prime}132}\ \overleftarrow{1^{\prime}3^{\prime}2^{\prime}1% ^{\prime}1}\\ \text{U}\bar{\text{U}}\text{16-3ab: }&\overrightarrow{2^{\prime}3^{\prime}1^{% \prime}1}\ \overleftarrow{3^{\prime}2^{\prime}1^{\prime}1[1^{\prime}123}\ % \overrightarrow{1^{\prime}132]}\\ \end{split}\right.\end{split}

(7.28)

The construction of the second-generation TDI paths above involves simply splicing two closed paths. From the splicing process, it can be seen that it essentially involves shifting along the time axis. The number of relatively independent paths obtained from splicing is related to the relative positions of the two closed paths along the time axis. This means that during the splicing process of two paths, one path remains stationary along the time axis while the other path translates along the time direction. Given the conditions for splicing, the number of combinations of relative positions along the time axis determines the number of independent interference paths. In practical terms, this is easy to understand: the time difference $dt(t)$ after one complete loop is a function of time. Calculation shows that when the same paths are traversed, the more beams that are involved, the smaller the computed time differences.

7.2.2.3 TDI with more than sixteen links

For TDI paths with more than 16 connections, their forms are complex and diverse, and the number of paths grows exponentially with the number of connections $n$ . Therefore, it’s impractical to enumerate all interference paths, and the previous splicing construction method cannot showcase all of them. In this section, we apply the methods previously used to demonstrate how to construct paths with connection numbers $n>16$ .

In the case of $n=16$ using self-splicing methods, we use Relay (U) and its reverse calculated path (denoted as $-U$ ) for splicing. As an example, illustrated in Figure 7.27, two interference paths are depicted with solid lines (bottom-left path for U and top-right for $-U$ ). The calculation directions are indicated with arrows, and the labels of spacecraft on both paths are marked.

To connect U’s $1b$ and $-U$ ’s $3\bar{c}$ using one path, when reaching $1b$ on U’s path, proceed via a bridging path to $3\bar{c}$ , circle $-U$ counterclockwise once, return to $3\bar{c}$ , then return via the added connection to $1b$ , and complete the remaining path along U. This process is expressed in Eq. (7.29) as U18-1b3c. Additionally, in Figure 7.27, it’s possible to connect U’s $2c$ and $-U$ ’s $3\bar{b}$ using two connections, thereby constructing an interference path with $n=20$ , expressed as U20-2c3a in Eq. (7.29). The newly added bridging paths in both interference paths are highlighted in bold in Eq. (7.29).

\begin{split}\text{U18-1b3c: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1% }\ \overleftarrow{3^{\prime}}\ \overrightarrow{\mathbf{2}}[\overleftarrow{1^{% \prime}3^{\prime}2^{\prime}}\ \overrightarrow{11^{\prime}2^{\prime}3^{\prime}}% \ \overleftarrow{1]\mathbf{2}2^{\prime}1^{\prime}1}\\ \text{U20-2c3a: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}\ % \overleftarrow{3^{\prime}2^{\prime}1^{\prime}\mathbf{3^{\prime}}}\ % \overrightarrow{\mathbf{2}[11^{\prime}2^{\prime}3^{\prime}}\ \overleftarrow{11% ^{\prime}3^{\prime}2^{\prime}]\mathbf{2}}\ \overrightarrow{\mathbf{3^{\prime}}% }\ \overleftarrow{1}\end{split}

(7.29)

Based on our previous analysis, U18-1b3c and U20-2c3a have the same relative positions along the time axis in their closed paths. Therefore, it’s inferred that their numerical computation results at the second-generation TDI are the same, which has been verified through numerical calculations. The numerical results for both are shown in Figure 7.28.

According to the above method, examples of interference paths constructed using Beacon (P) and $-P$ for $n=18$ and $n=20$ are shown in Figure 7.30. The path expressions are given in Eq. (7.30), and the computation results are illustrated in Figure 7.31.

\begin{split}\text{P18-2b1b: }&\overrightarrow{3^{\prime}1^{\prime}1\mathbf{3}% [211^{\prime}}\ \overleftarrow{2}\ \overrightarrow{3}\ \overleftarrow{11^{% \prime}3^{\prime}]\mathbf{3}3^{\prime}}\ \overrightarrow{2}\ \overleftarrow{1^% {\prime}12}\\ \text{P20-3b3a: }&\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{% \prime}}\ \overrightarrow{2\mathbf{11^{\prime}}[11^{\prime}}\ \overleftarrow{2% }\ \overrightarrow{3^{\prime}}\ \overleftarrow{11^{\prime}3^{\prime}}\ % \overrightarrow{2}]\overleftarrow{\mathbf{1^{\prime}1}1^{\prime}12}\\ \end{split}

(7.30)

Examples of interference paths constructed using Beacon (P) and Monitor (E) for $n=18$ and $n=20$ are shown in Figure 7.30. The path expressions are given in Eq. (7.31). The computation results are illustrated in Figure 7.31.

\begin{split}\text{PE18-3b2a: }&\overrightarrow{3^{\prime}1^{\prime}1}\ % \overleftarrow{3^{\prime}}\ \overrightarrow{2\mathbf{1}[1^{\prime}13}\ % \overleftarrow{2^{\prime}1^{\prime}1}\ \overrightarrow{2^{\prime}}\ % \overleftarrow{3]\mathbf{1}1^{\prime}12}\\ \text{PE20-2c3c: }&\overrightarrow{3^{\prime}1^{\prime}1}\ \overleftarrow{3^{% \prime}}\ \overrightarrow{2}\ \overleftarrow{1^{\prime}}\ \overrightarrow{% \mathbf{32}[2^{\prime}}\ \overleftarrow{3}\ \overrightarrow{1^{\prime}13}\ % \overleftarrow{2^{\prime}1^{\prime}1]\mathbf{23}12}\\ \end{split}

(7.31)

When the connection number $n\geqslant 18$ , it is possible to have combinations of two pairs of second-generation TDI paths. Some paths are combined in a way that results in paired cancellation, reducing the total number of connections $n$ below the sum of the original connection numbers. Below are two examples. Eq. (7.32) shows the connection of four first-generation interference paths using Relay types U, -U, V, and -V.

\begin{split}\left.\begin{split}\left.\begin{split}\text{V: }&\overrightarrow{% 1^{\prime}2^{\prime}3^{\prime}3}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime% }3}\\ \text{-V: }&\overrightarrow{33^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{3% 3^{\prime}2^{\prime}1^{\prime}}\\ \end{split}\right\}\overset{\text{\tiny{splicing}}}{\Rightarrow}\text{2nd-% generation TDI}&\\ \left.\begin{split}\text{U: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}% \ \overleftarrow{3^{\prime}2^{\prime}1^{\prime}1}\\ \text{-U: }&\overrightarrow{11^{\prime}2^{\prime}3^{\prime}}\ \overleftarrow{1% 1^{\prime}3^{\prime}2^{\prime}}\end{split}\right\}\overset{\text{\tiny{% splicing}}}{\Rightarrow}\text{2nd-generation TDI}&\\ \end{split}\right\}\overset{\text{\tiny{splicing}}}{\Rightarrow}\text{2nd-% generation TDI}\\ \end{split}

(7.32)

The connections of U and -U form a second-generation TDI path, and similarly, V and -V form another second-generation TDI path. Connecting these two second-generation TDI paths results in at least a second-generation TDI path (potentially achieving higher-order TDI paths). The construction process is shown in Eq. (7.33).

\begin{split}\left.\begin{split}\text{V: }&\overrightarrow{1^{\prime}2^{\prime% }3^{\prime}3}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}3}\\ \text{U: }&\overrightarrow{2^{\prime}3^{\prime}1^{\prime}1}\ \overleftarrow{3^% {\prime}2^{\prime}1^{\prime}1}\end{split}\right\}\Rightarrow\ &\overrightarrow% {2^{\prime}3^{\prime}1^{\prime}1[3}\ \overleftarrow{2^{\prime}1^{\prime}3^{% \prime}3}\ \overrightarrow{\mathbf{1^{\prime}2^{\prime}3^{\prime}}}]% \overleftarrow{\mathbf{3^{\prime}2^{\prime}1^{\prime}}1}\\ =\ &\overrightarrow{2^{\prime}3^{\prime}1^{\prime}13}\ \overleftarrow{2^{% \prime}1^{\prime}3^{\prime}31}\\ \left.\begin{split}\Rightarrow\ &\overrightarrow{2^{\prime}3^{\prime}1^{\prime% }13}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}31}\\ \text{-V: }&\overrightarrow{33^{\prime}1^{\prime}2^{\prime}}\ \overleftarrow{3% 3^{\prime}2^{\prime}1^{\prime}}\\ \end{split}\right\}\Rightarrow\ &\overrightarrow{2^{\prime}3^{\prime}[33^{% \prime}1^{\prime}2^{\prime}}\ \overleftarrow{33^{\prime}2^{\prime}\mathbf{1^{% \prime}}}]\overrightarrow{\mathbf{1^{\prime}}13}\ \overleftarrow{2^{\prime}1^{% \prime}3^{\prime}31}\\ =\ &\overrightarrow{2^{\prime}3^{\prime}33^{\prime}1^{\prime}2^{\prime}}\ % \overleftarrow{33^{\prime}2^{\prime}}\ \overrightarrow{13}\ \overleftarrow{2^{% \prime}1^{\prime}3^{\prime}31}\\ \left.\begin{split}\Rightarrow&\overrightarrow{2^{\prime}3^{\prime}33^{\prime}% 1^{\prime}2^{\prime}}\ \overleftarrow{33^{\prime}2^{\prime}}\ \overrightarrow{% 13}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}31}\\ \text{-U: }&\overrightarrow{11^{\prime}2^{\prime}3^{\prime}}\ \overleftarrow{1% 1^{\prime}3^{\prime}2^{\prime}}\\ \end{split}\right\}\Rightarrow\ &\overrightarrow{2^{\prime}3^{\prime}33^{% \prime}1^{\prime}[11^{\prime}2^{\prime}3^{\prime}}\ \overleftarrow{11^{\prime}% 3^{\prime}\mathbf{2^{\prime}}}]\overrightarrow{\mathbf{2^{\prime}}}\ % \overleftarrow{33^{\prime}2^{\prime}}\ \overrightarrow{13}\ \overleftarrow{2^{% \prime}1^{\prime}3^{\prime}31}\\ =\ &\overrightarrow{2^{\prime}3^{\prime}33^{\prime}1^{\prime}11^{\prime}2^{% \prime}3^{\prime}}\ \overleftarrow{11^{\prime}3^{\prime}33^{\prime}2^{\prime}}% \ \overrightarrow{13}\ \overleftarrow{2^{\prime}1^{\prime}3^{\prime}31}\end{split}

(7.33)

Connections that cancel each other are highlighted in bold. Finally, we obtain a interference path with $n=22$ . The numerical computation results are shown in Figure 7.34. The process of connection on the spacecraft’s world lines is illustrated more clearly in Figure 7.32, which depicts the cancellation process of arms in U, V, and -V connections: the left panel shows the initial situation before arm cancellation, while the right panel shows the situation after arm cancellation.

For the construction of four different types of first-generation TDI paths, such as the Michelson type X and -X, and the Relay type $\bar{\text{W}}$ and - $\bar{\text{W}}$ , the connection process is illustrated in Eq. (7.34). The construction process is shown in Eq. (7.35). Finally, we obtain an interference path with $n=18$ . The numerical results are shown in Figure 7.34.

\begin{split}\left.\begin{split}\left.\begin{split}\bar{\text{W}}:\ &% \overrightarrow{2^{\prime}213}\ \overleftarrow{2^{\prime}231}\\ -\bar{\text{W}}:\ &\overrightarrow{1322^{\prime}}\ \overleftarrow{3122^{\prime% }}\\ \end{split}\right\}\overset{\text{\tiny{splicing}}}{\Rightarrow}\text{2nd-% generation TDI}&\\ \left.\begin{split}\text{X}:\ &\overrightarrow{3^{\prime}322^{\prime}}\ % \overleftarrow{33^{\prime}2^{\prime}2}\\ -\text{X}:\ &\overrightarrow{22^{\prime}3^{\prime}3}\ \overleftarrow{2^{\prime% }233^{\prime}}\end{split}\right\}\overset{\text{\tiny{splicing}}}{\Rightarrow}% \text{2nd-generation TDI}&\\ \end{split}\right\}\overset{\text{\tiny{splicing}}}{\Rightarrow}\text{2nd-% generation TDI}\\ \end{split}

(7.34)

\begin{split}\left.\begin{split}\bar{\text{W}}:\ &\overrightarrow{2^{\prime}21% 3}\ \overleftarrow{2^{\prime}231}\\ \text{X: }&\overrightarrow{3^{\prime}322^{\prime}}\ \overleftarrow{33^{\prime}% 2^{\prime}2}\end{split}\right\}\Rightarrow\ &\overrightarrow{2^{\prime}21% \mathbf{3}}[\overleftarrow{\mathbf{3}3^{\prime}2^{\prime}2}\ \overrightarrow{3% ^{\prime}\mathbf{322^{\prime}}}]\overleftarrow{\mathbf{2^{\prime}23}1}\\ =\ &\overrightarrow{2^{\prime}21}\ \overleftarrow{3^{\prime}2^{\prime}2}\ % \overrightarrow{3^{\prime}}\ \overleftarrow{1}\\ \left.\begin{split}\Rightarrow\ &\overrightarrow{2^{\prime}21}\ \overleftarrow% {3^{\prime}2^{\prime}2}\ \overrightarrow{3^{\prime}}\ \overleftarrow{1}\\ -\bar{\text{W}}:\ &\overrightarrow{1322^{\prime}}\ \overleftarrow{3122^{\prime% }}\\ \end{split}\right\}\Rightarrow\ &\overrightarrow{2^{\prime}21}\ \overleftarrow% {3^{\prime}[3122^{\prime}}\ \overrightarrow{13\mathbf{22^{\prime}}}]% \overleftarrow{\mathbf{2^{\prime}2}}\ \overrightarrow{3^{\prime}}\ % \overleftarrow{1}\\ =\ &\overrightarrow{2^{\prime}21}\ \overleftarrow{3^{\prime}3122^{\prime}}\ % \overrightarrow{133^{\prime}}\ \overleftarrow{1}\\ \left.\begin{split}\Rightarrow&\overrightarrow{2^{\prime}21}\ \overleftarrow{3% ^{\prime}3122^{\prime}}\ \overrightarrow{133^{\prime}}\ \overleftarrow{1}\\ -\text{X}:\ &\overrightarrow{22^{\prime}3^{\prime}3}\ \overleftarrow{2^{\prime% }233^{\prime}}\\ \end{split}\right\}\Rightarrow\ &\overrightarrow{2^{\prime}21}\ \overleftarrow% {3^{\prime}31\mathbf{2}}[\overrightarrow{\mathbf{2}2^{\prime}3^{\prime}3}\ % \overleftarrow{2^{\prime}233^{\prime}]2^{\prime}}\ \overrightarrow{133^{\prime% }}\ \overleftarrow{1}\\ =\ &\overrightarrow{2^{\prime}21}\ \overleftarrow{3^{\prime}31}\ % \overrightarrow{2^{\prime}3^{\prime}3}\ \overleftarrow{2^{\prime}233^{\prime}2% ^{\prime}}\ \overrightarrow{133^{\prime}}\ \overleftarrow{1}\end{split}

(7.35)

Chapter 8 Conclusions and Outlooks

In the work, we discussed the principles of time-delay interferometry, the CGC2.7 planetary ephemeris framework, and the orbit selection and optimization for the ASTROD-GW mission. We focused on the numerical simulation for the TDI in space-based GW detection. We covered the paths of first-generation TDI and the process of constructing second-generation TDI paths using geometric methods based on existing first-generation paths. Based on these constructed interference paths, extensive numerical calculations were performed for TDI in the ASTROD-GW mission. Additionally, building upon previous work, a detailed geometric analysis of TDI paths was conducted.

TDI is a crucial component in space-based GW detection, directly impacting the sensitivity of GW detection. The numerical results obtained from our studies represent a preliminary results in this research. In future work, we can further enhance our understanding by simulating various GW signals and evaluating the sensitivity of different TDI paths for GW detection. This research can be complemented by experimental efforts to identify interferometry methods that best suit practical needs.

References

[1] W.-T. Ni. Gravitational Waves, Dark Energy and Inflation. Modern Physics Letters A, 25, 922-935 (2010).
[2] 倪维斗. \CJKfamilygbsn引力波, 暗能量与暴涨. 中文研究报告 (2010).
[3] H. Poincaré. C. R. Acad. Sci. 140, 1504 (1905), and references therein.
[4] A. Einstein. Does the Inertia of a Body Depend on its Energy Content. Ann. d. Phys. 17, 891 (1905).
[5] J. Stachel. “History of Relativity” in “Twentieth Century Physics”. eds. L. M. Brown, A. Pais, and B. Pippard, Philadelphia, PA.: Institute of Physics Pub., Vol. 1, 317-319, 1995.
[6] A. Pais. Subtle is the Lord: The Science and the Life of Albert Einstein. Oxford University Press, Oxford, 1982.
[7] A. Einstein. Preuss. Akad. Wiss. Berlin, Sitzungsber. 778, 799, 831, 844 (1915).
[8] A. Einstein. Preuss. Akad. Wiss. Berlin, Sitzungsber. 688 (1916), 154 (1918).
[9] J. M. Weisberg, J. H. Taylor. Relativistic Binary Pulsar B1913+16: Thirty Years of Observations and Analysis, in Binary Radio Pulsars. Proc. Aspen Conference, ASP Conf. Series, Editors: F. A. Rasio and I. H. Stairs, 2004; arXiv:astro-ph/0407149.
[10] 倪维斗. \CJKfamilygbsn2007-2008 天文学学科发展报告引力波部分. 中国天文学会, 北京: 中国科学技术出版社, 100-104, 2008.
[11] 倪维斗. \CJKfamilygbsnASTROD 空间引力波探测优化方案与实验室研究项目建议书 (中国科学院空间科学预先研究项目之项目建议书, 研究方向, 天体号脉). 2009.
[12] W.-T. Ni. Testing Relativistic Gravity and Detecting Gravitational Waves in Space. Proceedings of in the Ninth Asia-Pacific International Conference on Gravitation and Astrophysics (ICGA9), Wuhan, June 28-July 3, Editors: J. Luo, Z.-B. Zhou, H.-C. Yeh and J.-P. Hsu, World Scientific, 2010 [arXiv:1001.0213].
[13] 倪维斗等. \CJKfamilygbsnASTROD 空间引力波探测优化方案: ASTROD-GW. 中国宇航学会深空探测技术专业委员会第六届学术年会论文集, 122-128, 2009.
[14] TAMA300 Collaboration. http://tamago.mtk.nao.ac.jp/.
[15] The GEO600 Team. http://geo600.aei.mpg.de.
[16] The LIGO Scientific Collaboration. http://www.ligo.caltech.edu/.
[17] The Virgo Collaboration. http://www.virgo.infn.it/.
[18] AdLIGO. http://www.ligo.caltech.edu/advLIGO/.
[19] AdVirgo. http://www.cascina.virgo.infn.it/advirgo/.
[20] LCGT. http://gw.icrr.u-tokyo.ac.jp/lcgt/.
[21] ET (Einstein Telescope). http://www.et.gw.eu/.
[22] LISA Study Team. LISA (Laser Interferometer Space Antenna): A Cornerstone Mission for the Observation of Gravitational Waves. ESA System and Technology Study Report, ESA-SCI, 11 (2000).
[23] W.-T. Ni. ASTROD and ASTROD I - Overview and Progress. International Journal of Modern Physics D, 17, 921 (2008); and references therein.
[24] A. Bec-Borsenberger et al. Astrodynamical Space Test of Relativity using Optical Devices ASTROD—A Proposal Submitted to ESA in Response to Call for Mission Proposals for Two Flexi-Mission F2/F3. January 31, (2000).
[25] W.-T. Ni. Super-ASTROD: Probing Primordial Gravitational Waves and Mapping the Outer Solar System. Classical and Quantum Gravity, 26, 075021 (2009), and references therein.
[26] S. Kawamura et al. The Japanese Space Gravitational Wave Antenna: DECIGO. Classical and Quantum Gravity, 23, S125 (2006).
[27] J. Crowder, N. J. Cornish. Beyond LISA: Exploring Future Gravitational Wave Missions. Physical Review D, 72, 083005 (2005); and references therein.
[28] N. Seto. Correlation Analysis of Stochastic Gravitational Wave Background Around 0.1-1 Hz. Physical Review D, 73, 063001 (2006); and references therein.
[29] 廖安琪, 倪维斗, 施宙聪. \CJKfamilygbsn激光天文动力学弱光锁相的研究. 云南天文台台刊, 3, 88-100, (2002).
[30] G. J. Dick, M. Tu, D. Strekalov et al. Optimal Phase Lock at Femtowatt Power Levels for Coherent Optical Deep-Space Transponder. IPN Progress Report, 42-175 (2008).
[31] W.-T. Ni et al. Progress in Mission Concept Study and Laboratory Development for the ASTROD (Astrodynamical Space Test of Relativity using Optical Devices). Proceedings of SPIE 3116: Small Spacecraft, Space Environments, and Instrumentation Technologies, Editors: F. A. Allahdadi, E. K. Casani, and T. D. Maclay, 105-116 (SPIE, 1997).
[32] W.-T. Ni et al. Astrodynamical Space Test of Relativity using Optical Devices. Advances in Space Research, 32, 7, 1437-1441 (2003).
[33] M. Tinto, S. V. Dhurandhar. Time-Delay Interferometry. Living Reviews in Relativity, 8, 4 (2005).
[34] V. A. Brumberg. Essential Relativistic Celestial Mechanics. Adam Hilger, England, 1991.
[35] 刘林. \CJKfamilygbsn航天器轨道理论. 北京: 国防工业出版社, 2000.
[36] D.-W. Chiou, W.-T. Ni. ASTROD Orbit Simulation and Accuracy of Relativistic Parameter Determination. Advances in Space Research, 25, 6, 1259-1262 (2000).
[37] D.-W. Chiou, W.-T. Ni. Orbit Simulation for the Determination of Relativistic and Solar-System Parameters for the ASTROD Space Mission. presented to 33rd COSPAR Scientific Assembly. Warsaw: 16-23 July, (2000); arXiv:astro-ph/0407570.
[38] 汤键仁, 倪维斗. \CJKfamilygbsnCGC2 历表框架. 云南天文台台刊, 3, 21-32 (2002).
[39] C.-J. Tang, W.-T. Ni. Asteroid Perturbations and Mass Determination for the ASTROD Space Mission. presented to 33rd COSPAR Scientific Assembly, Warsaw, 16-23 July, (2000); arXiv:astro-ph/0407570.
[40] 李广宇, 倪维斗, 田兰兰. \CJKfamilygbsnPMOE2003 行星历表框架(I)数学模型. 紫金山天文台台刊, 22, 12-31 (2003).
[41] 夏一飞, 黄天衣. \CJKfamilygbsn球面天文学. 南京: 南京大学出版社, 1995.
[42] 紫金山天文台. \CJKfamilygbsn中国天文年历. 北京: 科学出版社, 2002.
[43] 刘林, 胡松杰等. \CJKfamilygbsn航天动力学引论. 南京: 南京大学出版社, 2006.
[44] E. Bowell. The Asteroid Orbital Elements Database,
ftp://ftp.lowell.edu/pub/elgb/astorb.html.
[45] 徐萃薇, 孙绳武. \CJKfamilygbsn计算方法引论(第二版). 北京: 高等教育出版社, 2002.
[46] DE405. ftp://ssd.jpl.nasa.gov/pub/eph/planets.
[47] 门金瑞, 倪维斗, 王刚. \CJKfamilygbsnASTROD-GW 轨道设计. 天文学报, 51, 2, 198-209 (2010).
[48] J.-R. Men, W.-T. Ni, G. Wang. Design of ASTROD-GW Orbit. Chinese Astronomy and Astrophysics, 34, 434 (2010).
[49] 门金瑞. \CJKfamilygbsnASTROD-GW 任务轨道设计[D]. 指导教师: 倪维斗, 南京: 紫金山天文台, 2010.
[50] 胡中为, 萧乃园. \CJKfamilygbsn天文学教程(上册). 北京：高等教育出版社, 2003.
[51] 李广宇, 田兰兰. \CJKfamilygbsnPMO2003 行星历表框架(V)历表文件的生成和使用. 紫金山天文台台刊, 23, 160-170 (2004).
[52] X. X. Newhall. Numerical Representation of Planetary Ephemerides. Celestial Mechanics, 45, 305-310 (1989).
[53] S. M. Kopeikin. Post-Newtonian Limitations on Measurement of the PPN Parameters Caused by motion of gravitiong bodies. Monthly Notices of the Royal Astronomical Society, 399, 1539-1552 (2009).
[54] S. V. Dhurandhar et al. Time-Delay Interferometry for LISA with One Arm Dysfunctional. Classical and Quantum Gravity, 27 (2010).
[55] M. Tinto, F. B. Estabrook, J. W. Armstrong. Time delay interferometry with moving spacecraft arrays. Physical Review D, 69, 082001 (2004).
[56] M. Vallisneri. Geometric time delay interferometry. Physical Review D, 72, 042003 (2005).
[57] G. Vine et al. Experimental Demonstration of Time-Delay Interferometry for the Laser Interferometer Space Antenna. Physical Review Letters, 104, 211103 (2010).