additions

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

Martin Staab  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), D-30167 Hanover, Germany  Leibniz Universität Hannover, D-30167
Hanover, Germany  Marc Lilley  SYRTE, Observatoire de Paris, Université PSL,
CNRS, Sorbonne Université, LNE, 61 avenue de l’Observatoire 75014 Paris,
France  Jean-Baptiste Bayle  University of Glasgow, Glasgow G12 8QQ, United
Kingdom  Olaf Hartwig  SYRTE, Observatoire de Paris, Université PSL, CNRS,
Sorbonne Université, LNE, 61 avenue de l’Observatoire 75014 Paris, France  Max
Planck Institute for Gravitational Physics (Albert Einstein Institute),
D-30167 Hanover, Germany  Leibniz Universität Hannover, D-30167 Hanover,
Germany

(August 28, 2024)

###### Abstract

Time-delay interferometry (TDI) is a crucial step in the on-ground data
processing pipeline of the 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 $\mathrm{MHz}$ down
to a few $\mathrm{Hz}$, requires careful design of the anti-aliasing 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. Secondly, the post-processing 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 residual in the TDI
combination 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.

##  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}\text{\,}\mathrm{Hz}$ to $1\text{\,}\mathrm{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, the laser phase noise is the primary noise source and is 8
orders of magnitude larger than the GW signals that one hopes to detect. Time-
delay interferometry (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 phase noise to levels sufficiently low such that GWs
become detectable [2, 3]. In TDI the measurements are time-delayed by
multiples of the LISA arm lengths and combined in a specific scheme to achieve
laser phase noise reduction. Second-generation TDI, which is the current
baseline laser phase 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 phase noise is strongly
suppressed and the residual ia fundamentally limited by the arm length
mismatch of the virtual interferometer [6].

There exist other approaches to perform laser phase noise suppression. In
TDI-$\infty$ [7], the observables that cancel laser phase 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 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-$\infty$ to
the full LISA constellation with time-evolving arms. Computationally,
TDI-$\infty$ has the drawback that it requires the storage and manipulation of
very large matrices.

In [9], starting from the interferometric measurements and for non-evolving
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 phase noise is significantly suppressed. This approach, dubbed
“automated Principal Component Interferometry”, or aPCI, is formulated in both
the time and frequency domain. 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 the laser phase noise.

While approaches such as TDI-$\infty$ 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
for example 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 ecent 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 frequency 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 in so as to have a more complete description of the post-processing
delay. We consider the effect of on-board processing (i.e. filtering and
decimation), and the influence of TDI which uses post-processing delay
operations that are subject to ranging and interpolation errors, see Figure 1.
We compute analytic formulae for all laser noise residuals induced by these
processing steps and compare those to 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 non-commutativity 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
models for the aliased laser noise due to decimation and interpolation errors
were 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 organised as follows. In section 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 section III we discuss the residual laser noise induced by the non-
commutativity of the on-board filtering and decimation operations with the
propagation delays. The interpolation and ranging residuals, which both arise
from the on-ground processing steps, are discussed in section IV. In section
V, 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 four appendices.

Figure 1: Illustration of the formation of interferometric beatnote phases and
successive on-board processing steps on the left and post-processing TDI on
the right. The propagation and post-processing delay are defined in section II
and section IV, respectively

##  II Interferometric measurements

LISA produces two main interferometric measurements per movable optical sub-
assembly (MOSA) relevant for laser noise reduction111The split interferometry
configuration involves a third interferometer, the test-mass interferometer,
which is not relevant for the purpose of this study. Those are the inter-
spacecraft and reference interferometers given by

$\displaystyle\mathrm{isi}_{ij}(t)$
$\displaystyle={\mathbf{S}}{\mathbf{F}}\big{(}{\mathbf{D}_{ij}}\phi_{ji}(t)-\phi_{ij}(t)\big{)},$
(1a) $\displaystyle\mathrm{rfi}_{ij}(t)$
$\displaystyle={\mathbf{S}}{\mathbf{F}}\big{(}\phi_{ik}(t)-\phi_{ij}(t)\big{)}.$
(1b)

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 inter-spacecraft
interferometer $\mathrm{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

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

Here, $d_{ij}(t)$ is the pseudo-range, 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
spacecraft222Each 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)., the reference interferometer $\mathrm{rfi}_{ij}$ combines the
local and the adjacent lasers.

It is useful to decompose the total laser phase $\phi_{ij}(t)$ or the
frequency $\nu_{ij}(t)=\dot{\phi}_{ij}(t)$ into two variables because, as we
shall see below, the instrumental and data pre-processing effects we will
describe (see the diagram of fig. 1) couple differently to the phase ramp
$\phi_{ij}^{o}(t)$ and any in-band fluctuations $\phi_{ij}^{\epsilon}(t)$. We
thus write

$\displaystyle\phi_{ij}(t)$
$\displaystyle=\phi_{ij}^{o}(t)+\phi_{ij}^{\epsilon}(t),$ (3)
$\displaystyle\nu_{ij}(t)$
$\displaystyle=\nu_{ij}^{o}(t)+\nu_{ij}^{\epsilon}(t),$ (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\text{\,}\mathrm{THz}$$ and
$\nu_{ij}^{\epsilon}(t)=\dot{\phi}_{ij}^{\epsilon}(t)$ accounts for any
rapidly varying random fluctuations. This in-band part is dominated by laser
frequency noise with an amplitude spectral density (ASD)
$\sqrt{S_{\dot{p}}}=$30\text{\,}\mathrm{Hz}\mathrm{/}\sqrt{\mathrm{Hz}}$$.

The beatnote phases of the interferometers are read out using a digital phase
locked loop running at $80\text{\,}\mathrm{MHz}$. Multiple decimation stages
reduce the sampling rate down to $4\text{\,}\mathrm{Hz}$ in order to produce
the final data streams telemetered to ground. Each decimation stage consists
of an anti-aliasing filter, ${\mathbf{F}}$, and a downsampling stage,
${\mathbf{S}}$, which downsamples the data by an integer factor. In this work,
we compare the analytic models that we derive with the most recent LISA
simulation codes, which run at rates that are much lower than the
$80\text{\,}\mathrm{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, 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}\operatorname{sinc}(f_{s}t-n)\cdot x_{n},$ (5)

which reconstructs the continuous time signal $x(t)$ from discrete samples
$x_{n}$.

Let us first describe the onboard 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}.$ (6)

We use eq. 5 to represent the output $y_{n}$ in continuous time and find

$y(t)=\int_{\mathbb{R}}\\!\underbrace{\sum_{m}h_{m}\delta(\tau-
mT_{s})}_{h_{\mathbf{F}}(\tau)}x(t-\tau)\,\mathrm{d}\tau={\mathbf{F}}x(t),$
(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_{\mathbf{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)=\left|\tilde{h}_{\mathbf{F}}(f)\right|^{2}S_{x}(f)=\tilde{\mathbf{F}}S_{x}(f)$
(8)

with $\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}\,.$ (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}\operatorname{sinc}(f_{s}t-n)\cdot
x_{n\cdot M}.$ (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}}\,S_{x}(f)=\operatorname{rect}\left(\frac{f}{f_{s}}\right)\sum_{n=0}^{M-1}S^{(n)}_{x}(f)\,.$
(11)

Here, $\tilde{\mathbf{S}}\,S_{x}(f)$ is a shorthand notation representing the
action of the decimation operator in Fourier space on the PSD of $x(t)$ and
the rectangular function 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^{\mathrm{th}}$ alias,
$S^{(n)}_{x}(f)$, on the right-hand-side of eq. 11, is given by

$S_{x}^{(n)}(f)=\begin{cases*}S_{x}(nf_{n}+f)&if $n$ is even,\\\
S_{x}((n+1)f_{n}-f)&if $n$ is odd.\end{cases*}$ (12)

Equations 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 Commutator Residuals

When computing TDI laser noise residuals, expressions of the form

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

appear, with operators $A$ and $B$ that act on the time series $\phi_{ij}(t)$.
We call two operators non-commutative if eq. 13 is non-vanishing. 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 section V).

Let us now outline how the filtering and sampling operators enter as
additional commutators in TDI.

The basic building block of every TDI combination is the set of intermediary
variables $\eta_{ij}$ (defined later in section V). 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
$\mathcal{\hat{D}}_{ij}=\mathbf{D}_{ij}$) the variable $\eta_{ij}$ simplifies
to the difference between a local laser phase with a distant laser phase
delayed by the light travel time between spacecraft $i$ and $j$. For the left-
handed MOSAs, one has

$\eta_{ij}(t)={\mathbf{D}_{ij}}\phi_{jk}(t)-\phi_{ij}(t).$ (14)

The expression for the right-handed MOSAs is similar but the indexing is
different. Under these conditions, i.e., with the specific algebraic form of
eq. 14, the fundamental laser noise limit is the the usual delay commutator
(see e.g. eq. 52).

If we instead insert the interferometric measurements of eq. 1 into the
definition of $\eta_{ij}$ given in eq. 38a, the resulting expression cannot be
recast into the algebraic form of eq. 14. 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{S}}{\mathbf{F}},{\mathbf{D}_{ij}}]\phi_{ji}(t)=({\mathbf{S}}{\mathbf{F}}){\mathbf{D}_{ij}}\phi_{ji}(t)-{\mathbf{D}_{ij}}({\mathbf{S}}{\mathbf{F}})\phi_{ji}(t),$
(15)

the delay operator ${\mathbf{D}_{ij}}$ switches places with the decimation
stage ${\mathbf{S}}{\mathbf{F}}$ such that eq. 1a becomes

$\mathrm{isi}_{ij}(t)={\mathbf{D}_{ij}}{\mathbf{S}}{\mathbf{F}}\phi_{ji}(t)-{\mathbf{S}}{\mathbf{F}}\phi_{ij}(t)+[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{ij}}]\phi_{ji}(t).$
(16)

This expression has the same algebraic form as eq. 14 if we make the
replacement $\phi_{ij}\rightarrow{\mathbf{S}}{\mathbf{F}}\phi_{ij}$, with the
exception that the inter-spacecraft interferometric measurement now also
contains the commutator
$[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{ij}}]\phi_{ji}$. This commutator comes
as an additional residual in the final TDI expressions (this is reminiscent of
the way readout noise enters the measurements, see [11]).

Using common commutator rules we can further split
$[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{ij}}]\phi_{ji}$ into a filter-delay
commutator and a decimation-delay commutator

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

and compute explicit analytic expressions for each contribution, see section
III.1 and section III.2.

###  III.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)\simeq\dot{d}(t){\mathbf{D}}{\mathbf{G}}\nu(t).$
(18)

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

Because of the orbital dynamics of the LISA constellation, $\dot{d}(t)$ will
vary on time scales 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. 18 reads

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

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.

###  III.2 Decimation-delay commutator

The second commutator appearing in eq. 17 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).$
(20)

Those operations do not commute due to the non-linear 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),$ (21)

where the modulating factor $c_{n}$ is given by

$c_{n}(d)=\begin{cases*}\sin[2](\pi f_{s}d\frac{n}{2})&for $n$ even,\\\
\sin[2](\pi f_{s}d\frac{n+1}{2})&for $n$ odd.\end{cases*}$ (22)

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)\leq
4\cdot\sum_{n=1}^{\infty}S_{\phi}^{(n)}(f).$ (23)

This bound is independent of the delay, and corresponds to the case of full
anti-correlation between ${\mathbf{S}}{\mathbf{D}}\phi(t)$ and
${\mathbf{D}}{\mathbf{S}}\phi(t)$ in eq. 20.

###  III.3 Comparison with numerical simulations

Figure 2: Commutator residuals from filtering and decimation. We compare the
numerically simulated data (yellow) against analytic models (dashed lines). In
purple we show the model for the coupling of the decimation-delay commutator
(aliasing) and in indigo the flexing-filtering effect.

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
precision333A beatnote frequency of several $\mathrm{MHz}$ results in a
rapidly increasing phase ramp. This would require a large number of
significant digits to deal with its dynamic range
($\sim${10}^{14}\text{\,}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{s}$$
after $1\text{\,}\mathrm{yr}$ at a required precision of at least
$\mathrm{\SIUnitSymbolMicro}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{s}$
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
$\mathrm{\SIUnitSymbolMicro}\mathrm{c}\mathrm{y}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{s}$
requirement with a knee at $2\text{\,}\mathrm{mHz}$ corresponds to a precision
of approximately $20\text{\,}\mathrm{nHz}$ which gives 15 orders of magnitude
when considering a $20\text{\,}\mathrm{MHz}$ beatnote frequency.. 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,

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

We can then easily rewrite the commutator given in eq. 17 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}},{{\mathbf{\dot{D}}}}]\nu(t)={\mathbf{S}}{\mathbf{F}}{{\mathbf{\dot{D}}}}\nu(t)-{{\mathbf{\dot{D}}}}{\mathbf{S}}{\mathbf{F}}\nu(t).$
(25)

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

$\displaystyle
S_{\delta\nu}^{{\mathbf{S}}[{\mathbf{F}},{{\mathbf{\dot{D}}}}]}(f)$
$\displaystyle\simeq\tilde{\mathbf{S}}\left((2\pi
f)^{2}\dot{d}^{2}\left|\frac{1}{2\pi}\frac{\mathrm{d}\tilde{h}_{\mathbf{F}}(f)}{\mathrm{d}f}\right|^{2}S_{\nu}(f)\right),$
(26a) $\displaystyle
S_{\delta\nu}^{[{\mathbf{S}},{{\mathbf{\dot{D}}}}]{\mathbf{F}}}(f)$
$\displaystyle\simeq
4\sum_{n=1}^{\infty}c_{n}(d)\left(\tilde{\mathbf{F}}S_{\nu}\right)^{(n)}\\!(f).$
(26b)

We note that for the full commutator in eq. 17, we need to apply
${\mathbf{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. 26, we
compare it to numerical simulations that generate time series corresponding to
eq. 25. We consider white laser frequency noise in $\nu(t)$ with an ASD of
$30\text{\,}\mathrm{Hz}\mathrm{/}\sqrt{\mathrm{Hz}}$444In reality, we expect
laser noise to increase towards 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. and neglect the central laser frequency of
$2.816\text{\times}{10}^{14}\text{\,}\mathrm{Hz}$ (its coupling to the
commutator is vanishing). The simulation is performed at a sampling rate of
$16\text{\,}\mathrm{Hz}$ and is then decimated down to
$4\text{\,}\mathrm{Hz}$. Anti-aliasing is performed using an FIR designed
according to appendix A. The delay operator is modeled by numerical
interpolation of the data using Lagrange polynomials (c.f. section IV).

In fig. 2 we compare the simulated data against the analytical model presented
in eq. 26. We use the Welch method from the SciPy package [26] to estimate the
PSD of the time series. The numerical result is explained by aliasing at low
frequencies and by flexing-filtering at high frequencies.

##  IV Post-processing residuals in TDI

The working principle of TDI is to time-shift the recorded beatnote phase
measurements and linearly combine them in order to reduce laser frequency
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 in on-ground data
processing as ${\hat{\mathcal{D}}}$. This operator acts on a discrete time
series and is therefore only a numerical approximation. Furthermore, time-
shifting is performed with imperfect knowledge of the delay. The discreteness
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{\mathcal{D}}}$ as

$({\hat{\mathcal{D}}}-{\mathbf{D}})\phi(t)=({\hat{\mathcal{D}}}-{\hat{\mathbf{D}}})\phi(t)+({\hat{\mathbf{D}}}-{\mathbf{D}})\phi(t),$
(27)

where the first term on the right-hand side represents the interpolation
residual and the second term the residual stemming from ranging errors. In the
following we study both effects in detail.

###  IV.1 Interpolation residual

We follow [27] and implement a post-processing delay 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 post-processing 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\leq
n\leq 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}.$
(28)

Before the convolution of eq. 28 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_{\mathcal{D}}(\tau)=\sum_{m=-N/2}^{N/2-1}k_{m}(\epsilon)\cdot\delta(\tau-(j+m)T_{s}).$
(29)

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

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

and use eq. 8 to derive the residual in terms of PSD. We find

$S_{\delta\phi}^{\mathcal{D}}(f)=\underbrace{\left|\tilde{h}_{\mathcal{D}}(f)e^{2\pi
ifd}-1\right|^{2}}_{\tilde{\bm{\Delta}}}S_{\phi}(f).$ (31)

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 non-stationary 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$555A rigorous proof is needed to validate this assumption which
was only found to be true empirically, in our work, and specifically for
Lagrange interpolation., we can derive an upper bound for the residual induced
by interpolation.

As suggested in [27], 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}\text{\,}\mathrm{Hz}1\text{\,}\mathrm{Hz}$) when using high
interpolation orders666Assuming a sampling rate of $4\text{\,}\mathrm{Hz}$,
typical filter lengths used for post-processing delays in TDI are in the range
32 to 66.. Alternative interpolation kernels are under study that use less
coefficients and optimize their performance over the entire band (and not only
at DC).

###  IV.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))\simeq{\mathbf{D}}\phi(t)-r(t){\mathbf{D}}\dot{\phi}(t),$
(32)

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)$
(33)

with $\nu(t)=\dot{\phi}(t)$. Similarly to laser phase or frequency, c.f. eq.
4, we decompose the ranging error $r(t)$ into an out-of-band component
$r^{o}(t)$ and an in-band component $r^{\epsilon}(t)$. Here, $r^{o}(t)$
absorbs ranging biases that are of the order of $3\text{\,}\mathrm{ns}$ [13]
and might be slowly drifting. The in-band component $r^{\epsilon}(t)$ has a
root-mean-squared value of $\sim$100\text{\,}\mathrm{fs}$$ (assuming the PSD
in eq. 71). Therefore, $r^{\epsilon}(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 $\mathrm{MHz}$ beatnote
frequency

$r(t){\mathbf{D}}\nu(t)\simeq
r^{o}(t){\mathbf{D}}\nu^{\epsilon}(t)+r^{\epsilon}(t){\mathbf{D}}\nu^{o}(t).$
(34)

For completeness, in appendix D 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. 33 as

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

###  IV.3 Comparison with numerical simulations

Figure 3: Processing residuals in application of the fractional delay filter
containing ranging errors. We compare the numerically simulated data (yellow)
against the analytical models for the ranging residual caused by ranging noise
(green) and a constant bias (wine). Additionally, we plot models for the
interpolation residual (teal) for different interpolation kernel lengths
$N=8,14,62$.

Once again, we check the validity of our analytic model using simulations
performed in units of frequency. Using eq. 24, we write down the expression
corresponding to the left-hand-side of eq. 27 for frequency data as

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

In order to simulate time series data corresponding to eq. 36, we first
generate a generic beatnote frequency with a constant offset of
$10\text{\,}\mathrm{MHz}$ and a white laser frequency noise component with an
ASD of $30\text{\,}\mathrm{Hz}\mathrm{/}\sqrt{\mathrm{Hz}}$ at a sampling
frequency of $4\text{\,}\mathrm{Hz}$. Both the post-processing delay
${\hat{\mathcal{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 post-processing delay. The ranging
error present in the post-processing delay is modeled by a bias
$B=${10}^{-8}\text{\,}\mathrm{s}$$, and ranging noise777We 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 section V. with an ASD of
$${10}^{-15}\text{\,}\mathrm{s}\mathrm{/}\sqrt{\mathrm{Hz}}$\left(\frac{$\mathrm{Hz}$}{f}\right)$
. The nominal value of the delay is taken to be equal to
$d=$8.125\text{\,}\mathrm{s}$$. This yields the worst case interpolation
error, with a fractional part $\epsilon=0.5$.

In fig. 3 we compare the PSD of the numerical time series corresponding to eq.
36 with the analytic expressions for each of the two components of eqs. 31 and
35 re-expressed for frequency data, i.e.,

$\displaystyle S_{\delta\nu}^{\mathcal{D}}(f)$
$\displaystyle=\tilde{\bm{\Delta}}S_{\nu}(f),$ (37a) $\displaystyle
S_{\delta\nu}^{\hat{\mathbf{D}}}(f)$ $\displaystyle=(2\pi f)^{2}\cdot
S_{\delta\phi}^{\hat{\mathbf{D}}}(f).$ (37b)

The interpolation residual (dashed teal) is strongly dependent on the length
$N$ of the interpolation kernel. For this reason, we show the interpolation
residual obtained for $N=8,\,14$, and $62$. The post-processing delay used in
the numerical implementation is performed with an interpolation order $N=14$.
As shown in the figure, the model of eq. 37 agrees with the data for all
frequencies.

##  V Second-generation Michelson combinations

Figure 4: Illustration of intermediary variables defined in eq. 38 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 build up from the previously defined ones such that the last
variable $\sigma$ is represented by the entire path. Additionally, the chain
of locked lasers for the locking configuration N1-12 is shown. The primary
laser is highlighted by the grey box.

In this section we present how the effects stemming from commutators and post-
processing appear in the second-generation TDI combination $X_{2}$. In order
to optimize numerical precision and to save computational cost, we calculate
$X_{2}$ in several stages using the following intermediary variables. To start
with, the variable $\eta$ is constructed from the inter-spacecraft 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
thus have

$\displaystyle\eta_{ij}$
$\displaystyle=\begin{cases*}\mathrm{isi}_{ij}-{\hat{\mathcal{D}}_{ij}}\frac{\mathrm{rfi}_{jk}-\mathrm{rfi}_{ji}}{2}&if
$\epsilon_{ijk}=1$,\\\
\mathrm{isi}_{ij}+\frac{\mathrm{rfi}_{ik}-\mathrm{rfi}_{ij}}{2}&if
$\epsilon_{ijk}=-1$,\end{cases*}$ (38a) $\displaystyle\pi_{ij}$
$\displaystyle=\eta_{ij}+{\hat{\mathcal{D}}_{ij}}\eta_{ji},$ (38b)
$\displaystyle\rho_{ij}$
$\displaystyle=\pi_{ij}+{\hat{\mathcal{D}}_{iji}}\pi_{ik},$ (38c)
$\displaystyle\sigma_{ij}$
$\displaystyle=\rho_{ij}+{\hat{\mathcal{D}}_{ijiki}}\rho_{ik}.$ (38d)

Figure 4 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 $\sigma$ we can finally express the second-generation Michelson
combination $X_{2}$ as

$X_{2}=\sigma_{13}-\sigma_{12}.$ (39)

The contracted post-processing delays ${\mathcal{D}_{i_{1}\cdots i_{N}}}$
appearing in eq. 38 are applied in two steps. First, the nested delay
$d_{i_{1}\cdots i_{N}}$ is calculated, and then, the fractional delay filter
is applied. The nested delay can be calculated using the recursive
operation888The 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.

$d_{i_{1}\cdots i_{N}}=d_{i_{1}i_{2}}+{\mathcal{D}_{i_{1}i_{2}}}d_{i_{2}\cdots
i_{N}}.$ (40)

Alternatively, contracted post-processing delays can be decomposed into atomic
delay operations

${\mathcal{D}_{i_{1}\cdots
i_{N}}}\rightarrow{\mathcal{D}_{i_{1}i_{2}}}\cdots{\mathcal{D}_{i_{N-1}i_{N}}}.$
(41)

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 below. Whether a single contracted delay operator or its decomposition in
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 over delaying a time series with a large one, since the
interpolation residual vanishes for an integer delay.

To obtain an explicit expression for the residual in $X_{2}$ we use the
factorization given in eqs. 38 and 39 together with the expressions for the
inter-spacecraft and reference interferometers in eqs. 16 and 1b. As the
couplings discussed in sections III and IV 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}=\dot{\mathrm{isi}}_{ij}^{o}$ and
$b_{i}=\dot{\mathrm{rfi}}_{ij}^{o}=-\dot{\mathrm{rfi}}_{ji}^{o}$999This
relation holds for $ij=12,23,31$ as adjacent reference interferometers track
the same beatnote up to a sign. as in [13]. The beatnote phase fluctuations
are given by

$\displaystyle\mathrm{isi}_{ij}^{\epsilon}$
$\displaystyle={\mathbf{D}_{ij}}{\mathbf{S}}{\mathbf{F}}p_{ji}-{\mathbf{S}}{\mathbf{F}}p_{ij}+[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{ij}}]p_{ji},$
(42a) $\displaystyle\mathrm{rfi}_{ij}^{\epsilon}$
$\displaystyle={\mathbf{S}}{\mathbf{F}}p_{ik}-{\mathbf{S}}{\mathbf{F}}p_{ij}.$
(42b)

Here, $p_{ij}$ denotes the phase noise of laser $ij$. In the following we
write as a short-hand notation $\bar{p}_{ij}={\mathbf{S}}{\mathbf{F}}p_{ij}$.

Next, we replace each occurrence of the post-processing delay
${\hat{\mathcal{D}}_{ij}}$ by the relation given in eq. 27. Here, we neglect
all second-order terms in the residuals. Doing so, we find that $X_{2}$ breaks
down into the following residuals

$X_{2}=\delta X_{2}^{[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}}]}+\delta
X_{2}^{\mathcal{D}}+\delta X_{2}^{{\hat{\mathbf{D}}}}+\delta X_{2},$ (43)

where the origin of the first three terms on the right-hand side is described
in sections III and IV. Additionally, the residuals are propagated through TDI
and therefore appear multiple times with different delays. The fourth term in
eq. 43 is the usual delay commutator describing the arm length mismatch in the
virtual interferometer.

In the following we discuss the constituents of eq. 43 individually. To
simplify the expressions we assume equal arms and denote the transfer that is
common to all residuals as

$\displaystyle{\mathbf{C}}$
$\displaystyle\equiv(1-{\mathbf{D}}^{4})(1-{\mathbf{D}}^{2}),$ (44)
$\displaystyle\tilde{\mathbf{C}}$ $\displaystyle\equiv 16\sin^{2}(2\pi
fd)\sin^{2}(4\pi fd).$ (45)

The residuals induced by filtering and decimation enter in the inter-
spacecraft interferometer (cf., eq. 16) and are thus only propagated through
TDI. This is consistent with [15], and they read

$\begin{split}\delta
X_{2}^{[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}}]}&={\mathbf{C}}\\{[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{13}}]p_{31}+{\mathbf{D}}[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{31}}]p_{13}\\\
&\quad\;\;-[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{12}}]p_{21}-{\mathbf{D}}[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{21}}]p_{12}\\}.\end{split}$
(46)

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

$\begin{split}\delta
X_{2}^{{\mathcal{D}}}&={\mathbf{C}}{\mathbf{D}}[\bm{\Delta}_{12}\bar{p}_{21}-{\mathbf{D}}\bm{\Delta}_{31}\bar{p}_{13}-\bm{\Delta}_{13}\bar{p}_{31}\\\
&\quad+{\mathbf{D}}(\bm{\Delta}_{13}\\!+\\!\bm{\Delta}_{31}\\!-\\!\bm{\Delta}_{12}\\!+\\!\bm{\Delta}_{121}\\!-\\!\bm{\Delta}_{131})\bar{p}_{12}].\end{split}$
(47)

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

$\begin{split}\delta
X_{2}^{{\mathcal{D}}}&={\mathbf{C}}{\mathbf{D}}\\{\bm{\Delta}_{12}\bar{p}_{21}+{\mathbf{D}}\bm{\Delta}_{21}\bar{p}_{12}\\\
&\qquad\qquad-{\mathbf{D}}\bm{\Delta}_{31}\bar{p}_{13}-\bm{\Delta}_{13}\bar{p}_{31}\\}.\end{split}$
(48)

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. In LISA, two additional modulations on the laser beam are
used to measure the inter-satellite range: the pseudo-random noise (PRN)
modulation (absolute ranging) and the clock sidebands. The ranging processing
described in appendix B, which is mostly adopted from [13], combines both
measurements. The resulting ranging estimates $\hat{d}_{ij}$ inherit a bias
$B_{ij}$ from the PRN measurements and a stochastic term from the sideband
measurements. Therefore, the ranging error reads

$r_{ij}(t)=B_{ij}-{\mathbf{D}_{ij}}M_{j}(t)+M_{i}(t),$ (49)

where $M_{i}$ denotes the modulation noise on left-handed MOSAs. The resulting
residual $\delta X_{2}^{{\hat{\mathbf{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

$\displaystyle\begin{split}\delta
X_{2}^{B}&={\mathbf{C}}\\{{\mathbf{D}}^{2}B_{31}\dot{\bar{p}}_{13}+{\mathbf{D}}B_{13}\dot{\bar{p}}_{31}\\\
&\qquad-{\mathbf{D}}^{2}B_{21}\dot{\bar{p}}_{12}-{\mathbf{D}}B_{12}\dot{\bar{p}}_{21}\\},\end{split}$
(50a) $\displaystyle\begin{split}\delta
X_{2}^{M}&={\mathbf{C}}\\{(a_{13}-a_{12}+(1-{\mathbf{D}}^{2})b_{1})M_{1}\\\
&\qquad-{\mathbf{D}}a_{21}M_{2}+{\mathbf{D}}a_{31}M_{3}\\},\end{split}$ (50b)

where $\dot{\bar{p}}_{ij}$ is the time derivative of the filtered and
decimated laser phase fluctuations.

Let us finally give the expression for 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 proportional to the arm length mismatch
$\Delta d$. For the second-generation Michelson variable $X_{2}$, the residual
reads

$\displaystyle\delta X_{2}$
$\displaystyle=[{\mathbf{D}_{13121}},{\mathbf{D}_{12131}}]{\mathbf{S}}{\mathbf{F}}\phi_{12}$
(51) $\displaystyle\simeq-\Delta
d_{X_{2}}{\mathbf{D}_{131212131}}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{S}}{\mathbf{F}}\phi_{12}.$
(52)

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

$\Delta d_{X_{2}}(t)=d_{131212131}(t)-d_{121313121}(t).$ (53)

Equation 52 can be further split up into a deterministic out-of-band drift and
an in-band component by plugging in eq. 4. As the deterministic part, we
recover101010As the phase ramp produced by the $\mathrm{THz}$ central laser
frequency is unaffected by filtering and decimation we can drop both
operations.

$\displaystyle\delta X_{2}^{o}=-\Delta
d_{X_{2}}{\mathbf{D}_{131212131}}\nu_{12}^{o},$ (54)

where the travel time difference can be efficiently computed from the time
delays $d_{ij}$ as described in appendix C. 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)=\left(\Delta d_{X_{2}}\right)^{2}(2\pi
f)^{2}\,\tilde{\mathbf{S}}\tilde{\mathbf{F}}S_{p}(f).$ (55)

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],

$\begin{split}\Delta
d_{X_{2}}&=(d_{131}\dot{d}_{121}-d_{121}\dot{d}_{131})(\dot{d}_{121}+\dot{d}_{131})\\\
&\qquad-(d_{131}\ddot{d}_{121}-d_{121}\ddot{d}_{131})(d_{121}+d_{131}),\end{split}$
(56)

which is of the order ${10}^{-12}\text{\,}\mathrm{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 echos 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 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-12111111Note that in this work, we shall use the
locking configuration naming convention introduced in [24]. introduced in
[23].

###  V.1 Six independent lasers

For the case of six lasers stabilized to their individual cavities, we assume
that their in-band phase noises $p_{ij}$ are independent. Thus, we can compute
the total PSD as the sum of the PSDs of each laser’s contribution.

Using eq. 19 and the upper bound given in eq. 23, it is easy to compute the
PSD of eq. 46 corresponding to the coupling of the filter-delay and
decimation-delay commutator. It reads

$\displaystyle S_{\delta X_{2}}^{{\mathbf{S}}[{\mathbf{F}},{\mathbf{D}}]}(f)$
$\displaystyle=4\tilde{\mathbf{C}}\tilde{\mathbf{S}}\left(\bar{\dot{d}}^{2}\left|\frac{1}{2\pi}\frac{\mathrm{d}\tilde{h}_{\mathbf{F}}(f)}{\mathrm{d}f}\right|^{2}S_{\dot{p}}(f)\right),$
(57a) $\displaystyle S_{\delta
X_{2}}^{[{\mathbf{S}},{\mathbf{D}}]{\mathbf{F}}}(f)$ $\displaystyle\leq
16\tilde{\mathbf{C}}\cdot\sum_{n=1}^{\infty}(\tilde{\mathbf{F}}S_{p})^{(n)}(f),$
(57b)

with the effective squared delay derivative

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

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 the
interpolation error. We use eqs. 31 and 48 and find

$S_{\delta X_{2}}^{{\mathcal{D}}}(f)\leq
4\tilde{\mathbf{C}}\tilde{\bm{\Delta}}\tilde{\mathbf{S}}\tilde{\mathbf{F}}S_{p}(f).$
(59)

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

Finally, using eqs. 35 and 50, the contribution from the ranging error yields
a PSD equal to [13]

$\begin{split}S_{\delta
X_{2}}^{{\hat{\mathbf{D}}}}(f)&=\tilde{\mathbf{C}}\Big{(}4\bar{B}^{2}(2\pi
f)^{2}\tilde{\mathbf{S}}\tilde{\mathbf{F}}S_{p}(f)\\\
&\qquad\qquad\qquad+A_{X_{2}}^{M}(f)S_{M}(f)\Big{)},\end{split}$ (60)

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

$\displaystyle\bar{B}^{2}$
$\displaystyle=\frac{B_{12}^{2}+B_{21}^{2}+B_{13}^{2}+B_{31}^{2}}{4},$ (61)
$\displaystyle\begin{split}A_{X_{2}}^{M}(f)&=(a_{12}-a_{13})^{2}+a_{21}^{2}+a_{31}^{2}\\\
&\qquad-4b_{1}(a_{12}-a_{13}-b_{1})\sin^{2}(2\pi fd).\end{split}$ (62)

###  V.2 Locked lasers

The baseline design of LISA foresees locked lasers to ensure that all beatnote
frequencies fall into the sensitive bandwidth of the photoreceivers onboard
the spacecraft ($5\text{\,}\mathrm{MHz}25\text{\,}\mathrm{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 offset
frequency $o_{ij}$. The so-called locking conditions for the inter-spacecraft
and reference locking interferometer are given by

$\displaystyle\dot{\mathrm{isi}}_{ij}$
$\displaystyle={{\mathbf{\dot{D}}}_{ij}}\nu_{ji}-\nu_{ij}=o_{ij}$ (63a)
$\displaystyle\dot{\mathrm{rfi}}_{ij}$
$\displaystyle=\nu_{ik}-\nu_{ij}=o_{ij}$ (63b)

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 the offset frequencies $o_{ij}$ 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. 4. 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{split}\begin{aligned}
p_{12}&=p,&p_{23}&={\mathbf{D}_{21}}p,&p_{31}&={\mathbf{D}_{31}}p,\\\
p_{13}&=p,&p_{32}&={\mathbf{D}_{31}}p,&p_{21}&={\mathbf{D}_{21}}p.\end{aligned}\end{split}$
(64)

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. 64 into the general expressions for laser noise related
residuals listed in section V. 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}}].$
(65)

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 indices as it is broken by laser
locking; this is easily verified from eq. 64.

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}$121212For example,
the inter-spacecraft measurements $\mathrm{isi}_{12}\sim\pi_{12}$ and
$\mathrm{isi}_{13}\sim\pi_{13}$.. 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 $[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}}]$ commutator we recover

$\delta
X_{2}^{[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}}]}={\mathbf{C}}([{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{121}}]-[{\mathbf{S}}{\mathbf{F}},{\mathbf{D}_{131}}])p.$
(66)

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

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

For the sample-delay commutator, the upper bound given in eq. 57b is still
valid. Furthermore, we note that both residuals vanish non-trivially 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 to laser noise for locked
lasers we use contracted delays (see eq. 47). We obtain

$\delta
X_{2}^{{\mathcal{D}}}={\mathbf{C}}{\mathbf{D}}^{2}(\bm{\Delta}_{131}-\bm{\Delta}_{121}){\mathbf{S}}{\mathbf{F}}p.$
(68)

Again, the upper bound given in eq. 59 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. 50b stays unchanged as the laser noise
component of the beatnote is not involved. The contribution of the ranging
biases in each arm, see eq. 50a, is given by

$\delta
X_{2}^{B}={\mathbf{C}}{\mathbf{D}}^{2}(B_{12}+B_{21}-B_{13}-B_{31}){\mathbf{S}}{\mathbf{F}}p.$
(69)

Therefore, the PSD of ranging error contributions is expressed as eq. 60 with
an effective squared bias equal to

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

###  V.3 Comparison with numerical simulations

Figure 5: ASDs of the second-generation Michelson variable $\dot{X}_{2}$ for
six independent lasers. The dashed lines correspond to the theoretical
predictions (upper bound for aliasing) for the different TDI residuals. Solid
lines are numerical estimates resulting from the three simulations described
in the text. The black vertical at ${10}^{-1}\text{\,}\mathrm{Hz}$ divides the
x-axis into logarithmically and linearly scaled. Note that the effects of
interpolation appears for frequencies larger than $1\text{\,}\mathrm{Hz}$ and
are shown in fig. 7. Figure 6: ASDs of the second-generation Michelson
variable $\dot{X}_{2}$ for locked lasers (N1-12). The dashed lines correspond
to the theoretical predictions (upper bound for aliasing) for the different
TDI residuals. The full lines are numerical estimates resulting from the three
simulations described in the text. The black vertical at
${10}^{-1}\text{\,}\mathrm{Hz}$ divides the x-axis into logarithmically and
linearly scaled. Note that the effects of interpolation appears for
frequencies larger than $1\text{\,}\mathrm{Hz}$ and are shown in fig. 7.
Figure 7: ASDs of residual noise in second-generation Michelson variable
$\dot{X}_{2}$ due to interpolation errors for six independent laser noises
(left) and locked lasers (right). The numerical ASDs using nested (cyan) and
flat delay operators (yellow) are compared against (exact) models derived from
eq. 47 and eq. 48, respectively.

We now compare the theoretical model described above to simulated LISA data,
obtained using LISA Instrument [29, 24] to generate the LISA measurements and
PyTDI [30] to calculate the second-generation Michelson variables. To compare
the couplings described above to the simulated data, individually, we run
three simulations with increasing complexity. In the first simulation, all
sources of noise are disabled. Any residual in the TDI variables can be
attributed to numerical effects that are due to rounding errors in the
simulated floating-point variables. In the second simulation, we add white
laser frequency noise with an ASD equal to
$30\text{\,}\mathrm{Hz}\mathrm{/}\sqrt{\mathrm{Hz}}$. This gives rise to the
residuals caused by the flexing-filtering effect, the coupling of the
sampling-delay commutator and interpolation errors. Finally, we introduce
modulation noise and ranging biases of the order of $30\text{\,}\mathrm{ns}$
(a factor 10 higher than in [13] in order to accentuate the coupling) 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\text{\times}{10}^{-15}\text{\,}\mathrm{s}\mathrm{/}\sqrt{Hz}$\left(\frac{f}{$\mathrm{Hz}$}\right)^{-2/3},$
(71)

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

All simulations span approximately 3 days in duration with the LISA
constellation following realistic heliocentric orbits provided by ESA. Physics
is simulated at $16\text{\,}\mathrm{Hz}$ and the filter presented in appendix
A is used for anti-aliasing before decimating to $4\text{\,}\mathrm{Hz}$.
Here, the filter design is not optimized using a trade-off between the number
of taps and the level of laser noise residual but rather chosen such that the
processing residuals are above numerical noises and can be validated against
analytical models.

LISA Instrument produces measurements in units of frequency in order to
circumvent numerical issues specific to phase units. As a result, the data
produced contains the total beatnote frequencies of the inter-spacecraft and
the reference interferometers, which have offsets between
$5\text{\,}\mathrm{MHz}25\text{\,}\mathrm{MHz}$. In addition to the carrier-
to-carrier beatnotes, sideband beatnotes that precisely track the delay
derivatives are also available.

Prior to any processing all measurements are converted 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 B). Then,
the second-generation Michelson variable $X_{2}$ is calculated using the
factorization expression in eq. 38. Here, we use Lagrange interpolation with
$n=62$ coefficients. To convert our expressions to frequency units, all delay
operators ${\hat{\mathcal{D}}_{ij}}$ are replaced by
${\dot{\hat{\mathcal{D}}}_{ij}}$, which are defined analogously to eq. 24. 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
C, and inserting it into the time derivative of eq. 54.

The resulting ASDs of the simulations for six independent and locked lasers
are presented in figs. 5, 6 and 7. To match the analytical models at high
frequencies (i.e. flexing-filtering effect and interpolation error) with the
numerical ASDs for locked lasers we had to relax the equal-arms assumption and
use the expressions for six unequal but constant arms. Figures 5 and 6 show
that for the laser-noise-only simulations, the residuals are well explained by
the commutator residuals. For both six independent and locked lasers, we plot
the same upper bound for aliased laser noise from eq. 57b. We note that the
locked case exhibits a slightly reduced noise level, which can be explained by
partial cancellation of the two contributions appearing in eq. 66. Adding
modulation noise and ranging bias to the simulation completely dominates the
aforementioned effects. Indeed, the residual in $X_{2}$ is now dominated by
the ranging and interpolation residuals, which are again well explained by
their analytical counterparts. At very low frequencies
(${10}^{-4}\text{\,}\mathrm{Hz}{10}^{-3}\text{\,}\mathrm{Hz}$) we observe
deviations from the models. Those can explained by the numerical noise floor
arising from rounding-errors in the simulation. We also observe this numerical
limit in the noiseless simulation shown in grey.

In fig. 7, we show the PSDs of interpolation errors for six independent and
locked lasers using either contracted or atomic delays, respectively. As a
side note, this residual only affects frequencies outside the LISA band
($>$1\text{\,}\mathrm{Hz}$$) and thus should be irrelevant from the point of
view of the performance of TDI. This being said, the focus of this paper is to
correctly model the residuals, which is why we include this computation. In
the left hand plot, we show the exact models (assuming six unequal but
constant arms) which can be seen to match the numerical results. In the right
hand plot, we show a similar plot but for locked lasers. On the one hand, we
find that, in this particular case, it is indeed advantageous to use
contracted delays to calculate $X_{2}$ for locked lasers in the N1-12 locking
configuration as it produces a residual that is smaller by approximately a
factor of two compared to using atomic delays. On the other hand, for six
independent lasers we virtually see no difference. As discussed earlier in
this section, there is no general rule for the optimal factorization of delays
as it is dependent on their particular value. However, in the worst case
scenario ($\epsilon=0.5$), atomic delays outperform contracted delays for six
lasers and vice versa for locked lasers (see eqs. 48 and 68).

##  VI 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, onboard processing steps, namely filtering and decimation, give rise to
additive noise in the inter-spacecraft interferometers due to non-commutation
with the delay operation. Secondly, the post-processing 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 and present all laser noise couplings in a consistent framework.
Furthermore, we investigate the impact of laser locking by discussing the
example of the locking configuration N1-12. 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
$\mathrm{mHz}$. This trend depends solely on the out-of-band delays due to
orbital dynamics and the $\mathrm{THz}$ frequency of the involved laser. It
can therefore be computed and removed by appropriately modeling the orbits and
the $\mathrm{THz}$ frequency evolution. This detrending step is a reversible
(the trend can always be added in again) part of pre-processing. 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. To be consistent with the existing literature, we first derive
analytical models for six independent laser. Then, we repeat the calculation
for locked lasers, more specifically the configuration N1-12. In this
configuration, the spacecraft 2 and 3 act as transponders directly sending
light back to spacecraft 1. This means the inter-spacecraft interferometer
beatnotes recorded on spacecraft 1 already represent the signal combinations
$\eta_{12}+{\mathbf{D}_{12}}\eta_{21}$ and
$\eta_{13}+{\mathbf{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, analytic models become more
complicated for locked lasers due to the introduced correlations. Furthermore,
we find that generally the worst case scenario in terms of laser noise
residuals can be larger for locked lasers than for six independent laser.
Hence, future worst case studies should account for all possible locking
configurations.

The level of additional noise due to onboard processing (filtering and
decimation) is strongly dependent on the design of the anti-aliasing filter.
In this study, we used a FIR filter with a transition band ranging from
$0.1\text{\,}\mathrm{Hz}2\text{\,}\mathrm{Hz}$ so as to relax requirements on
filter implementation (the number of taps is reduced to 103 compared to 145 in
the standard LISA Instrument implementation) and accentuate the flexing-
filtering coupling in the LISA band. In theory, the flexing-filtering coupling
can be mitigated by flattening the response of the filter in the pass-band on
ground. 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 post-processing
and has to be taken care of before decimation.

Laser noise residuals stemming from post-processing delays depend on the
interpolation method and the ranging performance. In this paper, we
extensively rely on Lagrange interpolation, which has a maximal flat response
at DC. Alternative interpolation kernels are given by the family of “windowed
sinc” kernels [27] or numerically optimized kernels. 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 also 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.

The results presented in this paper should be invariant under 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}\text{\,}$ 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, appropriate
compensation filter are crucial since sideband beatnotes are also subject to
anti-alias filtering. Any departure from unity in the filter’s 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. 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 50OQ1801 and FKZ 50OQ2301). 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.L.’s
and O.H.’s work was supported by the Programme National GRAM of CNRS/INSU with
INP and IN2P3 co-funded by CNES. O.H. and M.L. gratefully acknowledge support
from the Centre National d’Études Spatiales (CNES).

##  Appendix A Anti-aliasing filter design

The main objective of the anti-aliasing filter is to prevent folding of power
during decimation. In the numerical simulations performed in this work, we
simulate the physics at $16\text{\,}\mathrm{Hz}$ and decimate to
$4\text{\,}\mathrm{Hz}$. Therefore, we need to design the anti-aliasing filter
such that it cuts off at a Nyquist frequency of $2\text{\,}\mathrm{Hz}$. We
use the routines from the SciPy Python package to design a FIR filter with the
so called “window method”. Here, we assume a pass-band frequency of
$0.1\text{\,}\mathrm{Hz}$, a stop-band frequency of $2\text{\,}\mathrm{Hz}$
and a minimum attenuation of $181\text{\,}\mathrm{dB}$. The resulting transfer
function and its derivative which is of relevance for the flexing-filtering
coupling (see section III.1) are plotted in fig. 8.

Figure 8: Magnitude of filter transfer function (blue) and its frequency
derivative (red). The grey dashed line indicates the cutoff frequency at
$2\text{\,}\mathrm{Hz}$.

##  Appendix B 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 inter-spacecraft and reference
interferometer read

$\displaystyle\mathrm{isi}_{ij}^{\mathrm{sb}}$
$\displaystyle={\mathbf{D}_{ij}}\phi_{ji}^{\mathrm{m}}-\phi_{ij}^{\mathrm{m}}-(\nu_{ji}^{\mathrm{m}}-\nu_{ij}^{\mathrm{m}})t,$
(72a) $\displaystyle\mathrm{rfi}_{ij}^{\mathrm{sb}}$
$\displaystyle=\phi_{ik}^{\mathrm{m}}-\phi_{ij}^{\mathrm{m}}-(\nu_{ik}^{\mathrm{m}}-\nu_{ij}^{\mathrm{m}})t,$
(72b)

where the phase of the modulation is given by

$\phi_{ij}^{\mathrm{m}}(t)=\nu_{ij}^{\mathrm{m}}\cdot(t+M_{ij}(t)).$ (73)

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

$\displaystyle\mathrm{isi}_{ij}^{\mathrm{sb}}$
$\displaystyle=-\nu_{ji}^{\mathrm{m}}d_{ij}+{\mathbf{D}_{ij}}\nu_{ji}^{\mathrm{m}}M_{ji}-\nu_{ij}^{\mathrm{m}}M_{ij},$
(74a) $\displaystyle\mathrm{rfi}_{ij}^{\mathrm{sb}}$
$\displaystyle=\nu_{ik}^{\mathrm{m}}M_{ik}-\nu_{ij}^{\mathrm{m}}M_{ij},$ (74b)

which has a similar algebraic structure to the carrier phases (c.f. eq. 1)
where the product $\nu_{ij}^{\mathrm{m}}M_{ij}$ takes the place of the laser
phase $\phi_{ij}$. We now use the definition of the intermediary variable
$\eta_{ij}$ for the sidebands131313The post-processing delays required to
calculate $\eta_{ij}^{\mathrm{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. (see eq. 38a) to cancel modulation noise contributions
stemming from right-handed MOSAs,

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

Here, we use the shorthand notation $M_{1}=M_{12}$, $M_{2}=M_{23}$ and
$M_{3}=M_{31}$.

Finally, the variables $\eta_{ij}^{\mathrm{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}^{\mathrm{sb}}}{\nu_{ji}^{\mathrm{m}}}=d_{ij}-{\mathbf{D}_{ij}}\frac{\nu_{jk}^{\mathrm{m}}}{\nu_{ji}^{\mathrm{m}}}M_{j}+\frac{\nu_{ij}^{\mathrm{m}}}{\nu_{ji}^{\mathrm{m}}}M_{i}.$
(76)

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. 49.

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. 76
yields an accurate measurement of the delay derivative $\dot{\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].

##  Appendix C Differential Doppler shift

As explained in section V, 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. 54). 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,$ (77)

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}$
(c.f. eqs. 1 and 38). 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,$ (78)

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 $\dot{d}_{ij}$ directly and form
$\Delta\dot{d}=\dot{X}$ following the procedure explained in [6].

##  Appendix D Coupling of ranging noise to laser noise

In section IV, 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 $\mathrm{MHz}$ beatnote frequency. However, for the sake of
completeness and as it becomes relevant in processing pipelines where 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,

$\displaystyle\mathrm{isi}_{ij}(t)$
$\displaystyle={\mathbf{D}_{ij}}p_{ji}(t)-p_{ij}(t),$ (79a)
$\displaystyle\mathrm{rfi}_{ij}(t)$ $\displaystyle=p_{ik}(t)-p_{ij}(t).$ (79b)

For simplicity we omit the anti-aliasing filtering and decimation that was
considered in eq. 1.

Next, we insert eq. 79 into eq. 38a where the post-processing delay
${\hat{\mathbf{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}.$
(80)

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. 80 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. 38 and 39 to find the total laser noise residual in the
second-generation Michelson combination $X_{2}$. It consists of the residual
in eq. 80 that is propagated through TDI as well as the commutator of post-
processing delay operators

$\begin{split}[{\hat{\mathbf{D}}_{13121}},{\hat{\mathbf{D}}_{12131}}]&=[{\mathbf{D}_{13121}},{\mathbf{D}_{12131}}]\\\
&\quad+({\hat{\mathbf{D}}_{131212131}}-{\mathbf{D}_{131212131}})\\\
&\quad-({\hat{\mathbf{D}}_{121313121}}-{\mathbf{D}_{121313121}})\end{split}$
(81)

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

$\begin{split}\delta
X_{2}^{\hat{\mathbf{D}}}&={\mathbf{C}}\\{{\mathbf{D}}(r_{31}{\mathbf{D}}\dot{p}_{13})+(r_{13}{\mathbf{D}}\dot{p}_{31})\\\
&\qquad-{\mathbf{D}}(r_{21}{\mathbf{D}}\dot{p}_{12})-(r_{12}{\mathbf{D}}\dot{p}_{21})\\}\\\
&\quad+{\mathbf{C}}\\{r_{12}+{\mathbf{D}}r_{21}-r_{13}-{\mathbf{D}}r_{31}\\}\cdot{\mathbf{D}}^{8}\dot{p}_{1},\end{split}$
(82)

where we have used eq. 33 and have approximated nested ranging noise as

$\displaystyle r_{ijikikiji}$ $\displaystyle\simeq
r_{iji}+{\mathbf{D}}^{2}r_{iki}+{\mathbf{D}}^{4}r_{iki}+{\mathbf{D}}^{6}r_{iji},$
(83a) $\displaystyle r_{iji}$ $\displaystyle\simeq r_{ij}+{\mathbf{D}}r_{ji}.$
(83b)

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

$\begin{split}S_{\delta
X_{2}}^{{\hat{\mathbf{D}}}}(f)&=4\tilde{\mathbf{C}}(S_{r}(f^{\prime})*S_{\dot{p}}(f^{\prime}))(f)\\\
&\qquad+4(\tilde{\mathbf{C}}S_{r}(f^{\prime})*S_{\dot{p}}(f^{\prime}))(f),\end{split}$
(84)

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

## References

  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ , Laser interferometer space antenna (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Tinto and Armstrong [1999] 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).
  * Armstrong  _et al._ [1999] J. W. Armstrong, F. B. Estabrook, and M. Tinto, Time-delay interferometry for space-based gravitational wave searches, The Astrophysical Journal 527, 814 (1999).
  * Shaddock  _et al._ [2003] D. A. Shaddock, M. Tinto, F. B. Estabrook, and J. W. Armstrong, Data combinations accounting for LISA spacecraft motion, Phys. Rev. D 68, 061303 (2003), arXiv:gr-qc/0307080 .
  * Tinto  _et al._ [2004] M. Tinto, F. B. Estabrook, and J. W. Armstrong, Time delay interferometry with moving spacecraft arrays, Phys. Rev. D 69, 082001 (2004), arXiv:gr-qc/0310017 .
  * Bayle  _et al._ [2021] J.-B. Bayle, O. Hartwig, and M. Staab, Adapting time-delay interferometry for LISA data in frequency, Phys. Rev. D 104, 023006 (2021), arXiv:2103.06976 [gr-qc] .
  * Vallisneri  _et al._ [2021] M. Vallisneri, J.-B. Bayle, S. Babak, and A. Petiteau, Time-delay interferometry without delays, Phys. Rev. D 103, 082001 (2021), arXiv:2008.12343 [gr-qc] .
  * Houba  _et al._ [2023] N. Houba, S. Delchambre, G. Hechenblaikner, T. Ziegler, and W. Fichter, Time-delay interferometry infinity for tilt-to-length noise estimation in LISA, Class. Quant. Grav. 40, 107001 (2023).
  * Baghi  _et al._ [2021] Q. Baghi, J. Baker, J. Slutsky, and J. I. Thorpe, Model-independent time-delay interferometry based on principal component analysis, Physical Review D 104, 122001 (2021).
  * Baghi  _et al._ [2022] Q. Baghi, J. G. Baker, J. Slutsky, and J. I. Thorpe, Fully data-driven time-delay interferometry with time-varying delays, Annalen Phys. 2023, 2200447 (2022), arXiv:2209.10851 [gr-qc] .
  * Nam  _et al._ [2022] D. Q. Nam, Y. Lemiere, A. Petiteau, J.-B. Bayle, O. Hartwig, J. Martino, and M. Staab, TDI noises transfer functions for LISA (2022), arXiv:2211.02539 [gr-qc] .
  * Hartwig and Bayle [2021] O. Hartwig and J.-B. Bayle, Clock-jitter reduction in LISA time-delay interferometry combinations, Phys. Rev. D 103, 123027 (2021), arXiv:2005.02430 [astro-ph.IM] .
  * Hartwig  _et al._ [2022] O. Hartwig, J.-B. Bayle, M. Staab, A. Hees, M. Lilley, and P. Wolf, Time-delay interferometry without clock synchronization, Physical Review D 105, 122008 (2022).
  * Paczkowski  _et al._ [2022] 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).
  * Bayle  _et al._ [2019] 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), arXiv:1811.01575 [astro-ph.IM] .
  * Hartwig [2021] O. Hartwig,  _Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna_ , Ph.D. thesis, Leibniz U., Hannover (2021).
  * de Vine  _et al._ [2010] 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).
  * Mitryk  _et al._ [2010] S. J. Mitryk, V. Wand, and G. Mueller, Verification of time-delay interferometry techniques using the university of florida lisa interferometry simulator, Classical and Quantum Gravity 27, 084012 (2010).
  * Mitryk  _et al._ [2012] 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).
  * Yamamoto  _et al._ [2022] 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).
  * Gruning  _et al._ [2015] 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, Exper. Astron. 39, 281 (2015), arXiv:1309.1059 [astro-ph.IM] .
  * Vidal [2023] L. Vidal,  _Validation expérimentale des performances interférométriques de LISA_ , Ph.D. thesis, Université Paris Cité, Paris (2023).
  * Heinzel [2018] G. Heinzel,  _LISA Frequency planning_ , Tech. Rep. (Max-Planck Institute for Gravitational Physics, Albert Einstein Institute, 2018).
  * Bayle and Hartwig [2023] J.-B. Bayle and O. Hartwig, Unified model for the LISA measurements and instrument simulations, Phys. Rev. D 107, 083019 (2023), arXiv:2212.05351 [gr-qc] .
  * Jenkins and Watts [1968] G. M. Jenkins and D. G. Watts,  _Spectral analysis and its applications_ , Holden-Day series in time series analysis (Holden-Day, San Francisco, 1968).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python, Nature Meth. 17, 261 (2020), arXiv:1907.10121 [cs.MS] .
  * Shaddock  _et al._ [2004] D. A. Shaddock, B. Ware, R. E. Spero, and M. Vallisneri, Post-processed time-delay interferometry for LISA, Phys. Rev. D 70, 081101 (2004), arXiv:gr-qc/0406106 .
  * Cornish and Hellings [2003] N. J. Cornish and R. W. Hellings, The Effects of orbital motion on LISA time delay interferometry, Class. Quant. Grav. 20, 4851 (2003), arXiv:gr-qc/0306096 .
  * Bayle  _et al._ [2022] J.-B. Bayle, O. Hartwig, and M. Staab, Lisa instrument (2022).
  * Staab  _et al._ [2022] M. Staab, J.-B. Bayle, and O. Hartwig, Pytdi (2022).
  * Barke [2015] S. Barke,  _Inter-spacecraft frequency distribution for future gravitational wave observatories_ , Ph.D. thesis, Leibniz U., Hannover (2015).