# Laser noise residuals in LISA from onboard processing and time-delay interferometry

Martin Staab, Marc Lilley, Jean-Baptiste Bayle, Olaf Hartwig

► To cite this version:

Martin Staab, Marc Lilley, Jean-Baptiste Bayle, Olaf Hartwig. Laser noise residuals in LISA from onboard processing and time-delay interferometry. Physical Review D, 2024, 109 (4), pp.043040. 10.1103/PhysRevD.109.043040. hal-04154315

HAL Id: hal-04154315

https://hal.science/hal-04154315

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Laser noise residuals in LISA from onboard processing and time-delay interferometry

Martin Staab, Marc Lilley, Jean-Baptiste Bayle, Olaf Hartwig

► To cite this version:

Martin Staab, Marc Lilley, Jean-Baptiste Bayle, Olaf Hartwig. Laser noise residuals in LISA from onboard processing and time-delay interferometry. Physical Review D, 2024, 109 (4), pp.043040. 10.1103/PhysRevD.109.043040. hal-04154315

HAL Id: hal-04154315

https://hal.science/hal-04154315v1

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Laser noise residuals in LISA from on-board processing and time-delay interferometry

Martin Staab,¹,²,³,* Marc Lilley,³ Jean-Baptiste Bayle,⁴ and Olaf Hartwig³,¹,²,⁵

¹Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hanover, Germany

²Leibniz Universität Hannover, D-30167 Hanover, Germany

³SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université,
LNE, 61 avenue de l'Observatoire 75014 Paris, France

⁴University of Glasgow, Glasgow G12 8QQ, United Kingdom

⁵Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam, Germany

<figure><img src="image_2.png" /></figure>

(Received 22 June 2023; accepted 23 January 2024; published 22 February 2024)

Time-delay interferometry (TDI) is a crucial step in the on-ground data processing pipeline of the Laser Interferometer Space Antenna (LISA), as it reduces otherwise overwhelming laser noise and allows for the detection of gravitational waves (GWs). This being said, several laser noise couplings have been identified that limit the performance of TDI. First, on-board processing, which is used to decimate the sampling rate from tens of MHz down to a few Hz, requires careful design of the antialiasing filters to mitigate folding of laser noise power into the observation band. Furthermore, the flatness of those filters is important to limit the effect of the flexing-filtering coupling. Second, the postprocessing delays applied in TDI are subject to ranging and interpolation errors. All of these effects are partially described in the literature. In this paper, we present them in a unified framework and give a more complete description of aliased laser noise and the coupling of interpolation errors. Furthermore, for the first time, we discuss the impact of laser locking on laser noise residuals in the final TDI output. To verify the validity of the analytic power spectral density (PSD) models we derive, we run numerical simulations using LISA INSTRUMENT and calculate second-generation TDI variables with PyTDI. We consider a setup with six independent lasers and with locked lasers (locking configuration N1-12). We find that laser locking indeed affects the laser noise residuals in the TDI combinations as it introduces correlations among the six lasers inducing slight modulations of the PSDs compared to the case of six independent lasers. This implies further studies on laser noise residuals should consider the various locking configurations to produce accurate results.

DOI: 10.1103/PhysRevD.109.043040

I. INTRODUCTION

The Laser Interferometer Space Antenna (LISA) is a space
mission led by the European Space Agency (ESA), expected
to be launched in the 2030s. Its goal is to detect gravitational
waves (GWs) in a frequency band ranging from 10⁻⁴ Hz to
1 Hz [1]. High precision interferometric measurements will be
made via the exchange of laser beams among three spacecraft
orbiting the Sun and separated by 2.5 million kilometers, in
order to determine the variations in the distance between free-
falling test masses aboard each spacecraft to picometer
precision. In these measurements, laser noise is the primary

noise source and is 8 orders of magnitude larger than the GW
signals that one hopes to detect. TDI is a data processing
technique that combines the LISA measurements to construct
virtual equal-arm two-beam interferometers in order to reduce
the laser noise to levels sufficiently low such that GWs
become detectable [2,3]. In TDI the measurements are time-
shifted by multiples of the LISA arm lengths and combined in
a specific scheme to achieve laser noise reduction. Second-
generation TDI, which is the current baseline laser noise
reduction strategy for LISA, applies to the case in which the
arm lengths of the LISA constellation evolve slowly and
linearly in time [4,5]. In second generation TDI, laser noise is
strongly suppressed and the residual is fundamentally limited
by the arm length mismatch of the virtual interferometer [6].

There exist other approaches to perform laser noise
suppression. In TDI-∞ [7], the observables that cancel laser
noise are obtained numerically by solving for the null space
of the design matrix, i.e., the way the various noise sources
enter the interferometric measurements, for an arbitrary time
dependence of the arm lengths. The likelihood function that

<footnote>*Corresponding author: martin.staab@obspm.fr</footnote>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI. Open access publication funded by the Max Planck
Society.</footnote>
Laser noise residuals in LISA from on-board processing and time-delay interferometry

Martin Staab,¹,²,³,* Marc Lilley,³ Jean-Baptiste Bayle,⁴ and Olaf Hartwig³,¹,²,⁵

¹Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hanover, Germany

²Leibniz Universität Hannover, D-30167 Hanover, Germany

³SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université,
LNE, 61 avenue de l'Observatoire 75014 Paris, France

⁴University of Glasgow, Glasgow G12 8QQ, United Kingdom

⁵Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam, Germany

<figure><img src="image_2.png" /></figure>

(Received 22 June 2023; accepted 23 January 2024; published 22 February 2024)

Time-delay interferometry (TDI) is a crucial step in the on-ground data processing pipeline of the Laser Interferometer Space Antenna (LISA), as it reduces otherwise overwhelming laser noise and allows for the detection of gravitational waves (GWs). This being said, several laser noise couplings have been identified that limit the performance of TDI. First, on-board processing, which is used to decimate the sampling rate from tens of MHz down to a few Hz, requires careful design of the antialiasing filters to mitigate folding of laser noise power into the observation band. Furthermore, the flatness of those filters is important to limit the effect of the flexing-filtering coupling. Second, the postprocessing delays applied in TDI are subject to ranging and interpolation errors. All of these effects are partially described in the literature. In this paper, we present them in a unified framework and give a more complete description of aliased laser noise and the coupling of interpolation errors. Furthermore, for the first time, we discuss the impact of laser locking on laser noise residuals in the final TDI output. To verify the validity of the analytic power spectral density (PSD) models we derive, we run numerical simulations using LISA INSTRUMENT and calculate second-generation TDI variables with PyTDI. We consider a setup with six independent lasers and with locked lasers (locking configuration N1-12). We find that laser locking indeed affects the laser noise residuals in the TDI combinations as it introduces correlations among the six lasers inducing slight modulations of the PSDs compared to the case of six independent lasers. This implies further studies on laser noise residuals should consider the various locking configurations to produce accurate results.

DOI: 10.1103/PhysRevD.109.043040

# I. INTRODUCTION

The Laser Interferometer Space Antenna (LISA) is a space mission led by the European Space Agency (ESA), expected to be launched in the 2030s. Its goal is to detect gravitational waves (GWs) in a frequency band ranging from $10^{-4}$ Hz to 1 Hz [1]. High precision interferometric measurements will be made via the exchange of laser beams among three spacecraft orbiting the Sun and separated by 2.5 million kilometers, in order to determine the variations in the distance between free-falling test masses aboard each spacecraft to picometer precision. In these measurements, laser noise is the primary

noise source and is 8 orders of magnitude larger than the GW signals that one hopes to detect. TDI is a data processing technique that combines the LISA measurements to construct virtual equal-arm two-beam interferometers in order to reduce the laser noise to levels sufficiently low such that GWs become detectable [2,3]. In TDI the measurements are time-shifted by multiples of the LISA arm lengths and combined in a specific scheme to achieve laser noise reduction. Second-generation TDI, which is the current baseline laser noise reduction strategy for LISA, applies to the case in which the arm lengths of the LISA constellation evolve slowly and linearly in time [4,5]. In second generation TDI, laser noise is strongly suppressed and the residual is fundamentally limited by the arm length mismatch of the virtual interferometer [6].

There exist other approaches to perform laser noise suppression. In TDI-$\infty$ [7], the observables that cancel laser noise are obtained numerically by solving for the null space of the design matrix, i.e., the way the various noise sources enter the interferometric measurements, for an arbitrary time dependence of the arm lengths. The likelihood function that

<footnote>*Corresponding author: martin.staab@obspm.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.</footnote>
is used in GW source parameter estimation can then be written directly in the time domain in terms of the LISA interferometric measurements without having to reformulate the entire problem in terms of algebraically defined TDI variables. While the study in [7] was limited to an idealized toy model with a single Michelson interferometer, the authors of [8] applied TDI-∞ to the full LISA constellation with time-evolving arms. Computationally, TDI-∞ has the drawback that it requires the storage and manipulation of very large matrices.

In [9], starting from the interferometric measurements and for nonevolving LISA arms, the authors first form a matrix of integer-delayed measurements, which they decompose using principal component analysis (PCA) into high and low variance components. The latter correspond to the components for which the laser noise is significantly suppressed. This approach, dubbed “automated principal component interferometry”, or aPCI, is formulated in both the time and frequency domains. In [10], the same authors extend this approach to the case of time-evolving arms. Note that in [10], aPCI is shown not to perform as well as second-generation TDI in suppressing laser noise.

While approaches such as TDI-∞ and aPCI offer some interesting perspectives for a flexible data-driven formulation of TDI, “traditional” TDI can be formulated analytically. It is therefore tractable, better understood, and exact analytic transfer functions exist to describe the instrumental noise residuals present in the TDI variables. For instance, secondary noises such as test-mass acceleration noise and optical metrology system noise are dealt with in [11], clock noise is studied in [12,13], and tilt-to-length coupling in [14]. Note also that laser noise coupling residuals were discussed previously in [15,16]. It is also worth stressing that a good understanding of the noise content in the final TDI output is crucial to characterize the performance of LISA and guide the design of the instrument and that data analysis and parameter estimation will require accurate noise models in order to work reliably, making these studies particularly relevant for LISA.

In addition to the analytic and numerical studies available in the literature, there exist several hardware demonstrators that test various aspects of TDI experimentally. The LISA interferometry test-bed [17], while it could not reproduce the signal delays in a realistic way, did demonstrate for the first time that using first generation TDI both the laser and clock noise could be suppressed by 9 and 4–5 orders of magnitude, respectively. In UFLIS [18,19], using electronic phase delay units allowing for time-varying delays of the laser phases, the authors were able to demonstrate the efficacy of second-generation TDI. In more recent work [20], the Hexagon experiment demonstrated that clock synchronization can be achieved to sufficient accuracy to match the LISA requirements. Moreover, the authors find residual laser noise after TDI-like processing due to flexing-filtering, aliasing, and interpolation error. The LISA on table (LOT) experiment [21] (for recent progress see [22]) is an electro-optical setup aiming primarily at testing the laser noise suppression performance of TDI. In [22], the validity of second-generation TDI was demonstrated for linearly evolving LISA arms. It was also shown using an analytic model that the residual noise could be explained by the cascade integrator comb filtering and the decimation stages that are applied to the data.

In this paper, we study the coupling of laser noise residuals in standard TDI. We focus on the residual laser noise due to systematic effects and neglect most other noise sources. The one exception is noise in the ranging measurements which are used as delays in TDI, and which in principle couple to laser noise. Following [13], we assume this noise source will be strongly suppressed to the level of the highly-precise sideband interferometer readouts, such that its impact on the laser noise reduction is minor. We still include it so as to have a more complete description of the postprocessing delay. We consider the effect of on-board processing (i.e., filtering and decimation), and the influence of TDI which uses postprocessing delay operations that are subject to ranging and interpolation errors, see Fig. 1. We compute analytic formulas for all laser noise residuals induced by these processing steps and compare those to

<figure><img src="image_3.png" /><figcaption>FIG. 1. Flowchart illustrating beam propagation, online and offline processing for the interspacecraft interferometric beatnote phase. The propagation and postprocessing delay are defined in Sec. II and Sec. III B, respectively. Note that the online and offline processing of the reference interferometer, not shown here, is needed to compute the final output of this diagram, the TDI variables.</figcaption></figure>
is used in GW source parameter estimation can then be written directly in the time domain in terms of the LISA interferometric measurements without having to reformulate the entire problem in terms of algebraically defined TDI variables. While the study in [7] was limited to an idealized toy model with a single Michelson interferometer, the authors of [8] applied TDI-∞ to the full LISA constellation with time-evolving arms. Computationally, TDI-∞ has the drawback that it requires the storage and manipulation of very large matrices.

In [9], starting from the interferometric measurements and for nonevolving LISA arms, the authors first form a matrix of integer-delayed measurements, which they decompose using principal component analysis (PCA) into high and low variance components. The latter correspond to the components for which the laser noise is significantly suppressed. This approach, dubbed “automated principal component interferometry”, or aPCI, is formulated in both the time and frequency domains. In [10], the same authors extend this approach to the case of time-evolving arms. Note that in [10], aPCI is shown not to perform as well as second-generation TDI in suppressing laser noise.

While approaches such as TDI-∞ and aPCI offer some interesting perspectives for a flexible data-driven formulation of TDI, “traditional” TDI can be formulated analytically. It is therefore tractable, better understood, and exact analytic transfer functions exist to describe the instrumental noise residuals present in the TDI variables. For instance, secondary noises such as test-mass acceleration noise and optical metrology system noise are dealt with in [11], clock noise is studied in [12,13], and tilt-to-length coupling in [14]. Note also that laser noise coupling residuals were discussed previously in [15,16]. It is also worth stressing that a good understanding of the noise content in the final TDI output is crucial to characterize the performance of LISA and guide the design of the instrument and that data analysis and parameter estimation will require accurate noise models in order to work reliably, making these studies particularly relevant for LISA.

In addition to the analytic and numerical studies available in the literature, there exist several hardware demonstrators that test various aspects of TDI experimentally. The LISA interferometry test-bed [17], while it could not reproduce the signal delays in a realistic way, did demonstrate for the first time that using first generation TDI both the laser and clock noise could be suppressed by 9 and 4–5 orders of magnitude, respectively. In UFLIS [18,19], using electronic phase delay units allowing for time-varying delays of the laser phases, the authors were able to demonstrate the efficacy of second-generation TDI. In more recent work [20], the Hexagon experiment demonstrated that clock synchronization can be achieved to sufficient accuracy to match the LISA requirements. Moreover, the authors find residual laser noise after TDI-like processing due to flexing-filtering, aliasing, and interpolation error. The LISA on table (LOT) experiment [21] (for recent progress see [22]) is an electro-optical setup aiming primarily at testing the laser noise suppression performance of TDI. In [22], the validity of second-generation TDI was demonstrated for linearly evolving LISA arms. It was also shown using an analytic model that the residual noise could be explained by the cascade integrator comb filtering and the decimation stages that are applied to the data.

In this paper, we study the coupling of laser noise residuals in standard TDI. We focus on the residual laser noise due to systematic effects and neglect most other noise sources. The one exception is noise in the ranging measurements which are used as delays in TDI, and which in principle couple to laser noise. Following [13], we assume this noise source will be strongly suppressed to the level of the highly-precise sideband interferometer readouts, such that its impact on the laser noise reduction is minor. We still include it so as to have a more complete description of the postprocessing delay. We consider the effect of on-board processing (i.e., filtering and decimation), and the influence of TDI which uses postprocessing delay operations that are subject to ranging and interpolation errors, see Fig. 1. We compute analytic formulas for all laser noise residuals induced by these processing steps and compare those to

<figure><img src="image_3.png" /><figcaption>FIG. 1. Flowchart illustrating beam propagation, online and offline processing for the interspacecraft interferometric beatnote phase. The propagation and postprocessing delay are defined in Sec. II and Sec. III B, respectively. Note that the online and offline processing of the reference interferometer, not shown here, is needed to compute the final output of this diagram, the TDI variables.</figcaption></figure>
numerical simulations obtained using LISA INSTRUMENT and PyTDI. We do so for six independent lasers and also for locked lasers in the N1-12 locking scheme [23,24].

The laser noise residuals induced by the processing steps we derive below are already partially described in the literature. The flexing-filtering effect arising from the noncommutativity of the filtering and delay operation was previously discussed in [15]. The impact of ranging errors, i.e., modulation noise and ranging biases, in TDI was derived in [13], while a preliminary model for the aliased laser noise due to decimation and interpolation errors was derived in [16]. Here, we gather these results in a unified framework and check them against the most up-to-date instrumental setups. In addition, we include the effect of laser locking in our models, an important feature of LISA that was previously neglected. We show that laser locking can amplify the laser noise residual in the TDI combinations, so that future performance studies on laser noise coupling should consider the influence of the various locking configurations that are available. Finally, we correct the preliminary models presented in [16], which did not show perfect agreement with simulation results.

The paper is organized as follows. In Sec. II, we introduce the interferometric measurements available in LISA, the notion of discretely sampled time series and how filtering and decimation apply to discretized data. In Sec. III A we discuss the residual laser noise induced by the noncommutativity of the on-board filtering and decimation operations with the propagation delays. The interpolation and ranging residuals, which both arise in TDI, which is performed on ground, are discussed in Sec. III B. In Sec. IV, we apply the results of the previous sections to the case of the second-generation TDI Michelson variables for six independent lasers and for locked lasers. Additional details are given in six appendices.

## II. INTERFEROMETRIC MEASUREMENTS

LISA produces two main interferometric measurements per movable optical subassembly (MOSA) relevant for laser noise reduction.<fnref n="1" /> Those are the interspacecraft and reference interferometers given by

$$isi_{ij}(t) = \mathbf{SF}(\mathbf{D}_{ij}\phi_{ji}(t) - \phi_{ij}(t)), \quad (1a)$$

$$rfi_{ij}(t) = \mathbf{SF}(\phi_{ik}(t) - \phi_{ij}(t)). \quad (1b)$$

As illustrated in Fig. 1, each measurement represents the beatnote phase formed by two laser beams, whose phases are denoted by $\phi$ and labeled by the index pair $ij$. Here, we follow the conventions in [6], where $i$ denotes the hosting spacecraft and $j$ the spacecraft from which MOSA $ij$

receives light. The interspacecraft interferometer $isi_{ij}$ tracks the phase difference between the distant laser $\phi_{ji}$ that is propagated to the local MOSA and the local laser $\phi_{ij}$. Beam propagation is equivalent to applying the delay operator $\mathbf{D}_{ij}$ defined as

$$\mathbf{D}_{ij}\phi_{ji}(t) = \phi_{ji}(t - d_{ij}(t)). \quad (2)$$

Here, $d_{ij}(t)$ is the pseudorange, which includes the light travel time in some chosen global frame and any reference frame transformation accounting for the fact that the phases $\phi_{ji}$ are defined in their respective reference frame, $j$. In order to relate phases on the left and right MOSAs on each spacecraft,<fnref n="2" /> the reference interferometer $rfi_{ij}$ combines the local and the adjacent lasers.

It is useful to decompose the total laser phase $\phi_{ij}(t)$ (here expressed in units of cycles) or its derivative, the total frequency $\nu_{ij}(t) = \dot{\phi}_{ij}(t)$, into two variables. As we shall see below, laser and timing noise residuals due to on-board processing and TDI will couple differently to the phase ramp $\phi_{ij}^o(t)$ and to any in-band fluctuations $\phi_{ij}^e(t)$. We thus write

$$\phi_{ij}(t) = \phi_{ij}^o(t) + \phi_{ij}^e(t), \quad (3)$$

$$\nu_{ij}(t) = \nu_{ij}^o(t) + \nu_{ij}^e(t), \quad (4)$$

where $\nu_{ij}^o(t) = \dot{\phi}_{ij}^o(t)$ describes any slowly-varying drifts around the central laser frequency $\nu_0 = 281.6$ THz and $\nu_{ij}^e(t) = \dot{\phi}_{ij}^e(t)$ accounts for any rapidly varying random fluctuations. This in-band part is dominated by laser frequency noise with an ASD $\sqrt{S_p} = 30$ Hz/√Hz.

The beatnote phases of the interferometers are read out using a digital phase locked loop running at 80 MHz. Multiple decimation stages reduce the sampling rate down to 4 Hz in order to produce the final data streams telemetered to ground (see Fig. 1). Each decimation stage consists of an antialiasing filter, $F$, and a downsampling stage, $S$, which reduces the sampling rate of the data by an integer factor. In this work, we compare the analytic models that we derive to the most recent LISA simulation codes, which run at rates that are much lower than the 80 MHz quoted above, and thus only use a single decimation stage. This being said, the results obtained in this paper can easily be generalized to multiple decimation stages.

We shall express signals in continuous time so as to be compatible with the recent literature on TDI. However, the application of finite impulse response (FIR) filters,

<footnote n="1">¹The split interferometry configuration involves a third interferometer, the test-mass interferometer, which is not relevant for the purpose of this study.</footnote>

<footnote n="2">²Each spacecraft is equipped with two MOSAs. The left-handed MOSA on spacecraft $i$ refers to the one facing spacecraft $i+1$, while the right-handed MOSA refers to the one facing spacecraft $i-1$ (indices ranging from 1 to 3 cyclic).</footnote>
numerical simulations obtained using LISA INSTRUMENT and PyTDI. We do so for six independent lasers and also for locked lasers in the N1-12 locking scheme [23,24].

The laser noise residuals induced by the processing steps we derive below are already partially described in the literature. The flexing-filtering effect arising from the noncommutativity of the filtering and delay operation was previously discussed in [15]. The impact of ranging errors, i.e., modulation noise and ranging biases, in TDI was derived in [13], while a preliminary model for the aliased laser noise due to decimation and interpolation errors was derived in [16]. Here, we gather these results in a unified framework and check them against the most up-to-date instrumental setups. In addition, we include the effect of laser locking in our models, an important feature of LISA that was previously neglected. We show that laser locking can amplify the laser noise residual in the TDI combinations, so that future performance studies on laser noise coupling should consider the influence of the various locking configurations that are available. Finally, we correct the preliminary models presented in [16], which did not show perfect agreement with simulation results.

The paper is organized as follows. In Sec. II, we introduce the interferometric measurements available in LISA, the notion of discretely sampled time series and how filtering and decimation apply to discretized data. In Sec. III A we discuss the residual laser noise induced by the noncommutativity of the on-board filtering and decimation operations with the propagation delays. The interpolation and ranging residuals, which both arise in TDI, which is performed on ground, are discussed in Sec. III B. In Sec. IV, we apply the results of the previous sections to the case of the second-generation TDI Michelson variables for six independent lasers and for locked lasers. Additional details are given in six appendices.

## II. INTERFEROMETRIC MEASUREMENTS

LISA produces two main interferometric measurements per movable optical subassembly (MOSA) relevant for laser noise reduction.<fnref n="1" /> Those are the interspacecraft and reference interferometers given by

$$isi_{ij}(t) = \mathbf{SF}(\mathbf{D}_{ij}\phi_{ji}(t) - \phi_{ij}(t)), \quad (1a)$$

$$r\mathbf{fi}_{ij}(t) = \mathbf{SF}(\phi_{ik}(t) - \phi_{ij}(t)), \quad (1b)$$

As illustrated in Fig. 1, each measurement represents the beatnote phase formed by two laser beams, whose phases are denoted by $\phi$ and labeled by the index pair $ij$. Here, we follow the conventions in [6], where $i$ denotes the hosting spacecraft and $j$ the spacecraft from which MOSA $ij$

receives light. The interspacecraft interferometer $isi_{ij}$ tracks the phase difference between the distant laser $\phi_{ji}$ that is propagated to the local MOSA and the local laser $\phi_{ij}$. Beam propagation is equivalent to applying the delay operator $\mathbf{D}_{ij}$ defined as

$$\mathbf{D}_{ij}\phi_{ji}(t) = \phi_{ji}(t - d_{ij}(t)). \quad (2)$$

Here, $d_{ij}(t)$ is the pseudorange, which includes the light travel time in some chosen global frame and any reference frame transformation accounting for the fact that the phases $\phi_{ji}$ are defined in their respective reference frame, $j$. In order to relate phases on the left and right MOSAs on each spacecraft,<fnref n="2" /> the reference interferometer $r\mathbf{fi}_{ij}$ combines the local and the adjacent lasers.

It is useful to decompose the total laser phase $\phi_{ij}(t)$ (here expressed in units of cycles) or its derivative, the total frequency $\nu_{ij}(t) = \dot{\phi}_{ij}(t)$, into two variables. As we shall see below, laser and timing noise residuals due to on-board processing and TDI will couple differently to the phase ramp $\phi_{ij}^o(t)$ and to any in-band fluctuations $\phi_{ij}^e(t)$. We thus write

$$\phi_{ij}(t) = \phi_{ij}^o(t) + \phi_{ij}^e(t), \quad (3)$$

$$\nu_{ij}(t) = \nu_{ij}^o(t) + \nu_{ij}^e(t), \quad (4)$$

where $\nu_{ij}^o(t) = \dot{\phi}_{ij}^o(t)$ describes any slowly-varying drifts around the central laser frequency $\nu_0 = 281.6$ THz and $\nu_{ij}^e(t) = \dot{\phi}_{ij}^e(t)$ accounts for any rapidly varying random fluctuations. This in-band part is dominated by laser frequency noise with an ASD $\sqrt{S_p} = 30$ Hz/√Hz.

The beatnote phases of the interferometers are read out using a digital phase locked loop running at 80 MHz. Multiple decimation stages reduce the sampling rate down to 4 Hz in order to produce the final data streams telemetered to ground (see Fig. 1). Each decimation stage consists of an antialiasing filter, $F$, and a downsampling stage, $S$, which reduces the sampling rate of the data by an integer factor. In this work, we compare the analytic models that we derive to the most recent LISA simulation codes, which run at rates that are much lower than the 80 MHz quoted above, and thus only use a single decimation stage. This being said, the results obtained in this paper can easily be generalized to multiple decimation stages.

We shall express signals in continuous time so as to be compatible with the recent literature on TDI. However, the application of finite impulse response (FIR) filters,

<footnote n="1">¹The split interferometry configuration involves a third interferometer, the test-mass interferometer, which is not relevant for the purpose of this study.</footnote>

<footnote n="2">²Each spacecraft is equipped with two MOSAs. The left-handed MOSA on spacecraft $i$ refers to the one facing spacecraft $i+1$, while the right-handed MOSA refers to the one facing spacecraft $i-1$ (indices ranging from 1 to 3 cyclic).</footnote>
decimation, and interpolation requires some notion of discretely sampled time series. We will therefore make use of the Whittaker-Shannon interpolation formula [e.g., [25]],

$$x(t) = \sum_{n=-\infty}^{\infty} \text{sinc}(f_s t - n) \cdot x_n, \quad (5)$$

which reconstructs the continuous time signal $x(t)$ from discrete samples $x_n$.

Let us first describe the on-board processing, which consists of the application of a FIR filter and a decimation stage. A FIR filter is equivalent to a discrete convolution of the input time series $x_n$ with filter taps $h_m$,

$$y_n = \sum_m h_m \cdot x_{n-m}. \quad (6)$$

We use Eq. (5) to represent the output $y_n$ in continuous time and find

$$y(t) = \mathbf{F}x(t) = \underbrace{\int_{\mathbb{R}} \sum_m h_m \delta(\tau - mT_s) x(t - \tau) d\tau}_{h_F(\tau)}, \quad (7)$$

where we introduce the integral over the Dirac-delta distribution $\delta(t)$ to shift the time argument of $x(t)$. It follows that the application of a FIR filter is equivalent to a continuous time convolution with the filter kernel $h_F(\tau)$ defined above. Using the usual definition of the one-sided power spectral density (PSD), the PSD of the filtered process in Eq. (7) is given by

$$S_y(f) = |\tilde{\mathbf{F}}|^2 S_x(f) = |\tilde{h}_F(f)|^2 S_x(f) \quad (8)$$

with $\tilde{\mathbf{F}} = \tilde{h}(f)$ the Fourier transform of $h(t)$, while $S_x(f)$ is the PSD of $x(t)$.

Let us now discuss the decimation operator which we use to reduce the sampling rate by an integer factor $M$. On a discrete time grid, the resulting signal is given by

$$y_n = x_{n \cdot M}. \quad (9)$$

Again, we make use of Eq. (5) to find the corresponding continuous time representation,

$$y(t) = \mathbf{S}x(t) = \sum_{n=-\infty}^{\infty} \text{sinc}(f_s t - n) \cdot x_{n \cdot M}. \quad (10)$$

Here, $f_s$ denotes the sampling rate after decimation, and $\mathbf{S}$ symbolizes the action of the decimation operation in continuous time. The right-hand side of Eq. (10) is exactly equal to $x(t)$ if and only if it has a band limit that is less than the Nyquist rate after decimation, $f_n = f_s/2$. Otherwise aliasing occurs, which folds power from frequencies above

$f_n$ into the band $[0, f_n]$. This effect becomes apparent when looking at the corresponding one-sided PSD,

$$S_y(f) = |\tilde{\mathbf{S}}|^2 S_x(f) = \text{rect}\left(\frac{f}{f_s}\right) \sum_{n=0}^{M-1} S_x^{(n)}(f). \quad (11)$$

Here, we introduce the shorthand notation $|\tilde{\mathbf{S}}|^2$ which represents the action of the decimation operator on the PSD of $x(t)$. Note that decimation is not a linear operation and thus $|\tilde{\mathbf{S}}|^2$ is not a pure multiplicative factor as in Eq. (8).

The rectangular function introduced in Eq. (11) is defined to be equal to zero for $|f| > f_n$ and equal to one for $|f| < f_n$ such that the decimated signal is band-limited up to the new Nyquist rate. Finally, the $n$th alias, $S_x^{(n)}(f)$, on the right-hand-side of Eq. (11), is given by

$$S_x^{(n)}(f) = \begin{cases} S_x(n f_n + f) & \text{if } n \text{ is even,} \\ S_x((n+1)f_n - f) & \text{if } n \text{ is odd.} \end{cases} \quad (12)$$

Equation (11) and (12) highlight the typical folding into band of any spectral component that resides at frequencies higher than the new Nyquist rate (up to the highest frequency $Mf_n$, which corresponds to the Nyquist rate before decimation).

### III. SECOND-GENERATION MICHELSON COMBINATIONS

In this section we first explain how the numerical computation of the second-generation Michelson combination $X_2$ is implemented. Therefore, we need to introduce the postprocessing delay operator $\mathcal{D}$ which is used to time-shift the beatnote measurement in TDI. We then discuss how the two broad classes of residuals, namely the on-board processing residuals and TDI residuals, that contribute to $X_2$ arise. The detailed calculations for each residual, and the comparisons between model and simulation are given in appendices B and C.

In order to optimize numerical precision and to save computational cost, we calculate $X_2$ in several stages using the following intermediary variables. As shown below, the variable $\eta$ is constructed from the interspacecraft and reference interferometers. This step reduces the number of lasers from six to three. Then, the variables $\pi$, $\rho$, and $\sigma$ are constructed from $\eta$ by building round trip interferometers of increasing complexity. We have

$$\eta_{ij} = \begin{cases} \text{isi}_{ij} - \hat{\mathcal{D}}_{ij} \frac{\text{rfi}_{jk} - \text{rfi}_{ji}}{2} & \text{if } \epsilon_{ijk} = 1, \\ \text{isi}_{ij} + \frac{\text{rfi}_{ik} - \text{rfi}_{ji}}{2} & \text{if } \epsilon_{ijk} = -1, \end{cases} \quad (13a)$$

$$\pi_{ij} = \eta_{ij} + \hat{\mathcal{D}}_{ij} \eta_{ji}, \quad (13b)$$

$$\rho_{ij} = \pi_{ij} + \hat{\mathcal{D}}_{iji} \pi_{ik}, \quad (13c)$$

$$\sigma_{ij} = \rho_{ij} + \hat{\mathcal{D}}_{ijiki} \rho_{ik}. \quad (13d)$$
decimation, and interpolation requires some notion of discretely sampled time series. We will therefore make use of the Whittaker-Shannon interpolation formula [e.g., [25]],

$$x(t) = \sum_{n=-\infty}^{\infty} \text{sinc}(f_s t - n) \cdot x_n, \quad (5)$$

which reconstructs the continuous time signal $x(t)$ from discrete samples $x_n$.

Let us first describe the on-board processing, which consists of the application of a FIR filter and a decimation stage. A FIR filter is equivalent to a discrete convolution of the input time series $x_n$ with filter taps $h_m$,

$$y_n = \sum_m h_m \cdot x_{n-m}. \quad (6)$$

We use Eq. (5) to represent the output $y_n$ in continuous time and find

$$y(t) = \mathbf{F}x(t) = \underbrace{\int_{\mathbb{R}} \sum_m h_m \delta(\tau - mT_s) x(t - \tau) d\tau}_{h_F(\tau)}, \quad (7)$$

where we introduce the integral over the Dirac-delta distribution $\delta(t)$ to shift the time argument of $x(t)$. It follows that the application of a FIR filter is equivalent to a continuous time convolution with the filter kernel $h_F(\tau)$ defined above. Using the usual definition of the one-sided power spectral density (PSD), the PSD of the filtered process in Eq. (7) is given by

$$S_y(f) = |\tilde{\mathbf{F}}|^2 S_x(f) = |\tilde{h}_F(f)|^2 S_x(f) \quad (8)$$

with $\tilde{\mathbf{F}} = \tilde{h}(f)$ the Fourier transform of $h(t)$, while $S_x(f)$ is the PSD of $x(t)$.

Let us now discuss the decimation operator which we use to reduce the sampling rate by an integer factor $M$. On a discrete time grid, the resulting signal is given by

$$y_n = x_{n \cdot M}. \quad (9)$$

Again, we make use of Eq. (5) to find the corresponding continuous time representation,

$$y(t) = \mathbf{S}x(t) = \sum_{n=-\infty}^{\infty} \text{sinc}(f_s t - n) \cdot x_{n \cdot M}. \quad (10)$$

Here, $f_s$ denotes the sampling rate after decimation, and $\mathbf{S}$ symbolizes the action of the decimation operation in continuous time. The right-hand side of Eq. (10) is exactly equal to $x(t)$ if and only if it has a band limit that is less than the Nyquist rate after decimation, $f_n = f_s/2$. Otherwise aliasing occurs, which folds power from frequencies above

$f_n$ into the band $[0, f_n]$. This effect becomes apparent when looking at the corresponding one-sided PSD,

$$S_y(f) = |\tilde{\mathbf{S}}|^2 S_x(f) = \text{rect}\left(\frac{f}{f_s}\right) \sum_{n=0}^{M-1} S_x^{(n)}(f). \quad (11)$$

Here, we introduce the shorthand notation $|\tilde{\mathbf{S}}|^2$ which represents the action of the decimation operator on the PSD of $x(t)$. Note that decimation is not a linear operation and thus $|\tilde{\mathbf{S}}|^2$ is not a pure multiplicative factor as in Eq. (8).

The rectangular function introduced in Eq. (11) is defined to be equal to zero for $|f| > f_n$ and equal to one for $|f| < f_n$ such that the decimated signal is band-limited up to the new Nyquist rate. Finally, the $n$th alias, $S_x^{(n)}(f)$, on the right-hand-side of Eq. (11), is given by

$$S_x^{(n)}(f) = \begin{cases} S_x(n f_n + f) & \text{if } n \text{ is even,} \\ S_x((n+1)f_n - f) & \text{if } n \text{ is odd.} \end{cases} \quad (12)$$

Equation (11) and (12) highlight the typical folding into band of any spectral component that resides at frequencies higher than the new Nyquist rate (up to the highest frequency $Mf_n$, which corresponds to the Nyquist rate before decimation).

### III. SECOND-GENERATION MICHELSON COMBINATIONS

In this section we first explain how the numerical computation of the second-generation Michelson combination $X_2$ is implemented. Therefore, we need to introduce the postprocessing delay operator $\mathcal{D}$ which is used to time-shift the beatnote measurement in TDI. We then discuss how the two broad classes of residuals, namely the on-board processing residuals and TDI residuals, that contribute to $X_2$ arise. The detailed calculations for each residual, and the comparisons between model and simulation are given in appendices B and C.

In order to optimize numerical precision and to save computational cost, we calculate $X_2$ in several stages using the following intermediary variables. As shown below, the variable $\eta$ is constructed from the interspacecraft and reference interferometers. This step reduces the number of lasers from six to three. Then, the variables $\pi$, $\rho$, and $\sigma$ are constructed from $\eta$ by building round trip interferometers of increasing complexity. We have

$$\eta_{ij} = \begin{cases} \text{isi}_{ij} - \hat{\mathcal{D}}_{ij} \frac{\text{rfi}_{jk} - \text{rfi}_{ji}}{2} & \text{if } \epsilon_{ijk} = 1, \\ \text{isi}_{ij} + \frac{\text{rfi}_{ik} - \text{rfi}_{ji}}{2} & \text{if } \epsilon_{ijk} = -1, \end{cases} \quad (13a)$$

$$\pi_{ij} = \eta_{ij} + \hat{\mathcal{D}}_{ij} \eta_{ji}, \quad (13b)$$

$$\rho_{ij} = \pi_{ij} + \hat{\mathcal{D}}_{iji} \pi_{ik}, \quad (13c)$$

$$\sigma_{ij} = \rho_{ij} + \hat{\mathcal{D}}_{ijiki} \rho_{ik}. \quad (13d)$$
<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration of intermediary variables defined in Eq. (13) and locking configuration N1-12. The intermediary variables are depicted as an arrow representing the synthesized photon path of the long arm. They are incrementally built up from the previously defined ones such that the last variable σ is represented by the entire path. Additionally, the chain of locked lasers for the locking configuration N1-12 is shown (as discussed in Sec. IV B). The primary laser is highlighted by the gray box.</figcaption></figure>

Here, $\hat{D}$ denotes the postprocessing operator that also considers ranging errors as indicated by the hat.

Figure 2 provides an illustration of the intermediary variables as they synthesize two-beam interferometers with a long and a short arm. The long arm is depicted as an arrow propagating around the LISA constellation. Using the variable σ we can finally express the second-generation Michelson combination $X_2$ as

$$X_2 = \sigma_{13} - \sigma_{12}. \qquad (14)$$

The contracted postprocessing delay operators $D_{i_1 \cdots i_N}$ appearing in Eq. (13) are applied in two steps. First, the actual nested delay $d_{i_1 \cdots i_N}(t)$ is calculated, and then, the fractional delay filter is applied. The nested delay can be calculated using the recursive operation<fnref n="3" />

$$d_{i_1 \cdots i_N}(t) = d_{i_1 i_2}(t) + D_{i_1 i_2} d_{i_2 \cdots i_N}(t). \qquad (15)$$

Alternatively, contracted postprocessing delay operators can be decomposed into atomic delay operations

$$D_{i_1 \cdots i_N} \rightarrow D_{i_1 i_2} \cdots D_{i_{N-1} i_N}, \qquad (16)$$

which are applied from right to left to the time series. Note that the expressions on the left and right of the arrow are not

equivalent and will produce a different overall interpolation residual in TDI (see Sec. IV). Whether a single contracted delay operator or its decomposition into atomic delays performs better depends on the numerical values of the delays. As a general rule, delaying a time series by a small fractional delay is favorable, since the interpolation residual vanishes in the limit of an integer delay.

### A. On-board processing residuals

When computing TDI laser noise residuals, expressions of the form

$$[A, B]\phi_{ij}(t) = AB\phi_{ij}(t) - BA\phi_{ij}(t) \qquad (17)$$

appear, with operators A and B that act on the time series $\phi_{ij}(t)$. We call two operators noncommutative if Eq. (17) is nonvanishing. As an example, the fundamental limit for laser noise suppression in TDI can be described by a commutator of time-dependent delay operators, see e.g. [6]. This is discussed in more detail in Sec. IV.

Let us now outline how laser noise residuals due to filtering and sampling enter the interspacecraft interferometer via additional commutators.

The basic building block of every TDI combination is the set of intermediary variables $\eta_{ij}$. In the idealized case for which none of the data processing steps depicted in Fig. 1 are considered (i.e., $\mathbf{F} = \mathbf{S} = 1$ and $\hat{D}_{ij} = \mathbf{D}_{ij}$) the variable $\eta_{ij}$ simplifies to the difference between a local laser phase $\phi_i$ and a distant laser phase $\phi_j$ delayed by the light travel time between spacecraft i and j. Formally, one has

$$\eta_{ij}(t) = \mathbf{D}_{ij}\phi_j(t) - \phi_i(t). \qquad (18)$$

Here, we use the convention $\phi_1 = \phi_{12}$, $\phi_2 = \phi_{23}$ and $\phi_3 = \phi_{31}$. Under these conditions, i.e., with the specific algebraic form of Eq. (18), the fundamental laser noise limit is the usual delay commutator [see Eq. (33)].

If we instead insert the interferometric measurements of Eq. (1) into the definition of $\eta_{ij}$ given in Eq. (13a), the resulting expression cannot be recast into the algebraic form of Eq. (18). This is because of the order in which the filtering, decimation, and delay operations arise in Eq. (1a). If one introduces the following commutator into Eq. (1),

$$[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}(t) = (\mathbf{SF})\mathbf{D}_{ij}\phi_{ji}(t) - \mathbf{D}_{ij}(\mathbf{SF})\phi_{ji}(t), \qquad (19)$$

the delay operator $\mathbf{D}_{ij}$ switches places with the decimation stage SF such that Eq. (1a) becomes

$$\mathrm{isi}_{ij}(t) = \mathbf{D}_{ij}\mathbf{SF}\phi_{ji}(t) - \mathbf{SF}\phi_{ij}(t) + [\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}(t). \qquad (20)$$

This expression has the same algebraic form as Eq. (18) if we make the substitution $\phi_{ij} \rightarrow \mathbf{SF}\phi_{ij}$, with the exception that the interspacecraft interferometric measurement now

<footnote n="3">³The interpolation error arising in this operation might also be relevant as it contributes to the ranging error. To suppress it below the modulation noise level we use an appropriate interpolation orders (order 5 seems sufficient) in our numerical studies.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Illustration of intermediary variables defined in Eq. (13) and locking configuration N1-12. The intermediary variables are depicted as an arrow representing the synthesized photon path of the long arm. They are incrementally built up from the previously defined ones such that the last variable σ is represented by the entire path. Additionally, the chain of locked lasers for the locking configuration N1-12 is shown (as discussed in Sec. IV B). The primary laser is highlighted by the gray box.</figcaption></figure>

Here, $\hat{D}$ denotes the postprocessing operator that also considers ranging errors as indicated by the hat.

Figure 2 provides an illustration of the intermediary variables as they synthesize two-beam interferometers with a long and a short arm. The long arm is depicted as an arrow propagating around the LISA constellation. Using the variable σ we can finally express the second-generation Michelson combination $X_2$ as

$$X_2 = \sigma_{13} - \sigma_{12}. \qquad (14)$$

The contracted postprocessing delay operators $D_{i_1 \cdots i_N}$ appearing in Eq. (13) are applied in two steps. First, the actual nested delay $d_{i_1 \cdots i_N}(t)$ is calculated, and then, the fractional delay filter is applied. The nested delay can be calculated using the recursive operation<fnref n="3" />

$$d_{i_1 \cdots i_N}(t) = d_{i_1 i_2}(t) + D_{i_1 i_2} d_{i_2 \cdots i_N}(t). \qquad (15)$$

Alternatively, contracted postprocessing delay operators can be decomposed into atomic delay operations

$$D_{i_1 \cdots i_N} \rightarrow D_{i_1 i_2} \cdots D_{i_{N-1} i_N}, \qquad (16)$$

which are applied from right to left to the time series. Note that the expressions on the left and right of the arrow are not

equivalent and will produce a different overall interpolation residual in TDI (see Sec. IV). Whether a single contracted delay operator or its decomposition into atomic delays performs better depends on the numerical values of the delays. As a general rule, delaying a time series by a small fractional delay is favorable, since the interpolation residual vanishes in the limit of an integer delay.

### A. On-board processing residuals

When computing TDI laser noise residuals, expressions of the form

$$[A, B]\phi_{ij}(t) = AB\phi_{ij}(t) - BA\phi_{ij}(t) \qquad (17)$$

appear, with operators A and B that act on the time series $\phi_{ij}(t)$. We call two operators noncommutative if Eq. (17) is nonvanishing. As an example, the fundamental limit for laser noise suppression in TDI can be described by a commutator of time-dependent delay operators, see e.g. [6]. This is discussed in more detail in Sec. IV.

Let us now outline how laser noise residuals due to filtering and sampling enter the interspacecraft interferometer via additional commutators.

The basic building block of every TDI combination is the set of intermediary variables $\eta_{ij}$. In the idealized case for which none of the data processing steps depicted in Fig. 1 are considered (i.e., $F = S = 1$ and $\hat{D}_{ij} = D_{ij}$) the variable $\eta_{ij}$ simplifies to the difference between a local laser phase $\phi_i$ and a distant laser phase $\phi_j$ delayed by the light travel time between spacecraft i and j. Formally, one has

$$\eta_{ij}(t) = D_{ij}\phi_j(t) - \phi_i(t). \qquad (18)$$

Here, we use the convention $\phi_1 = \phi_{12}$, $\phi_2 = \phi_{23}$ and $\phi_3 = \phi_{31}$. Under these conditions, i.e., with the specific algebraic form of Eq. (18), the fundamental laser noise limit is the usual delay commutator [see Eq. (33)].

If we instead insert the interferometric measurements of Eq. (1) into the definition of $\eta_{ij}$ given in Eq. (13a), the resulting expression cannot be recast into the algebraic form of Eq. (18). This is because of the order in which the filtering, decimation, and delay operations arise in Eq. (1a). If one introduces the following commutator into Eq. (1),

$$[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}(t) = (\mathbf{SF})\mathbf{D}_{ij}\phi_{ji}(t) - \mathbf{D}_{ij}(\mathbf{SF})\phi_{ji}(t), \qquad (19)$$

the delay operator $\mathbf{D}_{ij}$ switches places with the decimation stage SF such that Eq. (1a) becomes

$$\mathrm{isi}_{ij}(t) = \mathbf{D}_{ij}\mathbf{SF}\phi_{ji}(t) - \mathbf{SF}\phi_{ij}(t) + [\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}(t). \qquad (20)$$

This expression has the same algebraic form as Eq. (18) if we make the substitution $\phi_{ij} \rightarrow \mathbf{SF}\phi_{ij}$, with the exception that the interspacecraft interferometric measurement now

<footnote n="3">³The interpolation error arising in this operation might also be relevant as it contributes to the ranging error. To suppress it below the modulation noise level we use an appropriate interpolation orders (order 5 seems sufficient) in our numerical studies.</footnote>
also contains the commutator $[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}$. This commutator comes as an additional residual in the final TDI expressions and enters via the usual TDI transfer function (this is reminiscent of the way readout noise enters the measurements, see [11]).

Using common commutator rules we can further split $[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}$ into a filter-delay commutator (flexing-filtering coupling) and a decimation-delay commutator

$$[\mathbf{SF}, \mathbf{D}]\phi(t) = \mathbf{S}[\mathbf{F}, \mathbf{D}]\phi(t) + [\mathbf{S}, \mathbf{D}]\mathbf{F}\phi(t), \quad (21)$$

and compute explicit analytic expressions for each contribution, see appendices B 1 and B 2.

## B. TDI residuals

The working principle of TDI is to time-shift the recorded beatnote phase measurements and linearly combine them in order to reduce laser noise in the resulting combinations. To achieve this we require an interpolation method and estimates of the aforementioned time shifts. In the following we denote a delay operation that is performed on-ground as $\hat{D}$. This operator acts on discrete time series and is therefore only a numerical approximation. Furthermore, time-shifting is performed with imperfect knowledge of the delay. Interpolation of the data and the error in ranging produce a residual with respect to the true propagation delay $\mathbf{D}$. To distinguish between these effects, we write the residual due to the imperfections of $\hat{D}$ as

$$(\hat{D} - \mathbf{D})\phi(t) = (\hat{D} - \hat{\mathbf{D}})\phi(t) + (\hat{\mathbf{D}} - \mathbf{D})\phi(t), \quad (22)$$

where the first term on the right-hand side represents the interpolation residual and the second term the residual stemming from ranging errors. We compute explicit expressions for each of these contributions in appendices C 1 and C 2, respectively.

## IV. RESIDUALS IN THE MICHELSON COMBINATIONS

Before considering either six independent lasers or the locking scheme N1-12 in Secs. IVA and IVB, let us first write down the general expression for the residuals in $X_2$. To do so, we use the factorization given in Eqs. (13) and (14) together with the expression of the reference interferometers in Eq. (1b) and that of the interspacecraft interferometers in Eq. (20). As the couplings discussed in Secs. III A and III B apply only to either the beatnote frequency offsets or the beatnote phase fluctuations, we split up the interferometers into those contributions. For brevity, we define $a_{ij} = \text{isi}_{ij}^o$ and $b_i = \text{rfi}_{ij}^o = -\text{rfi}_{ji}^o$ as in [13]. The beatnote phase fluctuations are given by

$$\text{isi}_{ij}^o = \mathbf{D}_{ij}\mathbf{SF}p_{ji} - \mathbf{SF}p_{ij} + [\mathbf{SF}, \mathbf{D}_{ij}]p_{ji}, \quad (23a)$$

$$\text{rfi}_{ij}^o = \mathbf{SF}p_{ik} - \mathbf{SF}p_{ij}. \quad (23b)$$

Here, $p_{ij}$ denotes the phase fluctuations of laser $ij$ that account for the natural phase noise of the laser and contributions that enter via laser locking (cf. [24] and Sec. IV B). In the following we write as a short-hand notation $\bar{p}_{ij} = \mathbf{SF}p_{ij}$.

Next, we replace each occurrence of the postprocessing delay $\hat{D}_{ij}$ by the relation given in Eq. (22). Here, we neglect all term that are second-order in the residuals (e.g., flexing-filtering coupling that enters via interpolation errors). Doing so, we find that $X_2$ breaks down into the following residuals

$$X_2 = \delta X_2^{[\mathbf{SF}, \mathbf{D}]} + \delta X_2^{\mathbf{D}} + \delta X_2^{\hat{\mathbf{D}}} + \delta X_2. \quad (24)$$

The first three terms on the right-hand side were described in Secs. III A and III B and the fourth term in Eq. (24) is the usual delay commutator describing the arm length mismatch in the virtual interferometer. This is discussed later in this section.

In the following we discuss each constituent of Eq. (24) individually. To simplify the expressions we assume equal arms and denote the multiplicative factor that is common to most residuals as

$$\mathbf{C} \equiv (1 - \mathbf{D}^4)(1 - \mathbf{D}^2), \quad (25)$$

$$|\tilde{\mathbf{C}}|^2 \equiv 16\sin^2(2\pi f d)\sin^2(4\pi f d). \quad (26)$$

The square root of $|\tilde{\mathbf{C}}|^2$ is commonly known as the “TDI transfer function” (or more precisely its magnitude).

The commutator residuals induced by filtering and decimation enter in the interspacecraft interferometer [cf., Eq. (20)] and are thus only propagated through TDI. This is consistent with [15], and they read

$$\delta X_2^{[\mathbf{SF}, \mathbf{D}]} = \mathbf{C}\{[\mathbf{SF}, \mathbf{D}_{13}]p_{31} + \mathbf{D}[\mathbf{SF}, \mathbf{D}_{31}]p_{13} \\ - [\mathbf{SF}, \mathbf{D}_{12}]p_{21} - \mathbf{D}[\mathbf{SF}, \mathbf{D}_{21}]p_{12}\}. \quad (27)$$

The interpolation residual depends on the factorization scheme used to compute the TDI variables. Using the factorization from Eq. (13) with “contracted delays,” we obtain

$$\delta X_2^{\mathbf{D}} = \mathbf{CD}[\mathbf{\Delta}_{12}\bar{p}_{21} - \mathbf{D}\mathbf{\Delta}_{31}\bar{p}_{13} - \mathbf{\Delta}_{13}\bar{p}_{31} \\ + \mathbf{D}(\mathbf{\Delta}_{13} + \mathbf{\Delta}_{31} - \mathbf{\Delta}_{12} + \mathbf{\Delta}_{121} - \mathbf{\Delta}_{131})\bar{p}_{12}]. \quad (28)$$

If one instead performs the interpolation with “atomic delays” (i.e., turning $\hat{D}_{i_1 \cdots i_N}$ into $\hat{D}_{i_1 i_2 \cdots \hat{D}_{i_{N-1} i_N}}$), one obtains

<footnote n="4">⁴This relation holds for $ij = 12, 23, 31$ as adjacent reference interferometers track the same beatnote up to a sign.</footnote>
also contains the commutator $[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}$. This commutator comes as an additional residual in the final TDI expressions and enters via the usual TDI transfer function (this is reminiscent of the way readout noise enters the measurements, see [11]).

Using common commutator rules we can further split $[\mathbf{SF}, \mathbf{D}_{ij}]\phi_{ji}$ into a filter-delay commutator (flexing-filtering coupling) and a decimation-delay commutator

$$[\mathbf{SF}, \mathbf{D}]\phi(t) = \mathbf{S}[\mathbf{F}, \mathbf{D}]\phi(t) + [\mathbf{S}, \mathbf{D}]\mathbf{F}\phi(t), \quad (21)$$

and compute explicit analytic expressions for each contribution, see appendices B 1 and B 2.

## B. TDI residuals

The working principle of TDI is to time-shift the recorded beatnote phase measurements and linearly combine them in order to reduce laser noise in the resulting combinations. To achieve this we require an interpolation method and estimates of the aforementioned time shifts. In the following we denote a delay operation that is performed on-ground as $\hat{D}$. This operator acts on discrete time series and is therefore only a numerical approximation. Furthermore, time-shifting is performed with imperfect knowledge of the delay. Interpolation of the data and the error in ranging produce a residual with respect to the true propagation delay $\mathbf{D}$. To distinguish between these effects, we write the residual due to the imperfections of $\hat{D}$ as

$$(\hat{D} - \mathbf{D})\phi(t) = (\hat{D} - \hat{\mathbf{D}})\phi(t) + (\hat{\mathbf{D}} - \mathbf{D})\phi(t), \quad (22)$$

where the first term on the right-hand side represents the interpolation residual and the second term the residual stemming from ranging errors. We compute explicit expressions for each of these contributions in appendices C 1 and C 2, respectively.

## IV. RESIDUALS IN THE MICHELSON COMBINATIONS

Before considering either six independent lasers or the locking scheme N1-12 in Secs. IVA and IVB, let us first write down the general expression for the residuals in $X_2$. To do so, we use the factorization given in Eqs. (13) and (14) together with the expression of the reference interferometers in Eq. (1b) and that of the interspacecraft interferometers in Eq. (20). As the couplings discussed in Secs. III A and III B apply only to either the beatnote frequency offsets or the beatnote phase fluctuations, we split up the interferometers into those contributions. For brevity, we define $a_{ij} = \text{isi}_{ij}^o$ and $b_i = \text{rfi}_{ij}^o = -\text{rfi}_{ji}^o$ as in [13]. The beatnote phase fluctuations are given by

$$\text{isi}_{ij}^o = \mathbf{D}_{ij}\mathbf{SF}p_{ji} - \mathbf{SF}p_{ij} + [\mathbf{SF}, \mathbf{D}_{ij}]p_{ji}, \quad (23a)$$

$$\text{rfi}_{ij}^o = \mathbf{SF}p_{ik} - \mathbf{SF}p_{ij}. \quad (23b)$$

Here, $p_{ij}$ denotes the phase fluctuations of laser $ij$ that account for the natural phase noise of the laser and contributions that enter via laser locking (cf. [24] and Sec. IV B). In the following we write as a short-hand notation $\bar{p}_{ij} = \mathbf{SF}p_{ij}$.

Next, we replace each occurrence of the postprocessing delay $\hat{D}_{ij}$ by the relation given in Eq. (22). Here, we neglect all term that are second-order in the residuals (e.g., flexing-filtering coupling that enters via interpolation errors). Doing so, we find that $X_2$ breaks down into the following residuals

$$X_2 = \delta X_2^{[\mathbf{SF}, \mathbf{D}]} + \delta X_2^{\mathbf{D}} + \delta X_2^{\hat{\mathbf{D}}} + \delta X_2. \quad (24)$$

The first three terms on the right-hand side were described in Secs. III A and III B and the fourth term in Eq. (24) is the usual delay commutator describing the arm length mismatch in the virtual interferometer. This is discussed later in this section.

In the following we discuss each constituent of Eq. (24) individually. To simplify the expressions we assume equal arms and denote the multiplicative factor that is common to most residuals as

$$\mathbf{C} \equiv (1 - \mathbf{D}^4)(1 - \mathbf{D}^2), \quad (25)$$

$$|\tilde{\mathbf{C}}|^2 \equiv 16\sin^2(2\pi f d)\sin^2(4\pi f d). \quad (26)$$

The square root of $|\tilde{\mathbf{C}}|^2$ is commonly known as the “TDI transfer function” (or more precisely its magnitude).

The commutator residuals induced by filtering and decimation enter in the interspacecraft interferometer [cf., Eq. (20)] and are thus only propagated through TDI. This is consistent with [15], and they read

$$\delta X_2^{[\mathbf{SF}, \mathbf{D}]} = \mathbf{C}\{[\mathbf{SF}, \mathbf{D}_{13}]p_{31} + \mathbf{D}[\mathbf{SF}, \mathbf{D}_{31}]p_{13} - [\mathbf{SF}, \mathbf{D}_{12}]p_{21} - \mathbf{D}[\mathbf{SF}, \mathbf{D}_{21}]p_{12}\}. \quad (27)$$

The interpolation residual depends on the factorization scheme used to compute the TDI variables. Using the factorization from Eq. (13) with “contracted delays,” we obtain

$$\delta X_2^{\mathbf{D}} = \mathbf{CD}[\mathbf{\Delta}_{12}\bar{p}_{21} - \mathbf{D}\mathbf{\Delta}_{31}\bar{p}_{13} - \mathbf{\Delta}_{13}\bar{p}_{31} + \mathbf{D}(\mathbf{\Delta}_{13} + \mathbf{\Delta}_{31} - \mathbf{\Delta}_{12} + \mathbf{\Delta}_{121} - \mathbf{\Delta}_{131})\bar{p}_{12}]. \quad (28)$$

If one instead performs the interpolation with “atomic delays” (i.e., turning $\hat{D}_{i_1 \cdots i_N}$ into $\hat{D}_{i_1 i_2 \cdots \hat{D}_{i_{N-1} i_N}}$), one obtains

<footnote n="4">⁴This relation holds for $ij = 12, 23, 31$ as adjacent reference interferometers track the same beatnote up to a sign.</footnote>
$$ \delta X_2^D = CD\{\Delta_{12}\bar{p}_{21} + D\Delta_{21}\bar{p}_{12} - D\Delta_{31}\bar{p}_{13} - \Delta_{13}\bar{p}_{31}\}. \qquad (29) $$

In the expressions above, $\Delta_{ij}$ denotes the effect of imperfect interpolation when time-shifting by the arm $ij$. This is discussed in more detail in Appendix C 1.

To compute the TDI residual induced by the ranging error we need a description of the ranging measurement and any additional processing employed to reduce this ranging error. An in-depth study of all ranging observables and how to best combine them is given in [26]. In our analysis we consider the onboard ranging methods that are realized by two additional modulations on the laser beams: the pseudo-random noise (PRN) modulation (absolute ranging) and the clock sidebands. The range estimate described in Appendix D, which is mostly adopted from [13], combines both measurements. The estimate, denoted $\hat{d}_{ij}$, inherits a bias $B_{ij}$ from the PRN measurements and a stochastic term from the sideband measurements. The resulting deviation from the true range $d_{ij}$, the ranging error, reads<fnref n="5" />

$$ r_{ij}(t) = B_{ij} - D_{ij}M_j(t) + M_i(t), \qquad (30) $$

where $M_i$ denotes the modulation noise on the left-handed MOSA on satellite $i$. The resulting residual $\delta X_2^D$ can be decomposed into the component originating from the ranging bias $\delta X_2^B$ and the modulation noise $\delta X_2^M$. These two contributions are consistent with [13] and read

$$ \delta X_2^B = C\{\mathbf{D}^2 B_{31} \dot{\hat{p}}_{13} + \mathbf{D} B_{13} \dot{\hat{p}}_{31} - \mathbf{D}^2 B_{21} \dot{\hat{p}}_{12} - \mathbf{D} B_{12} \dot{\hat{p}}_{21}\}, \qquad (31a) $$

$$ \delta X_2^M = C\{(a_{13} - a_{12} + (1 - \mathbf{D}^2)b_1)M_1 - \mathbf{D}a_{21}M_2 + \mathbf{D}a_{31}M_3\}, \qquad (31b) $$

where $\dot{\hat{p}}_{ij}$ is the time derivative of the filtered and decimated laser phase fluctuations. Note that we also present the subdominant coupling of the stochastic in-band ranging error to laser frequency noise in Appendix F.

Let us also give the expression for the fourth contribution in Eq. (24) corresponding to the usual TDI delay commutator that fundamentally limits laser noise reduction. Due to the flexing of the LISA constellation the round-trip times of any synthesized two-beam interferometer are not exactly identical. Hence, laser noise in the two beams does not cancel but enters into the TDI combination with a prefactor proportional to the arm length mismatch $\Delta d$. For the second-generation Michelson variable $X_2$, the residual reads

$$ \delta X_2 = [\mathbf{D}_{13121}, \mathbf{D}_{12131}] \mathbf{SF} \phi_{12} \qquad (32) $$

$$ \approx - \Delta d_{X_2} \mathbf{D}_{131212131} \frac{d}{dt}(\mathbf{SF} \phi_{12}). \qquad (33) $$

In the second line we use the property that a delay commutator acts like a derivative, as already described in [4,6,13,27]. The arm length mismatch in a second-generation Michelson interferometer is given by

$$ \Delta d_{X_2}(t) = d_{131212131}(t) - d_{121313121}(t). \qquad (34) $$

Equation (33) can be further split up into a deterministic out-of-band drift and an in-band component by plugging in Eq. (4) in the above expression. For the deterministic part, we recover<fnref n="6" />

$$ \delta X_2^o = -\Delta d_{X_2} \mathbf{D}_{131212131} \nu_{12}^o, \qquad (35) $$

where the travel time difference can be efficiently computed from the time delays $d_{ij}$ as described in Appendix E. Additionally, a model for the absolute laser frequency $\nu_{12}^o$ must be provided to subtract the trend. The stochastic in-band component of the delay commutator is characterized by its PSD

$$ S_{\delta X_2}(f) = (\Delta d_{X_2})^2 (2\pi f)^2 |\tilde{\mathbf{S}}|^2 |\tilde{\mathbf{F}}|^2 S_p(f). \qquad (36) $$

For this derivation we have assumed that $\Delta d_{X_2}$ is constant. In reality, the amplitude of residual laser noise is modulated due to orbital dynamics governing the motion of the three spacecraft. We expand the travel time difference $\Delta d_{X_2}$ up to second order in velocity $\dot{d}_{iji}$ and up to first order in acceleration $\ddot{d}_{iji}$ to derive a good approximation of $\Delta d_{X_2}$ [6],

$$ \Delta d_{X_2} = (d_{131} \dot{d}_{121} - d_{121} \dot{d}_{131})(\dot{d}_{121} + \dot{d}_{131}) - (d_{131} \ddot{d}_{121} - d_{121} \ddot{d}_{131})(d_{121} + d_{131}), \qquad (37) $$

which is of the order $10^{-12}$ s [13].

In previous studies [e.g., [11]], it was demonstrated that laser locking does not affect the coupling of path-length noises in TDI combinations. At first this seems counter-intuitive as locked lasers generate echoes of any path-length noise imprinted on the reference beam. However, those echoes are canceled out in TDI since any in-band component in all six laser phases is, by construction, strongly suppressed by the algorithm.

The couplings described above introduce residual laser phase in the TDI combinations. Thus, laser noise and path-length noises imprinted on the laser will enter the

<footnote n="5">⁵In the notation of Appendix C 2, where the PSD of the ranging residual commutator is derived, $B_{ij}$ is written as a generic bias, $r^o$, while the combination $-D_{ij}M_j(t) + M_i(t)$ is written as a generic stochastic term, $r^c(t)$.</footnote>

<footnote n="6">⁶Because the phase ramp produced by the THz central laser frequency is unaffected by filtering and decimation we can drop both **F** and **S** assuming an acausal filter implementation with vanishing group delay (cf. [15]).</footnote>
$$ \delta X_2^D = CD\{\Delta_{12}\bar{p}_{21} + D\Delta_{21}\bar{p}_{12} - D\Delta_{31}\bar{p}_{13} - \Delta_{13}\bar{p}_{31}\}. \qquad (29) $$

In the expressions above, $\Delta_{ij}$ denotes the effect of imperfect interpolation when time-shifting by the arm $ij$. This is discussed in more detail in Appendix C 1.

To compute the TDI residual induced by the ranging error we need a description of the ranging measurement and any additional processing employed to reduce this ranging error. An in-depth study of all ranging observables and how to best combine them is given in [26]. In our analysis we consider the onboard ranging methods that are realized by two additional modulations on the laser beams: the pseudo-random noise (PRN) modulation (absolute ranging) and the clock sidebands. The range estimate described in Appendix D, which is mostly adopted from [13], combines both measurements. The estimate, denoted $\hat{d}_{ij}$, inherits a bias $B_{ij}$ from the PRN measurements and a stochastic term from the sideband measurements. The resulting deviation from the true range $d_{ij}$, the ranging error, reads<fnref n="5" />

$$ r_{ij}(t) = B_{ij} - D_{ij}M_j(t) + M_i(t), \qquad (30) $$

where $M_i$ denotes the modulation noise on the left-handed MOSA on satellite $i$. The resulting residual $\delta X_2^D$ can be decomposed into the component originating from the ranging bias $\delta X_2^B$ and the modulation noise $\delta X_2^M$. These two contributions are consistent with [13] and read

$$ \delta X_2^B = C\{\mathbf{D}^2 B_{31} \dot{\hat{p}}_{13} + \mathbf{D} B_{13} \dot{\hat{p}}_{31} - \mathbf{D}^2 B_{21} \dot{\hat{p}}_{12} - \mathbf{D} B_{12} \dot{\hat{p}}_{21}\}, \qquad (31a) $$

$$ \delta X_2^M = C\{(a_{13} - a_{12} + (1 - \mathbf{D}^2)b_1)M_1 - \mathbf{D}a_{21}M_2 + \mathbf{D}a_{31}M_3\}, \qquad (31b) $$

where $\dot{\hat{p}}_{ij}$ is the time derivative of the filtered and decimated laser phase fluctuations. Note that we also present the subdominant coupling of the stochastic in-band ranging error to laser frequency noise in Appendix F.

Let us also give the expression for the fourth contribution in Eq. (24) corresponding to the usual TDI delay commutator that fundamentally limits laser noise reduction. Due to the flexing of the LISA constellation the round-trip times of any synthesized two-beam interferometer are not exactly identical. Hence, laser noise in the two beams does not cancel but enters into the TDI combination with a prefactor proportional to the arm length mismatch $\Delta d$. For the second-generation Michelson variable $X_2$, the residual reads

$$ \delta X_2 = [\mathbf{D}_{13121}, \mathbf{D}_{12131}] \mathbf{SF} \phi_{12} \qquad (32) $$

$$ \approx - \Delta d_{X_2} \mathbf{D}_{131212131} \frac{d}{dt}(\mathbf{SF} \phi_{12}). \qquad (33) $$

In the second line we use the property that a delay commutator acts like a derivative, as already described in [4,6,13,27]. The arm length mismatch in a second-generation Michelson interferometer is given by

$$ \Delta d_{X_2}(t) = d_{131212131}(t) - d_{121313121}(t). \qquad (34) $$

Equation (33) can be further split up into a deterministic out-of-band drift and an in-band component by plugging in Eq. (4) in the above expression. For the deterministic part, we recover<fnref n="6" />

$$ \delta X_2^o = -\Delta d_{X_2} \mathbf{D}_{131212131} \nu_{12}^o, \qquad (35) $$

where the travel time difference can be efficiently computed from the time delays $d_{ij}$ as described in Appendix E. Additionally, a model for the absolute laser frequency $\nu_{12}^o$ must be provided to subtract the trend. The stochastic in-band component of the delay commutator is characterized by its PSD

$$ S_{\delta X_2}(f) = (\Delta d_{X_2})^2 (2\pi f)^2 |\tilde{\mathbf{S}}|^2 |\tilde{\mathbf{F}}|^2 S_p(f). \qquad (36) $$

For this derivation we have assumed that $\Delta d_{X_2}$ is constant. In reality, the amplitude of residual laser noise is modulated due to orbital dynamics governing the motion of the three spacecraft. We expand the travel time difference $\Delta d_{X_2}$ up to second order in velocity $\dot{d}_{iji}$ and up to first order in acceleration $\ddot{d}_{iji}$ to derive a good approximation of $\Delta d_{X_2}$ [6],

$$ \Delta d_{X_2} = (d_{131} \dot{d}_{121} - d_{121} \dot{d}_{131})(\dot{d}_{121} + \dot{d}_{131}) - (d_{131} \ddot{d}_{121} - d_{121} \ddot{d}_{131})(d_{121} + d_{131}), \qquad (37) $$

which is of the order $10^{-12}$ s [13].

In previous studies [e.g., [11]], it was demonstrated that laser locking does not affect the coupling of path-length noises in TDI combinations. At first this seems counter-intuitive as locked lasers generate echoes of any path-length noise imprinted on the reference beam. However, those echoes are canceled out in TDI since any in-band component in all six laser phases is, by construction, strongly suppressed by the algorithm.

The couplings described above introduce residual laser phase in the TDI combinations. Thus, laser noise and path-length noises imprinted on the laser will enter the

<footnote n="5">⁵In the notation of Appendix C 2, where the PSD of the ranging residual commutator is derived, $B_{ij}$ is written as a generic bias, $r^o$, while the combination $-D_{ij}M_j(t) + M_i(t)$ is written as a generic stochastic term, $r^c(t)$.</footnote>

<footnote n="6">⁶Because the phase ramp produced by the THz central laser frequency is unaffected by filtering and decimation we can drop both **F** and **S** assuming an acausal filter implementation with vanishing group delay (cf. [15]).</footnote>
combination. However, the effect of the latter is subdominant, because laser noise dominates the residuals. The impact of locking is thus only relevant for laser noise. In the following sections we first describe the residuals assuming six independent lasers, i.e., when each laser is locked to an individual cavity; then, we derive the same residuals assuming the standard locking configuration N1-12⁷ introduced in [23].

### A. Six independent lasers

In the first step let us consider six lasers stabilized to their individual cavities. This is a theoretical setup, which cannot be used in practice, but worth studying for pedagogical reasons.

Using Eq. (B2) and the upper bound given in Eq. (B6), it is easy to compute the PSD of Eq. (27) corresponding to the coupling of the filter-delay and decimation-delay commutator. It reads

$$S_{\delta X_2}^{\mathbf{S}[\mathbf{F}, \mathbf{D}]}(f) = 4|\tilde{\mathbf{C}}|^2|\tilde{\mathbf{S}}|^2 \left( \tilde{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\hat{p}}(f) \right), \quad (38a)$$

$$S_{\delta X_2}^{[\mathbf{S}, \mathbf{D}] \mathbf{F}}(f) \le 16 |\tilde{\mathbf{C}}|^2 \cdot \sum_{n=1}^{\infty} (|\tilde{\mathbf{F}}|^2 S_p)^{(n)}(f), \quad (38b)$$

where

$$\tilde{d}^2 = \frac{\dot{d}_{12}^2 + \dot{d}_{21}^2 + \dot{d}_{13}^2 + \dot{d}_{31}^2}{4}. \quad (39)$$

Next, we investigate the PSD of the interpolation residual contribution. For six independent lasers we choose to use atomic delays over contracted delays because it results in much simpler couplings. Furthermore, we assume the worst-case interpolation error in order to derive an upper bound on its contribution. We use Eqs. (C4) and (29) and find

$$S_{\delta X_2}^{\mathbf{D}}(f) \le 4|\tilde{\mathbf{C}}|^2|\tilde{\mathbf{\Delta}}|^2|\tilde{\mathbf{S}}|^2|\tilde{\mathbf{F}}|^2 S_p(f). \quad (40)$$

Here, $\tilde{\mathbf{\Delta}}$ (without indices) represents the worst case interpolation error coupling.

Finally, using Eqs. (C8) and (31), the contribution from the ranging error yields a PSD equal to [13]

$$S_{\delta X_2}^{\hat{\mathbf{D}}}(f) = |\tilde{\mathbf{C}}|^2 (4\bar{B}^2(2\pi f)^2 |\tilde{\mathbf{S}}|^2 |\tilde{\mathbf{F}}|^2 S_p(f) + A_{X_2}^M(f)S_M(f)), \quad (41)$$

with an effective squared bias $\bar{B}^2$ and a modulating function $A_{X_2}^M(f)$ defined as

$$\bar{B}^2 = \frac{B_{12}^2 + B_{21}^2 + B_{13}^2 + B_{31}^2}{4}, \quad (42)$$

$$A_{X_2}^M(f) = (a_{12} - a_{13})^2 + a_{21}^2 + a_{31}^2 - 4b_1(a_{12} - a_{13} - b_1)\sin^2(2\pi fd). \quad (43)$$

### B. Locked lasers

The baseline design of LISA foresees locked lasers to ensure that all beatnote frequencies fall into the sensitive bandwidth of the photoreceivers on-board the spacecraft (5 MHz to 25 MHz). To achieve this the primary laser is stabilized using a cavity that serves as a frequency reference. The five remaining lasers are frequency offset locked in succession to the primary following a locking topology. Laser locking of one laser is achieved by adjusting the frequency of this laser source so that the beatnote frequency of the locking interferometer follows a predetermined time-varying offset frequency $o_{ij}(t)$. The so-called locking conditions for the interspacecraft and reference locking interferometer are given by

$$\dot{i}\ddot{s}_{ij} = \mathbf{D}_{ij}\nu_{ji} - \nu_{ij} = o_{ij}(t), \quad (44a)$$

$$\mathbf{r}\ddot{\mathbf{f}}_{ij} = \nu_{ik} - \nu_{ij} = o_{ij}(t), \quad (44b)$$

where $\nu_{ji}$ and $\nu_{ik}$, respectively, denote the frequencies of the reference lasers and $\nu_{ij}$ is the frequency of the laser that is controlled. As in [24], we assume the laser lock to be perfect with infinite gain, such that the locked beatnotes follow the programmed offset frequencies $o_{ij}(t)$ exactly. As the later only have out-of-band components, any locked laser is simply “echoing” the incoming phase noise of the reference laser. Therefore, laser noise becomes correlated among the six lasers.

In this section, we take as an example the N1-12 locking configuration as depicted in Fig. 2. Here, N1 specifies the locking topology and 12 the index of the primary laser (see [23,24] for an overview of the locking topologies). For this particular locking configuration the in-band phase noise of the six lasers is given by

$$\begin{align} p_{12} &= p, & p_{23} &= \mathbf{D}_{21}p, & p_{31} &= \mathbf{D}_{31}p, \tag{45} \\ p_{13} &= p, & p_{32} &= \mathbf{D}_{31}p, & p_{21} &= \mathbf{D}_{21}p. \end{align}$$

For better readability we drop the index on $p = p_{12}$ denoting the in-band phase fluctuations of the primary laser. We proceed by inserting the expressions of Eq. (45) into the general expressions for laser noise related residuals listed in Sec. IV. To simplify the expressions we make use of the commutator rule

$$[\mathbf{S}\mathbf{F}, \mathbf{D}_{ij}]\mathbf{D}_{jk} = [\mathbf{S}\mathbf{F}, \mathbf{D}_{ijk}] - \mathbf{D}_{ij}[\mathbf{S}\mathbf{F}, \mathbf{D}_{jk}]. \quad (46)$$

We note that the following results are very particular to the choice of locking configuration. Moreover, results for $X_2$, $Y_2$ and $Z_2$ no longer exhibit rotational symmetry in the

<footnote n="7">⁷Note that in this work, we shall use the locking configuration naming convention introduced in [24].</footnote>
combination. However, the effect of the latter is subdominant, because laser noise dominates the residuals. The impact of locking is thus only relevant for laser noise. In the following sections we first describe the residuals assuming six independent lasers, i.e., when each laser is locked to an individual cavity; then, we derive the same residuals assuming the standard locking configuration N1-12⁷ introduced in [23].

## A. Six independent lasers

In the first step let us consider six lasers stabilized to their individual cavities. This is a theoretical setup, which cannot be used in practice, but worth studying for pedagogical reasons.

Using Eq. (B2) and the upper bound given in Eq. (B6), it is easy to compute the PSD of Eq. (27) corresponding to the coupling of the filter-delay and decimation-delay commutator. It reads

$$S_{\delta X_2}^{\mathbf{S}[\mathbf{F}, \mathbf{D}]}(f) = 4|\tilde{\mathbf{C}}|^2|\tilde{\mathbf{S}}|^2 \left( \tilde{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\hat{p}}(f) \right), \quad (38a)$$

$$S_{\delta X_2}^{[\mathbf{S}, \mathbf{D}]\mathbf{F}}(f) \le 16|\tilde{\mathbf{C}}|^2 \cdot \sum_{n=1}^{\infty} (|\tilde{\mathbf{F}}|^2 S_p)^{(n)}(f), \quad (38b)$$

where

$$\tilde{d}^2 = \frac{\dot{d}_{12}^2 + \dot{d}_{21}^2 + \dot{d}_{13}^2 + \dot{d}_{31}^2}{4}. \quad (39)$$

Next, we investigate the PSD of the interpolation residual contribution. For six independent lasers we choose to use atomic delays over contracted delays because it results in much simpler couplings. Furthermore, we assume the worst-case interpolation error in order to derive an upper bound on its contribution. We use Eqs. (C4) and (29) and find

$$S_{\delta X_2}^{\mathbf{D}}(f) \le 4|\tilde{\mathbf{C}}|^2|\tilde{\mathbf{\Delta}}|^2|\tilde{\mathbf{S}}|^2|\tilde{\mathbf{F}}|^2 S_p(f). \quad (40)$$

Here, $\tilde{\mathbf{\Delta}}$ (without indices) represents the worst case interpolation error coupling.

Finally, using Eqs. (C8) and (31), the contribution from the ranging error yields a PSD equal to [13]

$$S_{\delta X_2}^{\hat{\mathbf{D}}}(f) = |\tilde{\mathbf{C}}|^2 (4\bar{B}^2(2\pi f)^2 |\tilde{\mathbf{S}}|^2 |\tilde{\mathbf{F}}|^2 S_p(f) + A_{X_2}^M(f)S_M(f)), \quad (41)$$

with an effective squared bias $\bar{B}^2$ and a modulating function $A_{X_2}^M(f)$ defined as

$$\bar{B}^2 = \frac{B_{12}^2 + B_{21}^2 + B_{13}^2 + B_{31}^2}{4}, \quad (42)$$

$$A_{X_2}^M(f) = (a_{12} - a_{13})^2 + a_{21}^2 + a_{31}^2 - 4b_1(a_{12} - a_{13} - b_1)\sin^2(2\pi fd). \quad (43)$$

## B. Locked lasers

The baseline design of LISA foresees locked lasers to ensure that all beatnote frequencies fall into the sensitive bandwidth of the photoreceivers on-board the spacecraft (5 MHz to 25 MHz). To achieve this the primary laser is stabilized using a cavity that serves as a frequency reference. The five remaining lasers are frequency offset locked in succession to the primary following a locking topology. Laser locking of one laser is achieved by adjusting the frequency of this laser source so that the beatnote frequency of the locking interferometer follows a predetermined time-varying offset frequency $o_{ij}(t)$. The so-called locking conditions for the interspacecraft and reference locking interferometer are given by

$$\dot{i}\ddot{s}_{ij} = \mathbf{D}_{ij}\nu_{ji} - \nu_{ij} = o_{ij}(t), \quad (44a)$$

$$\mathbf{r}\ddot{\mathbf{f}}_{ij} = \nu_{ik} - \nu_{ij} = o_{ij}(t), \quad (44b)$$

where $\nu_{ji}$ and $\nu_{ik}$, respectively, denote the frequencies of the reference lasers and $\nu_{ij}$ is the frequency of the laser that is controlled. As in [24], we assume the laser lock to be perfect with infinite gain, such that the locked beatnotes follow the programmed offset frequencies $o_{ij}(t)$ exactly. As the later only have out-of-band components, any locked laser is simply “echoing” the incoming phase noise of the reference laser. Therefore, laser noise becomes correlated among the six lasers.

In this section, we take as an example the N1-12 locking configuration as depicted in Fig. 2. Here, N1 specifies the locking topology and 12 the index of the primary laser (see [23,24] for an overview of the locking topologies). For this particular locking configuration the in-band phase noise of the six lasers is given by

$$\begin{align} p_{12} &= p, & p_{23} &= \mathbf{D}_{21}p, & p_{31} &= \mathbf{D}_{31}p, \tag{45} \\ p_{13} &= p, & p_{32} &= \mathbf{D}_{31}p, & p_{21} &= \mathbf{D}_{21}p. \nonumber \end{align}$$

For better readability we drop the index on $p = p_{12}$ denoting the in-band phase fluctuations of the primary laser. We proceed by inserting the expressions of Eq. (45) into the general expressions for laser noise related residuals listed in Sec. IV. To simplify the expressions we make use of the commutator rule

$$[\mathbf{S}\mathbf{F}, \mathbf{D}_{ij}]\mathbf{D}_{jk} = [\mathbf{S}\mathbf{F}, \mathbf{D}_{ijk}] - \mathbf{D}_{ij}[\mathbf{S}\mathbf{F}, \mathbf{D}_{jk}]. \quad (46)$$

We note that the following results are very particular to the choice of locking configuration. Moreover, results for $X_2$, $Y_2$ and $Z_2$ no longer exhibit rotational symmetry in the

<footnote n="7">⁷Note that in this work, we shall use the locking configuration naming convention introduced in [24].</footnote>
indices as it is broken by laser locking; this is easily verified from Eq. (45).

As an example, we derive the laser noise residuals for the second-generation Michelson $X_2$ variable. In general, we find that expressions simplify. This is due to the fact that the locking configuration N1-12 inherently generates round-trip measurements required to build $X_2$.<fnref n="8" /> For $Y_2$ and $Z_2$ we find more complicated expressions involving more terms. However, the TDI transfer function remains and additional factors only modulate the residual slightly. Those are caused by cross-spectral densities among correlated laser noise residuals.

For the [SF, D] commutator we recover

$$ \delta X_2^{[SF,D]} = C([SF, D_{121}] - [SF, D_{131}]) p. \quad (47) $$

We split it up further [as in Eq. (21)] into contributions from the filter-delay commutator and sample-delay commutator, respectively. The PSD of the flexing-filtering coupling is given by Eq. (38) with an effective squared delay derivative equal to

$$ \bar{d}^2 = \frac{(\dot{d}_{121} - \dot{d}_{131})^2}{4}. \quad (48) $$

For the sample-delay commutator, the upper bound given in Eq. (38b) is still valid. Furthermore, we note that both residuals vanish nontrivially if the round-trip delays $d_{121}$ and $d_{131}$ become equal (for the flexing-filtering coupling, their derivatives must also be equal).

To simplify the coupling of the interpolation error for locked lasers we use contracted delays [see Eq. (28)]. We obtain

$$ \delta X_2^D = CD^2(\Delta_{131} - \Delta_{121})SFp. \quad (49) $$

Again, the upper bound given in Eq. (40) still holds and the residual vanishes for equal round-trip times $d_{121}$ and $d_{131}$.

Finally, we discuss the residual caused by ranging errors. The contribution coming from modulation noise in Eq. (31b) stays unchanged as the laser noise component of the beatnote is not involved. The contribution of the ranging biases in each arm, see Eq. (31a), is given by

$$ \delta X_2^B = CD^2(B_{12} + B_{21} - B_{13} - B_{31})SFp. \quad (50) $$

Therefore, the PSD of ranging error contributions is expressed as Eq. (41) with an effective squared bias equal to

$$ \bar{B}^2 = \frac{(B_{12} + B_{21} - B_{13} - B_{31})^2}{4}. \quad (51) $$

We remind that these cancellations appear only in the Michelson combination $X$ centered on the spacecraft of the

primary laser, while $Y$ and $Z$ gain more complicated couplings which do not show the same cancellations.

### C. Comparison with numerical simulations

We now compare the theoretical models described above to simulated LISA data, obtained using LISA INSTRUMENT [24,28] to generate the measurements and PYTDI [29] to calculate the second-generation Michelson variables. In order to compare, one-by-one, the couplings described above to the simulated data, we run three different sets of simulations that single out the effects of (i) onboard processing, (ii) interpolation errors and (iii) ranging errors individually. As a reference we also run a simulation, where all noises are disabled, which serves as an indicator of the numerical limit due to rounding errors in the floating-point variables that we use to represent the physical quantities in the simulation.

To put the resulting residuals in context we compare them to the one picometer reference curve, which represents the coupling of a GW strain with an equivalent displacement amplitude of 1 pm/√Hz on all six links with relaxation at low frequencies, and is commonly used as a benchmark in LISA. Converted to frequency units and propagated through TDI it reads

$$ S_{\tilde{X}_2}^{\text{ref}}(f) = 4|\tilde{C}|^2(2\pi f)^2 \left(1 + \left(\frac{2 \text{ MHz}}{f}\right)^4\right) \left(\frac{1 \text{ pm}/\sqrt{\text{Hz}}}{1064 \text{ nm}}\right)^2. \quad (52) $$

Here, the factor 4 accounts for the fact that there is a total of four links involved in the computation of the Michelson variables. The relaxation below 2 mHz mimics the shape of the LISA noise budget. Overall, the reference curve above is about a factor 10 below the total optical metrology noise allocation in LISA.

Let us now describe the nominal setup of the simulator and the processing performed on the measurements. All simulations span approximately 3 days in duration with the LISA constellation following realistic heliocentric orbits provided by ESA. The beatnote frequencies<fnref n="9" /> of the inter-spacecraft and the reference interferometers, which range between 5 MHz to 25 MHz, are simulated at a sampling rate of 16 Hz by LISA INSTRUMENT. They are formed by

<footnote n="9">⁹Units of frequency are used to circumvent numerical issues specific to phase units. A beatnote frequency of several MHz results in a rapidly increasing phase ramp. This would require a large number of significant digits to deal with its dynamic range (~10¹⁴ cycles after 1 yr at a required precision of at least µ cycles which results in 20 significant digits). On the other hand, its derivative, the beatnote frequency, stays roughly constant and a double precision float (approximately 16 significant digits) is sufficient to represent it. Here, the picometer reference curve with a knee at 2 mHz corresponds to a precision of approximately 20 nHz which gives 15 orders of magnitude when considering a 20 MHz beatnote frequency.</footnote>

<footnote n="8">⁸For example, the interspacecraft measurements $\text{isi}_{12} \sim \pi_{12}$ and $\text{isi}_{13} \sim \pi_{13}$.</footnote>
indices as it is broken by laser locking; this is easily verified from Eq. (45).

As an example, we derive the laser noise residuals for the second-generation Michelson $X_2$ variable. In general, we find that expressions simplify. This is due to the fact that the locking configuration N1-12 inherently generates round-trip measurements required to build $X_2$.<fnref n="8" /> For $Y_2$ and $Z_2$ we find more complicated expressions involving more terms. However, the TDI transfer function remains and additional factors only modulate the residual slightly. Those are caused by cross-spectral densities among correlated laser noise residuals.

For the [SF, D] commutator we recover

$$ \delta X_2^{[SF,D]} = C([SF, D_{121}] - [SF, D_{131}]) p. \quad (47) $$

We split it up further [as in Eq. (21)] into contributions from the filter-delay commutator and sample-delay commutator, respectively. The PSD of the flexing-filtering coupling is given by Eq. (38) with an effective squared delay derivative equal to

$$ \bar{d}^2 = \frac{(\dot{d}_{121} - \dot{d}_{131})^2}{4}. \quad (48) $$

For the sample-delay commutator, the upper bound given in Eq. (38b) is still valid. Furthermore, we note that both residuals vanish nontrivially if the round-trip delays $d_{121}$ and $d_{131}$ become equal (for the flexing-filtering coupling, their derivatives must also be equal).

To simplify the coupling of the interpolation error for locked lasers we use contracted delays [see Eq. (28)]. We obtain

$$ \delta X_2^D = CD^2(\Delta_{131} - \Delta_{121})SFp. \quad (49) $$

Again, the upper bound given in Eq. (40) still holds and the residual vanishes for equal round-trip times $d_{121}$ and $d_{131}$.

Finally, we discuss the residual caused by ranging errors. The contribution coming from modulation noise in Eq. (31b) stays unchanged as the laser noise component of the beatnote is not involved. The contribution of the ranging biases in each arm, see Eq. (31a), is given by

$$ \delta X_2^B = CD^2(B_{12} + B_{21} - B_{13} - B_{31})SFp. \quad (50) $$

Therefore, the PSD of ranging error contributions is expressed as Eq. (41) with an effective squared bias equal to

$$ \bar{B}^2 = \frac{(B_{12} + B_{21} - B_{13} - B_{31})^2}{4}. \quad (51) $$

We remind that these cancellations appear only in the Michelson combination $X$ centered on the spacecraft of the

primary laser, while $Y$ and $Z$ gain more complicated couplings which do not show the same cancellations.

### C. Comparison with numerical simulations

We now compare the theoretical models described above to simulated LISA data, obtained using LISA INSTRUMENT [24,28] to generate the measurements and PYTDI [29] to calculate the second-generation Michelson variables. In order to compare, one-by-one, the couplings described above to the simulated data, we run three different sets of simulations that single out the effects of (i) onboard processing, (ii) interpolation errors and (iii) ranging errors individually. As a reference we also run a simulation, where all noises are disabled, which serves as an indicator of the numerical limit due to rounding errors in the floating-point variables that we use to represent the physical quantities in the simulation.

To put the resulting residuals in context we compare them to the one picometer reference curve, which represents the coupling of a GW strain with an equivalent displacement amplitude of 1 pm/√Hz on all six links with relaxation at low frequencies, and is commonly used as a benchmark in LISA. Converted to frequency units and propagated through TDI it reads

$$ S_{\tilde{X}_2}^{\text{ref}}(f) = 4|\tilde{C}|^2(2\pi f)^2 \left(1 + \left(\frac{2 \text{ MHz}}{f}\right)^4\right) \left(\frac{\text{pm}/\sqrt{\text{Hz}}}{1064 \text{ nm}}\right)^2. \quad (52) $$

Here, the factor 4 accounts for the fact that there is a total of four links involved in the computation of the Michelson variables. The relaxation below 2 mHz mimics the shape of the LISA noise budget. Overall, the reference curve above is about a factor 10 below the total optical metrology noise allocation in LISA.

Let us now describe the nominal setup of the simulator and the processing performed on the measurements. All simulations span approximately 3 days in duration with the LISA constellation following realistic heliocentric orbits provided by ESA. The beatnote frequencies<fnref n="9" /> of the inter-spacecraft and the reference interferometers, which range between 5 MHz to 25 MHz, are simulated at a sampling rate of 16 Hz by LISA INSTRUMENT. They are formed by

<footnote n="9">⁹Units of frequency are used to circumvent numerical issues specific to phase units. A beatnote frequency of several MHz results in a rapidly increasing phase ramp. This would require a large number of significant digits to deal with its dynamic range (~10¹⁴ cycles after 1 yr at a required precision of at least µ cycles which results in 20 significant digits). On the other hand, its derivative, the beatnote frequency, stays roughly constant and a double precision float (approximately 16 significant digits) is sufficient to represent it. Here, the picometer reference curve with a knee at 2 mHz corresponds to a precision of approximately 20 nHz which gives 15 orders of magnitude when considering a 20 MHz beatnote frequency.</footnote>

<footnote n="8">⁸For example, the interspacecraft measurements $\text{isi}_{12} \sim \pi_{12}$ and $\text{isi}_{13} \sim \pi_{13}$.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. ASDs of the residuals from onboard processing for the second-generation Michelson variable $\dot{X}_2$ using the antiaiasing filter design that is optimized to reduce the number of coefficients to save computational cost (see Appendix A). The dashed blue and red lines correspond to the theoretical predictions for the different TDI residuals. The solid lines are the numerical estimates resulting from the noise-free simulation (gray) and the simulation containing laser noise (blue). For reference we also show the picometer reference curve (dashed black). It is multiplied by the transfer function of the filter to account for the significant pass-band drop of the optimized filter design (see Fig. 6). The black vertical line at $10^{-1}$ Hz divides the x-axis into logarithmically and linearly scaled.</figcaption></figure>

pairwise interference of the six lasers which are assumed to be stabilized to their individual cavity. Thus, we assume independent laser frequency noises with an ASD equal to $30 \text{ Hz}/\sqrt{\text{Hz}}$.<fnref n="10" /> Finally, the beatnote frequencies are decimated down to 4 Hz using the default antiaiasing filter of the simulator. In addition to the carrier-to-carrier beatnotes, sideband beatnotes that precisely track the delay derivatives and PRN measurements (absolute ranging) are also available. The nominal setup excludes any errors in the sideband beatnotes and absolute ranging measurements.

Prior to any processing, all measurements are promoted from a 64-bit to an 80-bit floating-point variable. This makes sure that no additional numerical noise is introduced downstream. In the first processing step we extract the high-precision ranging information from the sidebands (see Appendix D). Then, the second-generation Michelson variable $X_2$ is calculated using the factorization expression

in Eq. (13). Here, we nominally use Lagrange interpolation with $N = 62$ coefficients to suppress the interpolation residual in-band. To convert the expressions to frequency units, all delay operators $\hat{D}_{ij}$ are replaced by $\hat{D}_{ij}$, which are defined analogously to Eq. (B7) [6]. In the last processing step, we subtract any out-of-band drifts from the Michelson combinations to reduce the effect of spectral leakage at DC. This is achieved by computing the differential Doppler shift, as explained in Appendix E, and inserting it into the time derivative of Eq. (35). To compare the numerical results with the models we use the Welch method from the SciPy package [30] to estimate the numerical ASDs of the detrended Michelson variables. To match the analytic models at high frequencies with the numerical ASDs we had to relax the equal-arms assumption and use the models for six unequal but constant arms.

In a first simulation, we focus on the effects of onboard processing. As discussed previously the data rate on board the spacecraft has to be reduced from tens of megahertz down to a few hertz. This will be performed in several stages and the choice of intermediate sampling rates and the exact design of the antiaiasing filters is still under study. The LISA simulators currently only model the final

<footnote n="10">¹⁰In reality, we expect laser noise to increase toward lower frequencies. In this manuscript, we assume that it is a white noise to easily study the shape of various residuals. However, the residual transfer functions derived here hold for any laser noise spectrum.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. ASDs of the residuals from onboard processing for the second-generation Michelson variable $\dot{X}_2$ using the antiaiasing filter design that is optimized to reduce the number of coefficients to save computational cost (see Appendix A). The dashed blue and red lines correspond to the theoretical predictions for the different TDI residuals. The solid lines are the numerical estimates resulting from the noise-free simulation (gray) and the simulation containing laser noise (blue). For reference we also show the picometer reference curve (dashed black). It is multiplied by the transfer function of the filter to account for the significant pass-band drop of the optimized filter design (see Fig. 6). The black vertical line at $10^{-1}$ Hz divides the x-axis into logarithmically and linearly scaled.</figcaption></figure>

pairwise interference of the six lasers which are assumed to be stabilized to their individual cavity. Thus, we assume independent laser frequency noises with an ASD equal to $30 \text{ Hz}/\sqrt{\text{Hz}}$.<fnref n="10" /> Finally, the beatnote frequencies are decimated down to 4 Hz using the default antiaiasing filter of the simulator. In addition to the carrier-to-carrier beatnotes, sideband beatnotes that precisely track the delay derivatives and PRN measurements (absolute ranging) are also available. The nominal setup excludes any errors in the sideband beatnotes and absolute ranging measurements.

Prior to any processing, all measurements are promoted from a 64-bit to an 80-bit floating-point variable. This makes sure that no additional numerical noise is introduced downstream. In the first processing step we extract the high-precision ranging information from the sidebands (see Appendix D). Then, the second-generation Michelson variable $X_2$ is calculated using the factorization expression

in Eq. (13). Here, we nominally use Lagrange interpolation with $N = 62$ coefficients to suppress the interpolation residual in-band. To convert the expressions to frequency units, all delay operators $\hat{D}_{ij}$ are replaced by $\hat{D}_{ij}$, which are defined analogously to Eq. (B7) [6]. In the last processing step, we subtract any out-of-band drifts from the Michelson combinations to reduce the effect of spectral leakage at DC. This is achieved by computing the differential Doppler shift, as explained in Appendix E, and inserting it into the time derivative of Eq. (35). To compare the numerical results with the models we use the Welch method from the SciPy package [30] to estimate the numerical ASDs of the detrended Michelson variables. To match the analytic models at high frequencies with the numerical ASDs we had to relax the equal-arms assumption and use the models for six unequal but constant arms.

In a first simulation, we focus on the effects of onboard processing. As discussed previously the data rate on board the spacecraft has to be reduced from tens of megahertz down to a few hertz. This will be performed in several stages and the choice of intermediate sampling rates and the exact design of the antiaiasing filters is still under study. The LISA simulators currently only model the final

<footnote n="10">¹⁰In reality, we expect laser noise to increase toward lower frequencies. In this manuscript, we assume that it is a white noise to easily study the shape of various residuals. However, the residual transfer functions derived here hold for any laser noise spectrum.</footnote>
<figure><img src="image_6.png" /><figcaption>FIG. 4. ASDs of residual laser noise in second-generation Michelson variable $\dot{X}_2$ due to interpolation errors using Lagrange interpolation with $N = 42$. The numerical ASD using “atomic” delay operators to perform TDI (solid blue) is compared against the model (dashed blue) and the picometer reference curve (dashed black). The vertical dashed gray line indicates the high frequency edge of the LISA band.</figcaption></figure>

decimation stage going from 16 Hz down to 4 Hz. The default implementation of the filter in LISA INSTRUMENT is quite conservative (see Fig. 6) and consequently suppresses the effect of aliasing and flexing-filtering coupling well below the numerical noise floor. With such a filter, one is thus unable to compare the analytic model describing aliasing and flexing-filtering to the results of the simulation. For this reason, we employ a filter design that uses a reduced number of coefficients (five-fold compared to the default implementation, for more details see Appendix A). This filter design reduces computational cost but amplifies the coupling of the aliasing and flexing-filtering effect. In Fig. 3 the numerical ASD (solid blue) of the residual in the Michelson variable $\dot{X}_2$ is compared to the analytical models for both effects. We observe that for our particular choice of orbits the flexing-filtering effect (dashed red) is subdominant and the entire residual is well explained by the aliasing effect (dashed blue). At low frequencies ($10^{-4}$ Hz to $10^{-3}$ Hz) the simulated data approaches the numerical noise floor given by the light gray curve that is calculated from the noise free simulation. We confirm that the residual

<figure><img src="image_7.png" /><figcaption>FIG. 5. ASDs of the laser and timing noise residuals caused by ranging errors for the second-generation Michelson variable $\dot{X}_2$. The blue dashed and dotted lines correspond to the theoretical predictions for the ranging bias and the modulation noise coupling, respectively. Solid lines are numerical estimates resulting from the noise free simulation (gray) and the simulation containing laser noise, modulation noise and ranging biases (blue). For reference we plot the picometer reference curve in dashed black. Here, the attenuation by the transfer function of the filter can be neglected as we use the default filter design of LISA INSTRUMENT which is extremely flat in the LISA band. The black vertical line at $10^{-1}$ Hz divides the x-axis into logarithmically and linearly scaled.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 4. ASDs of residual laser noise in second-generation Michelson variable $\dot{X}_2$ due to interpolation errors using Lagrange interpolation with $N = 42$. The numerical ASD using “atomic” delay operators to perform TDI (solid blue) is compared against the model (dashed blue) and the picometer reference curve (dashed black). The vertical dashed gray line indicates the high frequency edge of the LISA band.</figcaption></figure>

decimation stage going from 16 Hz down to 4 Hz. The default implementation of the filter in LISA INSTRUMENT is quite conservative (see Fig. 6) and consequently suppresses the effect of aliasing and flexing-filtering coupling well below the numerical noise floor. With such a filter, one is thus unable to compare the analytic model describing aliasing and flexing-filtering to the results of the simulation. For this reason, we employ a filter design that uses a reduced number of coefficients (five-fold compared to the default implementation, for more details see Appendix A). This filter design reduces computational cost but amplifies the coupling of the aliasing and flexing-filtering effect. In Fig. 3 the numerical ASD (solid blue) of the residual in the Michelson variable $\dot{X}_2$ is compared to the analytical models for both effects. We observe that for our particular choice of orbits the flexing-filtering effect (dashed red) is subdominant and the entire residual is well explained by the aliasing effect (dashed blue). At low frequencies ($10^{-4}$ Hz to $10^{-3}$ Hz) the simulated data approaches the numerical noise floor given by the light gray curve that is calculated from the noise free simulation. We confirm that the residual

<figure><img src="image_7.png" /><figcaption>FIG. 5. ASDs of the laser and timing noise residuals caused by ranging errors for the second-generation Michelson variable $\dot{X}_2$. The blue dashed and dotted lines correspond to the theoretical predictions for the ranging bias and the modulation noise coupling, respectively. Solid lines are numerical estimates resulting from the noise free simulation (gray) and the simulation containing laser noise, modulation noise and ranging biases (blue). For reference we plot the picometer reference curve in dashed black. Here, the attenuation by the transfer function of the filter can be neglected as we use the default filter design of LISA INSTRUMENT which is extremely flat in the LISA band. The black vertical line at $10^{-1}$ Hz divides the x-axis into logarithmically and linearly scaled.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Magnitude of transfer functions (solid lines) and their frequency derivatives (dashed lines) for the default implementation of the antialiasing filter in LISA INSTRUMENT (red) and the optimized design (blue). The gray and white frequency bands indicate the orientation of the subbands after folding where gray bands are mirrored before they are added to the new band (0 Hz to 2 Hz) and white bands are not.</figcaption></figure>

caused by onboard processing remains below the one picometer reference curve (dashed black).

In the second simulation, we use the default filter design of LISA INSTRUMENT. It strongly suppresses the commutator residuals such that it becomes possible to compare the interpolation residual to the corresponding theoretical model. To perform TDI we choose Lagrange interpolation with a reduced number of coefficients $N = 42$ and consider “atomic” delays instead of their contracted equivalent [see Eq. (16)]. Figure 4 shows the numerical result in solid blue and the model for the interpolation residual in dashed blue. We observe that the model explains the data well up to the numerical noise floor (light gray) and that the interpolation residual drops off quickly from the LISA band edge at 1 Hz toward lower frequencies. Furthermore, we confirm that $N = 42$ is sufficient to stay below the one picometer reference curve (dashed black) in the LISA band.

Finally, in the third simulation, we introduce modulation noise and ranging biases of the order of 30 ns (a factor 10 higher than in [13] in order to enhance the size of the coupling) so as to produce the couplings related to ranging errors. The ASD of the modulation noise on left-handed MOSAs is given by

$$ \sqrt{S_M(f)} = 8.3 \times 10^{-15} \text{ s/}\sqrt{\text{Hz}} \left(\frac{f}{\text{Hz}}\right)^{-2/3}, \quad (53) $$

while on the right-handed MOSAs it is a factor of 10 larger [31]. As already presented in [12,13], we can obtain an estimate of the pseudo-range $\hat{d}_{ij}$ that is free of the modulation errors of the right-handed MOSAs. To perform TDI we use the nominal number of interpolation coefficients

$N = 62$ and use the ranges that are contaminated by the ranging error. We present the resulting ASDs in Fig. 5. At low frequencies the numerical residual (solid blue) is explained by the coupling to modulation noise (dotted blue) which exceeds the one picometer reference curve at frequencies below $3 \times 10^{-2}$ Hz, highlighting that this is a significant contribution to the overall LISA noise budget. At high frequencies the effect of ranging biases (dashed blue) becomes significant but stays below the reference curve for the chosen numerical values of the biases.

Repeating the same set of simulation runs for locked lasers, we find that they also agree with the analytic models for the laser noise coupling residuals discussed in Sec. IV B. In addition, the plots we obtain differ only very slightly from the ones where six independent lasers were used. For this reason, we choose not to include them in the paper.

## V. CONCLUSION

In this paper, we have presented a comprehensive study of residual laser noise in the TDI variables for LISA. We have identified two categories of couplings. First, on-board processing steps, namely filtering and decimation, give rise to additive noise in the interspacecraft interferometers due to noncommutation with the delay operation. Second, the postprocessing delays employed to calculate TDI combinations only partially mitigate laser noise. This is because this offline computation relies on an interpolation method, which induces interpolation errors, and because the offline delays used for TDI include ranging errors. For both categories of laser noise couplings, we provide analytical models for the residuals in the second-generation TDI combination $X_2$; we validate those models using numerical simulations.

In the existing literature, the flexing-filtering effect [15] and the coupling of ranging errors [13] for six independent lasers are already described and preliminary models of laser noise residuals due to aliasing and interpolation errors are presented in [16]. In this study, we remedy the shortcomings of the latter two and also present all laser noise couplings in a consistent framework. Furthermore, we investigate the impact of laser locking by discussing the case of locked lasers. Finally, because we perform TDI in total frequency units, we explain the deterministic trend that is present in second-generation TDI combinations: differential Doppler shifts in the round-trip paths of the synthesized beams produce a beatnote of a few mHz. This trend depends solely on the out-of-band delays due to orbital dynamics and the THz frequency of the laser involved. It can therefore be computed and removed by appropriately modeling the orbits and the THz frequency evolution. This detrending step is a reversible (the trend can always be added in again) part of preprocessing. It occurs before parameter estimation and reduces spectral leakage in
<figure><img src="image_8.png" /><figcaption>FIG. 6. Magnitude of transfer functions (solid lines) and their frequency derivatives (dashed lines) for the default implementation of the antialiasing filter in LISA INSTRUMENT (red) and the optimized design (blue). The gray and white frequency bands indicate the orientation of the subbands after folding where gray bands are mirrored before they are added to the new band (0 Hz to 2 Hz) and white bands are not.</figcaption></figure>

caused by onboard processing remains below the one picometer reference curve (dashed black).

In the second simulation, we use the default filter design of LISA INSTRUMENT. It strongly suppresses the commutator residuals such that it becomes possible to compare the interpolation residual to the corresponding theoretical model. To perform TDI we choose Lagrange interpolation with a reduced number of coefficients $N = 42$ and consider “atomic” delays instead of their contracted equivalent [see Eq. (16)]. Figure 4 shows the numerical result in solid blue and the model for the interpolation residual in dashed blue. We observe that the model explains the data well up to the numerical noise floor (light gray) and that the interpolation residual drops off quickly from the LISA band edge at 1 Hz toward lower frequencies. Furthermore, we confirm that $N = 42$ is sufficient to stay below the one picometer reference curve (dashed black) in the LISA band.

Finally, in the third simulation, we introduce modulation noise and ranging biases of the order of 30 ns (a factor 10 higher than in [13] in order to enhance the size of the coupling) so as to produce the couplings related to ranging errors. The ASD of the modulation noise on left-handed MOSAs is given by

$$ \sqrt{S_M(f)} = 8.3 \times 10^{-15} \text{ s/}\sqrt{\text{Hz}} \left(\frac{f}{\text{Hz}}\right)^{-2/3}, \quad (53) $$

while on the right-handed MOSAs it is a factor of 10 larger [31]. As already presented in [12,13], we can obtain an estimate of the pseudo-range $\hat{d}_{ij}$ that is free of the modulation errors of the right-handed MOSAs. To perform TDI we use the nominal number of interpolation coefficients

$N = 62$ and use the ranges that are contaminated by the ranging error. We present the resulting ASDs in Fig. 5. At low frequencies the numerical residual (solid blue) is explained by the coupling to modulation noise (dotted blue) which exceeds the one picometer reference curve at frequencies below $3 \times 10^{-2}$ Hz, highlighting that this is a significant contribution to the overall LISA noise budget. At high frequencies the effect of ranging biases (dashed blue) becomes significant but stays below the reference curve for the chosen numerical values of the biases.

Repeating the same set of simulation runs for locked lasers, we find that they also agree with the analytic models for the laser noise coupling residuals discussed in Sec. IV B. In addition, the plots we obtain differ only very slightly from the ones where six independent lasers were used. For this reason, we choose not to include them in the paper.

## V. CONCLUSION

In this paper, we have presented a comprehensive study of residual laser noise in the TDI variables for LISA. We have identified two categories of couplings. First, on-board processing steps, namely filtering and decimation, give rise to additive noise in the interspacecraft interferometers due to noncommutation with the delay operation. Second, the postprocessing delays employed to calculate TDI combinations only partially mitigate laser noise. This is because this offline computation relies on an interpolation method, which induces interpolation errors, and because the offline delays used for TDI include ranging errors. For both categories of laser noise couplings, we provide analytical models for the residuals in the second-generation TDI combination $X_2$; we validate those models using numerical simulations.

In the existing literature, the flexing-filtering effect [15] and the coupling of ranging errors [13] for six independent lasers are already described and preliminary models of laser noise residuals due to aliasing and interpolation errors are presented in [16]. In this study, we remedy the shortcomings of the latter two and also present all laser noise couplings in a consistent framework. Furthermore, we investigate the impact of laser locking by discussing the case of locked lasers. Finally, because we perform TDI in total frequency units, we explain the deterministic trend that is present in second-generation TDI combinations: differential Doppler shifts in the round-trip paths of the synthesized beams produce a beatnote of a few mHz. This trend depends solely on the out-of-band delays due to orbital dynamics and the THz frequency of the laser involved. It can therefore be computed and removed by appropriately modeling the orbits and the THz frequency evolution. This detrending step is a reversible (the trend can always be added in again) part of preprocessing. It occurs before parameter estimation and reduces spectral leakage in
PSD estimates, a feature which is relevant for the present study.

Contrary to unsuppressed noises [11] (e.g., path length noises), the coupling of laser noise residuals is dependent on the underlying locking configuration. This can be explained by the fact that locked lasers follow the primary laser with configuration-specific time lags, which introduce correlations among all lasers. In order to connect with the existing literature, we first derive analytical models for six independent lasers. Then, we repeat the calculation for locked lasers, more specifically the configuration N1-12. In this configuration, spacecraft 2 and 3 act as transponders, directly sending light back to spacecraft 1. This means the interspacecraft interferometer beatnotes recorded on spacecraft 1 already represent the signal combinations $\eta_{12} + D_{12}\eta_{21}$ and $\eta_{13} + D_{13}\eta_{31}$ for laser noise, simplifying the expression for the Michelson $X_2$ variable. However, this simplification does not apply to $Y_2$ and $Z_2$. In general, because of the additional correlations, analytic models become more complicated for locked lasers. Furthermore, we find that the worst case scenario (i.e., the one that maximizes the size of the laser noise residuals) can be larger for locked lasers than for six independent lasers. This suggests that future investigations should account for all possible locking configurations.

The level of additional noise due to on-board processing (filtering and decimation) is strongly dependent on the design of the antialiasing filter. In this study, we propose an improved FIR filter design that reduces the number of filter taps by a factor of five compared to the default (more conservative) implementation in LISA INSTRUMENT. This relaxes the computational requirements of the onboard filter implementation and we show that the resulting residual due to the flexing-filtering effect and aliasing still respects the picometer reference curve in the LISA band. In theory, the flexing-filtering coupling can be mitigated by flattening the response of the antialiasing filter on ground by compensating for its pass-band droop. The appropriate design of compensation filters is the subject of on-going efforts in the LISA community. Indeed, any aliased noise due to insufficient attenuation in the stop-band cannot be reduced in postprocessing and has to be taken care of before decimation.

Laser noise residuals stemming from postprocessing delays depend on the interpolation method and the ability to make precise ranging measurements. In this paper, we extensively rely on Lagrange interpolation, which has a maximal flat response at DC. “Windowed sinc” kernels [32] or numerically optimized kernels provide alternative interpolation options. These are currently under study. Interpolation kernels of shorter length have smaller computational cost and result in less truncation at the boundaries (where the interpolation kernel does not completely overlap with the data). This problem becomes more critical in the presence of gaps. In this paper, we studied the impact of

contracting processing delays, i.e., combining nested delays first to form a single delay operation. We find that the best delay contraction strategy depends on the locking configuration and the particular numerical values of the delays.

As mentioned in the previous paragraph, the determination of accurate and precise intersatellite ranges is crucial for reducing laser noise in TDI. We show that a ranging bias of a few meters is sufficient to suppress laser noise below the picometer reference curve. Furthermore, the current estimate for the in-band component of the ranging error is modulation noise and exceeds the reference curve at low frequencies. This highlights that modulation noise is a more significant contribution in the overall LISA noise budget.

The results presented in this paper should be independent of the time reference frames the measurements are defined in. Therefore, the general findings should still be valid for measurements sampled according to realistic clocks that are processed using “Time-delay interferometry without clock synchronization” presented in [13]. Additionally to flexing arms due to orbital dynamics, clock drifts of the order of $10^{-7}$ become relevant for the flexing-filtering effect. Furthermore, extra care must be taken when extracting the delay estimates from the sideband measurements since measured pseudo-ranges have an in-band component. Here, equivalent filtering must be applied to the sideband beatnotes such that they are compatible with the clock noise contributions in the carrier beatnotes. Any difference in the filters’ transfer function will produce additional ranging noise in the delay estimates.

## ACKNOWLEDGMENTS

The authors would like to thank G. Heinzel for the valuable discussions and the feedback on the manuscript. Additionally, the authors thank the LISA Simulation Working Group and the LISA Simulation Expert Group for the lively discussions on all simulation-related activities. M. S. and O. H. acknowledge the support of the German Space Agency, DLR. The work is supported by the Federal Ministry for Economic Affairs and Climate Action based on a decision by the German Bundestag (FKZ 500Q1801 and FKZ 500Q2301). This work is also supported by the Max-Planck-Society within the LEGACY (“Low-Frequency Gravitational Wave Astronomy in Space”) collaboration (M.IF.A.QOP18098). J.-B. B. gratefully acknowledges support from the UK Space Agency via STFC [ST/W002825/1]. M. S., O. H. and M. L. gratefully acknowledge support from the Centre National d'Études Spatiales (CNES).

## APPENDIX A: ANTIALIASING FILTER DESIGN

The main objective of the antialiasing filter is to prevent power folding during decimation. In the numerical simulations performed in this work, we simulate the physics at 16 Hz and decimate to 4 Hz. As the LISA band only reaches up to 1 Hz we allow for aliasing between 1 Hz and
PSD estimates, a feature which is relevant for the present study.

Contrary to unsuppressed noises [11] (e.g., path length noises), the coupling of laser noise residuals is dependent on the underlying locking configuration. This can be explained by the fact that locked lasers follow the primary laser with configuration-specific time lags, which introduce correlations among all lasers. In order to connect with the existing literature, we first derive analytical models for six independent lasers. Then, we repeat the calculation for locked lasers, more specifically the configuration N1-12. In this configuration, spacecraft 2 and 3 act as transponders, directly sending light back to spacecraft 1. This means the interspacecraft interferometer beatnotes recorded on spacecraft 1 already represent the signal combinations $\eta_{12} + D_{12}\eta_{21}$ and $\eta_{13} + D_{13}\eta_{31}$ for laser noise, simplifying the expression for the Michelson $X_2$ variable. However, this simplification does not apply to $Y_2$ and $Z_2$. In general, because of the additional correlations, analytic models become more complicated for locked lasers. Furthermore, we find that the worst case scenario (i.e., the one that maximizes the size of the laser noise residuals) can be larger for locked lasers than for six independent lasers. This suggests that future investigations should account for all possible locking configurations.

The level of additional noise due to on-board processing (filtering and decimation) is strongly dependent on the design of the antialiasing filter. In this study, we propose an improved FIR filter design that reduces the number of filter taps by a factor of five compared to the default (more conservative) implementation in LISA INSTRUMENT. This relaxes the computational requirements of the onboard filter implementation and we show that the resulting residual due to the flexing-filtering effect and aliasing still respects the picometer reference curve in the LISA band. In theory, the flexing-filtering coupling can be mitigated by flattening the response of the antialiasing filter on ground by compensating for its pass-band droop. The appropriate design of compensation filters is the subject of on-going efforts in the LISA community. Indeed, any aliased noise due to insufficient attenuation in the stop-band cannot be reduced in postprocessing and has to be taken care of before decimation.

Laser noise residuals stemming from postprocessing delays depend on the interpolation method and the ability to make precise ranging measurements. In this paper, we extensively rely on Lagrange interpolation, which has a maximal flat response at DC. “Windowed sinc” kernels [32] or numerically optimized kernels provide alternative interpolation options. These are currently under study. Interpolation kernels of shorter length have smaller computational cost and result in less truncation at the boundaries (where the interpolation kernel does not completely overlap with the data). This problem becomes more critical in the presence of gaps. In this paper, we studied the impact of

contracting processing delays, i.e., combining nested delays first to form a single delay operation. We find that the best delay contraction strategy depends on the locking configuration and the particular numerical values of the delays.

As mentioned in the previous paragraph, the determination of accurate and precise intersatellite ranges is crucial for reducing laser noise in TDI. We show that a ranging bias of a few meters is sufficient to suppress laser noise below the picometer reference curve. Furthermore, the current estimate for the in-band component of the ranging error is modulation noise and exceeds the reference curve at low frequencies. This highlights that modulation noise is a more significant contribution in the overall LISA noise budget.

The results presented in this paper should be independent of the time reference frames the measurements are defined in. Therefore, the general findings should still be valid for measurements sampled according to realistic clocks that are processed using “Time-delay interferometry without clock synchronization” presented in [13]. Additionally to flexing arms due to orbital dynamics, clock drifts of the order of $10^{-7}$ become relevant for the flexing-filtering effect. Furthermore, extra care must be taken when extracting the delay estimates from the sideband measurements since measured pseudo-ranges have an in-band component. Here, equivalent filtering must be applied to the sideband beatnotes such that they are compatible with the clock noise contributions in the carrier beatnotes. Any difference in the filters’ transfer function will produce additional ranging noise in the delay estimates.

## ACKNOWLEDGMENTS

The authors would like to thank G. Heinzel for the valuable discussions and the feedback on the manuscript. Additionally, the authors thank the LISA Simulation Working Group and the LISA Simulation Expert Group for the lively discussions on all simulation-related activities. M. S. and O. H. acknowledge the support of the German Space Agency, DLR. The work is supported by the Federal Ministry for Economic Affairs and Climate Action based on a decision by the German Bundestag (FKZ 500Q1801 and FKZ 500Q2301). This work is also supported by the Max-Planck-Society within the LEGACY (“Low-Frequency Gravitational Wave Astronomy in Space”) collaboration (M.IF.A.QOP18098). J.-B. B. gratefully acknowledges support from the UK Space Agency via STFC [ST/W002825/1]. M. S., O. H. and M. L. gratefully acknowledge support from the Centre National d'Études Spatiales (CNES).

## APPENDIX A: ANTIALIASING FILTER DESIGN

The main objective of the antialiasing filter is to prevent power folding during decimation. In the numerical simulations performed in this work, we simulate the physics at 16 Hz and decimate to 4 Hz. As the LISA band only reaches up to 1 Hz we allow for aliasing between 1 Hz and
2 Hz, the Nyquist frequency after decimation. To save computational cost we aim to use a minimum number of filter coefficients. We optimize the filter design by adopting the Parks-McClellan algorithm [33,34] to account for the fact that the picometer reference curve has a f¹ shape in amplitude. The resulting transfer function and its derivative which is of relevance for the flexing-filtering coupling (see Appendix B 1) are plotted in Fig. 6. As described above the attenuation of the filter in the stop-band is relaxed between 2 Hz and 3 Hz, and 5 Hz and 7 Hz as those bands reside at 1 Hz to 2 HZ after folding. The final filter design that respects the picometer reference curve for aliasing discussed in Sec. IV and Appendix B has a total of 29 coefficients which represents a reduction by a factor of five compared to the default implementation of the LISA simulators which use 145 coefficients.

# APPENDIX B: DERIVATION OF ON-BOARD PROCESSING RESIDUALS

## 1. Flexing-filtering coupling

Let us derive the contribution coming from the commutator $[\mathbf{F}, \mathbf{D}]$, first described in [15] and dubbed “flexing-filtering coupling.” We assume that the delay $d(t)$ is slowly varying over the filter length and that its first derivative $\dot{d}(t)$ is small. We use Eq. (7), and expand $[\mathbf{F}, \mathbf{D}]\phi(t)$ at leading order in $\dot{d}(t)$ to find

$$[\mathbf{F}, \mathbf{D}]\phi(t) \approx \dot{d}(t)\mathbf{D}\mathbf{G}\nu(t). \quad (B1)$$

where $\mathbf{G}$ is a filter operator defined similarly to $\mathbf{F}$ in Eq. (7) with $h_G(\tau) = \tau \cdot h_F(\tau)$.

Because of the orbital dynamics of the LISA constellation, $\dot{d}(t)$ will vary on timescales of several months and, as a result, so will the level of the laser noise residual introduced by the flexing-filtering coupling. However, over sufficiently short observation times we can assume $\dot{d}$ to be constant. In such a case the PSD of Eq. (B1) reads

$$S_{\delta\phi}^{[\mathbf{F}, \mathbf{D}]}(f) = \dot{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\nu}(f). \quad (B2)$$

For longer observation times, one can use the maximum value of $\dot{d}$ as given by current predictions for the orbital dynamics of LISA in order to derive an upper bound for this PSD.

## 2. Decimation-delay commutator

The second commutator appearing in Eq. (21) is the commutator of the decimation and delay operations,

$$[\mathbf{S}, \mathbf{D}]\phi(t) = \mathbf{S}\mathbf{D}\phi(t) - \mathbf{D}\mathbf{S}\phi(t). \quad (B3)$$

Those operations do not commute due to the nonlinear nature of the decimation process. The PSD of this

expression can be derived using the definition in Eq. (10). The different aliases that are folded in band are modulated by a sine-squared factor. We obtain

$$S_{\delta\phi}^{[\mathbf{S}, \mathbf{D}]}(f) = 4 \cdot \sum_{n=1}^{M-1} c_n(d) \cdot S_{\phi}^{(n)}(f), \quad (B4)$$

where the modulating factor $c_n$ is given by

$$c_n(d) = \begin{cases} \sin^2(\pi f_s d^{\frac{n}{2}}) & \text{for } n \text{ even,} \\ \sin^2(\pi f_s d^{\frac{n+1}{2}}) & \text{for } n \text{ odd.} \end{cases} \quad (B5)$$

In the special case where the delay $d$ becomes an integer multiple of the sampling time $T_s = 1/f_s$ decimation and delay operations commute and the residual becomes zero. For a time-varying delay, this particular residual is non-stationary because its power is modulated as $d$ evolves in time. To account for this, we later consider only the upper bound obtained for $c_n = 1$ for all $n$,

$$S_{\delta\phi}^{[\mathbf{S}, \mathbf{D}]}(f) \le 4 \cdot \sum_{n=1}^{\infty} S_{\phi}^{(n)}(f). \quad (B6)$$

This bound is independent of the delay, and corresponds to the case of full anticorrelation between $\mathbf{S}\mathbf{D}\phi(t)$ and $\mathbf{D}\mathbf{S}\phi(t)$ in Eq. (B3).

## 3. Comparison with numerical simulations

The numerical simulations that are used in this work to validate the analytic models are performed in units of frequency in order to preserve numerical precision. As shown in [6], any delay operation on frequency data can be represented by the usual shift of the argument and a multiplicative Doppler factor,

$$\dot{\mathbf{D}}\nu(t) = (1 - \dot{d}(t)) \cdot \nu(t - d(t)). \quad (B7)$$

We can then easily rewrite the commutator given in Eq. (21) in terms of frequency data by replacing every occurrence of the delay operator $\mathbf{D}$ by its Doppler equivalent. It reads

$$[\mathbf{S}\mathbf{F}, \dot{\mathbf{D}}]\nu(t) = \mathbf{S}\mathbf{F}\dot{\mathbf{D}}\nu(t) - \dot{\mathbf{D}}\mathbf{S}\mathbf{F}\nu(t). \quad (B8)$$

Here, we need to account for the Doppler factor to cancel laser noise to first order. However, we find that it only has a negligible impact on the laser noise residual, and we can write for the PSDs

$$S_{\delta\nu}^{[\mathbf{S}, \dot{\mathbf{F}}]}(f) \approx |\tilde{\mathbf{S}}|^2 \left( (2\pi f)^2 \dot{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\nu}(f) \right), \quad (B9a)$$

$$S_{\delta\nu}^{[\mathbf{S}, \dot{\mathbf{D}}]\mathbf{F}}(f) \approx 4 \sum_{n=1}^{\infty} c_n(d) (|\tilde{\mathbf{F}}|^2 S_{\nu})^{(n)}(f). \quad (B9b)$$
2 Hz, the Nyquist frequency after decimation. To save computational cost we aim to use a minimum number of filter coefficients. We optimize the filter design by adopting the Parks-McClellan algorithm [33,34] to account for the fact that the picometer reference curve has a f¹ shape in amplitude. The resulting transfer function and its derivative which is of relevance for the flexing-filtering coupling (see Appendix B 1) are plotted in Fig. 6. As described above the attenuation of the filter in the stop-band is relaxed between 2 Hz and 3 Hz, and 5 Hz and 7 Hz as those bands reside at 1 Hz to 2 HZ after folding. The final filter design that respects the picometer reference curve for aliasing discussed in Sec. IV and Appendix B has a total of 29 coefficients which represents a reduction by a factor of five compared to the default implementation of the LISA simulators which use 145 coefficients.

# APPENDIX B: DERIVATION OF ON-BOARD PROCESSING RESIDUALS

## 1. Flexing-filtering coupling

Let us derive the contribution coming from the commutator $[\mathbf{F}, \mathbf{D}]$, first described in [15] and dubbed “flexing-filtering coupling.” We assume that the delay $d(t)$ is slowly varying over the filter length and that its first derivative $\dot{d}(t)$ is small. We use Eq. (7), and expand $[\mathbf{F}, \mathbf{D}]\phi(t)$ at leading order in $\dot{d}(t)$ to find

$$[\mathbf{F}, \mathbf{D}]\phi(t) \approx \dot{d}(t)\mathbf{D}\mathbf{G}\nu(t). \quad (B1)$$

where $\mathbf{G}$ is a filter operator defined similarly to $\mathbf{F}$ in Eq. (7) with $h_G(\tau) = \tau \cdot h_F(\tau)$.

Because of the orbital dynamics of the LISA constellation, $\dot{d}(t)$ will vary on timescales of several months and, as a result, so will the level of the laser noise residual introduced by the flexing-filtering coupling. However, over sufficiently short observation times we can assume $\dot{d}$ to be constant. In such a case the PSD of Eq. (B1) reads

$$S_{\delta\phi}^{[\mathbf{F}, \mathbf{D}]}(f) = \dot{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\nu}(f). \quad (B2)$$

For longer observation times, one can use the maximum value of $\dot{d}$ as given by current predictions for the orbital dynamics of LISA in order to derive an upper bound for this PSD.

## 2. Decimation-delay commutator

The second commutator appearing in Eq. (21) is the commutator of the decimation and delay operations,

$$[\mathbf{S}, \mathbf{D}]\phi(t) = \mathbf{S}\mathbf{D}\phi(t) - \mathbf{D}\mathbf{S}\phi(t). \quad (B3)$$

Those operations do not commute due to the nonlinear nature of the decimation process. The PSD of this

expression can be derived using the definition in Eq. (10). The different aliases that are folded in band are modulated by a sine-squared factor. We obtain

$$S_{\delta\phi}^{[\mathbf{S}, \mathbf{D}]}(f) = 4 \cdot \sum_{n=1}^{M-1} c_n(d) \cdot S_{\phi}^{(n)}(f), \quad (B4)$$

where the modulating factor $c_n$ is given by

$$c_n(d) = \begin{cases} \sin^2(\pi f_s d^{\frac{n}{2}}) & \text{for } n \text{ even,} \\ \sin^2(\pi f_s d^{\frac{n+1}{2}}) & \text{for } n \text{ odd.} \end{cases} \quad (B5)$$

In the special case where the delay $d$ becomes an integer multiple of the sampling time $T_s = 1/f_s$ decimation and delay operations commute and the residual becomes zero. For a time-varying delay, this particular residual is non-stationary because its power is modulated as $d$ evolves in time. To account for this, we later consider only the upper bound obtained for $c_n = 1$ for all $n$,

$$S_{\delta\phi}^{[\mathbf{S}, \mathbf{D}]}(f) \le 4 \cdot \sum_{n=1}^{\infty} S_{\phi}^{(n)}(f). \quad (B6)$$

This bound is independent of the delay, and corresponds to the case of full anticorrelation between $\mathbf{S}\mathbf{D}\phi(t)$ and $\mathbf{D}\mathbf{S}\phi(t)$ in Eq. (B3).

## 3. Comparison with numerical simulations

The numerical simulations that are used in this work to validate the analytic models are performed in units of frequency in order to preserve numerical precision. As shown in [6], any delay operation on frequency data can be represented by the usual shift of the argument and a multiplicative Doppler factor,

$$\dot{\mathbf{D}}\nu(t) = (1 - \dot{d}(t)) \cdot \nu(t - d(t)). \quad (B7)$$

We can then easily rewrite the commutator given in Eq. (21) in terms of frequency data by replacing every occurrence of the delay operator $\mathbf{D}$ by its Doppler equivalent. It reads

$$[\mathbf{S}\mathbf{F}, \dot{\mathbf{D}}]\nu(t) = \mathbf{S}\mathbf{F}\dot{\mathbf{D}}\nu(t) - \dot{\mathbf{D}}\mathbf{S}\mathbf{F}\nu(t). \quad (B8)$$

Here, we need to account for the Doppler factor to cancel laser noise to first order. However, we find that it only has a negligible impact on the laser noise residual, and we can write for the PSDs

$$S_{\delta\nu}^{[\mathbf{S}, \dot{\mathbf{F}}]}(f) \approx |\tilde{\mathbf{S}}|^2 \left( (2\pi f)^2 \dot{d}^2 \left| \frac{1}{2\pi} \frac{d\tilde{h}_{\mathbf{F}}(f)}{df} \right|^2 S_{\nu}(f) \right), \quad (B9a)$$

$$S_{\delta\nu}^{[\mathbf{S}, \dot{\mathbf{D}}]\mathbf{F}}(f) \approx 4 \sum_{n=1}^{\infty} c_n(d) (|\tilde{\mathbf{F}}|^2 S_{\nu})^{(n)}(f). \quad (B9b)$$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Commutator residuals from filtering and decimation. We compare the numerically simulated data (solid blue) against analytic models (dashed lines). In dashed blue we show the model for the coupling of the decimation-delay commutator (aliasing) and in dashed red the flexing-filtering effect. For reference we also plot the picometer reference curve in dashed black including the transfer function of the filter.</figcaption></figure>

We note that for the full commutator in Eq. (21), we need to apply S to the flexing-filtering contribution and account for the fact that laser noise is filtered prior to passing it through the decimation-delay commutator. In our numerical implementation, the effect of decimating the flexing-filtering residual was negligible compared to the in-band contribution. Nevertheless, we include it for completeness as it strongly depends on the filter design.

In order to test the validity of the theoretical model given by Eq. (B9), we compare it to numerical simulations that generate time series corresponding to Eq. (B8). We consider white laser frequency noise in $\nu(t)$ with an ASD of $30 \text{ Hz}/\sqrt{\text{Hz}}$ and neglect the central laser frequency of $2.816 \times 10^{14} \text{ Hz}$ (its coupling to the commutator is vanishing). The simulation is performed at a sampling rate of 16 Hz and is then decimated down to 4 Hz. Antialiasing is performed using an FIR designed according to Appendix A. The delay operator is modeled by numerical interpolation of the data using Lagrange polynomials (cf. Sec. III B).

In Fig. 7 we compare the simulated data against the analytical model presented in Eq. (B9). The numerical result is entirely explained by aliasing as the flexing-filtering remains subdominant across the band.

# APPENDIX C: DERIVATION OF TDI RESIDUALS

## 1. Interpolation residual

We follow [32] and implement the postprocessing delay operation in TDI as an interpolation with a fractional delay filter. Lagrange interpolation has demonstrated its suitability in this context, and we thus use it as the baseline in this work. In general, we can model postprocessing delays as

FIR filters. For convenience, we split the delay into an integer shift $j$ and a fractional shift $\epsilon$ ranging from 0 to 1. The latter defines the coefficients of the interpolation kernel $k_n(\epsilon)$ with $-N/2 \le n \le N/2 - 1$, which is convolved with the discrete data samples $\phi_n$,

$$ (\mathcal{D}\phi)_n = \sum_{m=-N/2}^{N/2-1} k_m(\epsilon) \cdot \phi_{n-j-m}. \quad (C1) $$

Before the convolution of Eq. (C1) is performed, the data samples are shifted by the integer shift $j$. The above formula holds for interpolation kernels of even length $N$. For odd $N$, a similar expression can be derived. Using Eq. (7), we find

$$ h_D(\tau) = \sum_{m=-N/2}^{N/2-1} k_m(\epsilon) \cdot \delta(\tau - (j+m)T_s). \quad (C2) $$

We define the additional phase residual caused by the interpolation error as

$$ \delta\phi^D(t) = (\mathcal{D} - \mathbf{D})\phi(t) = \underbrace{\mathbf{D}(\mathbf{D}^{-1}\mathcal{D} - 1)\phi(t)}_{\Delta} \quad (C3) $$

and use Eq. (8) to derive the residual in terms of PSD. We find

$$ S_{\delta\phi}^{D}(f) = \underbrace{|\tilde{h}_{D}(f)e^{2\pi ifd} - 1|^{2}}_{|\Delta|^{2}} S_{\phi}(f). \quad (C4) $$

In general, the interpolation kernel $k_m(\epsilon)$ has to be adjusted for every sample $n$ to account for time-dependent time shifts. Therefore, the flexing arms of LISA will produce a nonstationary interpolation residual as the fractional delay $\epsilon$ scans through different values. At $\epsilon = 0$ and $\epsilon = 1$, the delay is a pure integer shift and the residual vanishes. Assuming that the worst case is obtained for $\epsilon = 0.5$,<fnref n="11" /> we can derive an upper bound for the residual induced by interpolation.

As suggested in [32], and already mentioned at the beginning of this section, a suitable interpolation method is Lagrange interpolation, and we use in this work. The interpolation kernel $k_m$ is derived from fitting a Lagrange polynomial through a set of neighboring samples. This method is known for producing a maximally flat frequency response at DC and is therefore well suited for LISA data processing, as it performs well over the entire LISA band ($10^{-4}$ to 1 Hz) when using high interpolation orders.<fnref n="12" />

<footnote n="11">¹¹A rigorous proof is needed to validate this assumption which was only found to be true *empirically*, in our work, and specifically for Lagrange interpolation.</footnote>

<footnote n="12">¹²Assuming a sampling rate of 4 Hz, typical filter lengths used for postprocessing delays in TDI are in the range 32 to 66.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Commutator residuals from filtering and decimation. We compare the numerically simulated data (solid blue) against analytic models (dashed lines). In dashed blue we show the model for the coupling of the decimation-delay commutator (aliasing) and in dashed red the flexing-filtering effect. For reference we also plot the picometer reference curve in dashed black including the transfer function of the filter.</figcaption></figure>

We note that for the full commutator in Eq. (21), we need to apply S to the flexing-filtering contribution and account for the fact that laser noise is filtered prior to passing it through the decimation-delay commutator. In our numerical implementation, the effect of decimating the flexing-filtering residual was negligible compared to the in-band contribution. Nevertheless, we include it for completeness as it strongly depends on the filter design.

In order to test the validity of the theoretical model given by Eq. (B9), we compare it to numerical simulations that generate time series corresponding to Eq. (B8). We consider white laser frequency noise in $\nu(t)$ with an ASD of $30 \text{ Hz}/\sqrt{\text{Hz}}$ and neglect the central laser frequency of $2.816 \times 10^{14} \text{ Hz}$ (its coupling to the commutator is vanishing). The simulation is performed at a sampling rate of 16 Hz and is then decimated down to 4 Hz. Antialiasing is performed using an FIR designed according to Appendix A. The delay operator is modeled by numerical interpolation of the data using Lagrange polynomials (cf. Sec. III B).

In Fig. 7 we compare the simulated data against the analytical model presented in Eq. (B9). The numerical result is entirely explained by aliasing as the flexing-filtering remains subdominant across the band.

# APPENDIX C: DERIVATION OF TDI RESIDUALS

## 1. Interpolation residual

We follow [32] and implement the postprocessing delay operation in TDI as an interpolation with a fractional delay filter. Lagrange interpolation has demonstrated its suitability in this context, and we thus use it as the baseline in this work. In general, we can model postprocessing delays as

FIR filters. For convenience, we split the delay into an integer shift $j$ and a fractional shift $\epsilon$ ranging from 0 to 1. The latter defines the coefficients of the interpolation kernel $k_n(\epsilon)$ with $-N/2 \le n \le N/2 - 1$, which is convolved with the discrete data samples $\phi_n$,

$$ (\mathcal{D}\phi)_n = \sum_{m=-N/2}^{N/2-1} k_m(\epsilon) \cdot \phi_{n-j-m}. \quad (C1) $$

Before the convolution of Eq. (C1) is performed, the data samples are shifted by the integer shift $j$. The above formula holds for interpolation kernels of even length $N$. For odd $N$, a similar expression can be derived. Using Eq. (7), we find

$$ h_D(\tau) = \sum_{m=-N/2}^{N/2-1} k_m(\epsilon) \cdot \delta(\tau - (j+m)T_s). \quad (C2) $$

We define the additional phase residual caused by the interpolation error as

$$ \delta\phi^D(t) = (\mathcal{D} - \mathbf{D})\phi(t) = \underbrace{\mathbf{D}(\mathbf{D}^{-1}\mathcal{D} - 1)\phi(t)}_{\Delta} \quad (C3) $$

and use Eq. (8) to derive the residual in terms of PSD. We find

$$ S_{\delta\phi}^{D}(f) = \underbrace{|\tilde{h}_{D}(f)e^{2\pi ifd} - 1|^{2}}_{|\Delta|^{2}} S_{\phi}(f). \quad (C4) $$

In general, the interpolation kernel $k_m(\epsilon)$ has to be adjusted for every sample $n$ to account for time-dependent time shifts. Therefore, the flexing arms of LISA will produce a nonstationary interpolation residual as the fractional delay $\epsilon$ scans through different values. At $\epsilon = 0$ and $\epsilon = 1$, the delay is a pure integer shift and the residual vanishes. Assuming that the worst case is obtained for $\epsilon = 0.5$,<fnref n="11" /> we can derive an upper bound for the residual induced by interpolation.

As suggested in [32], and already mentioned at the beginning of this section, a suitable interpolation method is Lagrange interpolation, and we use in this work. The interpolation kernel $k_m$ is derived from fitting a Lagrange polynomial through a set of neighboring samples. This method is known for producing a maximally flat frequency response at DC and is therefore well suited for LISA data processing, as it performs well over the entire LISA band ($10^{-4}$ to 1 Hz) when using high interpolation orders.<fnref n="12" />

<footnote n="11">¹¹A rigorous proof is needed to validate this assumption which was only found to be true *empirically*, in our work, and specifically for Lagrange interpolation.</footnote>

<footnote n="12">¹²Assuming a sampling rate of 4 Hz, typical filter lengths used for postprocessing delays in TDI are in the range 32 to 66.</footnote>
Alternative interpolation kernels are under study that use less coefficients and optimize their performance over the entire band (and not only at DC).

## 2. Ranging residual

Any estimate of the delay $\hat{d}(t)$ differs from the true delay $d(t)$ by a ranging error $r(t)$. We define the corresponding delay operator as

$$ \hat{\mathbf{D}}\phi(t) = \phi(t-d(t)-r(t)) \approx \mathbf{D}\phi(t) - r(t)\mathbf{D}\dot{\phi}(t), \quad (C5) $$

where we have assumed $r(t)$ to be small and performed a series expansion to first order. The ranging residual is then given by

$$ \delta\phi^{\hat{\mathbf{D}}}(t) = (\hat{\mathbf{D}} - \mathbf{D})\phi(t) = -r(t)\mathbf{D}\nu(t) \quad (C6) $$

with $\nu(t) = \dot{\phi}(t)$. Similarly to laser phase or frequency, cf. Eq. (4), we decompose the ranging error $r(t)$ into an out-of-band component $r^o(t)$ and an in-band component $r^e(t)$. Here, $r^o(t)$ absorbs ranging biases that are of the order of 3 ns [13] and might be slowly drifting. The in-band component $r^e(t)$ has a root-mean-squared value of ~100 fs [assuming the PSD in Eq. (53)]. Therefore, $r^e(t) \ll r^o(t)$ and we find as the prominent in-band contributions the coupling of laser noise to the ranging bias and the coupling of ranging noise to the MHz beatnote frequency<fnref n="13" />

$$ r(t)\mathbf{D}\nu(t) \approx r^o(t)\mathbf{D}\nu^e(t) + r^e(t)\mathbf{D}\nu^o(t). \quad (C7) $$

For completeness, in Appendix F we present the coupling of the stochastic in-band ranging error to laser frequency noise.

To simplify the calculations we assume the out-of-band components of the ranging error and the beatnote frequency to be constant. We can readily write down the PSD of Eq. (C6) as

$$ S_{\delta\phi}^{\hat{\mathbf{D}}}(f) = (r^o)^2 S_{\nu}(f) + (\nu^o)^2 S_r(f). \quad (C8) $$

## 3. Comparison with numerical simulations

Once again, we check the validity of our analytic model using simulations performed in units of frequency. Using Eq. (B7), we write down the expression corresponding to the left-hand-side of Eq. (22) for frequency data as

$$ \delta\nu^{\hat{\mathbf{D}}}(t) = (1 - \dot{d}(t)) \cdot \hat{\mathbf{D}}\nu(t) - (1 - \dot{d}(t)) \cdot \mathbf{D}\nu(t). \quad (C9) $$

<figure><img src="image_10.png" /><figcaption>FIG. 8. Processing residuals in application of the fractional delay filter containing ranging errors. We compare the numerically simulated data (solid blue) against the analytical models for the ranging residual caused by ranging noise (dashed-dotted blue) and a constant bias (dashed blue). Additionally, we plot models for the interpolation residual (dashed red) for different interpolation kernel lengths $N = 14, 42$ alongside the picometer reference curve in dashed black for reference.</figcaption></figure>

In order to simulate time series data corresponding to Eq. (C9), we first generate a generic beatnote frequency with a constant offset of 10 MHz and a white laser frequency noise component with an ASD of $30 \text{ Hz}/\sqrt{\text{Hz}}$ at a sampling frequency of 4 Hz. Both the postprocessing delay $\hat{\mathbf{D}}$ and the propagation delay $\mathbf{D}$ are implemented as fractional delay filters. To simulate the latter, we use a very high interpolation order ($N = 502$), such that the interpolation error becomes negligible in comparison to that of the postprocessing delay. The ranging error present in the postprocessing delay is modeled by a bias $B = 10^{-8}$ s, and ranging noise<fnref n="14" /> with an ASD of $10^{-15}$ s/$\sqrt{\text{Hz}(\text{Hz}/f)}$. The nominal value of the delay is taken to be equal to $d = 8.125$ s. This yields the worst case interpolation error, with a fractional part $\epsilon = 0.5$.

In Fig. 8 we compare the PSD of the numerical time series corresponding to Eq. (C9) with the analytic expressions for each of the two components of Eqs. (C4) and (C8) reexpressed for frequency data, i.e.,

$$ S_{\delta\nu}^{\mathbf{D}}(f) = |\tilde{\mathbf{\Delta}}|^2 S_{\nu}(f), \quad (C10a) $$

$$ S_{\delta\nu}^{\hat{\mathbf{D}}}(f) = (2\pi f)^2 \cdot S_{\delta\phi}^{\hat{\mathbf{D}}}(f). \quad (C10b) $$

The interpolation residual (dashed teal) is strongly dependent on the length $N$ of the interpolation kernel.

<footnote n="13">¹³In the notation of Eq. (30), $r^o$ corresponds to $B_{ij}$ while $r^e(t)$ corresponds to the combination $-\mathbf{D}_{ij}M_j(t) + M_i(t)$.</footnote>

<footnote n="14">¹⁴We choose a red tilt for ranging noise to be easily distinguishable from coupling of the ranging bias to laser frequency noise. The level is comparable with realistic models of modulation noise discussed in Sec. III.</footnote>
Alternative interpolation kernels are under study that use less coefficients and optimize their performance over the entire band (and not only at DC).

## 2. Ranging residual

Any estimate of the delay $\hat{d}(t)$ differs from the true delay $d(t)$ by a ranging error $r(t)$. We define the corresponding delay operator as

$$ \hat{\mathbf{D}}\phi(t) = \phi(t-d(t)-r(t)) \approx \mathbf{D}\phi(t) - r(t)\mathbf{D}\dot{\phi}(t), \quad (C5) $$

where we have assumed $r(t)$ to be small and performed a series expansion to first order. The ranging residual is then given by

$$ \delta\phi^{\hat{\mathbf{D}}}(t) = (\hat{\mathbf{D}} - \mathbf{D})\phi(t) = -r(t)\mathbf{D}\nu(t) \quad (C6) $$

with $\nu(t) = \dot{\phi}(t)$. Similarly to laser phase or frequency, cf. Eq. (4), we decompose the ranging error $r(t)$ into an out-of-band component $r^o(t)$ and an in-band component $r^e(t)$. Here, $r^o(t)$ absorbs ranging biases that are of the order of 3 ns [13] and might be slowly drifting. The in-band component $r^e(t)$ has a root-mean-squared value of ~100 fs [assuming the PSD in Eq. (53)]. Therefore, $r^e(t) \ll r^o(t)$ and we find as the prominent in-band contributions the coupling of laser noise to the ranging bias and the coupling of ranging noise to the MHz beatnote frequency<fnref n="13" />

$$ r(t)\mathbf{D}\nu(t) \approx r^o(t)\mathbf{D}\nu^e(t) + r^e(t)\mathbf{D}\nu^o(t). \quad (C7) $$

For completeness, in Appendix F we present the coupling of the stochastic in-band ranging error to laser frequency noise.

To simplify the calculations we assume the out-of-band components of the ranging error and the beatnote frequency to be constant. We can readily write down the PSD of Eq. (C6) as

$$ S_{\delta\phi}^{\hat{\mathbf{D}}}(f) = (r^o)^2 S_{\nu}(f) + (\nu^o)^2 S_r(f). \quad (C8) $$

## 3. Comparison with numerical simulations

Once again, we check the validity of our analytic model using simulations performed in units of frequency. Using Eq. (B7), we write down the expression corresponding to the left-hand-side of Eq. (22) for frequency data as

$$ \delta\nu^{\hat{\mathbf{D}}}(t) = (1 - \dot{d}(t)) \cdot \hat{\mathbf{D}}\nu(t) - (1 - \dot{d}(t)) \cdot \mathbf{D}\nu(t). \quad (C9) $$

<figure><img src="image_10.png" /><figcaption>FIG. 8. Processing residuals in application of the fractional delay filter containing ranging errors. We compare the numerically simulated data (solid blue) against the analytical models for the ranging residual caused by ranging noise (dashed-dotted blue) and a constant bias (dashed blue). Additionally, we plot models for the interpolation residual (dashed red) for different interpolation kernel lengths $N = 14, 42$ alongside the picometer reference curve in dashed black for reference.</figcaption></figure>

In order to simulate time series data corresponding to Eq. (C9), we first generate a generic beatnote frequency with a constant offset of 10 MHz and a white laser frequency noise component with an ASD of $30 \text{ Hz}/\sqrt{\text{Hz}}$ at a sampling frequency of 4 Hz. Both the postprocessing delay $\hat{\mathbf{D}}$ and the propagation delay $\mathbf{D}$ are implemented as fractional delay filters. To simulate the latter, we use a very high interpolation order ($N = 502$), such that the interpolation error becomes negligible in comparison to that of the postprocessing delay. The ranging error present in the postprocessing delay is modeled by a bias $B = 10^{-8}$ s, and ranging noise<fnref n="14" /> with an ASD of $10^{-15}$ s/$\sqrt{\text{Hz}(\text{Hz}/f)}$. The nominal value of the delay is taken to be equal to $d = 8.125$ s. This yields the worst case interpolation error, with a fractional part $\epsilon = 0.5$.

In Fig. 8 we compare the PSD of the numerical time series corresponding to Eq. (C9) with the analytic expressions for each of the two components of Eqs. (C4) and (C8) reexpressed for frequency data, i.e.,

$$ S_{\delta\nu}^{\mathbf{D}}(f) = |\tilde{\mathbf{\Delta}}|^2 S_{\nu}(f), \quad (C10a) $$

$$ S_{\delta\nu}^{\hat{\mathbf{D}}}(f) = (2\pi f)^2 \cdot S_{\delta\phi}^{\hat{\mathbf{D}}}(f). \quad (C10b) $$

The interpolation residual (dashed teal) is strongly dependent on the length $N$ of the interpolation kernel.

<footnote n="13">¹³In the notation of Eq. (30), $r^o$ corresponds to $B_{ij}$ while $r^e(t)$ corresponds to the combination $-\mathbf{D}_{ij}M_j(t) + M_i(t)$.</footnote>

<footnote n="14">¹⁴We choose a red tilt for ranging noise to be easily distinguishable from coupling of the ranging bias to laser frequency noise. The level is comparable with realistic models of modulation noise discussed in Sec. III.</footnote>
For this reason, we show the interpolation residual obtained for N = 14, and 42. The postprocessing delay used in the numerical implementation is performed with an interpolation order N = 42. As shown in the figure, the model of Eq. (C10) agrees with the data over all frequencies.

## APPENDIX D: RANGING PROCESSING

The ranging processing discussed in this section is mostly adopted from [13] and we introduce only minor changes to the algorithm. The main differences are that we write it down in units of phase and reformulate the suppression of modulation noise contributions originating from right-handed optical benches. Additionally, we assume that the sideband phases are read out using a feed-forward scheme. This measure avoids tracking the carrier phase redundantly and also accounts for the difference in modulation frequencies on each MOSA. Therefore, the sideband phases of the interspacecraft and reference interferometer read

$$isi_{ij}^{\text{sb}} = \mathbf{D}_{ij}\phi_{ji}^{\text{m}} - \phi_{ij}^{\text{m}} - (\nu_{ji}^{\text{m}} - \nu_{ij}^{\text{m}})t, \quad (\text{D1a})$$

$$rfi_{ij}^{\text{sb}} = \phi_{ik}^{\text{m}} - \phi_{ij}^{\text{m}} - (\nu_{ik}^{\text{m}} - \nu_{ij}^{\text{m}})t, \quad (\text{D1b})$$

where the phase of the modulation is given by

$$\phi_{ij}^{\text{m}}(t) = \nu_{ij}^{\text{m}} \cdot (t + M_{ij}(t)). \quad (\text{D2})$$

Here, $M_{ij}(t)$ accounts for timing jitter due to modulation noise. Inserting this definition into Eq. (D1) yields

$$isi_{ij}^{\text{sb}} = -\nu_{ji}^{\text{m}}d_{ij} + \mathbf{D}_{ij}\nu_{ji}^{\text{m}}M_{ji} - \nu_{ij}^{\text{m}}M_{ij}, \quad (\text{D3a})$$

$$rfi_{ij}^{\text{sb}} = \nu_{ik}^{\text{m}}M_{ik} - \nu_{ij}^{\text{m}}M_{ij}, \quad (\text{D3b})$$

which has a similar algebraic structure to the carrier phases [cf. Eq. (1)] where the product $\nu_{ij}^{\text{m}}M_{ij}$ takes the place of the laser phase $\phi_{ij}$. We now use the definition of the intermed diary variable $\eta_{ij}$ for the sidebands<fnref n="15" /> [see Eq. (13a)] to cancel modulation noise contributions stemming from right-handed MOSAs,

$$\eta_{ij}^{\text{sb}} = -\nu_{ji}^{\text{m}}d_{ij} + \mathbf{D}_{ij}\nu_{jk}^{\text{m}}M_{j} - \nu_{ij}^{\text{m}}M_{i}. \quad (\text{D4})$$

Here, we use the shorthand notation $M_1 = M_{12}$, $M_2 = M_{23}$ and $M_3 = M_{31}$.

Finally, the variables $\eta_{ij}^{\text{sb}}$ have to be scaled by the respective modulation frequency and multiplied by -1 in order to yield a low noise estimate of the delay

$$\hat{d}_{ij} = -\frac{\eta_{ij}^{\text{sb}}}{\nu_{ji}^{\text{m}}} = d_{ij} - \mathbf{D}_{ij} \frac{\nu_{jk}^{\text{m}}}{\nu_{ji}^{\text{m}}} M_j + \frac{\nu_{ij}^{\text{m}}}{\nu_{ji}^{\text{m}}} M_i. \quad (\text{D5})$$

Assuming that modulation frequencies only differ by a fraction of a percent we can approximate the above expression by setting the ratios of frequencies to 1. Doing so, we recover the stochastic component in Eq. (30).

In practice, we process data in frequency units. Hence, the above procedure is rewritten in frequency by taking a global time derivative. By doing so, Eq. (D5) yields an accurate measurement of the delay derivative $\hat{d}_{ij}$, which then needs to be integrated to recover the delay itself, required for TDI. The integration constant is derived from the PRN ranging measurement [13] and is the origin of the appearance of the bias $B_{ij}$ in Eq. (30).

## APPENDIX E: DIFFERENTIAL DOPPLER SHIFT

As explained in Sec. III, any TDI combination representing a virtual two-beam interferometer does not cancel the laser phase perfectly for flexing arms but is limited by a residual given by the delay commutator. The origin of this residual is the travel time difference $\Delta d$ between the two virtual beams. The deterministic component of this residual can be calculated and subtracted from the TDI observable [see Eq. (35)]. Here, we present an efficient scheme to calculate $\Delta d$. First we recognize that the delay $d_{ij}$ can be written as

$$d_{ij}(t) = t - \mathbf{D}_{ij}t, \quad (\text{E1})$$

which has the same algebraic structure as $\eta_{ij}$ up to a sign. Here, the time argument $t$ takes the place of the laser phase $\phi_i$ and $\phi_j$ [cf. Eqs. (1) and (13)]. Using $\eta_{ij} = d_{ij}$ as inputs to TDI yields the travel time difference $\Delta d$ as

$$X = -[\mathbf{D}_A, \mathbf{D}_B]t = d_{AB} - d_{BA} = \Delta d, \quad (\text{E2})$$

where $AB$ and $BA$ denote the paths of the counter-propagating beams of an arbitrary TDI combination $X$ representing a two-beam interferometer.

As processing is performed in frequency units in this paper we are more interested in the derivative of $\Delta d$. To avoid numerical problems we thus operate on the delay derivatives $\hat{d}_{ij}$ directly and form $\Delta \hat{d} = \dot{X}$ following the procedure explained in [6].

## APPENDIX F: COUPLING OF RANGING NOISE TO LASER NOISE

In Sec. III B, we neglect the coupling of the stochastic component of the ranging error to laser noise as it appears to be much weaker compared with the coupling to the MHz beatnote frequency. However, for the sake of completeness and as it becomes relevant in processing pipelines where

<footnote n="15">¹⁵The postprocessing delays required to calculate $\eta_{ij}^{\text{sb}}$ have much more relaxed error requirements than is the case for laser noise cancellation. We can therefore use the delay estimates from the PRN.</footnote>
For this reason, we show the interpolation residual obtained for N = 14, and 42. The postprocessing delay used in the numerical implementation is performed with an interpolation order N = 42. As shown in the figure, the model of Eq. (C10) agrees with the data over all frequencies.

## APPENDIX D: RANGING PROCESSING

The ranging processing discussed in this section is mostly adopted from [13] and we introduce only minor changes to the algorithm. The main differences are that we write it down in units of phase and reformulate the suppression of modulation noise contributions originating from right-handed optical benches. Additionally, we assume that the sideband phases are read out using a feed-forward scheme. This measure avoids tracking the carrier phase redundantly and also accounts for the difference in modulation frequencies on each MOSA. Therefore, the sideband phases of the interspacecraft and reference interferometer read

$$isi_{ij}^{\text{sb}} = \mathbf{D}_{ij}\phi_{ji}^{\text{m}} - \phi_{ij}^{\text{m}} - (\nu_{ji}^{\text{m}} - \nu_{ij}^{\text{m}})t, \quad (\text{D1a})$$

$$rfi_{ij}^{\text{sb}} = \phi_{ik}^{\text{m}} - \phi_{ij}^{\text{m}} - (\nu_{ik}^{\text{m}} - \nu_{ij}^{\text{m}})t, \quad (\text{D1b})$$

where the phase of the modulation is given by

$$\phi_{ij}^{\text{m}}(t) = \nu_{ij}^{\text{m}} \cdot (t + M_{ij}(t)). \quad (\text{D2})$$

Here, $M_{ij}(t)$ accounts for timing jitter due to modulation noise. Inserting this definition into Eq. (D1) yields

$$isi_{ij}^{\text{sb}} = -\nu_{ji}^{\text{m}}d_{ij} + \mathbf{D}_{ij}\nu_{ji}^{\text{m}}M_{ji} - \nu_{ij}^{\text{m}}M_{ij}, \quad (\text{D3a})$$

$$rfi_{ij}^{\text{sb}} = \nu_{ik}^{\text{m}}M_{ik} - \nu_{ij}^{\text{m}}M_{ij}, \quad (\text{D3b})$$

which has a similar algebraic structure to the carrier phases [cf. Eq. (1)] where the product $\nu_{ij}^{\text{m}}M_{ij}$ takes the place of the laser phase $\phi_{ij}$. We now use the definition of the intermed diary variable $\eta_{ij}$ for the sidebands<fnref n="15" /> [see Eq. (13a)] to cancel modulation noise contributions stemming from right-handed MOSAs,

$$\eta_{ij}^{\text{sb}} = -\nu_{ji}^{\text{m}}d_{ij} + \mathbf{D}_{ij}\nu_{jk}^{\text{m}}M_{j} - \nu_{ij}^{\text{m}}M_{i}. \quad (\text{D4})$$

Here, we use the shorthand notation $M_1 = M_{12}$, $M_2 = M_{23}$ and $M_3 = M_{31}$.

Finally, the variables $\eta_{ij}^{\text{sb}}$ have to be scaled by the respective modulation frequency and multiplied by -1 in order to yield a low noise estimate of the delay

$$\hat{d}_{ij} = -\frac{\eta_{ij}^{\text{sb}}}{\nu_{ji}^{\text{m}}} = d_{ij} - \mathbf{D}_{ij} \frac{\nu_{jk}^{\text{m}}}{\nu_{ji}^{\text{m}}} M_j + \frac{\nu_{ij}^{\text{m}}}{\nu_{ji}^{\text{m}}} M_i. \quad (\text{D5})$$

Assuming that modulation frequencies only differ by a fraction of a percent we can approximate the above expression by setting the ratios of frequencies to 1. Doing so, we recover the stochastic component in Eq. (30).

In practice, we process data in frequency units. Hence, the above procedure is rewritten in frequency by taking a global time derivative. By doing so, Eq. (D5) yields an accurate measurement of the delay derivative $\hat{d}_{ij}$, which then needs to be integrated to recover the delay itself, required for TDI. The integration constant is derived from the PRN ranging measurement [13] and is the origin of the appearance of the bias $B_{ij}$ in Eq. (30).

## APPENDIX E: DIFFERENTIAL DOPPLER SHIFT

As explained in Sec. III, any TDI combination representing a virtual two-beam interferometer does not cancel the laser phase perfectly for flexing arms but is limited by a residual given by the delay commutator. The origin of this residual is the travel time difference $\Delta d$ between the two virtual beams. The deterministic component of this residual can be calculated and subtracted from the TDI observable [see Eq. (35)]. Here, we present an efficient scheme to calculate $\Delta d$. First we recognize that the delay $d_{ij}$ can be written as

$$d_{ij}(t) = t - \mathbf{D}_{ij}t, \quad (\text{E1})$$

which has the same algebraic structure as $\eta_{ij}$ up to a sign. Here, the time argument $t$ takes the place of the laser phase $\phi_i$ and $\phi_j$ [cf. Eqs. (1) and (13)]. Using $\eta_{ij} = d_{ij}$ as inputs to TDI yields the travel time difference $\Delta d$ as

$$X = -[\mathbf{D}_A, \mathbf{D}_B]t = d_{AB} - d_{BA} = \Delta d, \quad (\text{E2})$$

where $AB$ and $BA$ denote the paths of the counter-propagating beams of an arbitrary TDI combination $X$ representing a two-beam interferometer.

As processing is performed in frequency units in this paper we are more interested in the derivative of $\Delta d$. To avoid numerical problems we thus operate on the delay derivatives $\hat{d}_{ij}$ directly and form $\Delta \hat{d} = \dot{X}$ following the procedure explained in [6].

## APPENDIX F: COUPLING OF RANGING NOISE TO LASER NOISE

In Sec. III B, we neglect the coupling of the stochastic component of the ranging error to laser noise as it appears to be much weaker compared with the coupling to the MHz beatnote frequency. However, for the sake of completeness and as it becomes relevant in processing pipelines where

<footnote n="15">¹⁵The postprocessing delays required to calculate $\eta_{ij}^{\text{sb}}$ have much more relaxed error requirements than is the case for laser noise cancellation. We can therefore use the delay estimates from the PRN.</footnote>
one removes the phase ramp operates on the fluctuations directly we present the coupling mechanism below.

To suppress all other laser noise couplings in the final TDI combination we consider a setup where we have already removed the phase ramp from the interferometric measurements such that they only track the differential phase noise $p_{ij}$ of the six lasers,

$$isi_{ij}(t) = \mathbf{D}_{ij}p_{ji}(t) - p_{ij}(t), \quad (\text{F1a})$$

$$rfi_{ij}(t) = p_{ik}(t) - p_{ij}(t). \quad (\text{F1b})$$

For simplicity we omit the antialiasing filtering and decimation that was considered in Eq. (1).

Next, we insert Eq. (F1) into Eq. (13a) where the postprocessing delay $\hat{D}$ is used and only accounts for a stochastic ranging error $r(t)$. We can express $\eta_{ij}$ for the left and right-handed MOSAs as

$$\eta_{ij} = \hat{\mathbf{D}}_{ij} p_j - p_i - (\hat{\mathbf{D}}_{ij} - \mathbf{D}_{ij}) p_{ji}. \quad (\text{F2})$$

Here, we use the short-hand notation $p_1 = p_{12}$, $p_2 = p_{23}$ and $p_3 = p_{31}$. We recognize that the last term in Eq. (F2) is already a laser noise residual and we neglect higher order couplings in the following. Then, we use the intermediary variables, $\pi_{ij}$, $\rho_{ij}$ and $\sigma_{ij}$ defined in Eqs. (13) and (14) to find the total laser noise residual in the second-generation Michelson combination $X_2$. It consists of the residual in Eq. (F2) that is propagated through TDI as well as the commutator of postprocessing delay operators

$$[\hat{\mathbf{D}}_{13121}, \hat{\mathbf{D}}_{12131}] = [\mathbf{D}_{13121}, \mathbf{D}_{12131}] + (\hat{\mathbf{D}}_{131212131} - \mathbf{D}_{131212131}) - (\hat{\mathbf{D}}_{121313121} - \mathbf{D}_{121313121}) \quad (\text{F3})$$

applied to $p_1$. Here, we split the commutator into the “usual” delay commutator and two additional terms that produce further laser noise residuals. For equal arms, the full residual in $X_2$ reads

$$\delta X_2^\hat{D} = \mathbf{C}\{\mathbf{D}(r_{31}\mathbf{D}\dot{p}_{13}) + (r_{13}\mathbf{D}\dot{p}_{31}) - \mathbf{D}(r_{21}\mathbf{D}\dot{p}_{12}) - (r_{12}\mathbf{D}\dot{p}_{21})\} + \mathbf{C}\{r_{12} + \mathbf{D}r_{21} - r_{13} - \mathbf{D}r_{31}\} \cdot \mathbf{D}^8 \dot{p}_1, \quad (\text{F4})$$

where we have used Eq. (C6) and have approximated nested ranging noise as

$$r_{ijikiki} \approx r_{iji} + \mathbf{D}^2 r_{iki} + \mathbf{D}^4 r_{iki} + \mathbf{D}^6 r_{iji}, \quad (\text{F5a})$$

$$r_{iji} \approx r_{ij} + \mathbf{D} r_{ji}. \quad (\text{F5b})$$

Finally, we compute the PSD of Eq. (F4) by assuming that all laser and ranging noise terms are uncorrelated and have identical noise properties. We find

$$S_{\delta X_2}^{\hat{D}}(f) = 4|\tilde{\mathbf{C}}|^2 (S_r(f') * S_{\hat{p}}(f'))(f) + 4(|\tilde{\mathbf{C}}|^2 S_r(f') * S_{\hat{p}}(f'))(f), \quad (\text{F6})$$

where we have neglected any cross terms between the first two lines and the last line in Eq. (F4). The * sign denotes convolution (in frequency domain) which stems from the time domain products of ranging and laser noise contributions.

[1] P. Amaro-Seoane et al., Laser interferometer space antenna, arXiv:1702.00786.

[2] M. Tinto and J. W. Armstrong, Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation, Phys. Rev. D 59, 102003 (1999).

[3] J. W. Armstrong, F. B. Estabrook, and M. Tinto, Time-delay interferometry for space-based gravitational wave searches, Astrophys. J. 527, 814 (1999).

[4] D. A. Shaddock, M. Tinto, F. B. Estabrook, and J. W. Armstrong, Data combinations accounting for LISA spacecraft motion, Phys. Rev. D 68, 061303(R) (2003).

[5] M. Tinto, F. B. Estabrook, and J. W. Armstrong, Time delay interferometry with moving spacecraft arrays, Phys. Rev. D 69, 082001 (2004).

[6] J.-B. Bayle, O. Hartwig, and M. Staab, Adapting time-delay interferometry for LISA data in frequency, Phys. Rev. D 104, 023006 (2021).

[7] M. Vallisneri, J.-B. Bayle, S. Babak, and A. Petiteau, Time-delay interferometry without delays, Phys. Rev. D 103, 082001 (2021).

[8] N. Houba, S. Delchambre, G. Hechenblaikner, T. Ziegler, and W. Fichter, Time-delay interferometry infinity for tilt-to-length noise estimation in LISA, Classical Quantum Gravity 40, 107001 (2023).

[9] Q. Baghi, J. Baker, J. Slutsky, and J. I. Thorpe, Model-independent time-delay interferometry based on principal component analysis, Phys. Rev. D 104, 122001 (2021).

[10] Q. Baghi, J. G. Baker, J. Slutsky, and J. I. Thorpe, Fully data-driven time-delay interferometry with time-varying delays, Ann. Phys. (Amsterdam) 2023, 2200447 (2022).

[11] D. Q. Nam, Y. Lemiere, A. Petiteau, J.-B. Bayle, O. Hartwig, J. Martino, and M. Staab, TDI noises transfer functions for LISA, Phys. Rev. D 108, 082004 (2023).
one removes the phase ramp operates on the fluctuations directly we present the coupling mechanism below.

To suppress all other laser noise couplings in the final TDI combination we consider a setup where we have already removed the phase ramp from the interferometric measurements such that they only track the differential phase noise $p_{ij}$ of the six lasers,

$$isi_{ij}(t) = \mathbf{D}_{ij}p_{ji}(t) - p_{ij}(t), \quad (\text{F1a})$$

$$rfi_{ij}(t) = p_{ik}(t) - p_{ij}(t). \quad (\text{F1b})$$

For simplicity we omit the antialiasing filtering and decimation that was considered in Eq. (1).

Next, we insert Eq. (F1) into Eq. (13a) where the postprocessing delay $\hat{D}$ is used and only accounts for a stochastic ranging error $r(t)$. We can express $\eta_{ij}$ for the left and right-handed MOSAs as

$$\eta_{ij} = \hat{\mathbf{D}}_{ij} p_j - p_i - (\hat{\mathbf{D}}_{ij} - \mathbf{D}_{ij}) p_{ji}. \quad (\text{F2})$$

Here, we use the short-hand notation $p_1 = p_{12}$, $p_2 = p_{23}$ and $p_3 = p_{31}$. We recognize that the last term in Eq. (F2) is already a laser noise residual and we neglect higher order couplings in the following. Then, we use the intermediary variables, $\pi_{ij}$, $\rho_{ij}$ and $\sigma_{ij}$ defined in Eqs. (13) and (14) to find the total laser noise residual in the second-generation Michelson combination $X_2$. It consists of the residual in Eq. (F2) that is propagated through TDI as well as the commutator of postprocessing delay operators

$$[\hat{\mathbf{D}}_{13121}, \hat{\mathbf{D}}_{12131}] = [\mathbf{D}_{13121}, \mathbf{D}_{12131}] + (\hat{\mathbf{D}}_{131212131} - \mathbf{D}_{131212131}) - (\hat{\mathbf{D}}_{121313121} - \mathbf{D}_{121313121}) \quad (\text{F3})$$

applied to $p_1$. Here, we split the commutator into the “usual” delay commutator and two additional terms that produce further laser noise residuals. For equal arms, the full residual in $X_2$ reads

$$\delta X_2^\hat{D} = \mathbf{C}\{\mathbf{D}(r_{31}\mathbf{D}\dot{p}_{13}) + (r_{13}\mathbf{D}\dot{p}_{31}) - \mathbf{D}(r_{21}\mathbf{D}\dot{p}_{12}) - (r_{12}\mathbf{D}\dot{p}_{21})\} + \mathbf{C}\{r_{12} + \mathbf{D}r_{21} - r_{13} - \mathbf{D}r_{31}\} \cdot \mathbf{D}^8 \dot{p}_1, \quad (\text{F4})$$

where we have used Eq. (C6) and have approximated nested ranging noise as

$$r_{ijikiki} \approx r_{iji} + \mathbf{D}^2 r_{iki} + \mathbf{D}^4 r_{iki} + \mathbf{D}^6 r_{iji}, \quad (\text{F5a})$$

$$r_{iji} \approx r_{ij} + \mathbf{D} r_{ji}. \quad (\text{F5b})$$

Finally, we compute the PSD of Eq. (F4) by assuming that all laser and ranging noise terms are uncorrelated and have identical noise properties. We find

$$S_{\delta X_2}^{\hat{D}}(f) = 4|\tilde{\mathbf{C}}|^2 (S_r(f') * S_{\hat{p}}(f'))(f) + 4(|\tilde{\mathbf{C}}|^2 S_r(f') * S_{\hat{p}}(f'))(f), \quad (\text{F6})$$

where we have neglected any cross terms between the first two lines and the last line in Eq. (F4). The * sign denotes convolution (in frequency domain) which stems from the time domain products of ranging and laser noise contributions.

[1] P. Amaro-Seoane et al., Laser interferometer space antenna, arXiv:1702.00786.

[2] M. Tinto and J. W. Armstrong, Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation, Phys. Rev. D 59, 102003 (1999).

[3] J. W. Armstrong, F. B. Estabrook, and M. Tinto, Time-delay interferometry for space-based gravitational wave searches, Astrophys. J. 527, 814 (1999).

[4] D. A. Shaddock, M. Tinto, F. B. Estabrook, and J. W. Armstrong, Data combinations accounting for LISA spacecraft motion, Phys. Rev. D 68, 061303(R) (2003).

[5] M. Tinto, F. B. Estabrook, and J. W. Armstrong, Time delay interferometry with moving spacecraft arrays, Phys. Rev. D 69, 082001 (2004).

[6] J.-B. Bayle, O. Hartwig, and M. Staab, Adapting time-delay interferometry for LISA data in frequency, Phys. Rev. D 104, 023006 (2021).

[7] M. Vallisneri, J.-B. Bayle, S. Babak, and A. Petiteau, Time-delay interferometry without delays, Phys. Rev. D 103, 082001 (2021).

[8] N. Houba, S. Delchambre, G. Hechenblaikner, T. Ziegler, and W. Fichter, Time-delay interferometry infinity for tilt-to-length noise estimation in LISA, Classical Quantum Gravity 40, 107001 (2023).

[9] Q. Baghi, J. Baker, J. Slutsky, and J. I. Thorpe, Model-independent time-delay interferometry based on principal component analysis, Phys. Rev. D 104, 122001 (2021).

[10] Q. Baghi, J. G. Baker, J. Slutsky, and J. I. Thorpe, Fully data-driven time-delay interferometry with time-varying delays, Ann. Phys. (Amsterdam) 2023, 2200447 (2022).

[11] D. Q. Nam, Y. Lemiere, A. Petiteau, J.-B. Bayle, O. Hartwig, J. Martino, and M. Staab, TDI noises transfer functions for LISA, Phys. Rev. D 108, 082004 (2023).
[12] O. Hartwig and J.-B. Bayle, Clock-jitter reduction in LISA time-delay interferometry combinations, Phys. Rev. D **103**, 123027 (2021).

[13] O. Hartwig, J.-B. Bayle, M. Staab, A. Hees, M. Lilley, and P. Wolf, Time-delay interferometry without clock synchronization, Phys. Rev. D **105**, 122008 (2022).

[14] S. Paczkowski, R. Giusteri, M. Hewitson, N. Karnesis, E. D. Fitzsimons, G. Wanner, and G. Heinzel, Postprocessing subtraction of tilt-to-length noise in LISA, Phys. Rev. D **106**, 042005 (2022).

[15] J.-B. Bayle, M. Lilley, A. Petiteau, and H. Halloin, Effect of filters on the time-delay interferometry residual laser noise for LISA, Phys. Rev. D **99**, 084023 (2019).

[16] O. Hartwig, Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna, Ph.D. thesis, Leibniz University, Hannover, 2021.

[17] G. de Vine, B. Ware, K. McKenzie, R. E. Spero, W. M. Klipstein, and D. A. Shaddock, Experimental demonstration of time-delay interferometry for the Laser Interferometer Space Antenna, Phys. Rev. Lett. **104**, 211103 (2010).

[18] S. J. Mitryk, V. Wand, and G. Mueller, Verification of time-delay interferometry techniques using the University of Florida LISA interferometry simulator, Classical Quantum Gravity **27**, 084012 (2010).

[19] S. J. Mitryk, G. Mueller, and J. Sanjuan, Hardware-based demonstration of time-delay interferometry and TDI-ranging with spacecraft motion effects, Phys. Rev. D **86**, 122006 (2012).

[20] K. Yamamoto, C. Vorndamme, O. Hartwig, M. Staab, T. S. Schwarze, and G. Heinzel, Experimental verification of intersatellite clock synchronization at LISA performance levels, Phys. Rev. D **105**, 042009 (2022).

[21] P. Gruning, H. Halloin, P. Prat, S. Baron, J. Brossard, C. Buy, and A. Petiteau, Progress towards an electro-optical simulator for space based long arms interferometers, Exp. Astron. **39**, 281 (2015).

[22] L. Vidal, Validation expérimentale des performances interférométriques de LISA, Ph.D. thesis, Université Paris Cité, Paris, 2023.

[23] G. Heinzel, LISA frequency planning, Max-Planck Institute for Gravitational Physics, Albert Einstein Institute Technical Report No. LISA-AEI-INST-TN-002, 2018.

[24] J.-B. Bayle and O. Hartwig, Unified model for the LISA measurements and instrument simulations, Phys. Rev. D **107**, 083019 (2023).

[25] G. M. Jenkins and D. G. Watts, *Spectral Analysis and its Applications*, Holden-Day Series in Time Series Analysis (Holden-Day, San Francisco, 1968).

[26] J. N. Reinhardt, M. Staab, K. Yamamoto, J.-B. Bayle, A. Hees, O. Hartwig, K. Wiesner, and G. Heinzel, Ranging sensor fusion in LISA data processing: Treatment of ambiguities, noise, and on-board delays in LISA ranging observables, Phys. Rev. D **109**, 022004 (2024).

[27] N. J. Cornish and R. W. Hellings, The effects of orbital motion on LISA time delay interferometry, Classical Quantum Gravity **20**, 4851 (2003).

[28] J.-B. Bayle, O. Hartwig, and M. Staab, Lisa instrument (2022), https://zenodo.org/records/7071251.

[29] M. Staab, J.-B. Bayle, and O. Hartwig, Pytdi (2022), https://zenodo.org/records/6867012.

[30] P. Virtanen et al., SciPy 1.0-Fundamental algorithms for scientific computing in Python, Nat. Methods **17**, 261 (2020).

[31] S. Barke, Inter-spacecraft frequency distribution for future gravitational wave observatories, Ph.D. thesis, Leibniz University, Hannover, 2015.

[32] D. A. Shaddock, B. Ware, R. E. Spero, and M. Vallisneri, Post-processed time-delay interferometry for LISA, Phys. Rev. D **70**, 081101(R) (2004).

[33] T. Parks and J. McClellan, Chebyshev approximation for nonrecursive digital filters with linear phase, IEEE Trans. Circuit Theory **19**, 189 (1972).

[34] J. McClellan, T. Parks, and L. Rabiner, A computer program for designing optimum fir linear phase digital filters, IEEE Trans. Audio Electroacoust. **21**, 506 (1973).
[12] O. Hartwig and J.-B. Bayle, Clock-jitter reduction in LISA time-delay interferometry combinations, Phys. Rev. D **103**, 123027 (2021).

[13] O. Hartwig, J.-B. Bayle, M. Staab, A. Hees, M. Lilley, and P. Wolf, Time-delay interferometry without clock synchronization, Phys. Rev. D **105**, 122008 (2022).

[14] S. Paczkowski, R. Giusteri, M. Hewitson, N. Karnesis, E. D. Fitzsimons, G. Wanner, and G. Heinzel, Postprocessing subtraction of tilt-to-length noise in LISA, Phys. Rev. D **106**, 042005 (2022).

[15] J.-B. Bayle, M. Lilley, A. Petiteau, and H. Halloin, Effect of filters on the time-delay interferometry residual laser noise for LISA, Phys. Rev. D **99**, 084023 (2019).

[16] O. Hartwig, Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna, Ph.D. thesis, Leibniz University, Hannover, 2021.

[17] G. de Vine, B. Ware, K. McKenzie, R. E. Spero, W. M. Klipstein, and D. A. Shaddock, Experimental demonstration of time-delay interferometry for the Laser Interferometer Space Antenna, Phys. Rev. Lett. **104**, 211103 (2010).

[18] S. J. Mitryk, V. Wand, and G. Mueller, Verification of time-delay interferometry techniques using the University of Florida LISA interferometry simulator, Classical Quantum Gravity **27**, 084012 (2010).

[19] S. J. Mitryk, G. Mueller, and J. Sanjuan, Hardware-based demonstration of time-delay interferometry and TDI-ranging with spacecraft motion effects, Phys. Rev. D **86**, 122006 (2012).

[20] K. Yamamoto, C. Vorndamme, O. Hartwig, M. Staab, T. S. Schwarze, and G. Heinzel, Experimental verification of intersatellite clock synchronization at LISA performance levels, Phys. Rev. D **105**, 042009 (2022).

[21] P. Gruning, H. Halloin, P. Prat, S. Baron, J. Brossard, C. Buy, and A. Petiteau, Progress towards an electro-optical simulator for space based long arms interferometers, Exp. Astron. **39**, 281 (2015).

[22] L. Vidal, Validation expérimentale des performances interférométriques de LISA, Ph.D. thesis, Université Paris Cité, Paris, 2023.

[23] G. Heinzel, LISA frequency planning, Max-Planck Institute for Gravitational Physics, Albert Einstein Institute Technical Report No. LISA-AEI-INST-TN-002, 2018.

[24] J.-B. Bayle and O. Hartwig, Unified model for the LISA measurements and instrument simulations, Phys. Rev. D **107**, 083019 (2023).

[25] G. M. Jenkins and D. G. Watts, *Spectral Analysis and its Applications*, Holden-Day Series in Time Series Analysis (Holden-Day, San Francisco, 1968).

[26] J. N. Reinhardt, M. Staab, K. Yamamoto, J.-B. Bayle, A. Hees, O. Hartwig, K. Wiesner, and G. Heinzel, Ranging sensor fusion in LISA data processing: Treatment of ambiguities, noise, and on-board delays in LISA ranging observables, Phys. Rev. D **109**, 022004 (2024).

[27] N. J. Cornish and R. W. Hellings, The effects of orbital motion on LISA time delay interferometry, Classical Quantum Gravity **20**, 4851 (2003).

[28] J.-B. Bayle, O. Hartwig, and M. Staab, Lisa instrument (2022), https://zenodo.org/records/7071251.

[29] M. Staab, J.-B. Bayle, and O. Hartwig, Pytdi (2022), https://zenodo.org/records/6867012.

[30] P. Virtanen et al., SciPy 1.0-Fundamental algorithms for scientific computing in Python, Nat. Methods **17**, 261 (2020).

[31] S. Barke, Inter-spacecraft frequency distribution for future gravitational wave observatories, Ph.D. thesis, Leibniz University, Hannover, 2015.

[32] D. A. Shaddock, B. Ware, R. E. Spero, and M. Vallisneri, Post-processed time-delay interferometry for LISA, Phys. Rev. D **70**, 081101(R) (2004).

[33] T. Parks and J. McClellan, Chebyshev approximation for nonrecursive digital filters with linear phase, IEEE Trans. Circuit Theory **19**, 189 (1972).

[34] J. McClellan, T. Parks, and L. Rabiner, A computer program for designing optimum fir linear phase digital filters, IEEE Trans. Audio Electroacoust. **21**, 506 (1973).