1 Introduction

Breakthroughs in modern technology have made possible the construction of extremely large gravitational wave (GW) interferometer detectors both on the ground and in space. The past six decades’ heroic experimental efforts undertaken by physicists all over the world have finally culminated with the first direct observation of a GW signal announced by the Laser Interferometer Gravitational Wave Observatory (LIGO) project (Aasi et al. 2015; Abbott et al. 2016). On September 14, 2015, the two LIGO interferometers at Hanford (Washington) and Livingston (Louisiana), simultaneously measured and recorded strain data that indicated the presence of a GW signal emitted by a coalescing binary system containing two black-holes of masses \(M_1 = 36^{+5}_{-4} \, M_{\odot }\) and \(M_2 = 29^{+4}_{-4} \, M_{\odot }\) out to a luminosity distance of \(410^{+160}_{-180}\mathrm {\ Mpc}\) corresponding to a red-shift \(z = 0.09^{+0.03}_{-0.04}\) (the above uncertainties being at the 90 percent confidence level). Since the announcement of the first GW observation, more detections have been made by both LIGO and the VIRGO project, and it is expected that soon the KAGRA interferometer in Japan (Aso et al. 2013) will join its two western sisters in making astronomical observations.

Like Galileo Galilei (Favaro and Barbera 1966) in the year 1610 was looking for the first time at the marvels of the sky with one of the first-made optical telescopes, we are now just starting to explore the observational capabilities offered by GWs, which promise to unveil secrets of the Universe inaccessible by any other means (Thorne 1987).

Ground-based interferometers cover a frequency band ranging from a few tens of Hz to a few kHz, with the lower frequency cut-off determined by the Earth’s large seismic and gravity-gradient noises (Aasi et al. 2015; Accadia et al. 2012; Aso et al. 2013). To access lower regions of the GW spectrum, the European Space Agency (ESA) and the National Aeronautics and Space Administration (NASA) have been planning, since the late nineties, to jointly fly the Laser Interferometer Space Antenna (LISA). LISA, which is now expected to be launched in the year 2034 (Amaro-Seoane et al. 2017), is a space-based interferometer with three interplanetary spacecraft flying in an almost equilateral configuration and exchanging coherent laser beams along its three arms of 2.5 million km. It is designed to cover a bandwidth from \(10^{-4}\) to 1 Hz, unveiling a broad variety of GW sources unobservable by ground-based interferometers (Amaro-Seoane et al. 2017). Ground- and space-based detectors will complement each other in the observation of GWs in an essential way, analogous to the way optical, radio, X-ray, \(\gamma \)-ray, and other frequency-band observations have been doing for the electromagnetic spectrum.

The astrophysical sources observable in the mHz band include galactic binaries, extra-galactic super-massive black-hole binaries and coalescences, and stochastic GW background from the early Universe. Coalescing binaries are one of the important sources in this frequency region, as they include galactic and extra galactic stellar mass binaries, and massive and super-massive black-hole binaries. In addition, following the recent detections made by ground detectors, it has been estimated that a large population of small-mass (10–\(100 \, M_{\odot }\)) binary black-holes could also be observable by space-based interferometers. By taking advantage of the frequency evolution of the signals emitted by these systems, both ground- and space-based interferometers could observe the same signals in their observations frequency bands, thereby complementing their scientific inferences about these signals’ sources (Sesana 2016; Tinto and de Araujo 2016). Stellar mass binaries have also been shown by population synthesis studies to be present in large number in the frequency range below 2–3 mHz (Bender and Hils 1997; Nelemans et al. 2001). In the lower frequency range (\(\le 1\) mHz) there is a large number of such unresolvable sources in each of the frequency bins. These sources effectively form a stochastic GW background referred to as binary confusion noise.

Massive black-hole binaries are interesting sources both from the astrophysical and theoretical points of view. Coalescences of massive black holes from different galaxies after their merger during growth of the present galaxies would provide new and unique information on galaxy formation. Coalescence of binaries involving intermediate mass black holes could help to understand the formation and growth of massive black holes. The super-massive black-hole binaries are strong emitters of GWs and these spectacular events can be observable beyond red-shift of \(z=10\). These systems would help to determine the cosmological parameters independently. And, just as the cosmic microwave background is left over from the big bang, so too should there be a background of gravitational waves. Unlike electromagnetic waves, gravitational waves do not interact with matter after a few Planck times since the big bang, so they do not thermalize. Their spectrum today, therefore, is simply a red-shifted version of the spectrum they formed with, which would throw light on the physical conditions at the epoch of the early Universe.

Interferometric non-resonant detectors of gravitational radiation with frequency content \(f_{\mathrm {l}}< f < f_{\mathrm {u}}\) (\(f_{\mathrm {l}} , f_{\mathrm {u}}\) being respectively the lower and upper frequency cut-offs characterizing the detector’s operational bandwidth) use a coherent train of electromagnetic waves (of nominal frequency \(\nu _0 \gg f_{\mathrm {u}}\)) folded into several beams, and at one or more points where these intersect, monitor relative fluctuations of frequency or phase (homodyne detection). The observed low-frequency fluctuations are due to several causes:

  1. 1.

    frequency variations of the source of the electromagnetic signal about its nominal frequency \(\nu _0\),

  2. 2.

    relative motions of the electromagnetic source and the mirrors (or amplifying transponders) that do the folding,

  3. 3.

    temporal variations of the index of refraction along the beams, and, according to general relativity,

  4. 4.

    to any time-variable gravitational fields present, such as the transverse-traceless metric curvature of a passing plane gravitational-wave train.

To observe gravitational waves in this way, it is thus necessary to control, or monitor, the other sources of relative frequency fluctuations, and, in the data analysis, to use optimal algorithms based on the different characteristic interferometer responses to gravitational waves (the signal) and to the other sources (the noise) (Tinto and Estabrook 1995). By comparing phases of electromagnetic beams referenced to the same frequency generator and propagated along non-parallel equal-length arms, frequency fluctuations of the frequency reference can be removed, and gravitational-wave signals at levels many orders of magnitude lower can be detected.

In the present single-spacecraft Doppler tracking observations, for instance, many of the noise sources can be either reduced or calibrated by implementing appropriate microwave frequency links and by using specialized electronics (Armstrong 2006; Tinto 2002), so the fundamental limitation is imposed by the frequency (time-keeping) fluctuations inherent to the reference clock that controls the microwave system. Hydrogen maser clocks, currently used in Doppler tracking experiments, achieve their best performance at about 1000 sec. integration time, with a fractional frequency stability of a few parts in \(10^{-16}\). This is the reason why these one-arm interferometers in space—which have one Doppler readout and a “3-pulse” response to gravitational waves (Estabrook and Wahlquist 1975)—are most sensitive to mHz gravitational waves. This integration time is also comparable to the microwave propagation (or “storage”) time 2L/c (c being the speed-of-light in vacuum) to spacecraft en route to the outer solar system (for example \(L \simeq 5 \,\text{-- }\, 8 \mathrm {\ AU}\) for the Cassini spacecraft) (Armstrong 2006; Tinto 2002).

Low-frequency interferometric gravitational-wave detectors, such as the LISA mission and the Chinese TaiJi and TianQin projects (Amaro-Seoane et al. 2017; Hu and Wu 2017; Luo et al. 2016), have been proposed to achieve greater sensitivity to mHz gravitational waves. However, since the arm lengths of these space-based interferometers can differ by a few percent because of the effects of gravity on their three spacecraft, the direct recombination of the two beams at a photo detector will not effectively remove the laser frequency noise. This is because the frequency fluctuations of the laser will be delayed by different amounts within the two arms of unequal length. To cancel the laser frequency noise, the time-varying Doppler data must be recorded and post-processed to allow for arm-length differences (Tinto and Armstrong 1999). The data streams will have temporal structure, which can be described as due to many-pulse responses to \(\delta \)-function excitations, depending on time-of-flight delays in the response functions of the instrumental Doppler noises and in the response to incident plane-parallel, transverse, and traceless gravitational waves.

Although the theory of TDI can be used by any future space-based interferometer aiming to detect gravitational radiation, this article will focus on its implementation by the LISA mission (Amaro-Seoane et al. 2017).

The LISA design envisions a constellation of three spacecraft orbiting the Sun. Each spacecraft is to be equipped with two lasers sending beams to the other two (\(\sim 0.017\) AU away) while simultaneously measuring the beat frequencies between the local laser and the laser beams received from the other two spacecraft. The theory of TDI presented in this article will assume a successful prior removal of any first-order Doppler beat-notes due to spacecraft relative motions (Hellings 2001; Tinto et al. 2002b; Heinzel et al. 2011; Otto et al. 2012; Tinto and Hartwig 2018; Tinto and Yu 2015), giving six residual Doppler time series as the raw data of a time delay space interferometer. In one of the sections covering the experimental implementation of TDI we will return to the issue of the microwave beat frequency measurements and describe the two existing schemes for removing the clock noise from the one-way heterodyne phase measurements.

Following Tinto (1998), Armstrong et al. (1999) and Dhurandhar et al. (2002), we will regard LISA not as constituting one or more conventional Michelson interferometers, but rather, in a symmetrical way, a closed array of six one-arm delay lines between the test masses. In this way, during the course of the article, we will show that it is possible to synthesize new data combinations that cancel laser frequency noises, and estimate achievable sensitivities of these combinations in terms of the separate and relatively simple single arm responses both to gravitational wave and instrumental noise (cf. Tinto 1998; Armstrong et al. 1999; Dhurandhar et al. 2002).

In contrast to Earth-based interferometers, which operate in the long-wavelength limit (LWL) (arm lengths \(\ll \) gravitational wavelength \(\sim c/f_0\), where \(f_0\) is a characteristic frequency of the GW), LISA does not operate in the LWL over much of its frequency band. When the physical scale of a free mass optical interferometer intended to detect gravitational waves is comparable to or larger than the GW wavelength, time delays in the response of the instrument to the waves, and travel times along beams in the instrument, cannot be ignored and must be allowed for in computing the detector response used for data interpretation. It is convenient to formulate the instrumental responses in terms of observed differential frequency shifts—for short, Doppler shifts—rather than in terms of phase shifts usually used in interferometry, although of course these data, as functions of time, are inter-convertible through time-differentiation.

This third revision of our article on TDI is organized as follows. In Sect. 2 we provide an overview of the physical and historical motivations of TDI. In Sect. 3 we summarize the one-arm Doppler transfer functions of an optical beam between two carefully shielded test masses inside each spacecraft resulting from (i) frequency fluctuations of the lasers used in transmission and reception, (ii) fluctuations due to non-inertial motions of the spacecraft, and (iii) beam-pointing fluctuations and shot noise (Estabrook et al. 2000). Among these, the dominant noise is from the frequency fluctuations of the lasers and is several orders of magnitude (perhaps 7 or 8) above the other noises. This noise must be very precisely removed from the data in order to achieve the GW sensitivity at the level set by the remaining Doppler noise sources, which are at a much lower level and constitute the noise floor after the laser frequency noise is suppressed. We show that this can be accomplished by shifting and linearly combining the eighteen one-way Doppler data measured by LISA. The actual procedure can easily be understood in terms of properly defined time-delay operators that act on the one-way Doppler measurements. In Sect. 4 we develop a formalism involving the algebra of the time-delay operators which is based on the theory of rings and modules and computational commutative algebra. We show that the space of all possible interferometric combinations canceling the laser frequency noise is a module over the polynomial ring in which the time-delay operators play the role of the indeterminates (Dhurandhar et al. 2002). In the literature, this module is called the first module of syzygies (Becker and Weispfenning 1993; Kreuzer and Robbiano 2000). We show that the module can be generated from four generators(Armstrong et al. 1999; Dhurandhar et al. 2002), so that any data combination canceling the laser frequency noise is simply a linear combination formed from these generators. We would like to emphasize that this is the mathematical structure underlying TDI.

Also in Sect. 4 specific interferometric combinations are derived, and their physical interpretations are discussed. The expressions for the Sagnac interferometric combinations \((\alpha , \beta , \gamma , \zeta )\) are first obtained; in particular, the symmetric Sagnac combination \(\zeta \), for which each raw data set needs to be delayed by only a single arm transit time, distinguishes itself against all the other TDI combinations by having a higher order response to gravitational radiation in the LWL when the spacecraft separations are equal. We then express the unequal-arm Michelson combinations (XYZ) in terms of the \(\alpha \), \(\beta \), \(\gamma \), and \(\zeta \) combinations with further transit time delays. One of these interferometric data combinations would still be available if the links between one pair of spacecraft were lost. Other TDI combinations, which rely on only four of the possible six inter-spacecraft Doppler measurements (denoted P, E, and U) are also presented. They would of course be quite useful in case of potential loss of any two inter-spacecraft Doppler measurements.

TDI so formulated presumes the spacecraft-to-spacecraft light-travel-times to be constant in time, and independent from being up- or down-links. Reduction of data from moving interferometric laser arrays in solar orbit will in fact encounter non-symmetric up- and down-link light time differences that are significant, and need to be accounted for in order to exactly cancel the laser frequency fluctuations (Shaddock 2004; Cornish and Hellings 2003; Shaddock et al. 2003; Tinto et al. 2004; Rajesh Nayak and Vinet 2005; Dhurandhar 2009). In Sect. 5 we show that, by introducing a set of non-commuting time-delay operators, there exists a quite general procedure for deriving generalized TDI combinations that account for the effects of time-dependence of the arms. Using this approach it is possible to derive “flex-free” expression for the unequal-arm Michelson combinations \(X_1\), and obtain the generalized expressions for all the TDI combinations (Tinto et al. 2004). Alternatively, a rigorous mathematical formulation can be given in terms of rings and modules. But because of the non-commutativity of operators the polynomial ring is non-commutative. Thus the algebraic problem becomes extremely complex and a general solution seems difficult to obtain (Dhurandhar 2009). But we show that for the special case when one arm of LISA is dysfunctional a plethora of solutions can be found (Vallisneri et al. 2008; Dhurandhar et al. 2010). Such a possibility must be envisaged because of potential technical failures.

In Sect. 6 we address the question of maximization of the LISA signal-to-noise-ratio (SNR) to any gravitational-wave signal present in its data. This is done by treating the SNR as a functional over the space of all possible TDI combinations. As a simple application of the general formula we have derived, we apply our results to the case of signals randomly polarized and randomly distributed on the celestial sphere. We find that the standard LISA sensitivity figure derived for a single Michelson interferometer (Estabrook et al. 2000; Prince et al. 2002; Rajesh Nayak et al. 2003b) can be improved by a factor of \(\sqrt{2}\) in the low-part of the frequency band, and by more than \(\sqrt{3}\) in the remaining part of the accessible band. Further, we also show that if the location of the GW source is known, as the source appears to move in the LISA reference frame, it is possible then to optimally track the source by appropriately changing the data combinations during the course of its trajectory (Prince et al. 2002; Rajesh Nayak et al. 2003a). As an example of such type of source, we consider known binaries within our own galaxy.

In Sect. 7, we finally address aspects of TDI of more practical and experimental nature, and provide a list of references where more details about these topics can be found. It is worth mentioning that, as of today, TDI has already gone through several successful experimental tests (de Vine et al. 2010; Miller 2010; Spero et al. 2011; Mitryk et al. 2012; Grüning et al. 2015) and it has been endorsed by the LISA project as its baseline technique for achieving its required sensitivity to gravitational radiation (Amaro-Seoane et al. 2017).

We emphasize that, although this article will use LISA as baseline mission reference, the results here presented can easily be extended to other space mission concepts.

2 Physical and historical motivations of TDI

Equal-arm interferometer detectors of gravitational waves can observe gravitational radiation by canceling the laser frequency fluctuations affecting the light injected into their arms. This is done by comparing phases of split beams propagated along the equal (but non-parallel) arms of the detector. The laser frequency fluctuations affecting the two beams experience the same delay within the two equal-length arms and cancel out at the photodetector where relative phases are measured. This way gravitational-wave signals of dimensionless amplitude less than \(10^{-20}\) can be observed when using lasers whose frequency stability can be as large as roughly a few parts in \(10^{-13}\).

If the arms of the interferometer have different lengths, however, the exact cancellation of the laser frequency fluctuations, say C(t), will no longer take place at the photodetector. In fact, the larger the difference between the two arms, the larger will be the magnitude of the laser frequency fluctuations affecting the detector response. If \(L_1\) and \(L_2\) are the lengths of the two arms, it is easy to see that the amount of laser relative frequency fluctuations remaining in the response is equal to (units in which the speed of light \(c = 1\))

$$\begin{aligned} \Updelta C (t) = C(t - 2L_1) - C(t - 2L_2) . \end{aligned}$$
(1)

In the case of a space-based interferometer such as LISA, whose lasers are expected to display relative frequency fluctuations equal to about \(10^{-13}/\sqrt{\mathrm {Hz}}\) in the mHz band, and whose arms will differ by a few percent (Amaro-Seoane et al. 2017), Equation (1) implies the following expression for the amplitude of the Fourier components of the uncanceled laser frequency fluctuations (an over-imposed tilde denotes the operation of Fourier transform):

$$\begin{aligned} |{\widetilde{\Updelta C}} (f)| \simeq 4 \pi f \, |(L_1 - L_2)| \, |{{\widetilde{C}}} (f)| . \end{aligned}$$
(2)

At \(f = 10^{-3} \mathrm {\ Hz}\), for instance, and assuming \(|L_1 - L_2| \simeq 0.1 \mathrm {\ s}\), the uncanceled fluctuations from the laser are equal to \(\approx \, 1.3 \times 10^{-16}/\sqrt{\mathrm {Hz}}\). Since the LISA sensitivity goal is about \(10^{-20}/\sqrt{\mathrm {Hz}}\) in this part of the frequency band, it is clear that an alternative experimental approach for canceling the laser frequency fluctuations is needed.

A first attempt to solve this problem was presented by Faller and Bender (1984) and Faller et al. (1985, 1989), and the scheme proposed there can be understood through Fig. 1. In this idealized model the two beams exiting the two arms are not made to interfere at a common photodetector. Rather, each is made to interfere with the incoming light from the laser at a photodetector, decoupling in this way the phase fluctuations experienced by the two beams in the two arms. Now two Doppler measurements are available in digital form, and the problem now becomes one of identifying an algorithm for digitally canceling the laser frequency fluctuations from a resulting new data combination.

Fig. 1
figure 1

Light from a laser is split into two beams, each injected into an arm formed by pairs of free-falling mirrors. Since the length of the two arms, \(L_1\) and \(L_2\), are different, now the light beams from the two arms are not recombined at one photo detector. Instead each is separately made to interfere with the light that is injected into the arms. Two distinct photo detectors are now used, and phase (or frequency) fluctuations are then monitored and recorded there

The algorithm they first proposed, and refined subsequently in Giampieri et al. (1996), required processing the two Doppler measurements, say \(y_1 (t)\) and \(y_2 (t)\), in the Fourier domain. If we denote with \(h_1 (t)\), \(h_2 (t)\) the gravitational-wave signals entering into the Doppler data \(y_1\), \(y_2\), respectively, and with \(n_1\), \(n_2\) any other remaining noise affecting \(y_1\) and \(y_2\), respectively, then the expressions for the Doppler observables \(y_1\), \(y_2\) can be written in the following form:

$$\begin{aligned} y_1(t)& = C(t - 2L_1) - C(t) + h_1(t) + n_1(t) , \end{aligned}$$
(3)
$$\begin{aligned} y_2(t)& = C(t - 2L_2) - C(t) + h_2(t) + n_2(t) . \end{aligned}$$
(4)

From Eqs. (3) and (4) it is important to note the characteristic time signature of the random process C(t) in the Doppler responses \(y_1\), \(y_2\). The time signature of the noise C(t) in \(y_1(t)\), for instance, can be understood by observing that the frequency of the signal received at time t contains laser frequency fluctuations transmitted \(2 L_1 \mathrm {\ s}\) earlier. By subtracting from the frequency of the received signal the frequency of the signal transmitted at time t, we also subtract the frequency fluctuations C(t) with the net result shown in Eq. (3).

The algorithm for canceling the laser noise in the Fourier domain suggested in Faller and Bender (1984) works as follows. If we take an infinitely long Fourier transform of the data \(y_1\), the resulting expression of \(y_1\) in the Fourier domain becomes [see Eq. (3)]

$$\begin{aligned} {{\widetilde{y}}}_1 (f) = {{\widetilde{C}}} (f) \left[ e^{4 \pi i f L_1} - 1\right] + {{\widetilde{h}}}_1 (f) + {{\widetilde{n}}}_1 (f) . \end{aligned}$$
(5)

If the arm length \(L_1\) is known exactly, we can use the \({\widetilde{y}}_1\) data to estimate the laser frequency fluctuations \({{\widetilde{C}}} (f)\). This can be done by dividing \({{\widetilde{y}}}_1\) by the transfer function of the laser noise C into the observable \(y_1\) itself. By then further multiplying \({{\widetilde{y}}}_1/[e^{4 \pi i f L_1} - 1]\) by the transfer function of the laser noise into the other observable \({{\widetilde{y}}}_2\), i.e., \( [e^{4 \pi i f L_2} - 1]\), and then subtract the resulting expression from \({{\widetilde{y}}}_2\) one accomplishes the cancellation of the laser frequency fluctuations.

The problem with this procedure is the underlying assumption of being able to take an infinitely long Fourier transform of the data. Even if one neglects the variation in time of the LISA arms, by taking a finite-length Fourier transform of, say, \(y_1 (t)\) over a time interval 2T, the resulting transfer function of the laser noise C into \(y_1\) no longer will be equal to \([e^{4 \pi i f L_1} - 1]\). This can be seen by writing the expression of the finite length (2T) Fourier transform of \(y_1\) in the following way:

$$\begin{aligned} {\widetilde{y}}^{T}_1 \equiv \int _{-T}^{+T} y_1(t) \, e^{2 \pi i f t} \, dt = \int _{-\infty }^{+\infty } y_1(t) \, H(t) \, e^{2 \pi i f t} \, dt , \end{aligned}$$
(6)

where we have denoted with H(t) the function that is equal to 1 in the interval \([-T, +T]\), and zero everywhere else. Equation (6) implies that the finite-length Fourier transform, \({\widetilde{y}}^{T}_1\), of \(y_1(t)\) is equal to the convolution in the Fourier domain of the infinitely long Fourier transform of \(y_1 (t)\), \({{\widetilde{y}}}_1\), with the Fourier transform of H(t) (Jenkins and Watts 1969) (i.e., the “Sinc Function” of width 1/T). The key point here is that we can no longer use the transfer function \([e^{4 \pi i f L_i} - 1]\), \(i=1, 2\), for estimating the laser noise fluctuations from one of the measured Doppler data, without retaining residual laser noise into the combination of the two Doppler data \(y_1\), \(y_2\) valid in the case of infinite integration time. The amount of residual laser noise remaining in the Fourier-based combination described above, as a function of the integration time 2T and type of “window function” used, was derived in the appendix of Tinto and Armstrong (1999). There it was shown that, in order to suppress the residual laser noise below the LISA sensitivity level identified by secondary noises (such as proof-mass and optical-path noises) with the use of the Fourier-based algorithm an integration time of about six months was needed.

A solution to this problem was suggested in Tinto and Armstrong (1999), which works entirely in the time-domain. From Eqs. (3) and (4) we may notice that, by taking the difference of the two Doppler data \(y_1(t)\), \(y_2(t)\), the frequency fluctuations of the laser now enter into this new data set in the following way:

$$\begin{aligned} 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) . \end{aligned}$$
(7)

If we now compare how the laser frequency fluctuations enter into Eq. (7) against how they appear in Eqs. (3) and (4), we can further make the following observation. If we time-shift the data \(y_1(t)\) by the round trip light time in arm 2, \(y_1(t - 2L_2)\), and subtract from it the data \(y_2(t)\) after it has been time-shifted by the round trip light time in arm 1, \(y_2(t - 2L_1)\), we obtain the following data set:

$$\begin{aligned} 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) \\&\quad + n_1(t - 2L_2) - n_2(t - 2L_1) . \end{aligned}$$
(8)

In other words, the laser frequency fluctuations enter into \(y_1(t) - y_2(t)\) and \(y_1(t - 2L_2) - y_2(t - 2L_1)\) with the same time structure. This implies that, by subtracting Eq. (8) from Eq. (7) we can generate a new data set that does not contain the laser frequency fluctuations C(t),

$$\begin{aligned} X \equiv [y_1(t) - y_2(t)] - [y_1(t - 2L_2) - y_2(t - 2L_1)] . \end{aligned}$$
(9)

The expression above of the X combination shows that it is possible to cancel the laser frequency noise in the time domain by properly time-shifting and linearly combining Doppler measurements recorded by different Doppler readouts. This in essence is what TDI amounts to.

To gain a better physical understanding of how TDI works, let’s rewrite the above X combination in the following form

$$\begin{aligned} X = [y_1(t) + y_2(t - 2L_1)] - [y_2(t) + y_1(t - 2L_2)] , \end{aligned}$$
(10)

where we have simply rearranged the terms in Eq. (9) (Shaddock et al. 2003).

Fig. 2
figure 2

Schematic diagram for X, showing that it is a synthesized zero-area Sagnac interferometer. The optical path begins at an “x” and the measurement is made at an “o”

Equation (10) shows that X is the difference of two sums of relative frequency changes, each corresponding to a specific light path (the continuous and dashed lines in Fig. 2). The continuous line, corresponding to the first square-bracket term in Eq. (10), represents a light-beam transmitted from spacecraft 1 and made to bounce once at spacecraft 3 and 2 respectively. Since the other beam (dashed line) experiences the same overall delay as the first beam (although by bouncing off spacecraft 2 first and then spacecraft 3) when they are recombined they will cancel the laser phase fluctuations exactly, having both experienced the same total delays (assuming stationary spacecraft). For this reason the combination X can be regarded as a synthesized (via TDI) zero-area Sagnac interferometer, with each beam experiencing a delay equal to \((2L_1 + 2L_{2})\). In reality, there are only two beams in each arm (one in each direction) and the lines in Fig. 2 represent the paths of relative frequency changes rather than paths of distinct light beams.

In the following sections we will further elaborate and generalize TDI to the realistic LISA configuration.

3 Time-delay interferometry

The description of TDI for LISA is greatly simplified if we adopt the notation shown in Fig. 3, where the overall geometry of the LISA detector is defined. There are three spacecraft, six optical benches, six lasers, six proof masses, and eighteen photodetectors.Footnote 1 There are also six phase difference data going clock-wise and counter-clockwise around the LISA triangle. For the moment we will make the simplifying assumption that the array is stationary, i.e., the back and forth optical paths between pairs of spacecraft are simply equal to their relative distances (Shaddock 2004; Cornish and Hellings 2003; Shaddock et al. 2003; Tinto et al. 2004).

Several notations have been used in this context. The double index notation recently employed in Shaddock et al. (2003), where six quantities are involved, is self-evident. However, when algebraic manipulations are involved the following notation seems more convenient to use. The spacecraft are labeled 1, 2, 3 and their separating distances are denoted \(L_1\), \(L_2\), \(L_3\), with \(L_i\) being opposite spacecraft i. We orient the vertices 1, 2, 3 clockwise in Fig. 3. Unit vectors between spacecraft are \({\hat{n}}_i\), oriented as indicated in Fig. 3. We index the phase difference data to be analyzed as follows: The beam arriving at spacecraft i has subscript i and is primed or unprimed depending on whether the beam is traveling clockwise or counter-clockwise (the sense defined here with reference to Fig. 3) around the LISA triangle, respectively. Thus, as seen from the figure, \(s_{1}\) is the phase difference time series measured at reception at spacecraft 1 with transmission from spacecraft 2 (along \(L_3\)).

Fig. 3
figure 3

Schematic LISA configuration. The spacecraft are labeled 1, 2, and 3. The optical paths are denoted by \(L_i\), \(L_{i}'\) where the index i corresponds to the opposite spacecraft. The unit vectors \(\hat{\mathbf {n}}_i\) point between pairs of spacecraft, with the orientation indicated

Similarly, \(s_{1^{\prime}}\) is the phase difference series derived from reception at spacecraft 1 with transmission from spacecraft 3 (along \(L_{2}\)). The other four one-way phase difference time series from signals exchanged between the spacecraft are obtained by cyclic permutation of the indices: \(1 \rightarrow 2 \rightarrow 3 \rightarrow 1\). We also adopt a notation for delayed data streams, which will be convenient later for algebraic manipulations (Dhurandhar et al. 2002). We define the three time-delay operators \({\mathscr {D}}_i\), \(i = 1, 2, 3\), where for any data stream x(t)

$$\begin{aligned} {\mathscr {D}}_i x(t) \equiv x(t - L_i) , \end{aligned}$$
(11)

where \(L_i\), \(i = 1, 2, 3\), are the light travel times along the three arms of the LISA triangle (the speed of light c is assumed to be unity in this article). Thus, for example, \({\mathscr {D}}_2 s_{1} (t) = s_{1}(t - L_2)\), \({\mathscr {D}}_2 {\mathscr {D}}_3 s_{1}(t) = s_{1}(t - L_2 - L_3) = {\mathscr {D}}_3 {\mathscr {D}}_2 s_{1}(t)\), and so on. Note that here the operators commute. This is because the arm lengths have been assumed to be constant in time. If the \(L_i\) are functions of time then the operators no longer commute (Cornish and Hellings 2003; Tinto et al. 2004), as will be described in Sect. 4. The operator notation is very appropriate, because the delays can be written as products of the operators \({\mathscr {D}}_i\) and sums of those products and then we have to deal with polynomials in \({\mathscr {D}}_i\). These polynomials are subject to the usual rules of polynomial algebra. In case of time dependent arm lengths, the operators do not commute and due care must be taken in their manipulation respecting their order in products (Dhurandhar 2009; Dhurandhar et al. 2010).

In addition to the six inter-spacecraft one-way Doppler measurements, six phase difference series result from laser beams exchanged between adjacent optical benches within each spacecraft, and six more from measuring the displacements of the optical benches relative to the proof masses. These additional data sets are similarly indexed \(\tau _{i}\), \(\tau _{i'}\), and \(\varepsilon _{i}\), \(\varepsilon _{i'}\), \(i = 1, 2, 3\) respectively.

The proof-mass-plus-optical-bench assemblies for LISA spacecraft number 1 are shown schematically in Fig. 4. The photo receivers that generate the data \(s_{1}\), \(s_{1^{\prime}}\), \(\tau _{1}\), \(\tau _{1^{\prime}}\), \(\varepsilon _{1}\) and \(\varepsilon _{1^{\prime}}\) at spacecraft 1 are shown. The phase fluctuations from the six lasers, which need to be canceled, can be represented by six random processes \(p_{i}\), \(p_{i'}\), where \(p_{i}\), \(p_{i'}\) are the phases of the lasers in spacecraft i on the left and right optical benches, respectively, as shown in the figure. Note that this notation is in the same spirit as in Tinto et al. (2004) and Shaddock et al. (2003) in which moving spacecraft arrays have been analyzed.

We extend the cyclic terminology so that at vertex i, \(i = 1, 2, 3\), the random displacement vectors of the two proof masses are respectively denoted by \({\varvec{\delta }}_{i}(t)\), \({\varvec{\delta }}_{i'}(t)\), and the random displacements (perhaps several orders of magnitude greater) of their optical benches are correspondingly denoted by \({\varvec{\varDelta }}_{i}(t)\), \({\varvec{\varDelta }}_{i'} (t)\) where the primed and unprimed indices correspond to the right and left optical benches, respectively. As pointed out in Estabrook et al. (2000), the analysis does not assume that pairs of optical benches are rigidly connected, i.e., \({\varvec{\varDelta }}_{i} \ne {\varvec{\varDelta }}_{i'}\), in general. The present LISA design shows optical fibers transmitting signals both ways between adjacent benches. We ignore time-delay effects for these signals and will simply denote by \(\mu _i(t)\) the phase fluctuations upon transmission through the fibers of the laser beams with frequencies \(\nu _{i}\), and \(\nu _{i'}\). The \(\mu _i (t)\) phase shifts within a given spacecraft might not be the same for large frequency differences \(\nu _{i} - \nu _{i'}\). For the envisioned frequency differences (a few hundred MHz), however, the remaining fluctuations due to the optical fiber can be neglected (Estabrook et al. 2000). It is also assumed that the phase noise added by the fibers is independent of the direction of light propagation through them. For ease of presentation, in what follows we will assume the center frequencies of the lasers to be the same, and denote this frequency by \(\nu _0\).

The laser phase noise in \(s_{3^{\prime}}\) is therefore equal to \({\mathscr {D}}_{1^{\prime}} p_{2}(t) - p_{3^{\prime}}(t)\). Similarly, since \(s_{2}\) is the phase shift measured on arrival at spacecraft 2 along arm 1 of a signal transmitted from spacecraft 3, the laser phase noises enter into it with the following time signature: \({\mathscr {D}}_1 p_{3^{\prime}} (t) - p_{2}(t)\). Figure 4 endeavors to make the detailed light paths for these observations clear. An outgoing light beam transmitted to a distant spacecraft is routed from the laser on the local optical bench using mirrors and beam splitters. Note that inter-spacecraft light beams do not interact with the proof masses. As a beam is received it is routed to the photo receiver where it is mixed with light from the laser on that same optical bench. The inter-spacecraft phase data are denoted \(s_{1}\) and \(s_{1^{\prime}}\) in Fig. 4.

Fig. 4
figure 4

Simplified schematic diagram of the proof-mass and optical-bench assemblies for LISA spacecraft # 1. The random displacements of the two proof masses and two optical benches are indicated as lower case \({\varvec{\delta }}_{i} , {\varvec{\delta }}_{i'}\) and upper case \({\varvec{\varDelta }}_{i} , {\varvec{\varDelta }}_{i'}\) respectively. The left bench reads out a phase signal \(s_{1}\) from optical bench \(2'\) on board spacecraft # 2. The phase difference is measured by using the laser, the photo-detector on the left optical bench, and the phasemeter (not shown) where the base-band and digitization of the one-way measurements is performed. The motion of the optical bench relative to the proof mass is measured through internal metrology and results in the time series \(\varepsilon _{1}\). The relative phase fluctuations between the laser on the optical bench 1 and the laser on the optical bench \(1'\) are instead captured by the measurements \(\tau _{1}\) and \(\tau _{1^{\prime}}\) respectively

The expressions for the \(s_{i}\), \(s_{i'}\), \(\varepsilon _{i}\), \(\varepsilon _{i'}\), and \(\tau _{i}\), \(\tau _{i'}\) phase measurements can now be developed from Figs. 3 and 4, and they are for the particular LISA configuration in which all the lasers have the same nominal frequency \(\nu _0\), and the spacecraft are stationary with respect to each other. Consider for instance the \(s_{1^{\prime}} (t)\) process [Eq. (12) below]. The photo receiver on the right bench of spacecraft 1, which (in the spacecraft reference frame) experiences a time-varying displacement \({\varvec{\varDelta }}_{1^{\prime}}\), measures the phase difference \(s_{1^{\prime}}\) by first mixing the beam from the distant optical bench 3 in direction \({\hat{n}}_2\), and laser phase noise \(p_{3}\) and optical bench motion \({\varvec{\varDelta }}_{3}\) that have been delayed by \(L_{2}'\) seconds, with the local laser light (with phase noise \(p_{1^{\prime}}\)). Since for this simplified configuration no frequency offsets are present, there is of course no need for any heterodyne conversion (Tinto et al. 2002b). Similar considerations can be made for deriving the expressions for the \(\varepsilon \)- and \(\tau \)-measurements, and they are equal to (here units are such that \(2\pi \nu _0 = 1\)) (Otto et al. 2012):

$$\begin{aligned} s_{1^{\prime}}& = H_{1^{\prime}} + {\mathscr {D}}_{2^{\prime}} p_{3} - p_{1^{\prime}} - \mathbf{{n}}_{2} \cdot ({\mathscr {D}}_{2^{\prime}} {\varvec{\varDelta }}_{3} - {\varvec{\varDelta }}_{1^{\prime}}) + N_{1^{\prime}} , \end{aligned}$$
(12)
$$\begin{aligned} \varepsilon _{1^{\prime}}& = p_{1} - p_{1^{\prime}} + 2 \, \mathbf{{n}}_{2} \cdot ({\varvec{\delta }}_{1^{\prime}} - {\varvec{\varDelta }}_{1^{\prime}}) + \mu _1 + N_{1^{\prime}}^{\varepsilon } , \end{aligned}$$
(13)
$$\begin{aligned} \tau _{1^{\prime}}& = p_{1} - p_{1^{\prime}} + \mu _1 + N_{1^{\prime}}^{\tau } , \end{aligned}$$
(14)

while those from optical bench 1 are equal to

$$\begin{aligned} s_{1}& = H_{1} + {\mathscr {D}}_{3} p_{2^{\prime}} - p_{1} + \mathbf{{n}}_{3} \cdot ({\mathscr {D}}_{3} {\varvec{\varDelta }}_{2^{\prime}} - {\varvec{\varDelta }}_{1}) + N_{1} , \end{aligned}$$
(15)
$$\begin{aligned} \varepsilon _{1}& = p_{1^{\prime}} - p_{1} - 2 \, \mathbf{{n}}_{3} \cdot ({\varvec{\delta }}_{1} - {\varvec{\varDelta }}_{1}) + \mu _{1} + N_{1}^{\varepsilon } , \end{aligned}$$
(16)
$$\begin{aligned} \tau _{1}& = p_{1^{\prime}} - p_{1} + \mu _{1} + N_{1}^{\tau } . \end{aligned}$$
(17)

In Eqs. (1217) the H-terms are the contributions to the measured phase fluctuations due to a possibly present transverse-traceless gravitational wave signal; the p-terms represent the lasers’ phase noises; the N-terms are shot-noise phase fluctuations at the photo-detectors; the \(\mathbf{{n}}\)-terms are unit vectors along the directions of propagation of the laser beams; the \({\varvec{\varDelta }}\)-terms and \({\varvec{\delta }}\)-terms are vector random processes associated with the mechanical vibrations of the optical benches and proof masses with respect to the local inertial reference frame respectively; the \(\mu \)-terms are phase fluctuations due to the optical fibers linking the two optical benches and they can been assumed to be independent of the direction of propagation of the optical beams within them (see Otto et al. 2012 for a clear discussion about this point); finally the \({\mathscr {D}}_i\), \({\mathscr {D}}_{j'}\) are delay operators. Twelve other relations, for the readouts at vertices 2 and 3, are given by cyclic permutation of the indices in Eqs. (12, 1314, 151617). The gravitational-wave phase signal components \(H_{i} , H_{i'}\) \(i = 1, 2, 3\), in Eqs. (12) and (15) are given by integrating with respect to time the Eqs. (1) and (2) of Armstrong et al. (1999), which relate metric perturbations to optical frequency shifts. The optical path phase noise contributions \(N_{i'}\), \(N_{i}\), which include shot noise from the low SNR in the links between the distant spacecraft, can be derived from the corresponding terms given in Estabrook et al. (2000). The \(\varepsilon _{i}\), \(\varepsilon _{i'}\), \(\tau _{i}\), \(\tau _{i'}\) measurements will be made with high signal-to-noise ratios so that for them the shot noise is negligible.

4 Algebraic approach for canceling laser and optical bench noises

The arms of ground-based detectors are chosen to be of equal length so that the laser light experiences identical delay in each arm of the interferometer. This arrangement precisely cancels the laser phase (or frequency) noise at the photodetector. The required sensitivity of the instrument can thus only be achieved by near exact cancellation of the laser frequency noise. However, in LISA it is impossible to achieve equal distances between spacecraft, and the laser noise cannot be canceled in this way. It is, however, possible to combine the recorded data linearly with suitable time-delays corresponding to the three arm lengths of the giant triangular interferometer so that the laser phase noise is canceled. Here we present a systematic method based on modules over polynomial rings which guarantees all the data combinations to cancel both the laser phase and the optical bench motion noises.

We first consider the simpler case of a stationary LISA, neglect the optical-bench motion noise and focus only on the laser phase noise. We do this because the algebra is somewhat simpler and the method is easier to explain. The simplification amounts to physically considering each spacecraft rigidly carrying the assembly of lasers, beam-splitters, and photodetectors. The two lasers on each spacecraft could be considered to be locked, so effectively there would be only one laser on each spacecraft. This mathematically amounts to setting \({\varvec{\varDelta }}_i = {\varvec{\varDelta }}_{i'} = 0\) and \(p_i = p_{i'}\). The scheme we describe here for laser phase noise can be extended in a straight-forward way to include optical bench motion noise, which we address in the last part of this section.

The data combinations, when only the laser phase noise is considered, consist of the six suitably delayed data streams (inter-spacecraft), the delays being integer multiples of the light travel times between spacecraft, which can be conveniently expressed in terms of polynomials in the three delay operators \({\mathscr {D}}_1\), \({\mathscr {D}}_2\), \({\mathscr {D}}_3\). The laser noise cancellation condition puts three constraints on the six polynomials of the delay operators corresponding to the six data streams. The problem, therefore, consists of finding six-tuples of polynomials which satisfy the laser noise cancellation constraints. These polynomial tuples form a module.Footnote 2 called the first module of syzygies. There exist standard methods for obtaining the module, by which we mean methods for obtaining the generators of the module so that the linear combinations of the generators generate the entire module. The procedure first consists of obtaining a Gröbner basis for the ideal generated by the coefficients appearing in the constraints. This ideal is in the polynomial ring in the variables \({\mathscr {D}}_1\), \({\mathscr {D}}_2\), \({\mathscr {D}}_3\) over the domain of rational numbers (or integers if one gets rid of the denominators). To obtain the Gröbner basis for the ideal, one may use the Buchberger algorithm or use an application such as Mathematica (Wolfram 2014). From the Gröbner basis there is a standard way to obtain a generating set for the required module. This procedure has been described in the literature (Becker and Weispfenning 1993; Kreuzer and Robbiano 2000). We thus obtain seven generators for the module. However, the method does not guarantee a minimal set and we find that a generating set of 4 polynomial six-tuples suffice to generate the required module. Alternatively, we can obtain generating sets by using the software Macaulay 2 (Grayson and Stillman 2019).

The importance of obtaining more data combinations is evident: they provide the necessary redundancy—different data combinations produce different transfer functions for GWs and the system noises so specific data combinations could be optimal for given astrophysical source parameters in the context of maximizing SNR, detection probability, improving parameter estimates, and so on and so forth. In addition, as we will show later on, there exist TDI combinations that require a reduced number of inter-spacecraft measurements, covering for the eventuality of subsystems failures.

4.1 Cancellation of laser phase noise

We now only have six data streams \(s_i\) and \(s_{i'}\), where \(i = 1,2,3\). These can be regarded as 3 component vectors \({\mathbf {s}}\) and \({\mathbf {s}}'\), respectively. The six data streams with terms containing only the laser frequency noise are

$$\begin{aligned} s_1& = {\mathscr {D}}_3 p_2 - p_1 , \\ s_{1^{\prime}}& = {\mathscr {D}}_2 p_3 - p_1 \end{aligned}$$
(18)

and their cyclic permutations. Note that we have intentionally excluded from the data additional phase fluctuations due to the GW signal and the measurement noises. Since our immediate goal is to cancel the laser frequency noise we have only kept the relevant terms. Combining the streams for canceling the laser noise will introduce transfer functions for the other noises and the GW signal. This is important and will be discussed subsequently in the article.

The goal of the analysis is to add suitably delayed beams together so that the laser frequency noise terms add up to zero.

This amounts to seeking data combinations that cancel the laser frequency noise. In the notation/formalism that we have invoked, the delay is obtained by applying the operators \({\mathscr {D}}_k\) to the beams \(s_i\) and \(s_{i'}\). A delay of \(k_1 L_1 + k_2 L_2 + k_3 L_3\) is represented by the operator \({\mathscr {D}}_1^{k_1} {\mathscr {D}}_2^{k_2} {\mathscr {D}}_3^{k_3}\) acting on the data, where \(k_1\), \(k_2\), and \(k_3\) are integers. In general, a polynomial in \({\mathscr {D}}_k\), which is a polynomial in three variables, applied to, say, \(s_1\) combines the same data stream \(s_1(t)\) with different time-delays of the form \(k_1 L_1 + k_2 L_2 + k_3 L_3\). This notation conveniently rephrases the problem. One must find six polynomials say \(q_i ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3)\), \(q_{i'} ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3)\), \(i=1,2,3 \), such that

$$\begin{aligned} \sum _{i = 1}^3 q_i s_i + q_{i'} s_{i'} = 0 . \end{aligned}$$
(19)

The zero on the right-hand side of the above equation signifies zero laser phase noise.

It is useful to express Eq. (18) in matrix form. This allows us to obtain a matrix operator equation whose solutions are \({\mathbf {q}}\) and \({\mathbf {q}}'\), where \(q_i\) and \(q_{i'}\) are written as column vectors. We can similarly express \(s_i\), \(s_{i'}\), \(p_i\) as column vectors \({\mathbf {s}}\), \({\mathbf {s}}'\), \({\mathbf {p}}\), respectively. In matrix form Eq. (18) becomes

$$\begin{aligned} {\mathbf {s}}= {\mathbf {D}}^T \cdot {\mathbf {p}}, \qquad {\mathbf {s}}' = {\mathbf {D}}\cdot {\mathbf {p}}, \end{aligned}$$
(20)

where \({\mathbf {D}}\) is a \(3 \times 3\) matrix given by

$$\begin{aligned} {\mathbf {D}}= \left( \begin{array}{ccc} -1 &{} 0 &{} {\mathscr {D}}_2 \\ {\mathscr {D}}_3 &{} -1 &{} 0 \\ 0 &{} {\mathscr {D}}_1 &{} -1 \end{array} \right) . \end{aligned}$$
(21)

The exponent ’T’ represents the transpose of the matrix. Equation (19) becomes

$$\begin{aligned} {\mathbf {q}}^T \cdot {\mathbf {s}}+ {\mathbf {q}}'^T \cdot {\mathbf {s}}' = ({\mathbf {q}}^T \cdot {\mathbf {D}}^T + {\mathbf {q}}'^T \cdot {\mathbf {D}}) \cdot {\mathbf {p}}= 0 , \end{aligned}$$
(22)

where we have taken care to put \({\mathbf {p}}\) on the right-hand side of the operators. Since the above equation must be satisfied for an arbitrary vector \({\mathbf {p}}\), we obtain a matrix equation for the polynomials \(({\mathbf {q}}, {\mathbf {q}}')\):

$$\begin{aligned} {\mathbf {q}}^T \cdot {\mathbf {D}}^T + {\mathbf {q}}' \cdot {\mathbf {D}}= 0 . \end{aligned}$$
(23)

Note that since the \({\mathscr {D}}_k\) commute, the order in writing these operators is unimportant. In mathematical terms, the polynomials form a commutative ring.

4.2 Cancellation of laser phase noise in the unequal-arm interferometer

The use of commutative algebra is very conveniently illustrated with the help of the simpler example of the unequal-arm interferometer. Here there are only two arms instead of three as we have for LISA, and the mathematics is much simpler and so it is easy to see both physically and mathematically how commutative algebra can be applied to this problem of laser phase noise cancellation. The procedure is well known for the unequal-arm interferometer, but here we will describe the same method in terms of the delay operators that we have introduced.

Let \(\phi (t)\) denote the laser phase noise entering the arms as shown in Fig. 5. Consider the laser noise \(\phi (t)\) making a round trip around arm 1 whose length we take to be \(L_1\). If we interfere this returning light with the one entering into the arm we get the following two-way phase difference, \(\phi _1 (t)\)

$$\begin{aligned} \phi _1 (t) = \phi (t - 2 L_1) - \phi (t) \equiv ({\mathscr {D}}_1^2 - 1) \, \phi (t) . \end{aligned}$$
(24)

The second expression we have written in terms of the delay operators. This makes the procedure transparent as we shall see. We can do the same for the arm 2 to get the other two-way phase difference, \(\phi _2 (t)\)

$$\begin{aligned} \phi _2 (t) = \phi (t - 2 L_2) - \phi (t) \equiv ({\mathscr {D}}_2^2 - 1) \, \phi (t) . \end{aligned}$$
(25)

Clearly, if \(L_1 \ne L_2\), then the difference in phase \(\phi _2 (t) - \phi _1 (t)\) is not zero and the laser phase noise does not cancel out. However, if one further delays the phases \(\phi _1 (t)\) and \(\phi _2 (t)\) and constructs the following combination,

$$\begin{aligned} X(t) = [\phi _2 (t - 2 L_1) - \phi _2 (t)] - [\phi _1 (t - 2 L_2) - \phi _1 (t)] , \end{aligned}$$
(26)

then the laser phase noise does cancel out. We have already encountered this combination at the end of Sect. 2. It was first proposed by Tinto and Armstrong (1999).

Fig. 5
figure 5

Schematic diagram of the unequal-arm Michelson interferometer. The beam shown corresponds to the term \(({\mathscr {D}}_2^2 - 1)({\mathscr {D}}_1^2 - 1) \phi (t)\) in X(t) which is first sent around arm 1 followed by arm 2. The second beam (not shown) is first sent around arm 2 and then through arm 1. The difference in these two beams constitutes X(t)

The cancellation of laser frequency noise becomes obvious from the operator algebra in the following way. In the operator notation,

$$\begin{aligned} X(t)& = ({\mathscr {D}}_1^2 - 1) \, \phi _2 (t) - ({\mathscr {D}}_2^2 - 1) \, \phi _1 (t) \\& = [({\mathscr {D}}_1^2 - 1)({\mathscr {D}}_2^2 - 1) - ({\mathscr {D}}_2^2 - 1)({\mathscr {D}}_1^2 - 1)] \, \phi (t) \\& = 0 . \end{aligned}$$
(27)

From this expression one immediately sees that just the commutativity of the operators has been used to cancel the laser phase noise. The basic idea was to compute the lowest common multiple (L.C.M.) of the polynomials \({\mathscr {D}}_1^2 - 1\) and \({\mathscr {D}}_2^2 - 1\) (in this case the L.C.M. is just the product, because the polynomials are relatively prime) and use this fact to construct X(t) in which the laser phase noise is canceled. The operation is shown physically in Fig. 5.

The notions of commutativity of polynomials, L.C.M., and related operations belong to the field of commutative algebra. In fact we will be using the notion of a Gröbner basis, which is in a sense the generalization of the notion of the greatest common divisor (GCD). Since LISA has three spacecraft and six inter-spacecraft beams, the problem of the unequal-arm interferometer only gets technically more complex; in principle the problem is the same as in this simpler case. Thus, the simple operations that were performed here to obtain a laser noise free combination X(t) are not sufficient and more sophisticated methods need to be adopted from the field of commutative algebra. We address this problem in the forthcoming text.

4.3 The module of syzygies

Equation (23) has non-trivial solutions. Several solutions have been exhibited in Armstrong et al. (1999) and Estabrook et al. (2000). We merely mention these solutions here; in the forthcoming text we will discuss them in detail. The solution \(\zeta \) is given by \(-{\mathbf {q}}^T = {\mathbf {q}}'^T = ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3)\). The solution \(\alpha \) is described by \({\mathbf {q}}^T = -(1, {\mathscr {D}}_3, {\mathscr {D}}_1 {\mathscr {D}}_3)\) and \({\mathbf {q}}'^T = (1, {\mathscr {D}}_1 {\mathscr {D}}_2, {\mathscr {D}}_2)\). The solutions \(\beta \) and \(\gamma \) are obtained from \(\alpha \) by cyclically permuting the indices of \({\mathscr {D}}_k\), \({\mathbf {q}}\), and \({\mathbf {q}}'\). These solutions are important, because they consist of polynomials with lowest possible degrees and thus are simple. Other solutions containing higher degree polynomials can be generated conveniently from these solutions. Since the system of equations is linear, linear combinations of these solutions are also solutions to Eq. (23).

However, it is important to realize that we do not have a vector space here. Three independent constraints on a six-tuple do not produce a space which is necessarily generated by three basis elements. This conclusion would follow if the solutions formed a vector space but they do not. The polynomial six-tuple \({\mathbf {q}}\), \({\mathbf {q}}'\) can be multiplied by polynomials in \({\mathscr {D}}_1\), \({\mathscr {D}}_2\), \({\mathscr {D}}_3\) (scalars) which do not form a field—they form a ring. The multiplicative inverse in general does not exist in a ring. We have in fact a ring of polynomials in the operators \({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3\). We, therefore, have a module over this ring of polynomials and not a vector space. It is called the first module of syzygies. First we present the general methodology for obtaining the solutions to Eq. (23) and then apply it to Eq. (23).

There are three linear constraints on the polynomials given by Eq. (23). Since the equations are linear, the solutions space is a submodule of the module of six-tuples of polynomials. The module of six-tuples is a free module, i.e., it has six basis elements that not only generate the module but are linearly independent. A natural choice of the basis is \(f_m = (0, \dots , 1, \dots , 0)\) with 1 in the m-th place and 0 everywhere else; m runs from 1 to 6. The definitions of generation (spanning) and linear independence are the same as that for vector spaces. A free module is essentially like a vector space. But our interest lies in its submodule which need not be free and need not have just three generators as it would seem if we were dealing with vector spaces.

The problem at hand is of finding the generators of this submodule, i.e., any element of the submodule should be expressible as a linear combination of the generating set. In this way the generators are capable of spanning the full submodule or generating the submodule. In order to achieve our goal, we rewrite Eq. (23) explicitly component-wise:

$$\begin{aligned} q_1 + q_{1^{\prime}} - {\mathscr {D}}_3 q_{2^{\prime}} - {\mathscr {D}}_2 q_3& = 0 , \\ q_2 + q_{2^{\prime}} - {\mathscr {D}}_1 q_{3^{\prime}} - {\mathscr {D}}_3 q_1& = 0 , \\ q_3 + q_{3^{\prime}} - {\mathscr {D}}_2 q_{1^{\prime}} - {\mathscr {D}}_1 q_2& = 0 . \end{aligned}$$
(28)

The first step is to use Gaussian elimination to obtain \(q_1\) and \(q_2\) in terms of \(q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}}\),

$$\begin{aligned} q_1& = - q_{1^{\prime}} + {\mathscr {D}}_3 q_{2^{\prime}} + {\mathscr {D}}_2 q_3 , \\ q_2& = - q_{2^{\prime}} + {\mathscr {D}}_1 q_{3^{\prime}} + {\mathscr {D}}_3 q_1 , \\& = - {\mathscr {D}}_3 q_{1^{\prime}} - (1 - {\mathscr {D}}_3^2) q_{2^{\prime}} + {\mathscr {D}}_1 q_{3^{\prime}} + {\mathscr {D}}_2 {\mathscr {D}}_3 q_3 , \end{aligned}$$
(29)

and then substitute these values in the third equation to obtain a linear implicit relation between \(q_3\), \(q_{1^{\prime}}\), \(q_{2^{\prime}}\), \(q_{3^{\prime}}\). We then have:

$$\begin{aligned} (1 - {\mathscr {D}}_1 {\mathscr {D}}_2 {\mathscr {D}}_3) q_3 + ({\mathscr {D}}_1 {\mathscr {D}}_3 - {\mathscr {D}}_2) q_{1^{\prime}} + {\mathscr {D}}_1 (1 - {\mathscr {D}}_3^2) q_{2^{\prime}} + (1 - {\mathscr {D}}_1^2) q_{3^{\prime}} = 0 . \end{aligned}$$
(30)

Obtaining solutions to Eq. (30) amounts to solving the problem since the remaining polynomials \(q_1\), \(q_2\) have been expressed in terms of \(q_3\), \(q_{1^{\prime}}\), \(q_{2^{\prime}}\), \(q_{3^{\prime}}\) in Eq. (29). Note that we cannot carry on the Gaussian elimination process any further, because none of the polynomial coefficients appearing in Eq. (30) have an inverse in the ring.

We will assume that the polynomials have rational coefficients, i.e., the coefficients belong to \({\mathscr {Q}}\), the field of the rational numbers. The set of polynomials form a ring—the polynomial ring in three variables, which we denote by \({\mathscr {K}}={\mathscr {Q}}[{\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3]\). The polynomial vector \((q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}})\, \in \, {\mathscr {K}}^4\). The set of solutions to Eq. (30) is just the kernel of the homomorphism \(\varphi \) defined by,

$$\begin{aligned} \varphi : {\mathscr {K}}^4\longrightarrow & {} {\mathscr {K}}\, \\ (q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}})\longrightarrow & {} (1 - {\mathscr {D}}_1 {\mathscr {D}}_2 {\mathscr {D}}_3) q_3 + ({\mathscr {D}}_1 {\mathscr {D}}_3 - {\mathscr {D}}_2) q_{1^{\prime}} + {\mathscr {D}}_1 (1 - {\mathscr {D}}_3^2) q_{2^{\prime}} \, \\&\quad + (1 - {\mathscr {D}}_1^2) q_{3^{\prime}} \equiv \varPsi ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3) . \end{aligned}$$
(31)

The polynomial \(\varPsi ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3) \in {\mathscr {K}}\). If we set \(\varPsi \equiv 0\) we obtain the kernel of the homomorphism \(\varphi \) denoted by \( \ker \varphi \) which is all those 4-tuples \((q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}})\) which map to zero. Thus, the solution space is \( \ker \varphi \) and is also a submodule of \({\mathscr {K}}^4\). It is called the first module of syzygies. The physical significance of the kernel of a homomorphism is that the laser phase noise is mapped to zero. There are also second, third, …etc. modules of syzygies, but we will discuss them later in Sect. 4.6.

The generators of the first module of syzygies can be obtained from standard methods available in the literature. We briefly outline the method given in the books by Becker and Weispfenning (1993) and Kreuzer and Robbiano (2000) below.

4.4 Gröbner basis

The first step is to obtain the Gröbner basis for the ideal \({\mathscr {U}}\) generated by the coefficients in Eq. (30):

$$\begin{aligned} u_1 = 1 - {\mathscr {D}}_1 {\mathscr {D}}_2 {\mathscr {D}}_3 , \qquad u_2 = {\mathscr {D}}_1 {\mathscr {D}}_3 - {\mathscr {D}}_2 , \qquad u_3 = {\mathscr {D}}_1 (1 - {\mathscr {D}}_3^2) , \qquad u_4 = 1 - {\mathscr {D}}_1^2 . \end{aligned}$$
(32)

The ideal \({\mathscr {U}}\) consists of linear combinations of the form \(\sum v_i u_i\) where \(v_i\), \(i = 1, \dots , 4\) are polynomials in the ring \({\mathscr {K}}\). There can be several sets of generators for \({\mathscr {U}}\). A Gröbner basis is a set of generators which is ‘small’ in a specific sense.

There are several ways to look at the theory of Gröbner basis. One way is the following: Suppose we are given polynomials \(g_1, g_2, \dots , g_m\) in one variable over say \({\mathscr {Q}}\) and we would like to know whether another polynomial f belongs to the ideal generated by the g’s. A good way to decide the issue would be to first compute the GCD g of \(g_1\), \(g_2\), …, \(g_m\) and check whether f is a multiple of g. One can achieve this by doing the long division of f by g and checking whether the remainder is zero. All this is possible because \({\mathscr {Q}}[x]\) is a Euclidean domain and also a principle ideal domain (PID) wherein any ideal is generated by a single element. Therefore we have essentially just one polynomial—the GCD—which generates the ideal generated by \(g_1, g_2, \dots , g_m\). The ring of integers or the ring of polynomials in one variable over any field are examples of PIDs whose ideals are generated by single elements. However, when we consider more general rings (not PIDs) like the one we are dealing with here, we do not have a single GCD but a set of several polynomials that generate an ideal in general. A Gröbner basis of an ideal can be thought of as a generalization of the GCD. In the univariate case, the Gröbner basis reduces to the GCD.

Gröbner basis theory generalizes these ideas to multivariate polynomials which are neither Euclidean rings nor PIDs. Since there is in general not a single generator for an ideal, Gröbner basis theory comes up with the idea of dividing a polynomial with a set of polynomials, the set of generators of the ideal, so that by successive divisions by the polynomials in this generating set of the given polynomial, the remainder becomes zero. Clearly, every generating set of polynomials need not possess this property. Those special generating sets that do possess this property (and they exist!) are called Gröbner bases. For a division to be carried out in a sensible manner, an order must be put on the ring of polynomials, so that the final remainder after every division is strictly smaller than each of the divisors in the generating set. A natural order exists on the ring of integers or on the polynomial ring \({\mathscr {Q}}(x)\); the degree of the polynomial decides the order in \({\mathscr {Q}}(x)\). However, even for polynomials in two variables there is no natural order a priori (is \(x^2 + y\) greater or smaller than \(x + y^2\)?). But one can, by hand as it were, put an order on such a ring by saying \(x \gg y\), where \( \gg \) is an order, called the lexicographical order. We follow this type of order, \({\mathscr {D}}_1 \gg {\mathscr {D}}_2 \gg {\mathscr {D}}_3\) and ordering polynomials by considering their highest degree terms. It is possible to put different orderings on a given ring which then produce different Gröbner bases. Clearly, a Gröbner basis must have ‘small’ elements so that division is possible and every element of the ideal when divided by the Gröbner basis elements leaves zero remainder, i.e., every element modulo the Gröbner basis reduces to zero.

In the literature, there exists a well-known algorithm called the Buchberger algorithm, which may be used to obtain the Gröbner basis for a given set of polynomials in the ring. So a Gröbner basis of \({\mathscr {U}}\) can be obtained from the generators \(u_i\) given in Eq. (32) using this algorithm. The algorithm computes S-polynomials of polynomials in \({\mathscr {U}}\) pairwise, by canceling out the head terms and thus obtaining another polynomial which is ‘smaller’—we will explicitly demonstrate this procedure in Sect. 4.5 for a specific case. By repeating this procedure we obtain smaller polynomials until we reach the Gröbner basis. It is essentially again a generalization of the usual long division that we perform on univariate polynomials. More conveniently, we can use the well known application Mathematica (Wolfram 2014). The function GroebnerBasis in Mathematica yields a 3-element Gröbner basis \({\mathscr {G}}\) for \({\mathscr {U}}\):

$$\begin{aligned} {\mathscr {G}} = \{{\mathscr {D}}_3^2 - 1, {\mathscr {D}}_2^2 - 1, {\mathscr {D}}_1 - {\mathscr {D}}_2 {\mathscr {D}}_3 \} . \end{aligned}$$
(33)

One can easily check that all the \(u_i\) of Eq. (32) are linear combinations of the polynomials in \({\mathscr {G}}\) and hence \({\mathscr {G}}\) generates \({\mathscr {U}}\). One also observes that the elements look ‘small’ in the order mentioned above. However, one can satisfy oneself that \({\mathscr {G}}\) is a Gröbner basis by using the standard methods available in the literature. One method consists of computing the S-polynomials for all the pairs of the Gröbner basis elements and checking whether these reduce to zero modulo \({\mathscr {G}}\).

This Gröbner basis of the ideal \({\mathscr {U}}\) is then used to obtain the generators for the first module of syzygies. Note that although the Gröbner basis depends on the order we choose among the \({\mathscr {D}}_k\), the module itself is independent of the order (Becker and Weispfenning 1993).

4.5 Generating set for the first module of syzygies

A generating set for the module is obtained by further following the procedure in the literature (Becker and Weispfenning 1993; Kreuzer and Robbiano 2000). As we will show, we obtain seven generators for the module. These generators do not form a minimal set and there are relations among them; in fact this method does not guarantee a minimum set of generators. These generators can be expressed as linear combinations of \(\alpha \), \(\beta \), \(\gamma \), \(\zeta \) and also in terms of \(X^{(1)}\), \(X^{(2)}\), \(X^{(3)}\), \(X^{(4)}\) given below in Eq. (51). The importance in obtaining the seven generators is that the standard theorems guarantee that these seven generators do in fact generate the required module. Therefore, from this proven set of generators we can check whether a particular set is in fact a generating set. We present two important generating sets below.

We now follow the procedure and notation of Becker and Weispfenning (1993). We require the 4-tuple solutions \((q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}})\) to the equation:

$$\begin{aligned} (1 - xyz) \, q_3 + (xz - y) \, q_{1^{\prime}} + x(1 - z^2) \, q_{2^{\prime}} + (1 - x^2) \, q_{3^{\prime}} = 0 , \end{aligned}$$
(34)

where for convenience we have substituted \(x \equiv {\mathscr {D}}_1\), \(y \equiv {\mathscr {D}}_2\), \(z \equiv {\mathscr {D}}_3\). \(q_3\), \(q_{1^{\prime}}\), \(q_{2^{\prime}}\), \(q_{3^{\prime}}\) are polynomials in x, y, z with integer coefficients, i.e., in Z[xyz]. Consider the ideal in Z[xyz] (or \({\mathscr {Q}}[x,y,z]\) where \({\mathscr {Q}}\) denotes the field of rational numbers), formed by taking linear combinations of the coefficients in Eq. (34),

$$\begin{aligned} {\mathscr {F}}= \{f_1 = 1 - xyz, \,f_2 = xz - y, \,f_3 = x(1 - z^2), \,f_4 = 1 - x^2 \} . \end{aligned}$$
(35)

A Gröbner basis for this ideal is

$$\begin{aligned} {\mathscr {G}} = \{g_1 = z^2 - 1, g_2 = y^2 - 1, g_3 = x - yz\} . \end{aligned}$$
(36)

One can check that both the \(f_i\), \(i = 1,2,3,4\), and \(g_j\), \(j =1,2,3\), generate the same ideal because we can express one generating set in terms of the other and vice-versa:

$$\begin{aligned} f_i = d_{ij} g_j , \qquad g_j = c_{ji} f_i , \end{aligned}$$
(37)

where d and c are \(4 \times 3\) and \(3 \times 4\) polynomial matrices, respectively, and are given by

$$\begin{aligned} d = \left( \begin{array}{ccc} -1 &{} -z^2 &{} -yz \\ y &{} 0 &{} z \\ -x &{} 0 &{} 0 \\ -1 &{} - z^2 &{} -(x + yz) \end{array} \right) , \qquad c = \left( \begin{array}{cccc} 0 &{} 0 &{} -x &{} z^2 - 1 \\ -1 &{} -y &{} 0 &{} 0 \\ 0 &{} z &{} 1 &{} 0 \end{array} \right) . \end{aligned}$$
(38)

We now forge ahead to obtain the generators for the module. This is accomplished in two steps. First we obtain one set of generators which we denote by A. Then we obtain another set which we denote by \(B^*\). The full set of generators of the module is given by the set \(A \, \bigcup B^*\).

The matrices appearing in Eq. (38), or the relations in Eq. (37), can be used to obtain one set of generators, namely, A. We write:

$$\begin{aligned} f_i = d_{ij} g_j = d_{ij} c_{jk} f_k \equiv \delta _{ik} f_k , \end{aligned}$$
(39)

where, \(\delta _{ik}\) is the Kronecker delta or the unit matrix I. The last equality in the above Eq. (39) gives:

$$\begin{aligned} (\delta _{ik} - d_{ij} c_{jk}) f_k = 0 . \end{aligned}$$
(40)

Thus we have found solutions to Eq. (34), namely, \(a_{ik} = \delta _{ik} - d_{ij} c_{jk}\) for each i. In matrix form, A is the set of row vectors of the matrix \(I - d \cdot c \) where the dot denotes the matrix product and I is the identity matrix which is \(4 \times 4\) in our case. We have:

$$\begin{aligned} I - d \cdot c = \left( \begin{array}{cccc} 1 - z^2 &{} 0 &{} -x + yz &{} z^2 - 1 \\ 0 &{} 1 - z^2 &{} xy - z &{} y (1 - z^2) \\ 0 &{} 0 &{} 1 - x^2 &{} x (z^2 - 1) \\ - z^2 &{} xz &{} yz &{} z^2 \end{array} \right) . \end{aligned}$$
(41)

We list the 4 generators in the set A as the row vectors of the above matrix below:

$$\begin{aligned} a_1& = \left( 1 - z^2, 0, -x + yz, z^2 - 1 \right) , \\ a_2& = \left( 0, 1 - z^2, xy - z, y\left( 1 - z^2\right) \right) , \\ a_3& = \left( 0, 0, 1 - x^2, x\left( z^2 - 1\right) \right) , \\ a_4& = \left( -z^2, xz, yz, z^2\right) . \end{aligned}$$
(42)

Now we go on to compute \(B^*\).

The additional generators \((\in B^*)\) are obtained by first computing the S-polynomials of the Gröbner basis \({\mathscr {G}}\). The S-polynomial of two polynomials \(g_1, g_2\) is obtained by multiplying \(g_1\) and \(g_2\) by suitable terms and then adding, so that the highest terms cancel. For example in our case \(g_1 = z^2 - 1\) and \(g_2 = y^2 - 1\), and the highest terms are \(z^2\) for \(g_1\) and \(y^2\) for \(g_2\). Multiply \(g_1\) by \(y^2\) and \(g_2\) by \(z^2\) and subtract. Thus, the S-polynomial \(p_{12}\) of \(g_1\) and \(g_2\) is

$$\begin{aligned} p_{12} = y^2 g_1 - z^2 g_2 = z^2 - y^2 \equiv g_1 - g_2 . \end{aligned}$$
(43)

For the Gröbner basis of 3 elements we get 3 S-polynomials \(p_{12}\), \(p_{13}\), \(p_{23}\). In general, the \(p_{ij}\) must now be re-expressed in terms of the Gröbner basis \({\mathscr {G}}\). We thus have the equations:

$$\begin{aligned} p_{ij}& = s_{ij} g_i - s_{ji} g_j \, \\ p_{ij}& = q_{ijk} g_k . \end{aligned}$$
(44)

Subtracting the second equation from the first we get:

$$\begin{aligned} (s_{ij} - q_{iji}) g_i - (s_{ji} + q_{ijj}) g_j - \sum _{k \ne i,j} q_{ijk} g_k = 0 . \end{aligned}$$
(45)

For the S-polynomial \(p_{12}\), we have \(i = 1, j = 2\) and thus \(s_{12} = y^2, s_{21} = z^2\) from the first equation and \(q_{121} = - q_{122} = 1\) and \(q_{123} = 0\) from the second equation.

We now write Eq. (45) more compactly. We define:

$$\begin{aligned} r_{ijk} = {\left\{ \begin{array}{ll} s_{ij} - q_{iji} &{} k = i \, \\ - (s_{ji} + q_{ijj}) &{} k = j \, \\ - q_{ijk} &{} \mathrm {otherwise} , \end{array}\right. } \end{aligned}$$
(46)

and write Eq. (45) as,

$$\begin{aligned} r_{ijk} g_k = 0 . \end{aligned}$$
(47)

For each pair \(\{ij\}\) we have a row vector—a generator for the module of syzygies, say \(S_{\mathscr {G}}\) for the Gröbner basis elements. Taking the order \(p_{12}, p_{13}, p_{23}\) for the rows we get a matrix b which is given by:

$$\begin{aligned} b = \left( \begin{array}{ccc} y^2 - 1 &{} 1 - z^2 &{} 0 \\ x - yz &{} 0 &{} 1 - z^2 \\ 0 &{} x - yz &{} 1 - y^2 \end{array} \right) . \end{aligned}$$
(48)

Writing \(g_k = c_{km} f_m\) gives us the required generators for the set \({\mathscr {F}}\) because,

$$\begin{aligned} r_{ijk} c_{km} f_m = 0 . \end{aligned}$$
(49)

We have thus found more solutions to Eq. (34). This is conveniently achieved by multiplying b by the matrix c to obtain \(b^* = b \cdot c\). The row vectors \(b^*_i\), \(i = 1, 2, 3\), of \(b^*\) form the set \(B^*\):

$$\begin{aligned} b^*_1& = \left( z^2 - 1, y\left( z^2 - 1\right) , x\left( 1 - y^2\right) , \left( y^2 -1\right) \left( z^2 - 1\right) \right) , \\ b^*_2& = \left( 0, z\left( 1 - z^2\right) , 1 - z^2 - x\left( x - yz\right) , \left( x - yz\right) \left( z^2 - 1\right) \right) , \\ b^*_3& = \left( -x + yz, z - xy, 1 - y^2, 0\right) . \end{aligned}$$
(50)

This is the set \(B^*\) and the full set of generators is seven in number and is given by \(A \, \bigcup B^*\).

Alternatively, we may use a software package called Macaulay 2 (Grayson and Stillman 2019) that directly calculates the generators given the Eq. (28). Using Macaulay 2, we obtain six generators. Again, Macaulay’s algorithm does not yield a minimal set; we can express the last two generators in terms of the first four. Below we list this smaller set of four generators in the order \(X = (q_1, q_2, q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}})\):

$$\begin{aligned} X^{(1)}& = \left( y - xz, 0, 1 - z^2, 0, yz - x, z^2 - 1\right) , \\ X^{(2)}& = \left( - x, - y, - z, x, y, z \right) , \\ X^{(3)}& = \left( -1, - z, - xz, 1, xy, y \right) , \\ X^{(4)}& = \left( - xy, -1, - x, z, 1, yz \right) . \end{aligned}$$
(51)

Another set of generators are just \(\alpha \), \(\beta \), \(\gamma \), and \(\zeta \). This can be checked using Macaulay 2 (Grayson and Stillman 2019), or one can relate \(\alpha \), \(\beta \), \(\gamma \), and \(\zeta \) to the generators \(X^{(A)}\), \(A = 1, 2, 3, 4\), by polynomial matrices. We list them below for ready reference:

$$\begin{aligned} \alpha& = (-1, -z, -xz, 1, xy, y) , \\ \beta& = (-xy, -1, -x, z, 1, yz) , \\ \gamma& = (-y, -yz, -1, xz, x, 1) , \\ \zeta& = (-x, -y, -z, x, y, z) . \end{aligned}$$
(52)

In the generating set given in Eq. (51), the last three generators are just \(X^{(2)} = \zeta \), \(X^{(3)} =\alpha \), \(X^{(4)} = \beta \). An extra generator \(X^{(1)}\) is needed to generate all the solutions. In Appendix A, we express the seven generators we obtained following the literature, in terms of \(\alpha \), \(\beta \), \(\gamma \), and \(\zeta \). Also we express \(\alpha \), \(\beta \), \(\gamma \), and \(\zeta \) in terms of \(X^{(A)}\). This proves that all these sets generate the required module of syzygies.

The question now arises as to which set of generators we should choose that facilitates further analysis. The analysis is simplified if we choose a smaller number of generators. Also we would prefer low degree polynomials to appear in the generators so as to avoid cancellation of leading terms in the polynomials. By these two criteria we may choose \(X^{(A)}\) or \(\alpha \), \(\beta \), \(\gamma \), \(\zeta \). However, \(\alpha \), \(\beta \), \(\gamma \), \(\zeta \) possess the additional property that this set is left invariant under a cyclic permutation of indices 1, 2, 3. It is found that this set is more convenient to use because of this symmetry (Armstrong et al. 1999).

We would like to emphasise that the theorems we have used guarantee that all the TDIs are generated by taking the linear combinations of the generators—in particular the generating set \(\alpha , \beta , \gamma \) and \(\zeta \) generates all the TDI.

We remark that, the method described above, of Gröbner basis and the module of syzygies has been applied to a possible future space mission concept which involves six spacecraft arranged in the configuration of an octohedron (Wang et al. 2013). This constellation has 24 links and it is shown that these provide enough redundancy in the data to cancel not only the laser frequency noise but also the acceleration noise. The case of time-independent and equal arm-lengths has been treated and for this 7 generators have been found for the relevant module of syzygies. The seven generators generate the full set of noise cancelling data combinations, which is again guaranteed by the theorems, and this is achieved by taking their linear combinations.

An alternative and a totally different approach has been adopted in (Romano and Woan 2006) from the point of view of Bayesian statistical inference. The covariance matrix of the six elementary data streams sampled at integer multiples of the time-delay is first formed. Then a principal component analysis (PCA) of the covariance matrix is done and the space of TDIs is identified with the subspace formed by those eigenvectors belonging to small eigenvalues only; the large eigenvalues and the corresponding eigenvectors correspond to the laser frequency noise. As an example, they also obtain the TDI observable \(\alpha \) with their analysis. This analysis has been performed for the case of the LISA model with constant (time-independent) arm-lengths which fits into the discussion in this subsection as it deals with constant arm-lengths.

4.6 Relation to Hilbert’s syzygy theorem

In mathematics, the Hilbert’s syzygy theorem is one of the three fundamental theorems about polynomial rings over fields which lies at the roots of modern algebraic geometry. (The other two theorems are Hilbert’s basis theorem and Hilbert’s nullstellensatz—but they do not concern us here). Hilbert’s syzygy theorem is about the relations or syzygies (according to Hilbert) between generators of an ideal or more generally a module (an ideal is trivially a module). These syzygies form a module as we have seen before—the first module of syzygies. Now this module has generators, which we have explicitly calculated in the previous subsection. But there could be non-trivial relations between these generators, unless the generators are linearly independent. If it so happens that the generators are linearly independent, then the module has a basis and the module is said to be free. However, if there are non-trivial relations between the generators—the generators are not linearly independent, these relations again form a non-zero module—it is called the second module of syzygies. We may continue this process and obtain the \(k{\mathrm{th}}\) module of syzygies. Hilbert’s syzygy theorem states that if \({\mathscr {M}}\) is a finitely generated module over a polynomial ring \({\mathscr {F}}[x_1, x_2, \dots , x_n]\) in n variables over a field \({\mathscr {F}}\) then the \(k{\mathrm{th}}\) module of syzygies is free where \(k \le n\). The original ideal (or module) \({\mathscr {U}}\) is then said to be resolved after k steps. Hilbert’s theorem asserts that the module can be resolved in at most n steps, that is, \(k \le n\). Further, the number k does not depend on the choice of generating sets.

How does this concern us here? Let us go ahead and obtain the second module of syzygies—call it \({\mathscr {S}}\). We choose the generating set \(\alpha , \beta , \gamma , \zeta \) for the first module of syzygies. We then write the equation:

$$\begin{aligned} p_1 \alpha + p_2 \beta + p_3 \gamma + p_4 \zeta = 0 , \end{aligned}$$
(53)

where \(p_j \in {\mathscr {K}},\,\,j = 1, 2, 3, 4\), that is, each \(p_j\) is a polynomial in the variables xyz. We now look for the 4-tuple solutions \((p_1, p_2, p_3, p_4) \in {\mathscr {K}}^4\) satisfying Eq. (53). We find that the solutions are just multiples of the polynomial vector:

$$\begin{aligned} {\mathbf {e}}= (x - yz,\, y - xz,\, z - xy,\, xyz - 1) \in {\mathscr {K}}^4 , \end{aligned}$$
(54)

and the module \({\mathscr {S}}\) is given by

$$\begin{aligned} {\mathscr {S}}= \{ \,p(x, y, z)\,{\mathbf {e}}\,|\,\, p(x, y, z) \in {\mathscr {K}}\} , \end{aligned}$$
(55)

where p(xyz) is a polynomial in \({\mathscr {K}}\). \({\mathscr {S}}\) is in fact a free module of rank 1 because \({\mathbf {e}}\) generates \({\mathscr {S}}\) and is linearly independent. We find that the original module has resolved in only two steps which is less than three, thus respecting Hilbert’s theorem. So roughly speaking the original module is less ‘entangled’ because it did not require the maximum number of steps, namely, three in this case, to resolve. Also since the second module of syzygies essentially consists of just one non-trivial relation—it has only one generator—\(\alpha , \beta , \gamma \) and \(\zeta \) have essentially one non-trivial relation between them, they are close to being linearly independent.

We can also write the above modules and their homomorphisms as an exact sequence:

$$\begin{aligned} 0 \longrightarrow {\mathscr {K}}{\mathop {\longrightarrow }\limits ^{\chi }} {\mathscr {K}}^4 {\mathop {\longrightarrow }\limits ^{\psi }} {\mathscr {K}}^4 {\mathop {\longrightarrow }\limits ^{\varphi }} {\mathscr {U}}\longrightarrow 0 . \end{aligned}$$

The exact sequence implies that the kernel of a given homomorphism in the sequence is the image of the homomorphism on its left. This means the following in our case: \(\mathrm {Im}\ (\psi ) = \mathrm {ker}\ (\varphi )\) and \(\mathrm {Im}\ (\chi ) = \mathrm {ker}\ (\psi )\). Further, the left most module is the zero module since ker \((\chi ) = {0}\) or \(\chi \) is an isomorphism which maps the polynomial ring \({\mathscr {K}}\) to \({\mathscr {S}}\subset {\mathscr {K}}^4\). Similarly, the right most zero in the sequence implies that \(\varphi \) is onto.

We now define the homomorphism. There are three homomorphisms \(\chi , \psi , \varphi \) here. Starting from the right, we have the homomorphism \(\varphi \) which has already been defined earlier in Eq. (31). We now define the homomorphism \(\psi \). Let \(\{\varepsilon _1, \varepsilon _2, \varepsilon _3, \varepsilon _4 \}\) be the canonical basis of \({\mathscr {K}}^4\), that is, \(\varepsilon _1 = (1, 0, 0, 0)\) etc. \(\varepsilon _j\) has 1 in the \(j{\mathrm{th}}\) place and zeros everywhere else. Then \(\psi \) maps the \(\varepsilon _1, \varepsilon _2, \varepsilon _3, \varepsilon _4\) to the generators \(\alpha , \beta , \gamma , \zeta \) of the first module of syzygies \({\mathscr {M}}\) respectively. The homomorphism is then extended to whole of \({\mathscr {K}}^4\) by linearity and its image is \({\mathscr {M}}\subset {\mathscr {K}}^4\). The homomorphism \(\chi \) just takes a polynomial \(p (x, y, z) \in {\mathscr {K}}\) and maps it to \(p (x, y, z) {\mathbf {e}}\in {\mathscr {K}}^4\). As p(xyz) ranges over \({\mathscr {K}}\), the free module \({\mathscr {S}}\) is generated.

Writing the homomorphisms as an exact sequence expresses the concepts graphically and in a unified manner. Also that the original module is resolved in only two steps instead of three, and that the second module of syzygies is of rank one, may have important physical implications.

4.7 Canceling optical bench motion noise

There are eighteen Doppler data streams that have to be combined in an appropriate manner to cancel the noise from the laser as well as from the motion of the optical benches. From Eqs. (1217) and the other twelve data streams on spacecraft 2 and 3 obtained by cyclic permutations of them, we can now derive the expressions canceling the lasers and optical benches noises. We first note that the differences \(\varepsilon - \tau \) from each optical bench provide twice the displacement of the optical bench relative to the proof mass. Second, the laser phase-fluctuations with primed indices, \(p_{i'}\), can be expressed in terms of those with unprimed indices, \(p_i\), by taking suitable linear combinations of the \(s_i, s_{i'}, \varepsilon _i - \tau _i, \tau _i, \tau _{i'}\). So, let us first define the following linear combination

$$\begin{aligned} z_1 \equiv \frac{\tau _{1} - \tau _{1^{\prime}}}{2} = p_{1^{\prime}} - p_1, \end{aligned}$$
(56)

which only contains the two laser noises \(p_1\), \(p_{1^{\prime}}\).Footnote 3 By then defining the following variables \(\xi _1\), \(\xi _{1^{\prime}}\)

$$\begin{aligned} \xi _{1^{\prime}}\equiv & {} s_{1^{\prime}} + \frac{(\varepsilon _{1^{\prime}} - \tau _{1^{\prime}})}{2} + {\mathscr {D}}_{2^{\prime}} \frac{(\varepsilon _{3} - \tau _{3})}{2} , \end{aligned}$$
(57)
$$\begin{aligned} \xi _{1}\equiv & {} s_{1} + \frac{(\varepsilon _{1} - \tau _{1})}{2} + {\mathscr {D}}_{3}\frac{(\varepsilon _{2^{\prime}} - \tau _{2^{\prime}})}{2} , \end{aligned}$$
(58)

it is easy to see that they depend only on the laser noises \(p_i\) and \(p_{i'}\), and the proof-mass noises \({\varvec{\delta }}_{i} , {\varvec{\delta }}_{i'}\), as the optical bench noises have been removed. By then introducing the two observables, \(\eta _{1}\) and \(\eta _{1^{\prime}}\), defined as follows

$$\begin{aligned} \eta _{1^{\prime}}\equiv & {} \xi _{1^{\prime}} + z_1 \ = {\mathscr {D}}_{2^{\prime}} p_3 - p_1 , \end{aligned}$$
(59)
$$\begin{aligned} \eta _{1}\equiv & {} \xi _{1} - {\mathscr {D}}_{3} z_{2} \ = {\mathscr {D}}_{3} p_2 - p_1 , \end{aligned}$$
(60)

we have just reduced the problem of canceling the six laser noises and six optical bench noises to that of canceling the three laser noises \(p_i\). In other words, by synthesizing the six linear combinations \(\eta _i\), \(\eta _{i'}\) in terms of the eighteen measurements made by LISA we have reduce the problem to that for the simpler configuration with only three lasers, analyzed in the previous Sects. 4.14.4.

4.8 Physical interpretation of the TDI combinations

It is important to notice that the four interferometric combinations \((\alpha , \beta , \gamma , \zeta )\), which can be used as a basis for generating the entire TDI space, are actually synthesized Sagnac interferometers. This can be seen by rewriting the expression for \(\alpha \), for instance, in the following form,

$$\begin{aligned} \alpha = [\eta _{1^{\prime}} + {\mathscr {D}}_{2^{\prime}} \eta _{3^{\prime}} + {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_{1^{\prime}} \eta _{2^{\prime}}] - [\eta _{1} + {\mathscr {D}}_{3} \eta _{2} + {\mathscr {D}}_{3} {\mathscr {D}}_{1} \eta _{3}] , \end{aligned}$$
(61)

and noticing that the first square bracket on the right-hand side of Eq. (61) contains a combination of one-way measurements describing a light beam propagating clockwise around the array, while the other terms in the second square-bracket give the equivalent of another beam propagating counter-clockwise around the constellation.

Contrary to (\(\alpha \), \(\beta \), \(\gamma \)), \(\zeta \) can not be visualized as the difference (or interference) of two synthesized beams. However, it should still be regarded as a Sagnac combination since there exists a time-delay relationship between it and \(\alpha \), \(\beta \), and \(\gamma \) (Armstrong et al. 1999):

$$\begin{aligned} (I - {\mathscr {D}}_1 {\mathscr {D}}_2 {\mathscr {D}}_3) \zeta = ({\mathscr {D}}_1 - {\mathscr {D}}_2 {\mathscr {D}}_3) \alpha + ({\mathscr {D}}_2 - {\mathscr {D}}_3 {\mathscr {D}}_1) \beta + ({\mathscr {D}}_3 - {\mathscr {D}}_1 {\mathscr {D}}_2) \gamma . \end{aligned}$$
(62)

As a consequence of the time-structure of this relationship, \(\zeta \) has been called the symmetrized Sagnac combination.

By using the four generators, it is possible to construct several other interferometric combinations, such as the unequal-arm Michelson (XYZ), the Beacons (PQR), the Monitors (EFG), and the Relays (UVW). Contrary to the Sagnac combinations, these only use four of the six data combinations \(\eta _i\), \(\eta _{i'}\). For this reason they have obvious utility in the event of selected subsystem failures (Estabrook et al. 2000).

These observables can be written in terms of the Sagnac observables \((\alpha , \beta , \gamma , \zeta )\) in the following way,

$$\begin{aligned} {\mathscr {D}}_1 X& = {\mathscr {D}}_2 {\mathscr {D}}_3 \alpha - {\mathscr {D}}_2 \beta - D_3 \gamma + \zeta , \\ P& = \zeta - {\mathscr {D}}_1 \alpha , \\ E& = \alpha - {\mathscr {D}}_1 \zeta , \\ U& = {\mathscr {D}}_1 \gamma - \beta , \end{aligned}$$
(63)

as it is easy to verify by substituting the expressions for the Sagnac combinations into the above equations. Their physical interpretations are schematically shown in Fig. 6 by representing with arrows the one-way links they rely on.

Fig. 6
figure 6

Schematic diagrams of the unequal-arm Michelson, Monitor, Beacon, and Relay combinations. These TDI combinations rely only on four of the six one-way Doppler measurements, as illustrated here

In the case of the combination X, in particular, by writing it in the following form (Armstrong et al. 1999),

$$\begin{aligned} X = \left[ (\eta _{1^{\prime}} + {\mathscr {D}}_{2}\eta _3) + {\mathscr {D}}_{2} {\mathscr {D}}_2 (\eta _{1} + {\mathscr {D}}_3 \eta _{2})\right] - \left[ (\eta _{1} + {\mathscr {D}}_3 \eta _{2^{\prime}}) + {\mathscr {D}}_3 {\mathscr {D}}_{3} (\eta _{1^{\prime}} + {\mathscr {D}}_{2} \eta _{3})\right] , \end{aligned}$$
(64)

one can notice (as pointed out in Summers 2003; Shaddock et al. 2003) that this combination can be visualized as the difference of two sums of phase measurements, each corresponding to a specific light path from a laser on board spacecraft 1 having phase noise \(p_1\). The first square-bracket term in Eq. (64) represents a synthesized light-beam transmitted from spacecraft 1 and made to bounce once at spacecraft 2 and 3, respectively. The second square-bracket term instead corresponds to another beam also originating from the same laser, experiencing the same overall delay as the first beam, but bouncing off spacecraft 3 first and then spacecraft 2. When they are recombined they will cancel the laser phase fluctuations exactly, having both experienced the same total delay (assuming stationary spacecraft). The X combinations should therefore be regarded as the response of a zero-area Sagnac interferometer.

5 Time-delay interferometry with moving spacecraft

The rotational motion of the LISA array results in a difference of the light travel times in the two directions around a Sagnac circuit (Shaddock 2004; Cornish and Hellings 2003). Two time delays along each arm must be used, say \(L_{i}'\) and \(L_i\), for clockwise and counter-clockwise propagations respectively as they enter in any of the TDI combinations. Furthermore, since \(L_i\) and \(L_{i}'\) not only differ from one another but can be time dependent (they “flex”), it was shown that the “first generation” TDI combinations as derived so far do not suppress the laser phase noise (at least with present laser stability requirements) below a level determined by the secondary noises. For LISA, and assuming \(\dot{L}_i \simeq 10 \mathrm {\ m/s}\) (Amaro-Seoane et al. 2017), the estimated magnitude of the remaining frequency fluctuations from the laser can be about 30 times larger than the level set by the secondary noise sources in the center of the frequency band. To solve this potential problem, it has been shown that there exist new TDI combinations that are immune to first order shearing (flexing, or constant rate of change of delay times). These combinations can be derived by using the time-delay operators formalism introduced in the previous Sect. 4, although one has to keep in mind that now these operators no longer commute (Tinto et al. 2004).

To derive the new, “flex-free” TDI combinations we will start by taking specific combinations of the one-way data entering in each of the expressions derived in the previous Sect. 4. In what follows we focus only on the laser noises, and note that the expressions for the \(\eta \)-measurements now assume the following form as a result of the motion of the spacecraft

$$\begin{aligned} \eta _{1}& = {\mathscr {D}}_3 p_2 - p_1 , \qquad \eta _{1^{\prime}} = {\mathscr {D}}_{2^{\prime}} p_3 - p_1 , \end{aligned}$$
(65)
$$\begin{aligned} \eta _{2}& = {\mathscr {D}}_1 p_3 - p_2 , \qquad \eta _{2^{\prime}} = {\mathscr {D}}_{3^{\prime}} p_1 - p_2 , \end{aligned}$$
(66)
$$\begin{aligned} \eta _{3}& = {\mathscr {D}}_2 p_1 - p_3 , \qquad \eta _{3^{\prime}} = {\mathscr {D}}_{1^{\prime}} p_2 - p_3 . \end{aligned}$$
(67)

The combinations of the one-way measurements are chosen in such a way to retain only one of the three laser noises, if possible. In this way we can then implement an iterative procedure based on the use of these basic combinations and of time-delay operators, to cancel the laser noises after dropping terms that are quadratic in \({\dot{L}}/c\) or linear in the accelerations. This iterative time-delay method, to first order in the velocity, is illustrated abstractly as follows. Given a function of time \(\varPsi (t)\), time delay by \(L_i\) is now denoted either with the standard comma notation (Armstrong et al. 1999) or by applying the delay operator \({\mathscr {D}}_{i}\) introduced in the previous Sect. 4,

$$\begin{aligned} {\mathscr {D}}_{i} \varPsi \equiv \varPsi _{,i} \equiv \varPsi (t - L_i(t)) . \end{aligned}$$
(68)

We then impose a second time delay \(L_j(t)\):

$$\begin{aligned} {\mathscr {D}}_{j} {\mathscr {D}}_{i} \varPsi \equiv \varPsi _{;ij}& = \varPsi (t - L_j(t) - L_i(t - L_j(t))) \\\simeq & {} \varPsi (t - L_j(t) - L_i(t) + \dot{L}_i(t) L_j) \\\simeq & {} \varPsi _{,ij} + {\dot{\varPsi }}_{,ij} \dot{L}_i L_j . \end{aligned}$$
(69)

A third time delay \(L_k(t)\) gives

$$\begin{aligned} {\mathscr {D}}_{k} {\mathscr {D}}_{j} {\mathscr {D}}_{i} \varPsi = \varPsi _{;ijk}& = \varPsi (t - L_k(t) - L_j(t - L_k(t)) - L_i(t - L_k(t) - L_j(t - L_k(t)))) \\\simeq & {} \varPsi _{,ijk} + {\dot{\varPsi }}_{,ijk} \left[ \dot{L}_i (L_j + L_k) + \dot{L}_j L_k\right] , \end{aligned}$$
(70)

and so on, recursively; each delay generates a first-order correction proportional to its rate of change times the sum of all delays coming after it in the subscripts. Commas have now been replaced with semicolons (Shaddock et al. 2003), to remind us that we consider moving arrays. When the sum of these corrections to the terms of a data combination vanishes, the combination is called flex-free.

Also, note that each delay operator \({\mathscr {D}}_{i}\) has a unique inverse \({\mathscr {D}}^{-1}_{i}\), whose expression can be derived by requiring that \({\mathscr {D}}^{-1}_{i} {\mathscr {D}}_{i} = I\), and neglecting quadratic and higher order velocity terms. Its action on a time series \(\varPsi (t)\) is

$$\begin{aligned} {\mathscr {D}}^{-1}_{i} \varPsi (t) \equiv \varPsi (t + L_i (t + L_i)) . \end{aligned}$$
(71)

Note that this is not like an advance operator one might expect, since it advances not by \(L_i (t)\) but rather \(L_i (t + L_i)\).

5.1 The unequal-arm Michelson

The unequal-arm Michelson combination relies on the four measurements \(\eta _{1}\), \(\eta _{1^{\prime}}\), \(\eta _{2^{\prime}}\), and \(\eta _{3}\). Note that the two combinations \(\eta _{1} + \eta _{2',3}\), \(\eta _{1^{\prime}} + \eta _{3,2'}\) represent the two synthesized two-way data measured on board spacecraft 1, and can be written in the following form [see Eqs. (656667) for deriving the following synthesized two-way measurements]

$$\begin{aligned} \eta _{1} + \eta _{2',3}& = \left( {\mathscr {D}}_{3}{\mathscr {D}}_{3^{\prime}} - I \right) p_1 , \end{aligned}$$
(72)
$$\begin{aligned} \eta _{1^{\prime}} + \eta _{3,2'}& = \left( {\mathscr {D}}_{2^{\prime}}{\mathscr {D}}_{2} - I \right) p_1 , \end{aligned}$$
(73)

where I is the identity operator. Since in the stationary case any pairs of these operators commute, i.e., \({\mathscr {D}}_i {\mathscr {D}}_{j'} - {\mathscr {D}}_{j'} {\mathscr {D}}_i = 0\), from Eqs. (72, 73) it is easy to derive the following expression for the unequal-arm interferometric combination X which eliminates \(p_1\):

$$\begin{aligned} X = \left[ {\mathscr {D}}_{2^{\prime}}{\mathscr {D}}_{2} - I\right] (\eta _{1} + \eta _{2',3}) - \left[ {\mathscr {D}}_{3}{\mathscr {D}}_{3^{\prime}} - I \right] (\eta _{1^{\prime}} + \eta _{3,2'}) . \end{aligned}$$
(74)

If, on the other hand, the time-delays depend on time, the expression of the unequal-arm Michelson combination above no longer cancels \(p_1\). To derive the new expression for the unequal-arm interferometer that accounts for “flexing”, let us first consider the following two combinations of the one-way measurements entering into the X observable given in Eq. (74):

$$\begin{aligned} \left[ (\eta _{1^{\prime}} + \eta _{3;2'}) + (\eta _{1} + \eta _{2';3})_{;22'}\right]& = \left[ D_{2^{\prime}}D_{2}D_{3}D_{3^{\prime}} - I \right] p_1 , \end{aligned}$$
(75)
$$\begin{aligned} \left[ (\eta _{1} + \eta _{2';3}) + (\eta _{1^{\prime}} + \eta _{3;2'})_{;3'3}\right]& = \left[ D_{3}D_{3^{\prime}}D_{2^{\prime}}D_{2} - I \right] p_1 . \end{aligned}$$
(76)

Using Eqs. (7576), we can use the “delay technique” again to finally derive the following expression for the new unequal-arm Michelson combination \(X_1\) that accounts for the flexing effect:

$$\begin{aligned} X_1& = \left[ D_{2^{\prime}}D_{2}D_{3}D_{3^{\prime}} - I \right] \left[ (\eta _{1} + \eta _{2';3}) + (\eta _{1^{\prime}} + \eta _{3;2'})_{;3'3}\right] \\&\quad - \left[ D_{3}D_{3^{\prime}}D_{2^{\prime}}D_{2} - I \right] \left[ (\eta _{1^{\prime}} + \eta _{3;2'}) + (\eta _{1} + \eta _{2';3})_{;22'}\right] . \end{aligned}$$
(77)

As usual, \(X_2\) and \(X_3\) are obtained by cyclic permutation of the spacecraft indices. This expression is readily shown to be laser-noise-free to first order of spacecraft separation velocities \(\dot{L}_i\): it is “flex-free”.

5.2 The Sagnac combinations

In the above Sect. 5.1, we have used the same symbol X for the unequal-arm Michelson combination for both the rotating (i.e., constant delay times) and stationary cases. This emphasizes that, for this TDI combination (and, as we will see below, also for all the combinations including only four links) the forms of the equations do not change going from systems at rest to the rotating case. One needs only distinguish between the time-of-flight variations in the clockwise and counter-clockwise senses (primed and unprimed delays).

In the case of the Sagnac variables \((\alpha , \beta , \gamma , \zeta )\), however, this is not the case as it is easy to understand on simple physical grounds. In the case of \(\alpha \) for instance, light originating from spacecraft 1 is simultaneously sent around the array on clockwise and counter-clockwise loops, and the two returning beams are then recombined. Even if the array is only rigidly rotating, the two beams experience a different delay (the Sagnac effect), preventing the noise \(p_1\) from canceling in the \(\alpha \) combination.

To find the solution to this problem let us first rewrite \(\alpha \) in such a way to explicitly emphasize what it does: attempts to remove the same fluctuations affecting two beams that have been made to propagated clockwise and counter-clockwise around the array:

$$\begin{aligned} \alpha = [\eta _{1^{\prime}} + {\mathscr {D}}_{2^{\prime}} \eta _{3^{\prime}} + {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_{1^{\prime}} \eta _{2^{\prime}}] - [\eta _{1} + {\mathscr {D}}_{3} \eta _{2} + {\mathscr {D}}_{3} {\mathscr {D}}_{1} \eta _{3}] , \end{aligned}$$
(78)

where we have accounted for clockwise and counter-clockwise light delays. It is straight-forward to verify that this combination no longer cancels the laser noise. If, however, we expand the two terms inside the square-brackets on the right-hand side of Eq. (78) we find that they are equal to

$$\begin{aligned}&[ \eta _{1^{\prime}} + {\mathscr {D}}_{2^{\prime}} \eta _{3^{\prime}} + {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_{1^{\prime}} \eta _{2^{\prime}}] = [{\mathscr {D}}_{2^{\prime}}{\mathscr {D}}_{1^{\prime}}{\mathscr {D}}_{3^{\prime}} - I] \, p_1 , \end{aligned}$$
(79)
$$\begin{aligned}&[ \eta _{1} + {\mathscr {D}}_{3} \eta _{2} + {\mathscr {D}}_{3} {\mathscr {D}}_{1} \eta _{3}] = [{\mathscr {D}}_{3}{\mathscr {D}}_{1}{\mathscr {D}}_{2} - I] \, p_1 . \end{aligned}$$
(80)

If we now apply our iterative scheme to the combinations given in Eq. (80) we finally get the expression for the Sagnac combination \(\alpha _1\) that is unaffected by laser noise in presence of rotation,

$$\begin{aligned} \alpha _1& = [ {\mathscr {D}}_{3}{\mathscr {D}}_{1}{\mathscr {D}}_{2} - I] \, [\eta _{1^{\prime}} + {\mathscr {D}}_{2^{\prime}} \eta _{3^{\prime}} + {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_{1^{\prime}} \eta _{2^{\prime}}] \\&\quad -[ {\mathscr {D}}_{2^{\prime}}{\mathscr {D}}_{1^{\prime}}{\mathscr {D}}_{3^{\prime}} - I] \, [\eta _{1} + {\mathscr {D}}_{3} \eta _{2} + {\mathscr {D}}_{3} {\mathscr {D}}_{1} \eta _{3}] . \end{aligned}$$
(81)

If the delay-times are also time-dependent, we find that the residual laser noise remaining into the combination \(\alpha _1\) is actually equal to

$$\begin{aligned} \dot{p}_{1,1231'2'3'} \left[ \left( \dot{L}_1 + \dot{L}_2 + \dot{L}_3\right) \left( L_1^{\prime} + L_2^{\prime} + L_3^{\prime}\right) - \left( \dot{L}_1^{\prime} + \dot{L}_2^{\prime} + \dot{L}_3^{\prime}\right) \left( L_1 + L_2 + L_3\right) \right] . \end{aligned}$$
(82)

Fortunately, although first order in the relative velocities, the residual is small, as it involves the difference of the clockwise and counter-clockwise rates of change of the propagation delays on the same circuit. For LISA, the remaining laser phase noises in \(\alpha _i\), \(i = 1, 2, 3\), are several orders of magnitude below the secondary noises.

In the case of \(\zeta \), however, the rotation of the array breaks the symmetry and therefore its uniqueness. However, for a rigidly rotating array, there still exist three generalized TDI laser-noise-free data combinations that have properties very similar to \(\zeta \), and which can be used for the same scientific purposes (Tinto et al. 2001). These combinations, which we call \((\zeta _1, \zeta _2, \zeta _3)\), can be derived by applying again our time-delay operator approach.

Let us consider the following combination of the \(\eta _{i}\), \(\eta _{i'}\) measurements, each being delayed only once (Armstrong et al. 1999):

$$\begin{aligned} \eta _{3,3} - \eta _{3',3} + \eta _{1,1'}& = \left[ D_{3} D_{2} - D_{1^{\prime}}\right] p_1 , \end{aligned}$$
(83)
$$\begin{aligned} \eta _{1',1} - \eta _{2,2'} + \eta _{2',2'}& = \left[ D_{3^{\prime}} D_{2^{\prime}} - D_{1}\right] p_1 , \end{aligned}$$
(84)

where we have used the commutativity property of the delay operators in order to cancel the \(p_2\) and \(p_3\) terms as the array is rigidly rotating. Since both sides of the two equations above contain only the \(p_1\) noise, \(\zeta _1\) is found by the following expression:

$$\begin{aligned} \zeta _1 = \left[ D_{3^{\prime}}D_{2^{\prime}} - D_{1}\right] \left( \eta _{3,3} - \eta _{3',3} + \eta _{1,1'}\right) - \left[ D_{3}D_{2} - D_{1^{\prime}}\right] \left( \eta _{1',1} - \eta _{2,2'} + \eta _{2',2'}\right) . \end{aligned}$$
(85)

If the light-times in the arms are equal in the clockwise and counter-clockwise senses (e.g., no rotation) there is no distinction between primed and unprimed delay times. In this case, \(\zeta _1\) is related to our original symmetric Sagnac \(\zeta \) by \(\zeta _1 = \zeta _{,23} - \zeta _{,1}\). Thus, for the LISA case (arm length difference \(< 1\%\)), the SNR of \(\zeta _1\) will be the same as the SNR of \(\zeta \).

If the delay-times also change with time, the perfect cancellation of the laser noises is no longer achieved in the \((\zeta _1, \zeta _2, \zeta _3)\) combinations. However, it has been shown in Tinto et al. (2004) that the magnitude of the residual laser noises in these combinations are significantly smaller than the LISA secondary system noises, making their effects negligible.

The expressions for the Monitor, Beacon, and Relay combinations, accounting for the rotation and flexing of the LISA array, have been derived in the literature (Tinto et al. 2004) by applying the time-delay iterative procedure highlighted in this section. The interested reader is referred to that paper for details.

Second-generation TDI (or TDI-2) adequately suppress the laser noise below the levels identified by the secondary noises in LISA as well as in other missions currently under development (Luo et al. 2016; Hu and Wu 2017). The reader, however, might wonder whether TDI could successfully be implemented with arrays whose inter-spacecraft velocities are such as to make laser-noise terms quadratic in the inter-spacecraft velocities and the accelerations no longer negligible. Assuming such hypothetical missions could still perform heterodyne one-way Doppler measurements in the presence of extremely large laser beat-notes, one way to address this problem might be to further extend to higher orders the iterative procedure presented in this section. This approach would result in “TDI-3”, or in general, a “TDI-n” formulation of TDI.

Alternatively, two interesting techniques (Muratore et al. 2020; Vallisneri et al. 2020) have recently been proposed to address this problem. Both are numerically based and can in principle identify TDI combinations capable of canceling the laser noises in the presence of arbitrarily large inter-spacecraft velocities. The one proposed in Vallisneri et al. (2020) in particular follows from a Bayesian formulation of TDI (Romano and Woan 2006) in which the laser noises are characterized by having variances much larger than those associated with all other noises. The method relies on using the entire temporal duration of the one-way Doppler measurements and by relating them to the laser noises through an identified linear relationship. In so doing it is stated that data-gaps can be treated automatically within the method, and any gravitational wave filtering can be applied directly to the one-way Doppler data.The method has been termed as TDI-\(\infty \). However, here a proof of principle analysis has been described for two data streams, which in the case of LISA needs to be extended to six and this situation could technically be more complex. It would be desirable to perform a rigorous and extensive analysis of this method for future missions. The method looks promising and we are looking forward to see the future outcomes of this approach.

5.3 Algebraic approach to second-generation TDI

In this subsection we present a mathematical formulation of the “second-generation” TDI, which generalizes the one presented in Sect. 4 for stationary LISA. Although a full solution as in the case of stationary LISA seems difficult to obtain, significant progress can be made.

In the case of rigidly rotating array, in which only the Sagnac effect is considered (Cornish and Hellings 2003; Shaddock 2004), the up-down links are unequal and the delay-times remain constant. The mathematical formulation of Sect. 4 can be extended in a straight-forward way where now the six time-delays \({\mathscr {D}}_{i}\) and \({\mathscr {D}}_{i'}\) must be taken into account. The polynomial ring still remains commutative but it is now in six variables. The corresponding module of syzygies can be constructed over this larger polynomial ring (Rajesh Nayak and Vinet 2005).

When the arms are allowed to flex, that is, the operators themselves are functions of time, the operators no longer commute. One must then resort to non-commutative algebra. We outline the procedure below. Since lot of the discussion has been covered in the previous subsections we just describe the algebraic formulation. Equation (28) generalizes in two ways: (1) now we need to consider six operators \({\mathscr {D}}_i\) and \({\mathscr {D}}_{i'}\), and (2) we need to take into account the non-commutativity of the operators—the order of the operators is important. Accordingly Eq. (28) generalizes to,

$$\begin{aligned} q_1 + q_{1^{\prime}} - q_{2^{\prime}} {\mathscr {D}}_{3^{\prime}} - q_3 {\mathscr {D}}_2& = 0 , \\ q_2 + q_{2^{\prime}} - q_{3^{\prime}} {\mathscr {D}}_{1^{\prime}} - q_1 {\mathscr {D}}_3& = 0 , \\ q_3 + q_{3^{\prime}} - q_{1^{\prime}} {\mathscr {D}}_{2^{\prime}} - q_2 {\mathscr {D}}_1& = 0 . \end{aligned}$$
(86)

Since Eq. (86) are linear equations, they define a homomorphism of modules as follows: Eliminating \(q_1\) and \(q_2\) from the three Eqs. (86) while respecting the order of the variables we get,

$$\begin{aligned}&q_3 (1 - {\mathscr {D}}_2 {\mathscr {D}}_3 {\mathscr {D}}_1) + q_{1^{\prime}} ({\mathscr {D}}_{2^{\prime}} - {\mathscr {D}}_3 {\mathscr {D}}_1) + q_{2^{\prime}} ({\mathscr {D}}_{3^{\prime}} {\mathscr {D}}_3 - 1) {\mathscr {D}}_1 + q_{3^{\prime}} ({\mathscr {D}}_{1^{\prime}} {\mathscr {D}}_1 - 1) \, \\&\quad \equiv \psi ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3, {\mathscr {D}}_{1^{\prime}}, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_{3^{\prime}}) = 0 . \end{aligned}$$
(87)

Consider the polynomial ring \({\mathscr {Q}}({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3, {\mathscr {D}}_{1^{\prime}}, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_{3^{\prime}}) \equiv {\mathscr {K}}'\), in general non-commutative, of polynomials in the six variables \({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3, {\mathscr {D}}_{1^{\prime}}, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_{3^{\prime}}\) and coefficients in the rational field \({\mathscr {Q}}\). The operators \({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3, {\mathscr {D}}_{1^{\prime}}, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_{3^{\prime}}\) play the role of indeterminates. Equation (87) defines a homomorphism \(\varphi : {\mathscr {K}}'^4 \longrightarrow {\mathscr {K}}'\) where any polynomial vector \((q_3, q_{1^{\prime}}, q_{2^{\prime}}, q_{3^{\prime}}) \in {\mathscr {K}}'^4\) is mapped to the polynomial \(\psi ({\mathscr {D}}_1, {\mathscr {D}}_2, {\mathscr {D}}_3, {\mathscr {D}}_{1^{\prime}}, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_{3^{\prime}}) \in {\mathscr {K}}'\). The set of noise free TDI combinations is just the kernel of this homomorphism, namely, \(\varphi ^{-1} (0) \subset {\mathscr {K}}'^4\) which is a submodule of \({\mathscr {K}}'^4\). This is again a first module of syzygies over the polynomial ring \({\mathscr {K}}'\). This homomorphism can be extended to \({\mathscr {K}}'^6\) via the elimination equations for \(q_1\) and \(q_2\). Thus one obtains a module of noise-free TDI observables \({\mathscr {M}}\subset {\mathscr {K}}'^6\) which is isomorphic to \(\varphi ^{-1} (0)\). The polynomial vectors \((q_i, q_i')\) satisfying the above equations form a left module over \({\mathscr {K}}'\). A left module means that one can multiply a solution \((q_i, q_i')\) from the left by any polynomial in \({\mathscr {K}}'\), then it is also a solution to the Eqs. (86) and, therefore, in the module—the module of noise-free TDI observables. For details see Dhurandhar (2009).

When the operators do not commute, the algebraic problem is far more complex. If we follow on the lines of the commutative case, the first step would be to find a Gröbner basis for the ideal generated by the coefficients appearing in Eq. (87), namely, the set of polynomials:

$$\begin{aligned} \{1 - {\mathscr {D}}_2 {\mathscr {D}}_3 {\mathscr {D}}_1, {\mathscr {D}}_{2^{\prime}} - {\mathscr {D}}_3 {\mathscr {D}}_1, ({\mathscr {D}}_{3^{\prime}} {\mathscr {D}}_3 - 1) {\mathscr {D}}_1, {\mathscr {D}}_{1^{\prime}} {\mathscr {D}}_1 - 1 \}. \end{aligned}$$

Although we may be able to apply non-commutative Gröbner basis methods, the general solution seems quite difficult. However, simplifications are possible because of the inherent approximate symmetries in the problem and so the ring \({\mathscr {K}}'\), now non-commutative, can be quotiented by a certain ideal, simplifying the algebraic problem. One then needs to deal with a ‘smaller’ ring, which may be easier to deal with. We describe below how this can achieved with the help of certain commutators.

The level of non-commutativity can be found by computing commutators which occur in several of the well known TDI observables like the Michelson, Sagnac etc. We find that given our model of LISA, we require to go only up to the first order in \({\dot{L}}\). The reason is as follows. In Nayak et al. (2006) an optimal model of LISA was given which minimizes the flexing of the arms of LISA. The computations were done for the older model of LISA which had arm length of \(5 \times 10^6\) km. Since for the current design of LISA, the arm length has been reduced by a factor of 2, we apply scaling arguments to the results obtained therein. In Nayak et al. (2006) it was shown that the flexing of the arms of LISA can be derived from the geodesic deviation equation where one considers the Sun’s field in the Newtonian approximation. One must however, go to the second order in the separation vector in the geodesic deviation equation—the third derivative of the Newtonian potential of the Sun, expansion up-to the octupole term—in order to obtain the flexing. Since the separation vector comes in at the second order in the geodesic deviation equation, the \(\dot{L}\) and \({\ddot{L}}\) terms are scaled down by a factor of 4. We find that for the current design of LISA, \({\ddot{L}} \sim 3 \times 10^{-7}\mathrm {\ m/s}^2\) and \({\dot{L}}^2 / c^2 \sim 10^{-16}\). Since our aim is to cancel the laser frequency noise to better than 1 part in \(10^8\) a difference in path length of less than a few meters is sufficient to keep the level of the residual laser frequency noise below acceptable levels. Clearly, the \({\dot{L}}^2\) terms can be dropped, because their contribution is much smaller than 1 part in \(10^8\). In the case of the \({\ddot{L}}\), even if one considers say 12 successive optical paths, that is, about \(\Updelta t \sim 100\,\hbox {s}\) of light travel time, \(\Updelta t^2 {{\ddot{L}}} \sim 3 \times 10^{-3}\,\hbox {m}\). This is well below few meters and thus can be neglected in the residual laser noise computation. The calculations that follow neglect these terms.

Applying the operators twice in succession and dropping higher order terms as explained above,

$$\begin{aligned} k_2 k_1 \phi& = \phi (t - L_{k_1} (t - L_{k_2}) - L_{k_2}) , \\\approx & {} \phi (t -L_{k_1} - L_{k_2}) + L_{k_2} {\dot{L}}_{k_1} {\dot{\phi }} (t -L_{k_1} - L_{k_2}) , \end{aligned}$$
(88)

where \(k_1, k_2\) are any of the operators \({\mathscr {D}}_i, {\mathscr {D}}_{i'}\). The above formula can be easily generalized by induction to n operators. In several of the TDI observables, commutators play a major role. In general, a commutator of two operators xy is defined as the operator \([x, y] \equiv xy - yx\). In our situation x and y are strings of operators built up of the operators \(D_1, D_2, D_3, D_{1^{\prime}}, D_{2^{\prime}}, D_{3^{\prime}}\). For example, in the Sagnac combination the following commutator (Dhurandhar et al. 2008) occurs:

$$\begin{aligned}{}[D_1 D_2 D_3, D_{1^{\prime}} D_{2^{\prime}} D_{3^{\prime}}] \equiv D_1 D_2 D_3 D_{1^{\prime}} D_{2^{\prime}} D_{3^{\prime}} - D_{1^{\prime}} D_{2^{\prime}} D_{3^{\prime}} D_1 D_2 D_3 . \end{aligned}$$
(89)

Here \(x = D_1 D_2 D_3\) and \(y = D_{1^{\prime}} D_{2^{\prime}} D_{3^{\prime}}\). This commutator leads to the residual noise term given in Eq. (82) and which happens to be small. For the Sagnac combination, we obtain,

$$\begin{aligned} {[}D_1 D_2 D_3, D_{1^{\prime}} D_{2^{\prime}} D_{3^{\prime}}]& = (L_1 + L_2 + L_3)({\dot{L}}_1' + {\dot{L}}_2' + {\dot{L}}_3') \\&\quad - (L_1' + L_2' + L_3')({\dot{L}}_1 + {\dot{L}}_2 + {\dot{L}}_3) , \end{aligned}$$
(90)

where it is understood that the LHS acts on \(\phi \) while the RHS multiplies \({{\dot{\phi }}}\) at an appropriately delayed time. Note that the term in \(\phi \) cancels out on the RHS. The near vanishing of the above commutators implies that a vast simplification in the algebra is possible.

The result for the Sagnac combination can be generalized. To simplify notation we write \(x_k\) or \(y_m\) for the time-delay operators, where \(k, m = 1, 2, \dots , n\) and \(n \ge 2\), that is, \(x_k\) or \(y_m\) are any of the operators \(D_1, D_2, D_3, D_{1^{\prime}}, D_{2^{\prime}}, D_{3^{\prime}}\). Then a commutator is:

$$\begin{aligned}{}[x_1 x_2 \dots x_n, y_1 y_2 \dots y_n] = x_1 x_2 \dots x_n y_1 y_2 \dots y_n - y_1 y_2 \dots y_n x_1 x_2 \dots x_n . \end{aligned}$$
(91)

Up to the order of approximation we are working in, we compute the effect of the commutator on the phase \(\phi (t)\):

$$\begin{aligned}&[x_1 x_2 \dots x_n, y_1 y_2 \dots y_n] \,\phi (t) \, \\&\quad = \left( \sum _{k = 1}^n L_{x_k} \sum _{m = 1}^n {\dot{L}}_{y_m} - \sum _{m = 1}^n L_{y_m} \sum _{k = 1}^n {\dot{L}}_{x_k} \right) \, {{\dot{\phi }}} \left( t - \sum _{k = 1}^n L_{x_k} - \sum _{m = 1}^n L_{y_m} \right) . \end{aligned}$$
(92)

Note that the LHS acts on \(\phi (t)\), while the right-hand side multiplies \({{\dot{\phi }}}\) at an appropriately delayed time. Also the notation on the RHS is obvious: if for some k, we have, \(x_k = {\mathscr {D}}_{2^{\prime}}\) say, then \(L_{x_k} = L_2'\) and so on; the same holds for \(y_m\) for a given m. From this equation it immediately follows that if the operators \(y_1, y_2, \dots , y_n\) are a permutation of the operators \(x_1, x_2, \dots , x_n\), then the commutator,

$$\begin{aligned}{}[x_1 x_2 \dots x_n, y_1 y_2 \dots y_n] = 0 , \end{aligned}$$
(93)

up to the order we are working in. We can understand this by the following argument. If \(y_1 \dots y_n\) is a permutation of \(x_1 \dots x_n\) then both polynomials trace the same links, except that the nodes (spacecraft) of the links are taken in different orders. If the arm lengths were constant, the pathlengths would be identical and the commutator would be zero. But here, by neglecting \({{\ddot{L}}}\) terms and those of higher orders, we have effectively assumed that \({\dot{L}}\)s are constant, so the increments also cancel out, resulting in a vanishing commutator.

These vanishing commutators (in the approximation we are working in) can be used to simplify the algebra. We first construct the ideal \({\mathscr {U}}\) generated by the commutators such as those given by Eq. (93). Then we quotient the ring \({\mathscr {K}}'\) by \({\mathscr {U}}\), thereby constructing a smaller ring \({\mathscr {K}}'/{\mathscr {U}}\equiv {\bar{{\mathscr {K}}}}\). This ring is smaller because it has fewer distinct terms in a polynomial. Although, this reduces the complexity of the problem, a full solution to the TDI problem is still lacking.

In the following Sect. 5.4, we will consider the case where we have only two arms of LISA in operation, that is one arm is nonfunctional. The algebraic problem simplifies considerably and it turns out to be tractable.

5.4 Solutions with one arm nonfunctional

We must envisage the possibility that not all optical links of LISA can be operating at all times for various reasons like technical failure for instance or even the operating costs. An analysis covering the scientific capabilities achievable by LISA in the eventuality of loosing one and two links has been discussed in Vallisneri et al. (2008). Here we obtain the TDI combinations when one entire arm becomes dysfunctional. See Dhurandhar et al. (2010) for a full discussion. The results of this section are directly usable by the TaiJi (Hu and Wu 2017) and TianQin (Luo et al. 2016) missions.

We arbitrarily choose the non-functional arm to be the one connecting S/C 2 and S/C 3. This means from our labeling that the polynomials are now restricted to only the four operators \( {\mathscr {D}}_2, {\mathscr {D}}_{2^{\prime}}, {\mathscr {D}}_3, {\mathscr {D}}_{3^{\prime}}\) and we can set the polynomials \(q_2 = q_{3^{\prime}} = 0\). This simplifies the Eq. (86) considerably because the last two equations reduce to \(q_{2^{\prime}} = q_1 {\mathscr {D}}_3\) and \(q_3 = q_{1^{\prime}} {\mathscr {D}}_{2^{\prime}}\). Substituting these in the first equation gives just one equation:

$$\begin{aligned} q_1 (1 - {\mathscr {D}}_3 {\mathscr {D}}_{3^{\prime}}) + q_{1^{\prime}} (1 - {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_2) = 0 , \end{aligned}$$
(94)

If we can solve this equation for \(q_1, q_{1^{\prime}}\) then the full polynomial vector can be obtained because \(q_{2^{\prime}} = q_1 {\mathscr {D}}_3\) and \(q_3 = q_{1^{\prime}} {\mathscr {D}}_{2^{\prime}}\). It is clear that solutions are of the Michelson type. Also notice that the coefficients of this equation has the operators \(a = {\mathscr {D}}_3 {\mathscr {D}}_{3^{\prime}}\) and \(b = {\mathscr {D}}_{2^{\prime}} {\mathscr {D}}_2\) occurring in them. So the solutions \(q_1, q_{1^{\prime}}\) too will be in terms of a and b only. Physically, the operators a and b correspond to round trips.

One solution has already been given in the literature (Armstrong 2006; Vallisneri 2005). This solution in terms of ab is:

$$\begin{aligned} q_1& = 1 - b - ba + ab^2 , \\ q_{1^{\prime}}& = - (1 - a - ab + ba^2) . \end{aligned}$$
(95)

Writing,

$$\begin{aligned} \Updelta = q_1 (1 - a) + q_{1^{\prime}} (1 - b) , \end{aligned}$$
(96)

we get for Eq. (95), \(\Updelta = [ba, ab]\), which is a commutator that vanishes since ab is a permutation of ba. Thus, it is an element of \({\mathscr {U}}\) and Eq. (95) is a solution (over the quotient ring).

What we would like to emphasize is that there are more solutions of this type—in fact there are infinite number of such solutions. The solutions are based on vanishing of commutators. In Dhurandhar et al. (2010), such commutators are enumerated and for each such commutator there is a corresponding solution. Further an algorithm is given to construct such solutions. We briefly mention some results given in Dhurandhar et al. (2010). We start with the solutions \(q_1, q_{1^{\prime}}\) of Eq. (95). Note that these are of degree 3 in a and b. The commutator corresponding to this solution is \(\Updelta = [ba, ab] = ba^2b - a b^2 a\) and is of degree 4. There is only one such commutator at degree 4 and therefore one solution. The next higher degree solutions are found when the commutators have degree 8. The solutions \(q_1, q_{1^{\prime}}\) are of degree 7. There are three such commutators at degree 8. We demonstrate a systematic way to enumerate and list such solutions.

To obtain solutions with higher degree, we proceed as follows: We write,

$$\begin{aligned} r_1& = q_1 - b a^2 b , \\ r_1'& = q_{1^{\prime}} + a b^2 a , \end{aligned}$$
(97)

and check whether we get a commutator. We find,

$$\begin{aligned} \Updelta ' \equiv r_1 (1 - a) + r_1' (1 - b) = b a^2 ba - a b^2 ab . \end{aligned}$$
(98)

The lower degree terms cancel out to yield \(\Updelta '\) as above. However, it does not vanish—it is also not a commutator—and hence we conclude that \(r_1, r_1'\) is not a solution. The prime on \(\Updelta \) denotes the first iteration. Continuing further we write,

$$\begin{aligned} s_1& = r_1 + a b^2 a b , \\ s_1'& = r_1' - b a^2 b a, \end{aligned}$$
(99)

and similarly compute \(\Updelta '' \equiv s_1 (1 - a) + s_1' (1 - b) = b a^2 bab - a b^2 aba\) which also does not vanish, and again not a commutator. The double prime denotes the second iteration.

This process can be continued in a way so that the lower degree terms cancel out until we get a vanishing commutator, albeit, of higher degree. The next commutator is obtained by adding two more terms of successively higher degree in a similar way. The result is:

$$\begin{aligned} q_1& = 1 - b - ba + ab^2 - b a^2 b + a b^2 ab + ab^2 aba - ba^2 bab^2 , \\ q_{1^{\prime}}& = -(1 - a - ab + ba^2 - a b^2 a + b a^2 ba + ba^2 bab - ab^2 aba^2) . \end{aligned}$$
(100)

This results in the following commutator:

$$\begin{aligned} \Updelta = [a b^2 a, b a^2 b] , \end{aligned}$$
(101)

and hence (100) is a solution.

Another and the second solution is:

$$\begin{aligned} q_1& = 1 + a - b^2 - b^2 a - b^2 a^2 - b^2 a^3 + a^2 b^4 + a^2 b^4 a \\& = (1 - b^2 - b^2 a^2 + a^2 b^4)(1 + a) , \\ q_{1^{\prime}}& = -(1 + b - a^2 - a^2 b - a^2 b^2 - a^2 b^3 + b^2 a^4 + b^2 a^4 b) \\& = -(1 - a^2 - a^2 b^2 + b^2 a^4)(1 + b) , \end{aligned}$$
(102)

whose commutator is \(\Updelta = [b^2 a^2, a^2 b^2]\).

The third solution corresponds to the commutator [babaabab]. This solution is given by:

$$\begin{aligned} q_1& = 1 - b + ab - bab - baba + abab^2 - baba^2 b + abab^2 ab , \\ q_{1^{\prime}}& = -(1 - a + ba - aba - abab + baba^2 - abab^2 a + baba^2 ba) . \end{aligned}$$
(103)

Higher degree solutions can be constructed. An iterative algorithm has been described in Dhurandhar et al. (2010) for this purpose. The degrees of the commutators are in multiples of 4. If we call the degree of the commutators as 4n where \(n = 1, 2, \dots \), then the solutions \(q_1, q_{1^{\prime}}\) are of degree \(4n - 1\). The cases mentioned above, correspond to \(n = 1\) and \(n = 2\). The general formula for the number of commutators of degree 4n is \(^{2n -1} C_{n - 1}\). So at \(n = 3\) we have 10 commutators and so as many solutions \(q_1, q_{1^{\prime}}\) of degree 11.

These are the degrees of polynomials of \(q_1, q_{1^{\prime}}\) in the operators ab. But for the full polynomial vector, which has \(q_{2^{\prime}}\) and \(q_3\), we need to go over to the operators expressed in terms of \({\mathscr {D}}_i, {\mathscr {D}}_{i'}\). Then the degree of each of the \(q_1, q_{1^{\prime}}\) is doubled to \(8n - 2\), while \(q_{2^{\prime}}\) and \(q_3\) are each of degree \(8n - 1\). Thus, for a general value of n, the solution contains polynomials of maximum degree \(8n - 1\) in the time-delay operators.

From the mathematical point of view there is an infinite family of solutions. Note that no claim is made on exhaustive listing of solutions. However the family of solutions is sufficiently rich, because we can form linear combinations of these solutions and they also are solutions.

From the physical point of view, since terms in \({{\ddot{L}}}\) and \({\dot{L}}^2\) and higher orders have been neglected, a limit on the degree of the polynomial solutions arises. That is up to certain degree of the polynomials, we can safely assume the commutators to vanish. But as the degree of the polynomials increases it is not possible to neglect these higher order terms any longer and then such a limit becomes important. The limit is essentially set by \({{\ddot{L}}}\). We now investigate this limit and make a very rough estimate of it. From \({\ddot{L}} \sim 3 \times 10^{-3}\) m/s\(^2\), we compute the error in L, namely, \(\Updelta L \sim \frac{1}{2} \Updelta t^2 {{\ddot{L}}}\), where \(\Updelta t\) is the time the laser light takes to travel a specific path. If we allow the error \(\Updelta L\) to be no more than say 3 m, then we find \(\Updelta t \sim 4500\mathrm {\ s}\). Since each time-delay is about 8.3 s for LISA, the number of successive time-delays is about 540. This is the maximum degree of the polynomials. This means one can go up to \(n \lesssim 60\). If we set the limit more stringently at \(\Updelta L \sim 0.3 \mathrm {\ m}\), then the highest degree of the polynomial reduces to about 160, which means one can go up to \(n = 20\). Thus, there are a large number of TDI observables available to do the physics.

Some remarks are in order:

  • A geometric combinatorial approach was adopted in Vallisneri (2005) where several solutions were presented. Our approach is algebraic where the operations are algebraic operations on strings of operators. The algebraic approach has the advantage of easy manipulation of data streams while still displaying a geometrical interpretation in the derived TDI combinations.

  • Another important aspect is the GW response of such TDI observables. The GW response to a TDI observable may be calculated in the simplest way by assuming equal arms (the possible differences in lengths would be sensitive to frequencies outside the LISA bandwidth). This leads in the Fourier domain to polynomials in the same phase factor from which the signal to noise ratio can be found. A comprehensive and generic treatment of the responses of second-generation TDI observables can be found in Królak et al. (2004).

6 Optimal LISA sensitivity

All the above interferometric combinations have been shown to individually have rather different sensitivities (Estabrook et al. 2000), as a consequence of their different responses to gravitational radiation and system noises. Since LISA has the capability of simultaneously observing a gravitational-wave signal with many different interferometric combinations (all having different antenna patterns and noises), we should no longer regard LISA as a single detector system but rather as an array of gravitational-wave detectors working in coincidence. This suggests that the LISA sensitivity could be improved by optimally combining elements of the TDI space.

Before proceeding with this idea, however, let us consider again the so-called “second-generation” TDI Sagnac observables: \((\alpha _1, \alpha _2, \alpha _3)\). The expressions of the gravitational-wave signal and the secondary noise sources entering into \(\alpha _1\) will in general be different from those entering into \(\alpha \), the corresponding Sagnac observable derived under the assumption of a stationary LISA array (Armstrong et al. 1999; Estabrook et al. 2000). However, the other remaining secondary noises in LISA are so much smaller, and the rotation and systematic velocities in LISA are so intrinsically small, that index permutation may still be done for them (Tinto et al. 2004). It is therefore easy to derive the following relationship between the signal and secondary noises in \(\alpha _1\), and those entering into the stationary TDI combination \(\alpha \) (Shaddock et al. 2003; Tinto et al. 2004),

$$\begin{aligned} \alpha _1 (t) \simeq \alpha (t) - \alpha (t - L_1 - L_2 - L_3) , \end{aligned}$$
(104)

where \(L_i\), \(i=1, 2, 3\), are the unequal-arm lengths of the stationary LISA array. Equation (104) implies that any data analysis procedure and algorithm that will be implemented for the second-generation TDI combinations can actually be derived by considering the corresponding “first-generation” TDI combinations. For this reason, from now on we will focus our attention on the gravitational-wave responses of the first-generation TDI observables \((\alpha , \beta , \gamma , \zeta )\).

As a consequence of these considerations, we can still regard \((\alpha , \beta , \gamma , \zeta )\) as the generators of the TDI space, and write the most general expression for an element of the TDI space, \(\eta (f)\), as a linear combination of the Fourier transforms of the four generators \(({{\widetilde{\alpha }}}, {{\widetilde{\beta }}}, {\widetilde{\gamma }}, {{\widetilde{\zeta }}})\),

$$\begin{aligned} \eta (f) \equiv a_1 (f, {\varvec{\lambda }}) \, {{\widetilde{\alpha }}} (f) + a_2 (f, {\varvec{\lambda }}) \, {{\widetilde{\beta }}} (f) + a_3 (f, {\varvec{\lambda }}) \, {{\widetilde{\gamma }}} (f) + a_4 (f, {\varvec{\lambda }}) \, {{\widetilde{\zeta }}} (f) , \end{aligned}$$
(105)

where the \(\{a_i (f, {\varvec{\lambda }})\}^4_{i=1}\) are arbitrary complex functions of the Fourier frequency f, and of a vector \({\varvec{\lambda }}\) containing parameters characterizing the gravitational-wave signal (source location in the sky, waveform parameters, etc.) and the noises affecting the four responses (noise levels, their correlations, etc.). For a given choice of the four functions \(\{a_i \}^4_{i=1}\), \(\eta \) gives an element of the functional space of interferometric combinations generated by \((\alpha , \beta , \gamma , \zeta )\). Our goal is therefore to identify, for a given gravitational-wave signal, the four functions \(\{a_i \}^4_{i=1}\) that maximize the signal-to-noise ratio \(\mathrm {SNR}_{\eta }^2\) of the combination \(\eta \),

$$\begin{aligned} \mathrm {SNR}_{\eta }^2 = \int _{f_\mathrm {l}}^{f_\mathrm {u}} \frac{\left| a_1 \, {{\widetilde{\alpha }}_\mathrm {s}} + a_2 \, {{\widetilde{\beta }}_\mathrm {s}} + a_3 \, {{\widetilde{\gamma }}_\mathrm {s}} + a_4 {{\widetilde{\zeta }}_\mathrm {s}} \right| ^2}{\left\langle \left| a_1 \, {{\widetilde{\alpha }}_\mathrm {n}} + a_2 \, {{\widetilde{\beta }}_\mathrm {n}} + a_3 \, {{\widetilde{\gamma }}_\mathrm {n}} + a_4 \, {{\widetilde{\zeta }}_\mathrm {n}} \right| ^2 \right\rangle } \, df . \end{aligned}$$
(106)

In Eq. (106) the subscripts \(\mathrm {s}\) and \(\mathrm {n}\) refer to the signal and the noise parts of \(({\widetilde{\alpha }}, {\widetilde{\beta }}, {\widetilde{\gamma }}, {\widetilde{\zeta }})\), respectively, the angle brackets represent noise ensemble averages, and the interval of integration \((f_\mathrm {l}, f_\mathrm {u})\) corresponds to the frequency band accessible by LISA.

Before proceeding with the maximization of the \(\mathrm {SNR}_{\eta }^2\) we may notice from Eq. (62) that the Fourier transform of the totally symmetric Sagnac combination, \({\widetilde{\zeta }}\), multiplied by the transfer function \(1 - e^{2 \pi i f (L_1 + L_2 + L_3)}\) can be written as a linear combination of the Fourier transforms of the remaining three generators \(({\widetilde{\alpha }}, {\widetilde{\beta }}, {\widetilde{\gamma }})\). Since the signal-to-noise ratio of \(\eta \) and \((1 - e^{2 \pi i f (L_1 + L_2 + L_3)}) \eta \) are equal, we may conclude that the optimization of the signal-to-noise ratio of \(\eta \) can be performed only on the three observables \(\alpha , \beta , \gamma \). This implies the following redefined expression for \(\mathrm {SNR}_{\eta }^2\):

$$\begin{aligned} \mathrm {SNR}_{\eta }^2 = \int _{f_\mathrm {l}}^{f_\mathrm {u}} \frac{\left| a_1 \, {{\widetilde{\alpha }}_\mathrm {s}} + a_2 \, {{\widetilde{\beta }}_\mathrm {s}} + a_3 \, {{\widetilde{\gamma }}_\mathrm {s}} \right| ^2}{\left\langle \left| a_1 \, {{\widetilde{\alpha }}_\mathrm {n}} + a_2 \, {{\widetilde{\beta }}_\mathrm {n}} + a_3 \, {{\widetilde{\gamma }}_\mathrm {n}} \right| ^2 \right\rangle } \, df . \end{aligned}$$
(107)

The \(\mathrm {SNR}_{\eta }^2\) can be regarded as a functional over the space of the three complex functions \(\{ a_i \}^3_{i=1}\), and the particular set of complex functions that extremize it can of course be derived by solving the associated set of Euler–Lagrange equations.

To make the derivation of the optimal SNR easier, let us first denote by \(\mathbf{x}^{(\mathrm {s})}\) and \(\mathbf{x}^{(\mathrm {n})}\) the two vectors of the signals \(({\widetilde{\alpha }}_\mathrm {s}, {\widetilde{\beta }}_\mathrm {s}, {\widetilde{\gamma }}_\mathrm {s})\) and the noises \(({\widetilde{\alpha }}_\mathrm {n}, {\widetilde{\beta }}_\mathrm {n}, {\widetilde{\gamma }}_\mathrm {n})\), respectively. Let us also define \(\mathbf{a}\) to be the vector of the three functions \(\{a_i \}^3_{i=1}\), and denote with \(\mathbf{C}\) the Hermitian, non-singular, correlation matrix of the vector random process \(\mathbf{x}_\mathrm {n}\),

$$\begin{aligned} (\mathbf{C})_{rt} \equiv \left\langle \mathbf{x}^{(n)}_{r} \mathbf{x}^{(n)*}_{t} \right\rangle . \end{aligned}$$
(108)

If we finally define \((\mathbf{A})_{ij}\) to be the components of the Hermitian matrix \(\mathbf{x}^{(\mathrm {s})}_i \mathbf{x}^{(\mathrm {s})*}_j\), we can rewrite \(\mathrm {SNR}_{\eta }^2\) in the following form,

$$\begin{aligned} \mathrm {SNR}_{\eta }^2 = \int _{f_\mathrm {l}}^{f_\mathrm {u}} \frac{\mathbf{a}_i \mathbf{A}_{ij} \mathbf{a}^*_j}{\mathbf{a}_r \mathbf{C}_{rt} \mathbf{a}^*_t} \, df , \end{aligned}$$
(109)

where we have adopted the usual convention of summation over repeated indices. Since the noise correlation matrix \(\mathbf{C}\) is non-singular, and the integrand is positive definite or null, the stationary values of the signal-to-noise ratio will be attained at the stationary values of the integrand, which are given by solving the following set of equations (and their complex conjugated expressions):

$$\begin{aligned} \frac{\partial }{\partial \mathbf{a}_k} \left[ \frac{\mathbf{a}_i \mathbf{A}_{ij} \mathbf{a}^*_j}{\mathbf{a}_r \mathbf{C}_{rt} \mathbf{a}^*_t} \right] = 0 , \qquad k = 1, 2, 3 . \end{aligned}$$
(110)

After taking the partial derivatives, Eq. (110) can be rewritten in the following form,

$$\begin{aligned} (\mathbf{C}^{-1})_{ir} \, (\mathbf{A})_{rj} \, (\mathbf{a}^*)_j = \left[ \frac{\mathbf{a}_p \mathbf{A}_{pq} \mathbf{a}^*_q}{\mathbf{a}_l \mathbf{C}_{lm} \mathbf{a}^*_m} \right] (\mathbf{a}^*)_i , \qquad i = 1, 2, 3 , \end{aligned}$$
(111)

which tells us that the stationary values of the signal-to-noise ratio of \(\eta \) are equal to the eigenvalues of the the matrix \(\mathbf{C^{-1} \cdot A}\). The result in Eq. (110) is well known in the theory of quadratic forms, and it is called Rayleigh’s principle (Noble 1969; Selby 1964).

Now, to identify the eigenvalues of the matrix \(\mathbf{C^{-1} \cdot A}\), we first notice that the \(3 \times 3\) matrix \(\mathbf{A}\) has rank 1. This implies that the matrix \(\mathbf{C}^{-1} \cdot \mathbf{A}\) has also rank 1, as it is easy to verify. Therefore two of its three eigenvalues are equal to zero, while the remaining non-zero eigenvalue represents the solution we are looking for.

The analytic expression of the third eigenvalue can be obtained by using the property that the trace of the \(3 \times 3\) matrix \(\mathbf{C}^{-1} \cdot \mathbf{A}\) is equal to the sum of its three eigenvalues, and in our case to the eigenvalue we are looking for. From these considerations we derive the following expression for the optimized signal-to-noise ratio \(\mathrm {SNR}^2_{\eta \,\mathrm {opt}}\):

$$\begin{aligned} \mathrm {SNR}^2_{\eta \,\mathrm {opt}} = \int _{f_\mathrm {l}}^{f_\mathrm {u}} \mathbf{x}^{(\mathrm {s})*}_i \, (\mathbf{C}^{-1})_{ij} \, \mathbf{x}^{(\mathrm {s})}_j \, df . \end{aligned}$$
(112)

We can summarize the results derived in this section, which are given by Eqs. (107) and (112), in the following way:

  1. 1.

    Among all possible interferometric combinations LISA will be able to synthesize with its four generators \(\alpha \), \(\beta \), \(\gamma \), \(\zeta \), the particular combination giving maximum signal-to-noise ratio can be obtained by using only three of them, namely \((\alpha , \beta , \gamma ).\)

  2. 2.

    The expression of the optimal signal-to-noise ratio given by Eq. (112) implies that LISA should be regarded as a network of three interferometer detectors of gravitational radiation (of responses \((\alpha , \beta , \gamma )\)) working in coincidence (Finn 2001; Rajesh Nayak et al. 2003b).

6.1 General application

As an application of Eq. (112), here we calculate the sensitivity that LISA can reach when observing signals (i) from a specific direction and with a given polarization state as well as (ii) in the case of an ensemble of sources uniformly distributed on the celestial sphere and of random polarization.

To calculate the optimal signal-to-noise ratio we will also need to use a specific expression for the noise correlation matrix \(\mathbf{C}\). As a simplification, we will assume the LISA arm lengths to be equal to their nominal value \(L = 8.33 \mathrm {\ s}\), the optical-path noises to be equal and uncorrelated to each other and entering in the TDI responses with the same transfer function as the shot-noise, and finally the proof-mass noises to be also equal, uncorrelated to each other and to the optical-path noises Amaro-Seoane et al. (2017).Footnote 4 Under these assumptions the correlation matrix becomes real, its three diagonal elements are equal, and all the off-diagonal terms are equal to each other, as it is easy to verify by direct calculation (Estabrook et al. 2000). The noise correlation matrix \(\mathbf{C}\) is therefore uniquely identified by two real functions \(S_{\alpha }\) and \(S_{\alpha \beta }\) in the following way:

$$\begin{aligned} \mathbf{C} = \left( \begin{array}{lll} S_{\alpha } &{} S_{\alpha \beta } &{} S_{\alpha \beta } \\ S_{\alpha \beta } &{} S_{\alpha } &{} S_{\alpha \beta } \\ S_{\alpha \beta } &{} S_{\alpha \beta } &{} S_{\alpha } \end{array} \right) . \end{aligned}$$
(113)

The expression of the optimal signal-to-noise ratio assumes a rather simple form if we diagonalize this correlation matrix by properly “choosing a new basis”. There exists an orthogonal transformation of the generators \(({\widetilde{\alpha }}, {\widetilde{\beta }}, {\widetilde{\gamma }})\), which will transform the optimal signal-to-noise ratio into the sum of the signal-to-noise ratios of the “transformed” three interferometric combinations. The expressions of the three eigenvalues \(\{\mu _i\}^3_{i=1}\) (which are real) of the noise correlation matrix \(\mathbf{C}\) can easily be found by using the algebraic manipulator Mathematica (Wolfram 2014), and they are equal to

$$\begin{aligned} \mu _1 = \mu _2 = S_{\alpha } - S_{\alpha \beta } , \qquad \mu _3 = S_{\alpha } + 2 S_{\alpha \beta } . \end{aligned}$$
(114)

Note that two of the three real eigenvalues, (\(\mu _1\), \(\mu _2\)), are equal. This implies that the eigenvector associated to \(\mu _3\) is orthogonal to the two-dimensional space generated by the eigenvalue \(\mu _1\), while any chosen pair of eigenvectors corresponding to \(\mu _1\) will not necessarily be orthogonal. This inconvenience can be avoided by choosing an arbitrary set of vectors in this two-dimensional space, and by ortho-normalizing them. After some simple algebra, we have derived the following three ortho-normalized eigenvectors:

$$\begin{aligned} \mathbf{v_1} = \frac{1}{\sqrt{2}} \, (-1, 0, 1) \qquad \mathbf{v_2} = \frac{1}{\sqrt{6}} \, (1, -2, 1) \qquad \mathbf{v_3} = \frac{1}{\sqrt{3}} \, (1, 1, 1) . \end{aligned}$$
(115)

Equation (115) implies the following three linear combinations of the generators \(({\widetilde{\alpha }}, {\widetilde{\beta }}, {\widetilde{\gamma }})\):

$$\begin{aligned} A \equiv \frac{1}{\sqrt{2}} \left( {\widetilde{\gamma }} - {\widetilde{\alpha }}\right) , \ \qquad E \equiv \frac{1}{\sqrt{6}} \left( {\widetilde{\alpha }} - 2 {\widetilde{\beta }} + {\widetilde{\gamma }}\right) , \ \qquad T \equiv \frac{1}{\sqrt{3}} \left( {\widetilde{\alpha }} + {\widetilde{\beta }} + {\widetilde{\gamma }}\right) , \end{aligned}$$
(116)

where A, E, and T are italicized to indicate that these are “orthogonal modes”. Although the expressions for the modes A and E depend on our particular choice for the two eigenvectors (\(\mathbf{v_1}, \mathbf{v_2}\)), it is clear from our earlier considerations that the value of the optimal signal-to-noise ratio is unaffected by such a choice. From Eq. (116) it is also easy to verify that the noise correlation matrix of these three combinations is diagonal, and that its non-zero elements are indeed equal to the eigenvalues given in Eq. (114).

To calculate the sensitivity corresponding to the expression of the optimal signal-to-noise ratio averaged over signals randomly distributed on the sky and polarization states, we have proceeded similarly to what was done in Armstrong et al. (1999), Estabrook et al. (2000), Tinto et al. (2002a), and described in more detail in Prince et al. (2002). We assume an equal-arm LISA (\(L = 8.33 \mathrm {\ s}\)), and take the one-sided spectra of proof mass and aggregate optical-path-noises (on a single link), expressed as fractional frequency fluctuation spectra, to be respectively (see Amaro-Seoane et al. 2017)

$$\begin{aligned} S^{\, \mathrm {proof\, mass}}_y (f)& = 2.5 \times {10^{-48}} [f/1 \mathrm {\ Hz}]^{-2} \ \left[ 1 + \left( \frac{0.4 \ \mathrm {mHz}}{f} \right) ^2 \right] \ \\&\quad \times \left[ 1 + \left( \frac{f}{8 \ \mathrm {mHz}} \right) ^4 \right] \mathrm {\ Hz}^{-1} , \end{aligned}$$
(117)
$$\begin{aligned} S^{\, \mathrm {optical \, path}}_y (f)& = 4.4 \times 10^{-38} [f/1 \mathrm {\ Hz}]^2 \ \left[ 1 + \left( \frac{2 \ \mathrm {mHz}}{f} \right) ^4 \right] \mathrm {\ Hz}^{-1} . \end{aligned}$$
(118)

We also assume that aggregate optical path noise has the same transfer function as shot noise. To derive the expression of the LISA’s optimal sensitivity to a signal from a given direction and with a polarization state, let us rewrite the optimal SNR in terms of the combinations A, E, and T

$$\begin{aligned} \mathrm {SNR}^2_{\mathrm {opt}}& = \int _{f_\mathrm {l}}^{f_\mathrm {u}} |{{\widetilde{h}}} (f)|^2 \left[ \frac{|{{\widetilde{R}}}_A(\theta , \phi , \psi , \varepsilon , f, L)|^2}{S_A(f, L)} + \frac{|{{\widetilde{R}}}_E(\theta , \phi , \psi , \varepsilon , f, L)|^2}{S_E(f, L)} \right. \\&\quad\left. + \frac{|{{\widetilde{R}}}_T(\theta , \phi , \psi , \varepsilon , f, L)|^2}{S_T(f, L)} \right] , \end{aligned}$$
(119)

where \({{\widetilde{R}}}_A(\theta , \phi , \psi , \varepsilon , f, L)\), \({{\widetilde{R}}}_E(\theta , \phi , \psi , \varepsilon , f, L)\), \({{\widetilde{R}}}_T(\theta , \phi , \psi , \varepsilon , f, L)\), are the transfer functions of a gravitational wave signal into the combinations A, E, and T. Note we have made explicit the dependence of both the signal and the noise transfer functions on the array nominal arm length. In Eq. (119) \({{\widetilde{h}}} (f)\) is the amplitude of the gravitational wave at the Fourier frequency f, (\(\varepsilon , \psi \)) describe its polarization state, (\(\theta , \phi \)) are the Euler angles associated with the direction of propagation of the gravitational wave signal with respect to a selected coordinate system,Footnote 5 and \(S_A(f, L)\), \(S_E(f, L)\), \(S_T(f, L)\) are the spectral densities of the noises in the A, E, T combinations respectively.

From Eq. (119) it is easy to derive the following expression for the LISA optimal sensitivity, \(S_{h,\theta ,\phi ,\psi ,\varepsilon }^{opt} (\theta , \phi , \psi , \varepsilon , f, L)\), to a gravitational wave signal characterized by a specific direction of propagation (\(\theta , \phi \)) and polarization state (\(\epsilon , \psi \))

$$\begin{aligned} S_{h,\theta ,\phi ,\psi ,\varGamma }^{opt} (\theta , \phi , \psi , \varepsilon , f, L) \equiv \frac{S_A \, S_E \, S_T}{S_E \, S_T |{{\widetilde{R}}}_A|^2 + S_A \, S_T |{{\widetilde{R}}}_E|^2 + S_A \, S_E |{{\widetilde{R}}}_T|^2 } , \end{aligned}$$
(120)

where we have omitted the parameters and Fourier frequency dependence in the terms appearing on the right-hand-side of Eq. (120). In the case of an ensemble of sources randomly distributed over the celestial sphere emitting GW signals of random polarization states, the expression for the optimal sensitivity [Eq. (120)] assumes the following form

$$\begin{aligned} S_{h}^{\mathrm {opt}} (f) \equiv \frac{S_A (f,L) \, S_E(f,L) \, S_T(f,L)}{S_E(f,L) \, S_T(f,L) {{{\mathscr {R}}}}_A (f,L) + S_A \, S_T {{{\mathscr {R}}}}_E (f,L) + S_A \, S_E {{{\mathscr {R}}}}_E (f,L)}, \end{aligned}$$
(121)

where we have denoted with \({{{\mathscr {R}}}}_I (f,L) \equiv \langle |{{\widetilde{R}}}_A|^2 \rangle _{\theta ,\phi ,\psi ,\varepsilon } , \, I = A, E, T\), the average values, over signal direction and polarization states, of the module-squared transfer functions of the GW signals in the combinations A, E, and T respectively.

To quantify the sensitivity improvement of the averaged optimal data combination over that of a single LISA Michelson interferometer, we have analytically-numerically integrated the average of the module-squared transfer functions of the GW signals in the combinations A, E, and T, and then used the expression of the optimal averaged sensitivity derived in Eq. (121). From the expressions of \((\alpha , \beta , \gamma )\) and A, E, and T [Eqs. (61116)] one can then obtain the transfer functions of the Fourier components of an arbitrary gravitational wave signal and of the noises into the A, E, and T optimal combinations. The derivation of the averaged transfer functions was performed by first integrating analytically over the polarization parameters (\(\varepsilon , \psi \)) and then numerically over the Euler angles (\(\theta , \phi \)). To gain processing speed, the numerical integration algorithm relied on a Gauss–Legendre routine (Press et al. 2007). The integration was done for 10,000 Fourier frequencies in the \(\sim 10^{-4} \mathrm {\ Hz}\) to \(\sim 1 \mathrm {\ Hz}\) LISA band after deriving the analytic expressions of the Fourier transforms of the gravitational-wave response of \((\alpha , \beta , \gamma )\) from the formulas in Armstrong et al. (1999), Tinto et al. (2002a). From the Fourier transforms of the \((\alpha , \beta , \gamma )\) responses at each frequency, we construct the Fourier transforms of (AET). We then square and average to compute the mean-squared responses of (AET) at that frequency.

The noise spectra of (AET) are determined from the raw spectra of proof-mass and optical-path noises, and the transfer functions of these noises to (AET). Using the transfer functions given in Estabrook et al. (2000), the resulting spectra are equal to

$$\begin{aligned} S_{A}(f)& = S_{E}(f) = 16 \sin ^2(\pi f L) \, [3 + 2 \cos (2 \pi f L) + \cos (4 \pi f L)] \, S^{\, \mathrm {proof\ mass}}_y (f) \\&\quad + 8 \sin ^2(\pi f L) \, [2 + \cos (2 \pi f L)] \, S^{\, \mathrm {optical\, path}}_y (f) , \end{aligned}$$
(122)
$$\begin{aligned} S_{T}(f)& = 8 \sin ^2 (3 \pi f L) \, S^{\,\mathrm {proof\, mass}}_y + 2 [1 + 2 \cos (2 \pi f L)]^2 \, S^{\, \mathrm {optical\, path}}_y (f) . \end{aligned}$$
(123)

In Fig. 7 we show the sensitivity curve for the LISA equal-arm Michelson response as a function of the Fourier frequency, and the sensitivity curve from the optimum weighting of the data. Note that at frequencies where the LISA Michelson combination has best sensitivity, the improvement in signal-to-noise ratio provided by the optimal observable is slightly larger than \(\sqrt{2}\).

Fig. 7
figure 7

The LISA unequal-arm Michelson X sensitivity curve and the sensitivity curve for the optimal combination of the data, both as a function of Fourier frequency. LISA is assumed to have a nominal arm length \(L = 8.33\,\hbox {s}\)

Fig. 8
figure 8

The ratio \(\rho (f)\) between the sensitivity of the unequal-arm Michelson combination X and the optimal sensitivity, as a function of the Fourier frequency f. The sensitivity gain in the low-frequency band is equal to \(\sqrt{2}\), while it can get larger than 2 at selected frequencies in the high-frequency region of the accessible band. The proof mass and optical path noise spectra are those defined by the project (Amaro-Seoane et al. 2017). See the main body of the paper for a quantitative discussion of this point

Fig. 9
figure 9

The sensitivities of the three combinations (AET) and the resulting optimal sensitivity as functions of the Fourier frequency f. The sensitivities of A and E are equal over the entire frequency band. The sensitivity of T is significantly smaller than the other two in the low part of the frequency band, while is comparable to (and at times larger than) the sensitivities of the other two in the high-frequency region. See text for a complete discussion

In Fig. 8 we plot the ratio between the sensitivity of a single Michelson interferometer and the optimal sensitivity. For Fourier frequencies smaller than 1/L, the sensitivity improvement is \(\sqrt{2}\). For frequencies greater than or about equal to 1/L, the sensitivity improvement is larger and varies with the frequency, showing an average value of about \(\sqrt{3}\). In particular, for bands of frequencies centered on integer multiples of 1/L, the sensitivity of the combination T contributes strongly and the aggregate sensitivity in these bands can be greater than 2.

To better understand the contribution from the three different combinations to the optimal combination of the three generators, in Fig. 9 we plot the sensitivities of (AET) as well as the optimal sensitivity. The sensitivities of the three modes are plotted versus frequency. For the equal-arm case computed here, the sensitivity of A and E are equal across the band. In the long wavelength region of the band, modes A and E have SNRs much greater than mode T, where its contribution to the total sensitivity is negligible. At higher frequencies, however, the T combination has sensitivity greater than or comparable to the other modes and can dominate the sensitivity improvement at selected frequencies. Some of these results have also been obtained in Rajesh Nayak et al. (2003b).

6.2 Optimization of SNR for binaries with known direction but with unknown orientation of the orbital plane

Binaries are important sources for LISA and therefore the analysis of such sources is of major importance. One such class is of massive or super-massive binaries whose individual masses could range from \(10^3 M_{\odot }\) to \(10^8 M_{\odot }\) and which could be up to a few Gpc away. Another class of interest are known binaries within our own galaxy whose individual masses are of the order of a solar mass but are just at a distance of a few kpc or less. Here the focus will be on this latter class of binaries. It is assumed that the direction of the source is known, which is so for known binaries in our galaxy. However, even for such binaries, the inclination angle of the plane of the orbit of the binary is either poorly estimated or unknown. The optimization problem is now posed differently: the SNR is optimized after averaging over the orientations of the orbital planes of the binary systems, so the results obtained are optimal on the average, that is, the source is tracked with an observable which is optimal on the average (Rajesh Nayak et al. 2003b). For computing the average, a uniform distribution for the direction of the orbital angular momentum of the binary is assumed.

When the binary masses are of the order of a solar mass and the signal typically has a frequency of a few mHz, the GW frequency of the binary may be taken to be constant over the period of observation, which is typically taken to be of the order of an year. A complete calculation of the signal matrix and the optimization procedure of SNR is given in Rajesh Nayak et al. (2003a). Here we briefly mention the main points and the final results.

We consider a point source fixed in the Solar System Barycentric reference frame. But as the LISA constellation moves along its heliocentric orbit, the apparent direction of the source in the LISA frame changes as a function of time. Since there are two frames involved, we label the direction angles explicitly with corresponding subscripts. The direction angles of the source in barycentric frame are labelled as \((\theta _\mathrm {B}, \phi _\mathrm {B})\) while those in the LISA frame as \((\theta _\mathrm {L}. \phi _\mathrm {L})\). The corresponding LISA reference frame is denoted by \((x_\mathrm {L}, y_\mathrm {L}, z_\mathrm {L})\) whose orientation changes with time with respect to the barycentric frame. The LISA reference frame \((x_\mathrm {L}, y_\mathrm {L}, z_\mathrm {L})\) has been defined in Rajesh Nayak et al. (2003a) as follows: The origin lies at the center of the LISA triangle and the plane of LISA coincides with the \((x_\mathrm {L}, y_\mathrm {L})\) plane with spacecraft 2 lying on the \(x_\mathrm {L}\) axis. Figure 10 displays this apparent motion for a source lying in the ecliptic plane, that is with \(\theta _\mathrm {B} = 90^{\circ }\) and \(\phi _\mathrm {B} = 0^{\circ }\). The source in the LISA reference frame describes a figure of 8. Optimizing the SNR amounts to tracking the source with an optimal observable as the source apparently moves in the LISA reference frame.

Fig. 10
figure 10

Apparent position of the source in the sky as seen from LISA frame for \((\theta _\mathrm {B}= 90^{\circ }, \phi _\mathrm {B}=0^{\circ })\). The track of the source for a period of 1 year is shown on the unit sphere in the LISA reference frame

Since an average has been taken over the orientation of the orbital plane of the binary, the signal matrix is now different from the matrix \(\mathbf{A}\) obtained in Sect. 6.1.

The mutually orthogonal data combinations A, E, T are convenient in carrying out the computations because in this case as well, they simultaneously diagonalize the signal and the noise covariance matrix. For convenience we label them with index \(I = 1, 2, 3\) respectively. The basic GW wave strains of the binary in the Fourier domain are given by,

$$\begin{aligned} h_{+}(\varOmega )& = H \left[ \frac{1+\cos ^{2}\varepsilon }{2}\cos 2\psi -i\,\cos \varepsilon \,\sin 2\psi \right] , \\ h_{\times }(\varOmega )& = H \left[ -\frac{1+\cos ^{2}\varepsilon }{2}\sin 2\psi -i\,\cos \varepsilon \,\cos 2\psi \right] . \end{aligned}$$
(124)

where H is an overall dimensionless amplitude depending on the masses, the distance and the GW frequency of the binary system and the angles \((\varepsilon , \psi )\) describe the orientation of the binary orbit (\(\varepsilon , \psi \) could specify the direction the orbital angular momentum vector). The GW response for a generator is,

$$\begin{aligned} h^{(I)}(\varOmega )= F_{+}^{(I)}(\varOmega )h_{+}(\varOmega ) + F_{\times }^{(I)}(\varOmega )h_{\times }(\varOmega ), \end{aligned}$$
(125)

where the \(F_{+}^{(I)}\) and \(F_{\times }^{(I)}\) are the antenna pattern functions corresponding to the generator I where I is any of the generators AET. From these GW strains we define the signal covariance matrix,

$$\begin{aligned} h^{(I)}_{(J)}& = h^{(I)} h_{(J)}^{*}, \\& = \left( F_{+}^{(I)}h_{+}+F_{\times }^{(I)}h_{\times }\right) \left( F_{+ (J)}h_{+}+ F_{\times (J)}h_{\times }\right) ^{*}. \end{aligned}$$
(126)

In general we may not have any knowledge of the orientation of the orbital plane of the binary source. We therefore average over the direction of the orbital angular momentum of the binary uniformly distributed over the sphere. The orientation of the binary (its orbital angular momentum vector) has been described in terms of the angles \(\varepsilon \) and \(\psi \) in Eq. (124). We carry out the averaging of \(h^{(I)}_{(J)}\) over \((\varepsilon , \psi )\) which results in an overall factor of 2/5. The averaged matrix we denote by \({\mathscr {H}}^{(I)}_{(J)}\). In the Fourier domain it is given by,

$$\begin{aligned} {\mathscr {H}}^{(I)}_{(J)}(\varOmega ) = H^2 \left( \frac{2}{5} \right) \left( F_{+}^{(I)} F_{+ (J)}^{*} + F_{\times }^{(I)} F_{\times (J)}^{*}\right) (\varOmega ) . \end{aligned}$$
(127)

The signal matrix so averaged has the following properties:

  • \({\mathscr {H}}\) is the sum of outer products of two vectors: \(F_{+}^{(I)}\) with its complex conjugate and \(F_{\times }^{(I)}\) with its complex conjugate. Thus the natural basis for expressing \({\mathscr {H}}\) consists of the two vectors \(F_{+}^{(I)}\) and \(F_{\times }^{(I)}\).

  • \({\mathscr {H}}\) is of rank two, everywhere except on the \(\theta =\frac{\pi }{2}\) plane where it is one when \(F_{\times }^{(I)}\) goes identically to zero. In Sect. 6.1 since optimization is performed first before averaging, the signal matrix is constituted from a single vector and thus has rank one.

For a generic data combination with \(\alpha _{(I)}\) as polynomial coefficients in the delay operators, the SNR is given by:

$$\begin{aligned} \mathrm {SNR}^{2}=\frac{\alpha _{(I)}\alpha ^{(J) *}{\mathscr {H}}^{(I)}_{(J)}}{\alpha _{(I)}\alpha ^{(J) *}N^{(I)}_{(J)}} , \end{aligned}$$
(128)

where \(N^{(I)}_{(J)}\) is the noise matrix obtained by taking the outer product of the noise vectors for each I (see Rajesh Nayak et al. 2003b for details). Also since we are dealing with essentially monochromatic sources, the quantities \(\alpha _{(I)}\) just reduce to complex numbers. It is shown in Rajesh Nayak et al. (2003a) that the optimization problem reduces to an eigenvalue problem with the eigenvalues being the squares of the SNRs. There are two eigen-vectors which are labeled as \(\mathbf {v}_{+, \times }\) belonging to two non-zero eigenvalues. The two SNRs are labeled as \(\mathrm {SNR}_+\) and \(\mathrm {SNR}_{\times }\), corresponding to the two orthogonal (thus statistically independent) eigenvectors \(\mathbf {v}_{+, \times }\). As was done in Sect. 6.1, the two SNRs can be squared and added to yield a network SNR, which is defined through the equation

$$\begin{aligned} \mathrm {SNR}_\mathrm {network}^2 = \mathrm {SNR}_+^2 + \mathrm {SNR}_\times ^2 . \end{aligned}$$
(129)

The corresponding observable is called the network observable. The third eigenvalue is zero and the corresponding eigenvector orthogonal to \(\mathbf {v}_{+}\) and \(\mathbf {v}_{\times }\) gives zero signal. The eigenvectors and the SNRs are functions of the apparent source direction parameters \((\theta _\mathrm {L}, \phi _\mathrm {L})\) in the LISA reference frame, which in turn are functions of time. The eigenvectors optimally track the source as it moves in the LISA reference frame. Assuming an observation period of an year, the SNRs may be integrated over this period of time.

In the low frequency regime, we can obtain analytical expressions for the optimal SNRs. A large fraction of the sources for LISA fall into this category, for example, massive/supermassive blackhole binaries, several galactic and extragalactic binaries which contribute to the ‘confusion noise’. Here we do not give the detailed calculations but refer the reader to Rajesh Nayak et al. (2003a).

Integration over a period of \(T = T_{\odot }\) leads to the following results:

$$\begin{aligned} \mathrm {SNR}_{+}(\theta _B, \phi _B)& = \mathrm {SNR}_0 \, g_{+} (\cos \theta _B) , \\ \mathrm {SNR}_{\times } (\theta _B, \phi _B)& = \mathrm {SNR}_0 \, g_{\times } (\cos \theta _B), \end{aligned}$$
(130)

where,

$$\begin{aligned} \mathrm {SNR}_0 = \frac{3}{\sqrt{5}} \frac{H}{n_{A}} (\varOmega L)^2 \sqrt{T_{\odot }} , \end{aligned}$$
(131)

and

$$\begin{aligned} g_{+}^2 (x)& = \frac{1}{T_{\odot }} \int _0^{T_{\odot }} \left( \frac{1 + \cos ^2 \theta _L}{2} \right) ^2 dt \\& = \frac{1}{4} \left( 1 + \frac{x^2}{4}\right) ^2 \\&\quad + \frac{3}{16} (1 - x^2) \left( 1 + \frac{3}{4} x^2 + \frac{9}{32} \left( 1 - x^2 \right) \right) , \\ g_{\times }^2 (x)& = \frac{1}{T_{\odot }} \int _0^{T_{\odot }} \cos ^2 \theta _L dt \\& = \frac{1}{4} \left[ x^2 + \frac{3}{2} \left( 1 - x^2 \right) \right] . \end{aligned}$$
(132)

The quantity \(n_{A}\) is the PSD of the noise in the data combination A. We take the integration time \(T_{\odot }\) to be of the order of a year. In order to carry out the integrals one needs the relation between \(\theta _L\) and \((\theta _B, \phi _B)\). For the sake of completeness, we give the expressions for both \(\theta _L\) and \(\phi _L\) in terms of \(\theta _B\) and \(\phi _B\):

$$\begin{aligned} \cos \theta _{L}& = \frac{1}{2}\cos \theta _{B} -\frac{\sqrt{3}}{2}\sin \theta _{B}\sin \left( \phi _{B}-\omega t \right) , \\ \tan \phi _L& = \frac{w^{(2)}}{w^{(1)}} , \end{aligned}$$
(133)

where,

$$\begin{aligned} w^{(1)}& = \cos \omega t \sin \theta _B \cos (\phi _B - \omega t) - \frac{1}{2}\sin \omega t \sin \theta _B \sin (\phi _B - \omega t) \, \\&\quad -\, \frac{\sqrt{3}}{2} \, \sin \omega t \cos \theta _B , \\ w^{(2)}& = \sin \omega t \sin \theta _B \cos (\phi _B - \omega t) + \frac{1}{2}\cos \omega t \sin \theta _B \sin (\phi _B - \omega t) \\&\quad +\, \frac{\sqrt{3}}{2} \, \cos \omega t \cos \theta _B . \end{aligned}$$
(134)

Here, \(\omega = 2 \pi /T_{\odot }\) (to simplify the calculation, one could put \(\phi _B = 0\)). This relation is used in the above integrals and a simple calculation leads to the expressions for \(g_{+} (x)\) and \(g_{\times } (x)\) where \(x = \cos \theta _B\).

We have purposely not ‘simplified’ the formulae in powers of \(x^2\) because in this form it is easy to see the limits \(x = \pm 1, 0\) corresponding to \(\theta _B = 0, \pi , \pi /2\) respectively. In fact, since only \(x^2\) occurs in the expressions of \(g_{+, \times }\), there is symmetry about the ecliptic plane. The network SNR is just the root mean square of the two SNRs:

$$\begin{aligned} \mathrm {SNR}_{\mathrm {network}} (\theta _B, \phi _B) = \mathrm {SNR}_0 \, g_{\mathrm {network}} (\cos \theta _B), \end{aligned}$$
(135)

where, \(g_{\mathrm {network}}^2 (x) = g_{+}^2 (x) + g_{\times }^2 (x)\). The factors g are of the order of unity and the \(\hbox {SNR}_0\) gives essentially the integrated SNR of a GW source. The maximum integrated SNR is obtained for sources lying in the ecliptic plane \((\theta _B = 90^{\circ })\). This can be readily explained from Fig. 10 where the trajectory of such a source is plotted. One observes that a large fraction of the orbit in the LISA frame is away from the LISA plane \((\theta _L = 90^{\circ })\) where most of the contribution to the SNR comes from.

7 Experimental aspects of TDI

It is clear that the suppression of the laser phase fluctuations by more than seven orders of magnitude with the use of TDI is a very challenging experimental task. It requires developing and building subsystems capable of unprecedented accuracy and precision levels, and test their end-to-end performance in a laboratory environment that naturally precludes the availability of \(2.5 \times 10^6\) km delay lines! In what follows we will address some aspects related to the experimental implementation of TDI, and derive the performance specifications for the subsystems involved. We will not describe, however, any of the experimental efforts related to the verification of TDI in a laboratory environment. For that, we refer the interested reader to de Vine et al. (2010), Miller (2010), Spero et al. (2011) and Mitryk et al. (2012).

From simple physical grounds, it is easy to see that a successful implementation of TDI requires:

  1. 1.

    accurate knowledge of the time shifts, \(L_{i} (t), L_{i} '(t) \ i = 1, 2, 3\), to be applied to the heterodyne measurements \(s_i(t), s_{i'}(t), \varepsilon _i(t), \varepsilon _{i'}(t), \tau _i(t), \tau _{i'}(t) \ i = 1, 2, 3\);

  2. 2.

    accurate synchronization among the three clocksFootnote 6 on board the three spacecraft as these are used for time-stamping the recorded heterodyne phase measurements;

  3. 3.

    sampling time stability (i.e., clock stability) for successfully suppressing the laser noise to the desired level;

  4. 4.

    an accurate reconstruction algorithm of the phase measurements corresponding to the required time delays as these in general will not be equal to integer multiples of the sampling time;

  5. 5.

    a phase meter capable of a very large dynamic range in order to suppress the laser noise to the required level while still preserving the phase fluctuations induced by a gravitational-wave signal in the TDI combinations;

  6. 6.

    a precise interferometric design architecture and an accurate procedure for calibrating the phase fluctuations of the onboard clocks out of the down-converted one-way heterodyne Doppler measurements.

In addition, it has recently been pointed out (Bayle et al. 2019) that the low-pass filters applied to the one-way measurements (to comply with the stringent data rates limitations faced by LISA) do not commute with the delay operators of time-changing delays. Although this fact prevents the exact cancellation of the filtered laser noises in the second-generation TDI combinations, it is still possible to achieve adequate suppression of the laser noise by applying TDI to the filtered one-way data (Bayle et al. 2019).

In the following subsections, we will quantitatively address the issues listed above, and provide the reader with a related list of references where more details can be found.

7.1 Time-delays accuracies

The TDI combinations described in the previous sections (whether of the first- or second-generation) rely on the assumption of knowing the time-delays with infinite accuracy to exactly cancel the laser noise. Since the six delays will in fact be known only within the accuracies \(\delta L_i , \delta L_{i}' \ i=1, 2, 3\), the cancellation of the laser frequency fluctuations in, for instance, the combinations (\(\alpha , \beta , \gamma , \zeta \)) will no longer be exact. To estimate the magnitude of the laser fluctuations remaining in these data sets, let us define \({\hat{L}}_i , \hat{L_{i}'} \ i=1, 2, 3\) to be the estimated time-delays. They are related to the true delays \(L_i , L_{i}' \ i=1, 2, 3\), and the accuracies \(\delta L_i \ , \delta L_{i}' \ i=1, 2, 3\) through the following expressions

$$\begin{aligned} {\hat{L}}_i = L_i + \delta L_i , \qquad i=1, 2, 3 , \end{aligned}$$
(136)

and similarly for the primed delays. In what follows we will limit our derivation of the time-delays accuracies to only the first-generation TDI combinations and treat the three common delays \(L_i , \ i=1, 2, 3\) as constants equal to 8.33 light-seconds. We will also assume to know with infinite accuracies and precisions all the remaining physical quantities (listed at the beginning of Sect. 7) that are needed to successfully synthesize the TDI generators.

If we now substitute Eq. (136) into the expression for the TDI combination \(\alpha \), for instance, [Eq. (61)] and expand it to first order in \(\delta L_i\), it is easy to derive the following approximate expression for \({{\hat{\alpha }}} (t)\), which now will show a non-zero contribution from the laser noises

$$\begin{aligned} {{\hat{\alpha }}} (t) \simeq \alpha (t) + [{{\dot{\phi }}}_{2,12} - {{\dot{\phi }}}_{3,13}] \ \delta L_1 + [{{\dot{\phi }}}_{3,2} - {{\dot{\phi }}}_{1,123}] \ \delta L_2 + [{{\dot{\phi }}}_{1,123} - {{\dot{\phi }}}_{2,3}] \ \delta L_3 , \end{aligned}$$
(137)

where the symbol \(\dot{}\) denotes time derivative. Time-delay interferometry can be considered effective if the magnitude of the remaining fluctuations from the lasers are much smaller than the fluctuations due to the secondary (proof mass and optical path) noises entering \(\alpha (t)\). This requirement implies a limit in the accuracies of the measured delays.

Let us assume the laser phase fluctuations to be uncorrelated to each other, their one-sided power spectral densities to be equal, the three arm lengths to differ by a few percent, and the three arm length accuracies also to be equal. By requiring the magnitude of the remaining laser noises to be smaller than the secondary noise sources, it is straightforward to derive, from Eq. (137) and the expressions for the noise spectrum of the \(\alpha \) TDI combination given in Estabrook et al. (2000), the following constraint on the common arm length accuracy \(| \delta L_{\alpha }|\)

$$\begin{aligned}&|\delta L_{\alpha }| \ll \varDelta L_{\alpha } (f) \\&\quad \equiv \frac{1}{\sqrt{32} \pi f} \sqrt{\frac{[8 \sin ^2(3 \pi f L) + 16 \sin ^2(\pi f L) ] S^{\mathrm {\,proof \ mass}} + S^{\mathrm {\,optical \ path}} }{S_{\phi }}} , \end{aligned}$$
(138)

with similar inequalities also holding for \(\beta \) and \(\gamma \). Here \(S_{\phi }\), \(S^{\mathrm {\,proof \ mass}}\), \(S^{\mathrm {\,optical \ path}}\) are the one-sided power spectral densities of the relative frequency fluctuations of a stabilized laser, a single proof mass, and a single-link optical path respectively. If we take them to be equal to the following functions of the Fourier frequency f (Tinto and Armstrong 1999; Amaro-Seoane et al. 2017)

$$\begin{aligned} S_{\phi } (f)& = 10^{-28} \ f^{-2/3} + 6.3 \ \times \ 10^{-37} \ f^{-3.4} \mathrm {\ Hz}^{-1} , \end{aligned}$$
(139)
$$\begin{aligned} S^{\, \mathrm {proof\, mass}}_y (f)& = 2.5 \times {10^{-48}} [f/1 \mathrm {\ Hz}]^{-2} \ \left[ 1 + \left( \frac{0.4 \ \mathrm {mHz}}{f} \right) ^2 \right] \ \\&\quad \times \left[ 1 + \left( \frac{f}{8 \ \mathrm {mHz}} \right) ^4 \right] \mathrm {\ Hz}^{-1} , \end{aligned}$$
(140)
$$\begin{aligned} S^{\, \mathrm {optical \, path}}_y (f)& = 4.4 \times 10^{-38} [f/1 \mathrm {\ Hz}]^2 \ \left[ 1 + \left( \frac{2 \ \mathrm {mHz}}{f} \right) ^4 \right] \mathrm {\ Hz}^{-1} . \end{aligned}$$
(141)

(where f is in Hz), we find that \(\varDelta L_{\alpha } (f)\) [the right-hand side of the inequality given by Eq. (138)] reaches its minimum of about 62.5 m at the Fourier frequency \(f_{\min } = 3.4 \times 10^{-3}\mathrm {\ Hz}\), over the assumed \((10^{-4}, 1)\mathrm {\ Hz}\) LISA band. This implies that, if the arm length knowledge \(| \delta L_{\alpha } |\) can be made much smaller than 62.5 m, the magnitude of the residual laser noise affecting the \(\alpha \) combination can be regarded as negligible over the entire frequency band. This reflects the fact that \(\varDelta L_{\alpha } (f)\) is a “U-shaped” function of the Fourier frequency (see Fig. 11 in the following subsection, where we plot both the delays as well as the clocks timing-synchronization accuracies).

A perturbation analysis similar to the one described above can be performed for \(\zeta \), resulting into the following inequality for the required delay accuracy, \(| \delta L_{\zeta } |\)

$$\begin{aligned} |\delta L_{\zeta }| \ll \ \varDelta L_{\zeta } (f) \equiv \frac{1}{2 \pi f} \ \sqrt{\frac{4 \ \sin ^2(\pi f L) \ S^{\, \mathrm {proof\, mass}} + S^{\, \mathrm {optical\, path}} }{S_{\phi } }} . \end{aligned}$$
(142)

The inequality in Eq. (142) implies a minimum of the function \(\varDelta L_{\zeta } (f)\) equal to about 155 m at the Fourier frequency \(f_{\min } = 2.8 \times 10^{-3}\mathrm {\ Hz}\). Here again the right-hand-side of Eq. (142) is a “U-shaped” function of the Fourier frequency.

Arm length accuracies at the centimeters level have already been demonstrated in the laboratory (Esteban et al. 2011; Sutton et al. 2010; Wang et al. 2014; Heinzel et al. 2011), making us confident that the required level of time-delays accuracy will be available.

In relation to the accuracies derived above, it is interesting to calculate the time scales during which the arm lengths will change by an amount equal to the accuracies themselves. This identifies the minimum time required before updating the arm length values in the TDI combinations.

It has been calculated by Folkner et al. (1997) that the relative longitudinal speeds between the three pairs of spacecraft, during approximately the first year of the LISA mission, can be written in the following approximate form

$$\begin{aligned} V_{i,j} (t) = V^{(0)}_{i,j} \ \sin \left( {{2 \pi t} \over T_{i,j}}\right) \qquad (i,j) = (1,2) \ ; \ (1,3) \ ; \ (2,3) , \end{aligned}$$
(143)

where we have denoted with (1, 2), (1, 3), (2, 3) the three possible spacecraft pairs, \(V^{(0)}_{i,j}\) is a constant velocity, and \(T_{i,j}\) is the period for the pair (ij). In Folkner et al. (1997), it has also been shown that the LISA trajectory can be selected in such a way that two of the three arms’ rates of change are essentially equal during the first year of the mission. Following Folkner et al. (1997), but scaling the velocities to values corresponding to the current LISA configuration of \(2.5 \ \times 10^6\) km, we will assume \(V^{(0)}_{1,2} = V^{(0)}_{1,3} \ne V^{(0)}_{2,3}\), with \(V^{(0)}_{1,2} = 0.25 \mathrm {\ m/s}\), \(V^{(0)}_{2,3} = 3.1\mathrm {\ m/s}\), \(T_{1,2} = T_{1,3} \approx 4\) months, and \(T_{2,3} \approx 1\) year. From Eq. (143) it is easy to derive the variation of each arm length, for example \(\Updelta L_3 (t)\), as a function of time t and the time scale \(\delta t\) during which it takes place

$$\begin{aligned} \Updelta L_3 (t) = V^{(0)}_{1,2} \ \sin \left( {{2 \pi t} \over T_{1,2}}\right) \delta t . \end{aligned}$$
(144)

Equation (144) implies that a variation in arm length \(\Updelta L_3 \approx 10\mathrm {\ m}\) can take place during different time scales, depending on when during the mission this change takes place. For instance, if \(t \ll T_{1,2}\) we find that the arm length \(L_3\) changes by more than its accuracy (\(\approx 10\mathrm {\ m}\)) after a time \(\delta t = 9.2 \times 10^3\mathrm {\ s}\). If however \(t \simeq T_{1,2}/4\), the arm length will change by the same amount after only \( \delta t \simeq 40\mathrm {\ s}\) instead. As this value is comparable to the nominal one-way-light-time (8.33 sec), one might argue that the measured time-delay might not represent well enough the delay that needs to be applied in the TDI combinations at that particular time.

One way to address this problem is to treat the delays in the TDI combinations as parameters to be determined by a non-linear least-squares procedure. This is because the minimum value of the minimizer is achieved at the correct delays for which the laser noises exactly cancel in the TDI combinations. Such a technique, which was named Time-Delay Interferometric Ranging (TDIR) (Tinto et al. 2005), requires a starting point in the time-delays space in order to implement the minimization, and it works quite effectively jointly with the ranging data available onboard. At the time of this writing the LISA project has decided to apply TDIR to the mission data.

7.2 Clocks synchronization

The effectiveness of the TDI data combinations requires the clocks on board the three spacecraft to be synchronized. In what follows we will identify the minimum level of off-synchronization among the clocks that can be tolerated. To proceed with our analysis we will treat one of the three clocks (say the clock on board spacecraft 1) as the master clock defining the “LISA-time”, while the other two to be synchronized to it.

The relativistic (Sagnac) time-delay effect due to the fact that the LISA trajectory is a combination of two rotations, each with a period of one year, will have to be accounted for in the synchronization procedure and, as has already been discussed earlier, will be accounted for within the second-generation formulation of TDI.

Here, for simplicity, we will analyze an idealized non-rotating constellation in order to get a sense of the required level of clocks synchronization. Let us denote by \( \delta t_2\), \( \delta t_3\), the time accuracies (time-offsets) for the clocks on board spacecraft 2 and 3 respectively. If t is the time on board spacecraft 1, then what is believed to be time t on board spacecraft 2 and 3 is actually equal to the following times

$$\begin{aligned} {{\hat{t}}_2}& = t + \delta t_2 , \end{aligned}$$
(145)
$$\begin{aligned} {{\hat{t}}_3}& = t + \delta t_3 . \end{aligned}$$
(146)

If we now substitute Eqs. (145 and 146) into the TDI combination \(\zeta \), for instance, and expand it to first order in \(\delta t_i , \ i=2, 3\), it is easy to derive the following approximate expression for \({{\hat{\zeta }}} (t)\), which shows the following non-zero contribution from the laser noises

$$\begin{aligned} {{\hat{\zeta }}} (t) \simeq \zeta (t) + [{{\dot{\phi }}}_{1,23} - {{\dot{\phi }}}_{3,12} + {{\dot{\phi }}}_{2,13} - {\dot{\phi }}_{2,13}] \ \delta t_2 + [{{\dot{\phi }}}_{2,13} - {{\dot{\phi }}}_{1,23} + {{\dot{\phi }}}_{3,12} - {\dot{\phi }}_{3,12}] \ \delta t_3 . \end{aligned}$$
(147)

By requiring again the magnitude of the remaining fluctuations from the lasers to be smaller than the fluctuations due to the other (secondary) noise sources affecting \(\zeta (t)\), it is possible to derive an upper limit for the accuracies of the synchronization of the clocks. If we assume again the three laser phase fluctuations to be uncorrelated to each other, their one-sided power spectral densities to be equal, the three arm lengths to differ by a few percent, and the two time-offset’s magnitudes to be equal, by requiring the magnitude of the remaining laser noises to be smaller than the secondary noise sources it is easy to derive the following constraint on the time synchronization accuracy \(| \delta t_{\zeta }|\)

$$\begin{aligned} |\delta t_{\zeta }| \ll \varDelta t_{\zeta } (f) \equiv \ {1 \over {2 \pi f}} \ \sqrt{ {12 \ \sin ^2(\pi f L) \ S^{\mathrm {\,proof \ mass}} + 3 \ S^{\mathrm {\,optical \ path}} } \over {4 \ S_{\phi } } } , \end{aligned}$$
(148)

with \(S_{\phi }\), \(S^{\mathrm {\,proof \ mass}}\), \(S^{\mathrm {\,optical \ path}}\) again as given in Eqs. (139 – 141).

We find that the function \(\varDelta t_{\zeta } (f)\) defined in Eq. (148) reaches its minimum of about 448 nanoseconds at the Fourier frequency \(f_{\min } = 2.8 \times 10^{-3}\mathrm {\ Hz}\). This means that clocks synchronized at a level of accuracy significantly better than 448 nanoseconds will result into a residual laser noise that is much smaller than the secondary noise sources entering into the \(\zeta \) combination.

An analysis similar to the one described above can be performed for the remaining generators (\(\alpha , \beta , \gamma \)). For them we find that the corresponding inequality for the accuracy in the synchronization of the clocks is now equal to

$$\begin{aligned}&|\delta t_{\alpha }| \ll \varDelta t_{\alpha } (f) \\&\quad \equiv \ {1 \over {2 \pi f}} \ \sqrt{ {[4 \sin ^2(3 \pi f L) + 8 \sin ^2(\pi f L)] S^{\mathrm {\,proof \ mass}} + 3 S^{\mathrm {\,optical \ path}} } \over {4 \ S_{\phi }} } , \end{aligned}$$
(149)

with equal expressions holding also for \(\beta \) and \(\gamma \). The function \(\varDelta t_{\alpha } (f)\), given on the right-hand side of Eq. (149), has a minimum equal to 457 nanoseconds at the Fourier frequency \(f_{\min } = 2.9 \times 10^{-3}\mathrm {\ Hz}\). As for the arm length accuracies, also the timing accuracy requirements become less stringent at different frequencies because of their “U-shaped” dependence on the Fourier frequency.

Fig. 11
figure 11

The delays and clocks timing-synchronization functions \(\varDelta L_{\alpha } (f)\), \(\varDelta L_{\zeta } (f)\), \(\varDelta t_{\alpha } (f)\), \(\varDelta t_{\zeta } (f)\), as defined by the expressions appearing on the right-hand-sides of Eqs. (138,  142148149) respectively. The minimum values of these “U-shaped” functions correspond to the minimum accuracies required for the delays and clocks timing-synchronizations so that the laser noises are suppressed to the level identified by the secondary noises. The accuracies implemented by LISA should naturally be better than these values to guarantee a residual laser noise that is much smaller than the level identified by the secondary noises in the TDI combinations

In Fig. 11 we plot both delays and clocks timing-synchronization functions \(\varDelta L_{\alpha } (f)\), \(\varDelta L_{\zeta } (f)\), \(\varDelta t_{\alpha } (f)\), \(\varDelta t_{\zeta } (f)\), as defined by the expressions appearing on the right-hand-sides of Eqs. (138142148, 149) respectively. The minima of these curves identify the upper limits for the delays and clocks timing-synchronization accuracies for which the residual laser noise is smaller than the secondary noises in the TDI combinations. By requiring, for instance, the residual laser noise to be a factor of 10 smaller than the secondary noises over the entire LISA band, we would need to achieve levels of accuracies that are a factor of 10 smaller than those quoted above for the \(\alpha \) and \(\zeta \) combinations. A required clock synchronization accuracy of about 40 ns, for instance, translates into a ranging accuracy of 12 m, which has been experimentally shown to be easily achievable (Esteban et al. 2011; Sutton et al. 2010; Wang et al. 2014; Heinzel et al. 2011).

7.3 Clocks timing jitter

The sampling times of all the measurements needed for synthesizing the TDI combinations will not be constant, due to the intrinsic timing jitters of the onboard measuring system. Among all the subsystems involved in the data measuring process, the onboard clock is expected to be the dominant source of time jitter in the sampled data. Presently existing space qualified clocks can achieve an Allan standard deviation of about \(10^{-13}\) for integration times from 1 to 10,000 s. This timing stability translates into a time jitter of about \(10^{-13}\,\hbox {s}\) over a period of 1 s. A perturbation analysis including the three sampling time jitters due to the three clocks shows that any laser phase fluctuations remaining in the four TDI generators will also be proportional to the sampling time jitters. Since the latter are approximately five orders of magnitude smaller than the arm length and clocks synchronization accuracies derived earlier, we conclude that the magnitude of laser noise residual into the TDI combinations due to the sampling time jitters is negligible.

7.4 Sampling reconstruction algorithm

The derivations of the time-delays and clocks synchronization accuracies highlighted earlier presumed the availability of the phase measurement samples at the required time-delays. Since this condition will not be true in general, as the time-delays used by the TDI combinations will not be equal to integer-multiples of the sampling time, with a sampling rate of, let us say, 10 Hz, the time delays could be off their correct values by a tenth of a second, way more than the 10 nanoseconds time-delays and clocks synchronization accuracies estimated above.

Earlier suggestions (Hellings 2001) for addressing this problem envisioned sampling the data at very-high rates (perhaps of the order of hundreds of MHz), so reducing the additional error to the estimated time-delays to a few nanoseconds. Although in principle such a solution would allow us to suppress the residual laser noise to the required level, it would create an insurmountable problem for transmitting the science data to the ground due to the limited space-to-ground data rates.

An alternate scheme for obtaining the phase measurement points needed by TDI (Tinto et al. 2003) envisioned sampling the phase measurements at the required delayed times. This scheme naturally requires knowledge of the time-delays and synchronization of the clocks at the required accuracy levels during data acquisition. Although such a procedure could be feasible in principle, it would still leave open the possibility of irreversible corruption of the TDI combinations in the eventuality of performance degradation in the ranging and clock synchronization procedures.

Given that the data will need to be sampled at a rate of 10 Hz, an alternative options is to implement an interpolation scheme for reconstructing the required data points from the sampled measurements. An analysis (Tinto et al. 2003) based on the implementation of the truncated (Shannon 1998) formula, however, showed that several months of data were required in order to reconstruct the phase samples at the estimated time-delays with a sufficiently high accuracy. This conclusion implied that several months (at the beginning and end) of the entire data records measured by LISA would be of no use, resulting into a significant mission science degradation.

Although the truncated Shannon formula was proved to be impracticable (Tinto et al. 2003) for reconstructing phase samples at the required time-delays, it was then recognized that (Shaddock et al. 2004) a more efficient and accurate interpolation technique (Laakso et al. 1996) could be adopted. In what follows, we provide a brief account of this data processing technique, which is known as “fractional-delay filtering” (FDF).

To understand how FDF works, let’s write the truncated Shannon formula for the delayed sample, \(y_N(n - D)\), which we want to construct by filtering the sampled data y(n)

$$\begin{aligned} y_N(n - D) = \sum _{j = -N}^{N} y(n + j) \ \mathrm {sinc}(D - j) , \end{aligned}$$
(150)

where, as usual, \(\mathrm {sinc}(x) \equiv \sin (x)/x\), and N is an integer at which the Shannon formula has been truncated to. As pointed out in Shaddock et al. (2004), although the truncated Shannon formula is optimal in the least-squares sense, the sinc-function that appears in it is far from being ideal in reconstructing the transfer function \(e^{2 \pi i f D/f_s}\), where \(f_s\) is the sampling frequency. In fact, over the LISA observational band the sinc-function displays significant ringing, which can only be suppressed by taking N very large (as the error, \(\varepsilon \), decays slowly as 1/N). It was estimated that, in order to achieve an \(\varepsilon < 10^{-8}\), an \(N \ge 10^8\) is needed.

If, however, we give up on the requirement of minimizing the error in the least-squares sense and replace it with a mini-max criterion error applied to the absolute value of the difference between the ideal transfer function (i.e., \(e^{2 \pi i f D/f_s}\)) and a modified sinc-function, we will be able to achieve a rapid convergence while suppressing the ringing effects associated with the sinc-function.

One way to achieve this result is to modify the Shannon formula by multiplying the sinc-function by a window-function, W(j), in the following way

$$\begin{aligned} y_N(n - D) = \sum _{j = -N}^{N} y(n + j) \ W(j) \ \mathrm {sinc}(D - j) , \end{aligned}$$
(151)

where W(j) smoothly decays to zero at \(j = \pm N\). In Shaddock et al. (2004) several windows were tested, and the resulting values of N needed to accurately reconstruct the desired delayed samples were estimated, both on theoretical and numerical grounds. It was found that, with windows belonging to the family of Lagrange polynomials (Shaddock et al. 2004) a delayed sample could be reconstructed by using \(N \simeq 25\) samples while achieving a mini-max error \(\varepsilon < 10^{-12}\) between the ideal transfer function \(e^{2 \pi i f D/f_s}\), and the kernel of the modified truncated Shannon formula.

7.5 Data digitization and bit-accuracy requirement

It has been shown (Tinto et al. 2003) that the maximum of the ratio between the amplitudes of the laser and the secondary phase fluctuations occurs at the lower end of the LISA bandwidth (i.e., 0.1 mHz) and it is equal to about \(10^{10}\). This corresponds to the minimum dynamic range for the phasemeters to correctly measure the laser fluctuations and the weaker (gravitational-wave) signals simultaneously. An additional safety factor of \(\approx 10\) should be sufficient to avoid saturation if the noises are well described by Gaussian statistics. In terms of requirements on the digital signal processing subsystem, this dynamic range implies that approximately 36 bits are needed when combining the signals in TDI, only to bridge the gap between laser frequency noise and the other noises and gravitational-wave signals. More bits might be necessary to provide enough information to efficiently filter the data when extracting weak gravitational-wave signals embedded into noise.

The phasemeters will be the onboard instrument that will perform the phase measurements containing the gravitational signals. They will also need to simultaneously measure the time-delays to be applied to the TDI combinations via ranging tones over-imposed on the laser beams exchanged by the spacecraft. And they will need to have the capability of simultaneously measuring additional side-band tones that are required for the calibration of the onboard Ultra-Stable Oscillator used in the down-conversion of heterodyned carrier signal (Hellings 2001; Tinto et al. 2002b).

Work toward the realization of a phasemeter capable of meeting these very stringent performance and operational requirements has aggressively been performed both in the United States and in Europe (Shaddock et al. 2006; Gerberding et al. 2013a, b; Burnett 2010; Wang 2013), and we refer the reader interested in the technical details associated with the development studies of such device to the above references and references therein.

7.6 USO noise calibration

The inter-spacecraft heterodyne one-way measurements are made by interfering the received and outgoing laser light. Since the received and receiving frequencies of the laser beams can be different by tens to perhaps hundreds of MHz (a consequence of the Doppler effect from the relative inter-spacecraft velocities and the intrinsic frequency differences of the lasers), one needs to remove these large beat notes present in the heterodyne measurements before they are digitized. This is done by relying on the use of a microwave signal generated by an onboard clock, usually referred to as an Ultra-Stable Oscillator (USO). The magnitude of the frequency fluctuations introduced by the USOs into the heterodyne measurements depends linearly on the USOs’ noises themselves and the heterodyne beat-note frequencies determined by the inter-spacecraft relative velocities. Space-qualified, state-of-the-art clocks are oven-stabilized crystals characterized by an Allan standard deviation (Barnes et al. 1971) of \({\sigma }_A ( \tau ) \approx 10^{-13}\) for averaging times \( \tau \in [1 - 10{,}000]\,\mathrm {s}\), covering the frequency band of interest to space-based interferometers (Amaro-Seoane et al. 2017; Hu and Wu 2017; Luo et al. 2016; Tinto et al. 2013). In the case of the LISA mission, in particular, it has been estimated (Tinto et al. 2002b) that the magnitude of the power spectral density of the USOs’ relative frequency fluctuations appearing, for instance, in the first-generation TDI combination X may be 3 orders of magnitude larger than those due to the residual (optical path and proof mass) noise sources.

The technique currently baselined by the LISA project for removing the USO noise from the TDI combinations requires the modulation of the laser beams exchanged by the spacecraft, and the further measurement of six more inter-spacecraft relative phases by comparing the side-bands of the received beam against side-bands of the transmitted beam (Hellings 2001; Tinto et al. 2002b, 2007; Tinto and Hartwig 2018; Hartwig and Bayle 2020). The physical reason behind the use of modulated beams for calibrating the USOs’ noises is to exchange the USOs’ phase fluctuations with the same time delays as those experienced by the lasers among the three spacecraft by performing side-bands/side-bands measurements. In so doing, six additional phase measurements are generated that allow one to calibrate out the USOs’ phase fluctuations from the TDI combinations while preserving the gravitational wave signal in the resulting USO-calibrated TDI data (see Otto et al. 2012; Heinzel et al. 2011; Tinto and Hartwig 2018; Hartwig and Bayle 2020 for more details on the USO calibration method).

An alternative technique to the one highlighted above for calibrating the USOs noises has recently been proposed (Tinto and Yu 2015). It is based on the observation that, by coherently transferring the laser phase fluctuations to the microwave signal used in the heterodyne measurements, one would then need to cancel only the laser noise (the noise due to the microwave signal is now proportional to it!) by deriving some suitably modified TDI combinations. Coherently linking optical laser frequencies to microwave frequencies was thought to be impossible because of the inability to directly count the optical frequency of a laser. However, with the advent of the self-referenced octave-span optical frequency comb (OFC) scheme—for which Hall and Hänsch received the 2005 physics Nobel Prize (Hall 2006; Hänsch 2006)—it is possible to generate a microwave frequency signal that is coherent to the frequency of the laser at a level significantly better than the frequency stability required of a USO (Ma et al. 2004). Given the ongoing progress made towards space-qualified, self-referenced OFCs  (Wilken et al. 2013; Lee et al. 2014; Lezius et al. 2016), it is likely they will become commercially available within the next few years. Although this might be too late for incorporating them in the LISA design, the resulting hardware simplification and improved system reliability should seriously be considered by future space-based interferometer projects.

8 Concluding remarks

In this third revision of our article on TDI, we have summarized its ability to cancel the laser phase noise from heterodyne phase measurements performed by a constellation of three spacecraft tracking each other along arms of unequal length. Underlying the TDI technique is the mathematical structure of the theory of Gröbner basis and the algebra of modules over polynomial rings. These methods have been motivated and illustrated with the simple example of an unequal-arm interferometer in order to give a physical insight of TDI. Here, these techniques have been rigorously applied to the idealized case of a stationary interferometer such as the LISA mission. This allowed us to derive the generators of the module from which the entire TDI data set can be obtained. The formulation of TDI presented in this review article has been base-lined by other mission concepts (Luo et al. 2016; Hu and Wu 2017) besides the LISA mission (Ni 2016). A straight-forward generalization of it has also been extended to configurations with more than three satellites so as to cancel the spacecraft acceleration noises. Although an octahedron configuration (Wang et al. 2013) can in principle cancel simultaneously the six laser noises and the eighteen components of the acceleration noises by applying TDI to the twenty four one-way measurements made by such an array, the requirements on the trajectories of the six spacecraft become a major road-block for successfully flying such an interesting mission concept. The stationary LISA case was used as a propaedeutical introduction to the physical motivation of TDI, and for further extending it to the realistic LISA configuration of free-falling spacecraft orbiting around the Sun. The TDI data combinations canceling laser phase noise in this general case are referred to as second-generation TDI, and they contain twice as many terms as their corresponding first-generation combinations valid for the stationary configuration.

As a data analysis application we have shown that it is possible to identify specific TDI combinations that will allow LISA to achieve optimal sensitivity to gravitational radiation (Prince et al. 2002; Rajesh Nayak et al. 2003a, b). The resulting improvement in sensitivity over that of an unequal-arm Michelson interferometer, in the case of signals randomly distributed over the celestial sphere and of random polarization, is non-negligible. We have found this to be equal to a factor of \(\sqrt{2}\) in the low-part of the LISA frequency band, and slightly more than \(\sqrt{3}\) in the high-part of it. The SNR for binaries whose location in the sky is known, but their polarization is not, can also be optimized, and the degree of improvement depends on the location of the source in the sky.

We also addressed several experimental aspects of TDI, and emphasized that it has already been successfully tested experimentally (de Vine et al. 2010; Miller 2010; Spero et al. 2011; Mitryk et al. 2012).

As of the writing of this third edition of our Living Review article, it is very gratifying to see how much TDI has matured since the publishing of its first version, and how many scientists have been and are currently contributing to its theory and experimental implementation by the forthcoming LISA mission. The purpose of this third edition review of TDI was to provide the basic mathematical tools and knowledge of the current experimental results needed for working on future TDI projects. We hope to have accomplished this goal, and that the next generation of gravitational-wave physicists working on future space-based GW missions will further expand this fascinating field of research.