SATDI: Simulation and Analysis for Time-Delay Interferometry

LISA is designed to employ three S/C to form a triangular constellation with an arm length of $2.5\times 10^{6}$ km, trailing the Earth by 20^∘. TAIJI is another space GW mission utilizing a LISA-like orbit with a triangular interferometer of $3\times 10^{6}$ km, leading/trailing the Earth by 20^∘ [52, 53, 54, 55, 56]. These two orbital formations could be formulated using the Clohessy-Wiltshire frame [57, 58]. In previous works, we have developed a workflow to optimize the LISA-like orbits using an ephemeris framework and calculate the mismatches of TDI combinations in dynamic cases [59, 60, 61]. In addition to LISA and TAIJI, the orbits of ASTROD-GW [62, 63] and AMIGO are calculated for free-falling trajectories [64, 65, 66, 67, 68]. All these missions are proposed to utilize TDI and their numerical orbits serve as input for our framework.

The ephemeris framework serves as a fundamental tool for orbital simulation. To calculate the orbit with sufficient accuracy, various interactions are considered, including 1) Newtonian and first-order post-Newtonian interactions among the Sun, major planets, Pluto, Moon, Ceres, Pallas, and Vesta; 2) figure interactions between the Sun, Earth, Moon, and other bodies treated as point masses; 3) perturbations from a selected set of 340 asteroids; and 4) tidal effects on the Moon caused by the Earth [69, 70]. Once the initial conditions of the S/C and celestial bodies are set at starting mission time, the numerical integration is carried out under the gravitational forces in the solar system. The coordinate time of orbit calculations is based on barycentric dynamical time (TDB), and the time difference between the proper time on the S/C and coordinate time could be calculated by including the relativistic effects as specified in [71].

These pre-calculated orbit data are essential components of simulation and analysis. While simulation can utilize the orbital data directly assuming the positions and velocities are true values for the missions, the analysis should be based on the observed orbital data with errors in principle. The orbit positioning for a deep space mission with a separation of $\sim 5\times 10^{8}$ km from Earth may have an uncertainty ranging from less than a kilometer to tens of kilometers depending on ground tracking data and inter-S/C measurements [72, 73]. Our tests show that this level of orbital positioning error has a trivial impact on the current analysis. Therefore, we temporarily ignore the errors of orbit determination during the analysis process. In this work, we select the LISA orbit previously obtained in [61], which meets the requirement of relative velocities between S/C less than 5 m/s for $2.5\times 10^{6}$ km arm length [1].

TDI constructs equivalent equal arm interferometry by combining multiple link measurements. In the current payload designs of LISA, two optical benches are deployed on each S/C, and the three measurements are gathered on each optical bench [74, 24, 18, 19, and references therein]. The interferometric measurements on S/C$i$ are illustrated in Fig. 1.

In Fig. 1, the $s_{ji}$, $\varepsilon_{ij}$, and $\tau_{ij}$ refer to the science interferometer, test-mass interferometer, and reference interferometer on each optical bench, respectively. With only considering laser noise, acceleration noise, and optical metrology noise, three interferometric measurements are expressed as

where $p_{ij}$ denotes laser noise on the optical bench of S/C$i$ pointing to S/C$j$, $n^{\mathrm{op}}_{ij}$ represents optical path noise on the S/C$i$ facing $j$, $n^{\mathrm{acc}}_{ij}$ denotes the acceleration noise from test mass on the S/C$i$ pointing to $j$, and $\mathcal{D}_{ij}$ is a time-delay operator, $\mathcal{D}_{ij}y=y(t-L_{ij})$ where $L_{ij}$ is the arm length or propagation time from S/C$i$ to $j$.

The expressions of TDI have been formulated in previous literature [4, 5, 6, 7, 9, 15, 13, 12]. As the fiducial TDI configuration, the expression of the first-generation Michelson-X could be written as

where $\eta_{ji}$ is a combination of interferometric measurements from S/C$j$ to S/C$i$ as defined in [74, 24, 18]. For the links in the counterclockwise directions ($j=$S/C2 $\rightarrow$ $i=$ S/C1, S/C3$\rightarrow$S/C2 and S/C1$\rightarrow$S/C3), the $\eta_{ji}$ are

For the links in the clockwise directions, their expressions are slightly different,

The corresponding second-generation TDI channels Michelson-X1 could be written as

Their geometries are illustrated in Fig. 2, which will guide time-delay calculations in the following subsection.

For each TDI configuration, whether first-generation or second-generation, three ordinary channels could be constructed by cyclically permuting the indexes of three S/C. Then three (quasi-)orthogonal or optimal observables, denoted as (A, E, T), can be composed from ordinary observables $(a,b,c)$ using a linear combination [75, 76],

The A and E channels are expected to be science data streams capable of effectively detecting GW, while the T channel is expected to be GW-insensitive data and to characterize instrumental noises.

Time delays are key quantities in TDI used to operate interferometric measurements as shown in Eqs. (1)-(5). Their calculations are implemented along the geometry of a TDI combination as illustrated in Fig. 2. Three vertical lines in Fig. 2 represent the trajectories of three S/C in the time direction. The horizontal separations between vertical lines indicate the spatial distances between S/C. The green and magenta lines show the two beams in the TDI link combination. The ticks on the y-axis represent the time delay with respect to the interferometric time $\tau=0$. We select the Michelson-X observable to briefly explain the procedures of calculation, and a more detailed description can be found in [77]. The calculation starts from S/C1 at the earliest time point $\tau\simeq-4L$ moving toward S/C2 at $\tau\simeq-3L$ along the green line. The propagation time and position of laser arrival are determined at S/C2. The second step continues from S/C2 to S/C1 along the path of the green line to obtain the laser travel time and position of S/C1. When the calculation along green line reaches the S/C1 at the latest time point $\tau=0$, the calculation continues along the magenta line in the backward time direction. The calculation stops at the initial sending S/C after completing a loop. Afterward, the time delays with respect to the latest time points are resolved from the result from each step, as well as the S/C positions.

During the propagation time calculation, additional relativistic time delays caused by gravitational fields are considered [78]. The time delay from the sending time $T^{s}$ at $\bm{r}_{s}$ to the receiving time $T^{r}$ at $\bm{r}_{r}$ is obtained by [79, 80],

where $R$ is the coordinate distance between the sender and receiver, $c$ is the speed of light, and $\Delta T_{\text{PN}}$ is the relativistic time delay caused by gravitational field,

where $G$ is the gravitational constant, $M$ is the mass of the celestial body, $\bf{N}_{s}$ and $\bf{N}_{r}$ are unit vectors from the celestial body to the $\bm{r}_{s}$ and $\bm{r}_{r}$ positions, and $R_{s}$ and $R_{r}$ represent the radial distances of sender and receiver from the gravitating body, respectively. In the current calculation, the leading order of relativistic time delay caused by the gravitational field of the Sun is included, the effects from other planets are disregarded. Additionally, due to the orbital motion of the S/C, an iteration process is utilized to determine the arrival time and position of the receiver,

The Chebyshev polynomial interpolation is employed to precisely determine the position of S/C at any given moment [81, 82].

In this calculation, the propagation time between S/C are computed numerically based on their orbit data. In a realistic scenario, the ranging time from sender to receiver would be estimated from measurements taken on S/C [83, 84, 85, 86, references therein]. In the simulation of the next section, the propagation time with a ranging error is taken into account to evaluate the residual laser noise in the TDI channels. For analysis of the GW signal from MBHB in Section IV, the uncertainties of arm lengths are expected to be trivial and negligible.

The diagram illustrating simulation and analysis process is presented in Fig. 3. The initial step involves selecting a mission with noise budgets and setting targeting source(s) is(are) set. In the second step, a TDI configuration is chosen, and the geometry of the link combination is generated accordingly. Subsequently, the time delays and positions of S/C are calculated by implementing algorithms in Section II. Based on the noise budgets, time-domain noises are generated for each component with a specified sample rate (In this work, we choose 4 Hz), assuming that the different noises are uncorrelated. TDI is then implemented by synthesizing these noises, and the data at each time point is obtained by using Lagrange interpolation [87, 88]. After TDI, the noise data are yielded with laser noise suppression.

On the other hand, a GW signal from MBHB is simulated using a set of chosen parameters. The GW waveform in the time-domain is generated from the waveform model SEOBNRv4HM [89]. Subsequently, TDI is applied to the waveform in the SSB coordinate with TDB, and the responded waveform is computed considering the orbital motion of detector. The waveform is then injected into the noise data and yields the simulated data containing a GW signal. In the next step, the analysis is performed in the frequency domain. The frequency-domain waveforms are generated from a targeted parameter space using the reduced order model (ROM) of SEOBNRv4HM [90]. The waveforms multiplied by the response function of TDI are used for matched filtering to identify GW signs and infer the source parameters.

There are various types of observation noises for LISA [91], and our current focus is on laser frequency noise, acceleration noise, and optical metrology noise. The dominant laser frequency noise is generated by Nd:YAG lasers with a stability of $30\ \mathrm{Hz}/\sqrt{\mathrm{Hz}}$, corresponding to an amplitude spectral density (ASD) of $\tilde{p}(f)\simeq 1\times 10^{-13}/\sqrt{\rm Hz}$ [1]. In this investigation, the first-generation and second-generation TDI Michelson observables are examined for laser noise suppressions. A ranging error between S/C is considered which includes a bias of 3 ns (0.9 m) and noise with a spectrum of $1\times 10^{-15}/f/\sqrt{\rm Hz}$ [91, 92, 93].

Clock noise is assumed to be sufficiently canceled out by using the algorithm in [94] and is therefore ignored in current simulation. Both the acceleration noise $\mathrm{S_{n,acc}}$ and optical path noise $\mathrm{S_{n,op}}$ are included with spectra

where $A_{\mathrm{acc},ij}$ and $A_{\mathrm{op},ij}$ are tunable amplitudes of noise for each optical bench, with $A_{\mathrm{acc}}=3$ and $A_{\mathrm{op}}=10$ assigned to all components in current setups.

The examinations of laser noise suppressions in Michelson configurations are performed as tests, and the PSDs of residual laser noise are presented in Fig. 4. The results for the first-generation TDI channels, (X, Y, Z), are displayed in the upper plot, and the results for the second-generation TDI channels (X1, Y1, Z1) are depicted in the lower panel. Compared to the secondary noises including acceleration noise and optical metrology noise, residual laser noises in first-generation TDI channels are moderately lower at lower frequencies and slightly higher at characteristic frequencies in the current dataset. This level of laser noise is mainly limited by the path mismatches of two beams due to the orbital dynamics. For the second-generation TDI channels, the residual laser noise is significantly lower than secondary noises, and this level of suppression is subject to the ranging errors between S/C.

The conventions for waveform simulation in SSB coordinates are illustrated in Fig. 5 [96, 95]. For a source of MBHB, its location is described by the ecliptic longitude $\lambda$ and ecliptic latitude $\beta$. The propagation vector of GW will be

The orthogonal vectors in the longitude and latitude directions are defined as

By defining the polarization angle $\psi$, the principal directions $[\mathbf{p},\mathbf{q}]$ are rotated from $[\mathbf{u},\mathbf{v}]$,

The tensors for $+$ and $\times$ polarizations can be generated,

Therefore, the GW responded in a laser link from the sender at $\mathbf{p}_{s}$ to the receiver at $\mathbf{p}_{r}$ could be expressed as

where $\mathbf{n}_{sr}$ is the unit vector from sender to receiver, and $H_{+,\times}$ represents the incoming GW signal from a source. The arrival time $t_{r}$ with the influence of GW could be obtained from

where $L_{sr}$ is propagation time in the solar systems dynamics obtained from Eqs. (7)-(9), and the last term on the right side is the effects of GW. Approximating the wave propagation time as $t(\lambda)\approx t_{s}+\lambda/c$ and position as $\mathbf{p}(\lambda)=\mathbf{p}_{s}(t_{s})+\lambda\mathbf{n}_{sr}(t_{s})$, the GW at a given time point can be expressed as

The GW response in relative frequency deviation from sender to receiver could be written as [95]

With the responded GW in a single link, the time-domain waveform in a TDI channel will be

where the first term on the right side is cumulative effects of links from one (adding) beam, and the second term denotes the combination of another beam which will be subtracted from the first term. The response in the frequency domain could be expressed as

The frequency-domain waveform of TDI will be

The frequency evolution with time is calculated to obtain the instantaneous location of the detector and its response,

In contrast to Eq. (21), an approximation is usually applied $\mathbf{p}_{s}(t_{s})\approx\mathbf{p}_{s}(t_{r})$, and the GW in the single link is written as

For the LISA mission, the velocity of the S/C in SSB coordinates will be $\sim$30 km/s, and the position of sender will deviate from the original location by $\sim$250 km during a transmission time of 8.33 s. Therefore, the approximation in Eq. (26) may introduce some inaccuracy in the calculation. We choose one source for simulation to examine the deviations between two algorithms given by Eqs. (21) and (26). The parameters for waveform simulation are shown in Table. 1.

The simulated waveforms of the second-generation TDI Michelson configuration are shown in Fig. 6. In the upper panel, the time-domain waveforms in three optimal channels are present for the last $\sim$600 before the coalescence. The x-axis is the time in TDB relative to the arrival time of the merger signal at the SSB center. And the coalescence reaches the detector earlier than at the SSB center in current case. As depicted by the waveforms, the amplitudes and phases of waveforms are modulated by the TDI. The amplitudes of waveforms in A and E are significantly suppressed around $t\simeq-350$ s, corresponding to their lowest characteristic frequency, $f_{c}=1/(4L)\simeq 0.03$ Hz. For higher frequencies, the wavelength of GW is shorter than $4L$, causing the mixed phase to cross one cycle in TDI. The different amplitudes of A and E also verify that two observables have different antenna patterns [99], despite two channels having identical averaged-sensitivities. The waveform in the T channel is noticeably canceled in low frequencies, and the amplitude is amplified with the increase of frequencies.

In the middle plot of Fig. 6, Fourier transforms are applied to the time-domain waveforms, as well as the differences between two waveforms ($\delta h_{\mathrm{TDI}}=h_{\mathrm{TDI}}-h_{\mathrm{TDI,approx}}$), generated from Eqs (21) and (26). These differences are depicted by the dotted lines, and the noise ASDs of TDI channels are also shown for comparison. The ASDs of A (blue dashed-dotted line) and E (orange dashed-dotted line) are equal and overlapped, and the T channel exhibits lower than the other two channels at lower frequencies and becomes comparable at high frequencies. For current source, the differences between two algorithms vary among three channels. In the T channel, the discrepancy could exceed its noise level within the frequency range of [1 mHz, 10 mHz]. The residual waveform in the A channel may be comparable to the noise level at certain frequencies, whereas the waveform difference in the E channel is one order of magnitude lower than its noise. The divergence levels of GW signals may vary with different sources. On the other hand, the waveforms in the frequency domain clearly show that GW signals are significantly suppressed by TDI around their characteristic frequencies ($f_{c}=\frac{n}{4L}$ for $n=1,2,3...$).

In the lower panel of Fig. 6, in addition to the predominant $(l=2,m=2)$ harmonic mode of the GW, extra modes of waveform are simulated and plotted for the A channel. The dashed lines indicate the waveforms from time-domain simulation, and the solid lines represent the waveforms obtained by using the frequency-domain model. As shown by the curves, the waveforms exhibit good agreements in terms of amplitudes. However, the waveforms from time-domain and frequency-domain approximants are best fitted at different phases/times. As a result, their match decreases when matched filtering is implemented with integrated multiple harmonics. As a compromise, only the dominant $(l=2,m=2)$ mode is injected into the noise data to perform the analysis.

The fiducial strategy for identifying GW chirp signals involves filtering the data with a pre-built template bank [100, 101, 102, 40, and references therein]. The template bank should adequately cover the targeting parameter space with a required density [103, 104, 105, 106, and references therein]. Machine learning techniques have been applied to identify the chirp signals in advanced LIGO/Virgo data [107, 108]. For ground-based GW searches, chirp signals typically last seconds to dozens of seconds in the sensitive band of detectors, and the motion of detector with the Earth could be ignored. The template bank is built to cover the intrinsic parameters of the binary sources, such as masses and spins. However, for space-borne GW detectors, chirp signals could persist for days to months in their sensitive band, and the motion of interferometer could alter the response function to GWs. Moreover, TDI introduces the additional modulation into the observed signals. Therefore, the template bank may need to include extrinsic parameters to sufficiently match the observed signals. Weaving et al. [40] generated a template bank searched for MBHB in the LISA mock data, and verified the feasibility of the signal identification strategy. Given that the template bank searches are computing expensive, we only introduce the matched filtering algorithm here.

The matched filter is an optimal algorithm for identifying a GW signal with a known waveform. The output of the filtering can be expressed as [100]

where $\tilde{d}$ is the data in the frequency domain, $\tilde{h}_{\mathrm{TDI}}$ is a waveform template, $S_{n}$ is the noise PSD of the data stream, and $t_{0}$ is a reference time which refers to coalescence time here. A normalization factor can be obtained for the template,

Then, the signal-to-noise ratio (SNR) of the matched filtering is given by

Using a frequency-domain waveforms modeled with fitted parameters in Section III.3, a matched filter is implemented to the optimal channels of the second-generation TDI Michelson, and the results are presented in the upper panel of Fig. 7. Among the three optimal channels, the E channel exhibits the highest SNR ($\rho_{E}=356$), followed by the A channel with a lower SNR ($\rho_{A}=243$). The T channel shows a significantly lower SNR ($\rho_{T}=45$) compared to the other two channels. Additionally, the SNR profile of the T channel is broader than that of other two channels, and this should be due to that its SNR is mainly contributed from lower frequencies.

An optimal SNR is defined for the template fully matching the signal in data,

The ratio $\frac{2|\tilde{h}|\sqrt{f}}{\sqrt{S}}$ between the waveform and noise ASD partially reflect the SNR contribution across a log-scaled frequency range. The ratios for Michelson configuration from both first-generation and second-generation TDI are depicted in the lower panel of Fig. 7. As expected, the E channels exhibit the highest values among three observables, and the A channel is slightly lower than E for this instance. The A/E observables from both first-generation and second-generation TDI configurations are (nearly) identical, except at characteristic frequencies. Their largest SNR contributions are expected to occur in the frequency band around 5 mHz. However, noticeable differences are observed in the T channels between two different generations. As we have previously analyzed in [99, 51], benefiting from unequal arm interferometry, the T channels from Michelson TDI are more sensitive compared to the equal-arm case. The T channel may contain GW signals resulting from imperfect cancellation of the three ordinary channels (X, Y, Z) or (X1, Y1, Z1). The second-generation observables process more recurrent links than the first-generation, amplifying the effect of unequal arm length. Consequently, the SNR in the T channel of the second-generation is higher than that of the first-generation. On the other side, in these observables, especially for the second-generation, the ratios exhibit discontinuities at their characteristic frequencies, caused by the overestimated PSDs at these frequencies due to averaging.

Source parameters are inferred from the observed data by using Bayesian analysis. The posteriors of parameters are calculated from

where $p(\vec{\theta})$ is prior probabilities of parameters. $\mathcal{L}(d|\vec{\theta})$ is the likelihood function of a set of parameters [109],

and $\mathbf{C}$ is the correlation matrix of noises in selected TDI channels,

where $T_{\mathrm{obs}}$ is the observation time which is 30 days in this investigation, $S_{xx}$/$S_{xy}$ are PSD or cross-spectral density of data streams which estimated by using algorithm detailed in [100]. For networks with multiple detectors, the joint log-likelihood is given by

The prior setups are shown in Table 2. Instead of inferring individual component masses ($m_{1},m_{2}$), the redshifted chirp mass $\mathcal{M}_{c}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ and mass ratio $q=m_{2}/m_{1}$ are estimated. The prior of chirp mass is adjustable for different mass binaries, and the prior of mass ratio is set to be a range of [0.1, 0.2]. The prior ranges of these two parameters are set to be relatively smaller to reduce computing time, and it will not affect the final results. The luminosity distance $d_{L}$ uses a distribution in the comoving volume by using the cosmological parameters from Planck 2018 [97, 98]. The inclination $\iota$ follows a sinusoidal distribution, the sky location ($\lambda,\beta$) is uniformly distributed in a spherical surface, the reference phase $\phi_{c}$ is evenly distributed in $[0,2\pi]$, and the coalescence time $t_{c}$ is uniformly sampled within $\pm$180 seconds with respect to the injected time. The Nested sampler $\mathsf{MultiNest}$ [110, 111] is employed to perform the parameter inference.

Three cases are executed for the first comparison with different TDI datasets: 1) three optimal channels (A, E, T) of first-generation Michelson, 2) three optimal channels from second-generation Michelson, and 3) two science channels from second-generation Michelson (A, E). The inferred distributions of parameters are shown in Fig. 8. The parameters inferred from three scenarios are moderately different. The distributions of the mass parameters are almost identical, while the distributions for extrinsic parameters vary with datasets. As aforementioned, the detectability of the T channel is sensitive to arm differences during the orbital motion, and the time-varying information could be averaged and discarded in the frequency-domain analysis. As a result, accurately inferring the signal in the T channel could be tricky or otherwise cause bias in estimation. As shown in the result, since the T channel from second-generation Michelson accumulates more SNR than it does in first-generation, the distributions of extrinsic parameters inferred from second-generation channels diverge from the injected values more than the first-generation, especially for the luminosity distance $d_{L}$ and longitude $\lambda$. After removing the T data stream, two science channels from the second-generation can constrain all parameters in sensible regions. As we also can see from the corner plot, the $\mathcal{M}_{c}$ and $q$ are correlated as expected. The inclination $\iota$ not only degenerates with the distance $d_{L}$ but also the ecliptic latitude $\beta$, and the longitude is correlated with arrival time $t_{c}$. Both polarization angle $\psi$ and reference phase $\phi_{c}$ are inferred with bimodal distributions.

In another case, for a LISA-like mission, the orbital motion of interferometer will modulate the amplitude and phase of waveform. Hence, the effects of orbit dynamics should be taken into account when modeling TDI waveforms in the frequency-domain, especially for a long-duration signal. To access this impact, an inference with static constellation is performed using the first-generation Michelson optimal channels. The results are presented in Fig. 9. Compared to the dynamic model, the inferred values from the static model diverge from the true values, particularly for the intrinsic parameters, chirp mass $\mathcal{M}_{c}$ and mass ratio $q$, while the extrinsic parameters are still within the $3\sigma$ credible regions. Furthermore, the extrinsic parameters estimated from dynamic model exhibit better constraints with smaller uncertainties. One reason for this could be the better matched waveforms, and another reason could be that the motion of detector helps to resolve the location of source by demodulated response function.

In principle, simulating and analyzing a variety of sources would be necessary to validate the analysis performance. However, this process would be computationally costly. To assess the analysis performance more rigorously, a more challenging case could be utilized to examine whether the parameters can be accurately inferred. As proposed in [55], the inclusion of the TAIJIm observatory in the network with LISA can significantly promote the determinations of parameters for MBHBs. By jointly analyzing simulated data for TAIJIm along with LISA, the result, shown by the blue regions in Fig. 10, indicates that all parameters can be inferred accurately with much smaller uncertainties. This suggests the correctness and robustness of analysis. Additionally, the joint data analysis also mitigates the inaccurate influence from the Michelson T channels and produces sensible distributions.