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

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. 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 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 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 PSD 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 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, 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-∞ [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-∞ 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 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-∞ and aPCI offer some interesting perspectives for a flexible data-driven formu-lation 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 secondgeneration 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.

II. INTERFEROMETRIC MEASUREMENTS
LISA produces two main interferometric measurements per movable optical sub-assembly (MOSA) relevant for laser noise reduction1 .Those are the inter-spacecraft and reference interferometers given by Each measurement represents the beatnote phase formed by two laser beams, whose phases are denoted by ϕ 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 isi ij tracks the phase difference between the distant laser ϕ ji that is propagated to the local MOSA and the local laser ϕ ij .Beam propagation is equivalent to applying the delay operator D ij defined as 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 ϕ ji are defined in their respective reference frame, j.In order to relate phases on the left and right MOSAs on each spacecraft2 , the reference interferometer rfi ij combines the local and the adjacent lasers.
It is useful to decompose the total laser phase ϕ ij (t) or the frequency ν ij (t) = φij (t) into two variables because, as we shall see below, the instrumental and data preprocessing effects we will describe (see the diagram of fig. 1) couple differently to the phase ramp ϕ o ij (t) and any in-band fluctuations ϕ ϵ ij (t).We thus write where ν o ij (t) = φo ij (t) describes any slowly-varying drifts around the central laser frequency ν 0 = 281.6THz and ν ϵ ij (t) = φϵ ij (t) accounts for any rapidly varying random fluctuations.This in-band part is dominated by laser frequency noise with an amplitude spectral density (ASD) S ṗ = 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.Each decimation stage consists of an anti-aliasing filter, F, and a downsampling stage, 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 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], 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 , We use eq.( 5) to represent the output y n in continuous time and find where we introduce the integral over the Dirac-delta distribution δ(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 (τ ) 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 with 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 Again, we make use of eq. ( 5) to find the corresponding continuous time representation, Here, f s denotes the sampling rate after decimation, and 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, 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 onesided PSD, Here, 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 th alias, S x (f ), on the right-hand-side of eq.(11), is given by 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 M f n , which corresponds to the Nyquist rate before decimation).

III. COMMUTATOR RESIDUALS
When computing TDI laser noise residuals, expressions of the form appear, with operators A and B that act on the time series ϕ 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 η 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., F = S = 1 and Dij = D ij ) the variable η 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 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 η 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) , (15) the delay operator D ij switches places with the decimation stage SF such that eq.(1a) becomes This expression has the same algebraic form as eq.( 14) if we make the replacement ϕ ij → SFϕ ij , with the exception that the inter-spacecraft interferometric measurement now also contains the commutator [SF, D ij ]ϕ 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 [SF, D ij ]ϕ ji into a filter-delay commutator and a decimation-delay commutator and compute explicit analytic expressions for each contribution, see section III A and section III B.

A. Flexing-filtering coupling
Let us derive the contribution coming from the commutator [F, D], first described in [15] and dubbed "flexingfiltering coupling".We assume that the delay d(t) is slowly varying over the filter length and that its first derivative ḋ(t) is small.We use eq.( 7), and expand [F, D]ϕ(t) at leading order in ḋ(t) to find where G is a filter operator defined similarly to F in eq. ( 7) with h Because of the orbital dynamics of the LISA constellation, ḋ(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 ḋ to be constant.In such a case the PSD of eq. ( 18) reads For longer observation times, one can use the maximum value of ḋ as given by current predictions for the orbital dynamics of LISA in order to derive an upper bound for this PSD.

B. Decimation-delay commutator
The second commutator appearing in eq. ( 17) is the commutator of the decimation and delay operations, 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 where the modulating factor c n is given by for n even, 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 nonstationary 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, This bound is independent of the delay, and corresponds to the case of full anti-correlation between SDϕ(t) and DSϕ(t) in eq. ( 20).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.

C. 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 precision3 .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, We can then easily rewrite the commutator given in eq. ( 17) in terms of frequency data by replacing every occurrence of the delay operator D by its Doppler equivalent.It reads 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 We note that for the full commutator in eq. ( 17), 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. ( 26), we compare it to numerical simulations that generate time series corresponding to eq. ( 25).We consider white laser frequency noise in ν(t) with an ASD of 30 Hz/ √ Hz4 and neglect the central laser frequency of 2.816 × 10 14 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 (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 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 D. To distinguish between these effects, we write the residual due to the imperfections of D as 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.

A. 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 ϵ ranging from 0 to 1.The latter defines the coefficients of the interpolation kernel k n (ϵ) with −N/2 ≤ n ≤ N/2 − 1, which is convolved with the discrete data samples ϕ n , 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 We define the additional phase residual caused by the interpolation error as and use eq.( 8) to derive the residual in terms of PSD.We find In general, the interpolation kernel k m (ϵ) 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 ϵ scans through different values.At ϵ = 0 and ϵ = 1, the delay is a pure integer shift and the residual vanishes.Assuming that the worst case is obtained for ϵ = 0.5 5 , 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 Hz to 1 Hz) when using high interpolation orders 6 .Alternative interpolation kernels are under study that use less coefficients and optimize their performance over the entire band (and not only at DC).

B. Ranging residual
Any estimate of the delay d(t) differs from the true delay d(t) by a ranging error r(t).We define the corresponding delay operator as where we have assumed r(t) to be small and performed a series expansion to first order.The ranging residual is then given by with ν(t) = φ(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 ϵ (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 ϵ (t) has a root-mean-squared value of ∼ 100 fs (assuming the PSD in eq. ( 71)).Therefore, r ϵ (t) ≪ 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 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 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.

C. Comparison with numerical simulations
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 In order to simulate time series data corresponding to eq. (36), 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 Hz/ √ Hz at a sampling frequency of 4 Hz.Both the post-processing delay D and the propagation delay 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 s, and ranging noise7 with an ASD of 10 −15 s/ √ Hz 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 ϵ = 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)  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 σ 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.
re-expressed for frequency data, i.e., 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
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 η is constructed from the inter-spacecraft and reference interferometers.This step reduces the number of lasers from six to three.Then, the variables π, ρ and σ are constructed from η by building round trip interferome-ters of increasing complexity.We thus have 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 σ we can finally express the secondgeneration Michelson combination X 2 as The contracted post-processing delays D i1•••i N appearing in eq. ( 38) are applied in two steps.First, the nested delay d i1•••i N is calculated, and then, the fractional delay filter is applied.The nested delay can be calculated using the recursive operation8 Alternatively, contracted post-processing delays can be decomposed into atomic delay operations 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.(1b) and (16).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 9 as in [13].The beatnote phase fluctuations are given by Here, p ij denotes the phase noise of laser ij.In the following we write as a short-hand notation pij = SFp ij .Next, we replace each occurrence of the post-processing delay Dij 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 where the origin of the first three terms on the righthand 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 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 δX 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 If one instead performs the interpolation with "atomic delays" (i.e., turning 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 dij inherit a bias B ij from the PRN measurements and a stochastic term from the sideband measurements.Therefore, the ranging error reads where M i denotes the modulation noise on left-handed MOSAs.The resulting residual δX D 2 can be decomposed into the component originating from the ranging bias δX B 2 and the modulation noise δX M 2 .These two contributions are consistent with [13] and read where ṗ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 ∆d.For the second-generation Michelson variable X 2 , the residual reads In the second line we use the property that a delay commutator acts like a derivative, as already described in [4,6,13,28].The arm length mismatch in a secondgeneration Michelson interferometer is given by 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 recover10 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 ν o 12 must be provided to subtract the trend.The stochastic in-band component of the delay commutator is characterized by its PSD For this derivation we have assumed that ∆d X2 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 ∆d X2 up to second order in velocity ḋiji and up to first order in acceleration diji to derive a good approximation of ∆d X2 [6], 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 counterintuitive as locked lasers generate echoes of any pathlength 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-12 11 introduced in [23].

A. 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 with the effective squared delay derivative 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 Here, ∆ (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] with an effective squared bias B2 and a modulating function A M X2 (f ) defined as

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 onboard 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 offset frequency o ij .The so-called locking conditions for the inter-spacecraft and reference locking interferometer are given by where ν ji and ν ik , respectively, denote the frequencies of the reference lasers and ν 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 topolo-gies).For this particular locking configuration the inband phase noise of the six lasers is given by 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 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 12 .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 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 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 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 Therefore, the PSD of ranging error contributions is expressed as eq.( 60) with an effective squared bias equal to

C. Comparison with numerical simulations
We now compare the theoretical model described above to simulated LISA data, obtained using LISA Instrument [24,29] 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 Hz/ √ 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 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 while on the right-handed MOSAs it a factor of 10 larger [31].As already presented in [12,13], we can obtain an estimate of the pseudo-range dij 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 Hz and the filter presented in appendix A is used for antialiasing before decimating to 4 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 MHz to 25 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 Dij are replaced by Ḋij , which are defined analogously to eq. ( 24).In the last processing step, we subtract any outof-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 to 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 Hz to 10 −3 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 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 (ϵ = 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 secondgeneration 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 involved laser.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 pre-processing.It occurs before parameter estimation and reduces spectral leakage in PSD esti-mates, 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 η 12 + D 12 η 21 and η 13 + D 13 η 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 Hz to 2 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 become relevant for the flexingfiltering effect.Furthermore, extra care must be taken when extracting the delay estimates from the sideband measurements since measured pseudo-ranges have an inband 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.
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 ∆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 ∆d.First we recognize that the delay d ij can be written as which has the same algebraic structure as η ij up to a sign.
Here, the time argument t takes the place of the laser phase ϕ i and ϕ j (c.f.eqs.( 1) and ( 38)).Using η ij = d ij as inputs to TDI yields the travel time difference ∆d as where AB and BA denote the paths of the counterpropagating 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 ∆d.To avoid numerical problems we thus operate on the delay derivatives ḋij directly and form ∆ ḋ = Ẋ 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 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, isi ij (t) = D ij p ji (t) − p ij (t), (D1a) rfi ij (t) = p ik (t) − p ij (t). (D1b) For simplicity we omit the anti-aliasing filtering and decimation that was considered in eq. ( 1).Next, we insert eq.(D1) into eq.(38a) where the post-processing delay D is used and only accounts for a stochastic ranging error r(t).We can express η ij for the left and right-handed MOSAs as 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. ( D2) is already a laser noise residual and we neglect higher order couplings in the following.Then, we use the intermediary variables, π ij , ρ ij and σ 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.(D2) that is propagated through TDI as well as the commutator of post-processing delay operators (D5b) Finally, we compute the PSD of eq.(D4) by assuming that all laser and ranging noise terms are uncorrelated and have identical noise properties.We find where we have neglected any cross terms between the first two lines and the last line in eq.(D4).The * sign denotes convolution (in frequency domain) which stems from the time domain products of ranging and laser noise contributions.

Figure 1 .
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

Figure 2 .
Figure2.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.

Figure 4 .
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 σ 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.

Figure 5 .
Figure 5. ASDs of the second-generation Michelson variable Ẋ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 Hz divides the x-axis into logarithmically and linearly scaled.Note that the effects of interpolation appears for frequencies larger than 1 Hz and are shown in fig.7.

Figure 6 .
Figure 6.ASDs of the second-generation Michelson variable Ẋ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 Hz divides the x-axis into logarithmically and linearly scaled.Note that the effects of interpolation appears for frequencies larger than 1 Hz and are shown in fig.7.

Figure 7 .
Figure 7. ASDs of residual noise in second-generation Michelson variable Ẋ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.