# Unified model for the LISA measurements and instrument simulations

Jean-Baptiste Bayle, Olaf Hartwig

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

To cite this version:

Jean-Baptiste Bayle, Olaf Hartwig. Unified model for the LISA measurements and instrument simulations. Physical Review D, 2023, 107 (8), pp.083019. 10.1103/PhysRevD.107.083019. hal-03921038

HAL Id: hal-03921038

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

Submitted on 20 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Unified model for the LISA measurements and instrument simulations

Jean-Baptiste Bayle<sup>1,*</sup> and Olaf Hartwig<sup>2,3</sup>

<sup>1</sup>University of Glasgow, Glasgow G12 8QQ, United Kingdom

<sup>2</sup>SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université,
LNE, 61 avenue de l'Observatoire, 75014 Paris, France

<sup>3</sup>Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut),
Callinstraße 38, 30167 Hannover, Germany

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

(Received 13 December 2022; accepted 17 March 2023; published 11 April 2023)

LISA is a space-based mHz gravitational-wave observatory, with a planned launch in 2034. It is expected to be the first detector of its kind, and will present unique challenges in instrumentation and data analysis. An accurate preflight simulation of LISA data is a vital part of the development of both the instrument and the analysis methods. The simulation must include a detailed model of the full measurement and analysis chain, capturing the main features that affect the instrument performance and processing algorithms. Here, we propose a new model that includes, for the first time, proper relativistic treatment of reference frames with realistic orbits, a model for onboard clocks and clock synchronization measurements, proper modeling of total laser frequencies (including laser locking), frequency planning and Doppler shifts, better treatment of onboard processing, and updated noise models. We then introduce two implementations of this model, LISANode and LISA Instrument. We demonstrate that TDI processing successfully recovers gravitational-wave signals from the significantly more realistic and complex simulated data. LISANode and LISA Instrument are already widely used by the LISA community and, for example, currently provide the mock data for the LISA data challenges.

DOI: 10.1103/PhysRevD.107.083019

I. INTRODUCTION

Following the opening of the gravitational Universe by
the many observations of ground-based gravitational-wave
detectors [1–15], the European Space Agency (ESA) has
selected the Laser Interferometer Space Antenna (LISA) as
the L3 mission. LISA is a space-borne gravitational-wave
observatory sensitive to gravitational signals between
0.1 mHz and 1 Hz, where we expect a large diversity of
sources, ranging from supermassive black-hole binaries,
quasimonochromatic Galactic binaries, extreme mass-ratio
inspirals, and stellar-mass binaries [16]. In addition to these
expected sources, a number of potential signals might be
detected, including stochastic gravitational-wave signals
from the early Universe, cusps and kinks of cosmic strings
and other unmodeled burst sources. Precise measurements
of the source parameters will help answer many astro-
physical and cosmological questions, as well as constrain
models beyond the general theory of relativity.

Achieving these outstanding science objectives will present challenges in both instrumentation and data analysis. Contrary to ground-based gravitational-wave observatories, LISA is expected to be signal dominated, with tens of thousands of sources of different kind present in the LISA band at all times. Telling all of these sources apart and estimating their parameters requires novel approaches to data analysis (explored in the context of the LISA data challenges), the development and testing of which necessitates realistic simulated data. In addition, LISA will make use of sophisticated noise-reduction algorithms to reject the most dominant instrumental-noise sources. The core of these algorithms is known as time-delay interferometry (TDI), in which multiple data streams are combined with appropriate time shifts to generate virtual equal-arm interferometers in postprocessing. Understanding how different noise sources couple into the data is crucial to guide the development of such noise-reduction pipelines. Finally, with a planned launch in 2034, the LISA mission is currently preparing for adoption. The development of a simulation model is needed to support these activities, validate the instrument design and ensure that the science objectives can be achieved.

To fulfill these objectives, one needs to capture in the
simulation model the main features that affect the instrument
performance and processing algorithms. The simulated data

<footnote>j2b.bayle@gmail.com</footnote>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI.</footnote>
should be representative of the time series we will receive from the real instrument. Therefore, we focus in this paper on a time-domain instrument model. In addition, we must be able to simulate several years of data in a reasonable time to evaluate different instrument configurations for full mission duration, currently planned as 4 years [16]; this makes a detailed engineering-level simulation unfeasible.

This instrument model builds on a legacy of previous constellation-level LISA simulators. LISA Simulator was developed to quickly generate measurement data [17,18]. The simulator worked exclusively in the frequency domain and was based on transfer functions for a simple instrumental model. Synthetic LISA was a Python-based simulator that worked in the time domain and used an idealized (and now out-of-date) instrumental configuration to study the performance of noise reduction algorithms for a constellation with time-varying arm lengths [19]. TDISim was a prototype TDI simulation tool programmed in MATLAB. The simulation fully operated in the time domain and performed both data generation and TDI, including for the first time the updated split-interferometry optical bench design and a simplified state-space model for the motion of the test mass and the spacecraft [20].

We based our simulation efforts on LISACode, which was initiated with the similar ambition to include most of the ingredients that were thought to influence LISA's performance at the time [21]. Since then, developments in the instrument and mission design revealed new important effects that must be included in the simulations. The model that we propose in this paper is an attempt to extend LISACode's model to capture those effects.

Section II introduces the conventions we use, and, for the first time, a description of the time frames relevant for LISA Instrument simulations. In Sec. III, we describe the optical simulation model, which includes the up-to-date split-interferometry optical bench design. Contrary to previous simulators, we properly model the total laser frequencies, as well as realistic orbits and any Doppler effects arising from differential spacecraft motion. We also account for the sideband modulations used to correct for clock errors. Then, we describe in Sec. IV the readout of the interferometric beatnotes and how it is affected by imperfections of onboard clocks. Our treatment of the onboard processing is presented in Sec. V; here, we also give the equations for the final phasemeter readouts. In Sec. VI, we describe how we model laser locking, and its impact on the measurements. Lastly, in Sec. VII, we give a high-level model of the pseudoranging measurements that are used to estimate the arm lengths. We then introduce LISANode and LISA Instrument, two implementations of this simulation model, and discuss their performances in Sec. VIII. Finally, in Sec. IX, we show simulation results and highlight the main features that differ from previously simulated data. We demonstrate that despite the added complexity, we can recover

gravitational signals using the latest noise-reduction algorithms. We conclude in Sec. X.

## II. FRAMEWORK AND CONVENTIONS

### A. Constellation overview

LISA is an almost equilateral triangle, composed of three identical spacecraft, which we label 1, 2, 3 clockwise when looking down at their solar panels. These spacecraft exchange laser beams, which are combined on optical benches inside movable optical subassemblies (MOSAs).

To uniquely label these MOSAs, we use two indices. The first one is that of the spacecraft the MOSA is mounted, while the second index is that of the spacecraft the MOSA is pointing to. Most components of interest (such as the optical benches, test masses, etc.) can be uniquely associated to one of the MOSAs, in which case we use the same two indices. Elements that exist only once onboard a spacecraft, such as the ultrastable oscillators (USOs), are indexed by that spacecraft index. These labeling conventions, which are largely based on the proposed unified conventions of the LISA consortium [22], are summarized in Fig. 1.

Quantities that describe a process that involves the propagation between two spacecraft will be interpreted as being associated with the spacecraft in which the quantity is measured. For example, the gravitational-wave signal observed in the interferometer on MOSA 12 will be indexed with the same indices 12. The same convention applies to the propagation delay of a beam arriving on spacecraft 1 from spacecraft 2, which will be labeled by the indices 12.

In this paper, we often derive equations for a specific spacecraft or MOSA. The expressions for the other 2 spacecraft or the other 5 MOSAs can then be deduced by

<figure><img src="image_4.png" /><figcaption>FIG. 1. Labeling conventions used for spacecraft, light travel times, lasers, MOSAs, and interferometric measurements (from [23]).</figcaption></figure>
combining cyclic permutations {1 → 2, 2 → 3, 3 → 1}, and reflections {1 → 1, 2 → 3, 3 → 2}.

## B. Time-coordinate frames

The instrumental simulation mostly concerns itself with the physics inside a spacecraft (e.g., the evolution of laser beam phases and their interferometric beatnotes), which is best modeled in the three spacecraft proper time (TPS). These time frames are defined as the times shown by perfect clocks comoving with the spacecraft centers of mass. We denote them with $\tau_1$, $\tau_2$, and $\tau_3$.

The spacecraft proper times (TPSs) are idealized time-scales, which cannot be realized in practice. All measurements instead refer to an imperfect on-board timer, which represents an approximation of the associated TPS. We denote these three onboard clock time frames as $\hat{\tau}_1$, $\hat{\tau}_2$, and $\hat{\tau}_3$.

Finally, processes on the Solar-System scale are modeled according to a global time frame, such as the barycentric coordinate time (TCB), denoted $t$. This is the case for the spacecraft orbits or the gravitational waveforms. Our instrumental simulation does not make a direct use of the barycentric coordinate time. Instead, we rely on external tools, such as LISA Orbits [24], to directly compute quantities expressed in the TPSs.

In general, signals are expressed in their natural time coordinate. E.g., laser beam phases and beatnotes are expressed in the TPS of the spacecraft housing the laser. It is sometimes useful to express a signal in a different time coordinate. To prevent confusion, we will use the same symbol but add a superscript denoting the time coordinate. For example, a phase $\phi$ could be expressed as a function of the TPS 1, writing $\phi^{\tau_1}(x)$, or as a function of the clock time of that spacecraft, writing $\phi^{\hat{\tau}_1}(x)$. Note that the symbol used for the function argument is arbitrary, and does not specify the reference frame. We will often use $\tau$ without subscripts as a generic function argument.

Conversions between time coordinates can easily be expressed with these conventions. For example, $\tau_1^{\hat{\tau}_1}(\tau)$ is the TPS as a function of the clock time onboard spacecraft 1. Trivially,

$$t^l(\tau) = \tau_1^{\tau_1}(\tau) = \hat{\tau}_1^{\tau_1}(\tau) = \tau. \qquad (1)$$

It is often useful to model the deviation of the onboard clock time with respect to the associated TPS. We adopt the notation

$$\tau_1^{\hat{\tau}_1}(\tau) = \tau + \delta\tau_1^{\hat{\tau}_1}(\tau), \qquad (2a)$$

$$\hat{\tau}_1^{\tau_1}(\tau) = \tau + \delta\hat{\tau}_1^{\tau_1}(\tau). \qquad (2b)$$

One important class of signals we study are phases $\phi$ of electromagnetic waves. As scalar quantities, these are

invariant under coordinate transformations, such that they transform from one time frame to another using a simple time shift,

$$\phi^{\tau_1}(\tau) = \phi^{\hat{\tau}_1}(\hat{\tau}_1^{\tau_1}(\tau)). \qquad (3)$$

## III. OPTICAL SIMULATION

In this section, we derive the model for the generation and propagation of the laser beams, as well as their interference at the photodiodes.

### A. Optical-bench design

As illustrated in Fig. 1, each spacecraft hosts two optical benches. We usually refer to one optical bench as the *local* optical bench; the other optical bench hosted by the same spacecraft as the *adjacent* optical bench; we call the *distant* optical bench the one situated on the spacecraft exchanging light with the local optical bench. Each optical bench is associated with a laser source, a gravitational reference sensor (GRS) containing a free-falling test mass, and telescope to send and collect light to and from distant spacecraft.

Laser beams are combined in three different heterodyne interferometers. The interspacecraft interferometer (ISI) mixes the local beam with the distant beam (coming from the distant optical bench) and the test-mass interferometer (TMI) mixes the local and adjacent beams, after it has bounced on the local test mass; and the reference interferometer (RFI) mixes the local and adjacent beams without interaction with the test mass. Fig. 2 gives an overview of the optical bench 12.

In reality, each single interferometer output is implemented using redundant balanced detection with four quadrant photodiodes (QPDs). We do not simulate balanced detection, and only consider a single data stream for each interferometer. Additional readouts related to the laser

<figure><img src="image_5.png" /><figcaption>FIG. 2. Schematics of the optical design implemented in the simulation, along with notations for the laser beam and beatnote total phases, here for MOSA 12.</figcaption></figure>
beams alignment, such as the differential wave front sensing (DWS), are not included in the model presented here. We are currently working to implement them in the simulation by propagating additional independent variables representing the different beam tilts. We plan to describe this model in more detail in a follow-up paper.

## B. Laser-beam model

### 1. Simple laser beam

We use a number of assumptions to model the information carried by the electromagnetic field of a laser beam (in all generality, these are two 3-vector fields).

We work in the plane-wave approximation, and assume that any effects due to wave front imperfections can be modeled as equivalent longitudinal path length variations. In addition, we neglect effects related to the fields polarization, and assume that the waves propagate in a perfect vacuum, such that we only model the scalar electric field amplitude (we do not model the magnetic field amplitude, as it is completely determined by the electrical field amplitude [25]).

At any fixed point inside a spacecraft, the complex amplitude of the electromagnetic field associated with a laser beam can be written as

$$E(\tau) = E_0(\tau)e^{i2\pi\Phi(\tau)}, \quad (4)$$

where $\Phi(\tau)$ is the total phase in units of cycles.

LISA ultimately measures phase differences, such that we do not simulate the field amplitude term $E_0$, but only the phase $\Phi(\tau)$. We expect couplings between the field amplitude and the measured phase difference (e.g., the relative intensity noise [26]). We currently do not model these effects, but assume that they can be modeled as equivalent phase noise in the final readout.

### 2. Phase or frequency?

The optical frequency of the lasers is around $\nu_0 = 281.6$ THz, such that the total phase increases quickly with time. This makes using it challenging for numerical simulations, as any variable representing the total phase will either numerically overflow when using fixed-point arithmetic, or eventually suffer an unacceptable loss of precision when using floating-point arithmetic.<fnref n="1" />

To avoid these issues, we simulate frequencies instead of phase (given by $\nu = \dot{\Phi}$, since we express the phase in units of cycles), which are controlled to remain at the same order of magnitude during the whole mission duration. However, modeling the propagation of laser beams is often easier in

phase. Therefore, we will derive most of the equations of this paper both in units of phase and frequency.

### 3. Two-variable decomposition

In LISA, effects on the laser beams come into play at completely different timescales and dynamic ranges. On the one hand, some effects modulate the frequency of our beams on a time scale of the orbital revolution around the Sun, which lies well outside our measurement frequency band (below $10^{-4}$ Hz). These effects tend to have large dynamic ranges; for instance, the Doppler shifts caused by the relative spacecraft motion can fluctuate by several megahertz over the mission duration.

On the other hand, we want to track small phase or frequency fluctuations within our measurement band (from $10^{-4}$ Hz and up to 1 Hz), caused by gravitational-wave signals and instrumental noises. These fluctuations have a much smaller amplitude. The laser noise being the dominant effect, causing the heterodyne beatnote frequency to shift by about a few tens of Hertz, while gravitational waves typically cause frequency shifts of a few hundreds of nanohertz.

To address this problem, we model these different effects independently. We decompose the laser beam frequency into one constant and two variables,

$$\nu(\tau) = \nu_0 + \nu^o(\tau) + \nu^e(\tau). \quad (5)$$

The large frequency offsets $\nu^o(\tau)$ are used to represent frequency-plan offsets and Doppler shifts, both on the order of megahertz, as well as the gigahertz sideband frequency offsets. The small frequency fluctuations $\nu^e(\tau)$, on the other hand, are used to describe gravitational signals and noises. A simple laser beam would therefore be entirely represented by the couple $\{\nu^o(\tau), \nu^e(\tau)\}$.

Alternatively, we can express Eq. (5) in phase units by writing the total phase as

$$\Phi(\tau) = \nu_0\tau + \phi^o(\tau) + \phi^e(\tau) + \phi^0, \quad (6)$$

where the definitions of large phase drifts $\phi^o(\tau)$ and small phase fluctuations $\phi^e(\tau)$ follow from Eq. (5),

$$\nu^o(\tau) = \dot{\phi}^o(\tau) \quad \text{and} \quad \nu^e(\tau) = \dot{\phi}^e(\tau). \quad (7)$$

As we simulate frequencies, we do not track the initial phase of the laser beam $\phi^0 \in [0, 2\pi]$ in the following.

Let us stress that this decomposition is entirely artificial. In reality, we will only have access to the total phase or frequency. Therefore, to produce data representative of the real instrument telemetry, we always compute the total phase or frequency as the final simulation output.

<footnote n="1">¹For a precision better than a micro-cycle, a 64-bit integer representing the total phase will overflow every 0.07 s.</footnote>
### 4. Modulated beams

In LISA, laser beams are phase-modulated using a gigahertz signal derived from the local clock. The electric field reads

$$E(\tau) = E_0 e^{i2\pi\Phi_c(\tau)} e^{im\cos(2\pi\Phi_m(\tau))}, \quad (8)$$

where m is the modulation depth; $\Phi_c(\tau)$ is the total phase of the carrier, and $\Phi_m(\tau)$ is the total phase of the modulating signal, both expressed in cycles.

The Jacobi-Anger expansion lets us write the previous expression using Bessel functions. Because the modulation depth $m \approx 0.15$ is small [20], we can further expand the result to first order in m and write the complex field amplitude as the sum

$$E(\tau) \approx E_0 \left( e^{i2\pi\Phi_c(\tau)} + i \frac{m}{2} \left[ e^{i2\pi\Phi_{sb^+}(\tau)} + e^{i2\pi\Phi_{sb^-}(\tau)} \right] \right), \quad (9)$$

where we have defined the upper and lower sideband phases,

$$\Phi_{sb^+}(\tau) = \Phi_c(\tau) + \Phi_m(\tau), \quad (10a)$$

$$\Phi_{sb^-}(\tau) = \Phi_c(\tau) - \Phi_m(\tau). \quad (10b)$$

The modulated laser beam can then be written as the superposition of carrier, upper sideband, and lower sideband,

$$E(\tau) \approx E_c(\tau) + E_{sb^+}(\tau) + E_{sb^-}(\tau). \quad (11)$$

For the purpose of our simulation, the information content of the upper and lower sidebands are almost identical (one difference is that they lie at a different frequencies, and are thus affected differently by Doppler shifts). We make the assumption that they can be combined in such a way that we can treat them as one signal. Therefore, we only simulate the upper sideband. For clarity, we drop the sign in all sideband indices and simply use sb when we refer to the upper sideband.

We apply the same two-variable decomposition to the sideband total frequency. Ultimately, each modulated laser beam is then implemented using four quantities,

$$\nu(\tau) \equiv \{\nu_c^o(\tau), \nu_c^e(\tau), \nu_{sb}^o(\tau), \nu_{sb}^e(\tau)\}, \quad (12)$$

where $\nu_c^o$ and $\nu_{sb}^o$ are the carrier and sideband frequency offsets, respectively, and $\nu_c^e$ and $\nu_{sb}^e$ are the carrier and sideband frequency fluctuations.

### C. Local beams

#### 1. Local beam at laser source

As illustrated in Fig. 2, optical bench 12 has an associated laser source. We call *local beam* the modulated beam produced by this laser source. We denote the total phase and frequency of the carrier as $\Phi_{12,c}(\tau)$ and $\nu_{12,c}(\tau)$, respectively. Similarly, the sideband total phase and frequency are denoted as $\Phi_{12,sb}(\tau)$ and $\nu_{12,sb}(\tau)$. All these signals are functions of the spacecraft proper time (TPS) $\tau_1$.

The total phase $\Phi_{12,c}(\tau) = \nu_0\tau + \phi_{12,c}^o(\tau) + \phi_{12,c}^e(\tau)$ of the carrier is decomposed in terms of drifts and fluctuations, with

$$\phi_{12,c}^o(\tau) = \int_{\tau_{1,0}}^{\tau} O_{12}(\tau')d\tau', \quad (13a)$$

$$\phi_{12,c}^e(\tau) = p_{12}(\tau), \quad (13b)$$

where $O_{12}(\tau)$ is the carrier frequency offset for this laser source with respect to the central frequency $\nu_0$, and $p_{12}(\tau)$ is the laser source phase fluctuations expressed in cycles.

As explained in Sec. VI, $p_{ij}(\tau)$ can either describe the noise $N_{ij}^p(\tau)$ of a cavity-stabilized laser (cf. Appendix B) or the fluctuations resulting from an offset frequency lock. Likewise, $O_{12}(\tau)$ is either set as an offset from the nominal frequency,<fnref n="2" /> or computed based on the locking conditions. In terms of frequency, we simply have

$$\nu_{12,c}^o(\tau) = O_{12}(\tau), \quad (14a)$$

$$\nu_{12,c}^e(\tau) = \dot{p}_{12}(\tau). \quad (14b)$$

Let us now look at the sideband, which is derived from the local clock. As described in detail in Sec. IV, the modulating signal inherits any USO timing errors $q_1$, such that we have

$$\Phi_{12,m}(\tau) = \nu_{12}^m \cdot (\tau + q_1^o(\tau) + q_1^e(\tau) + M_{12}(\tau)) \quad (15)$$

for the total phase of the modulating signal. Here, $\nu_{12}^m = 2.4$ GHz is the constant nominal frequency of the modulating signal on optical bench 12. We use the same modulation frequency for all optical benches indexed cyclically (12, 23, 31), while the remaining ones (13, 32, 21) are instead at 2.401 GHz. The modulation noise term $M_{12}(\tau)$ accounts for any additional imperfections (either in the electrical frequency conversion to 2.4 GHz or the optical modulation).

<footnote n="2">²In the current mission baseline, there is no way to measure or set the absolute laser frequency with high precision. Therefore, the values set in the simulation cannot be accessed in reality.</footnote>
The total phase of the modulating signal can then be decomposed into

$$ \phi_{12}^{o}(\tau) = \nu_{12}^{m} \cdot (\tau + q_{1}^{o}(\tau)), \quad (16a) $$

$$ \phi_{12,m}^{\epsilon}(\tau) = \nu_{12}^{m} \cdot (q_{1}^{\epsilon}(\tau) + M_{12}(\tau)). \quad (16b) $$

Inserting these terms in Eq. (10), we get the phase and frequency offsets and fluctuations for the local sideband,

$$ \phi_{12,sb}^{o}(\tau) = \int_{\tau_{1,0}}^{\tau} O_{12}(\tau')d\tau' + \nu_{12}^{m}(\tau + q_{1}^{o}(\tau)), \quad (17a) $$

$$ \phi_{12,sb}^{\epsilon}(\tau) = p_{12}(\tau) + \nu_{12}^{m}(q_{1}^{\epsilon}(\tau) + M_{12}(\tau)), \quad (17b) $$

and

$$ \nu_{12,sb}^{o}(\tau) = O_{12}(\tau) + \nu_{12}^{m}(1 + \dot{q}_{1}^{o}(\tau)), \quad (18a) $$

$$ \nu_{12,sb}^{\epsilon}(\tau) = \dot{p}_{12}(\tau) + \nu_{12}^{m}(\dot{q}_{1}(\tau) + \dot{M}_{12}(\tau)). \quad (18b) $$

Note that there is only one clock per spacecraft, such that we use the same $q_1$ for sideband beams on both optical benches on spacecraft 1.

## 2. Local beams at the interspacecraft and reference interferometer photodiodes

As shown in Fig. 2, local beams propagate in the local optical bench 12 and interfere at the ISI, TMI, and RFI photodiodes. In our simulations, we neglect any phase term due to the propagation time. However, all beams pick up a generic optical path length noise term $N^{ob}(\tau)$ (different for each interferometer), which models all optical path length variations due to, e.g., jitters of optical components in the path of the laser beams. By convention, we choose that a positive value of the optical path length noise term corresponds to a decrease in the actual optical path length, which in turn corresponds to a positive shift in phase or frequency.

Therefore, we write the phase drifts and fluctuations of the local beams at the ISI and RFI photodiodes (valid for both carriers and sidebands) as

$$ \phi_{\text{isi/rfi}_{12} \leftarrow 12}^{o}(\tau) = \phi_{12}^{o}(\tau), \quad (19a) $$

$$ \phi_{\text{isi/rfi}_{12} \leftarrow 12}^{\epsilon}(\tau) = \phi_{12}^{\epsilon}(\tau) + \frac{\nu_0}{c} N_{\text{isi/rfi}_{12} \leftarrow 12}^{\text{ob}}(\tau). \quad (19b) $$

Equivalently, the frequency offsets and fluctuations of the same beams read

$$ \nu_{\text{isi/rfi}_{12} \leftarrow 12}^{o}(\tau) = \nu_{12}^{o}(\tau), \quad (20a) $$

$$ \nu_{\text{isi/rfi}_{12} \leftarrow 12}^{\epsilon}(\tau) = \nu_{12}^{\epsilon}(\tau) + \frac{\nu_0}{c} \dot{N}_{\text{isi/rfi}_{12} \leftarrow 12}^{\text{ob}}(\tau). \quad (20b) $$

<figure><img src="image_6.png" /><figcaption>FIG. 3. Definition of line-of-sight displacement of the test mass. Positive values of $N_{12}^{\delta}$ indicate a motion towards the measuring optical bench.</figcaption></figure>

## 3. Local beam at the test-mass interferometer photodiode

The local beam reflects off the test mass before impinging on the TMI photodiode. As a consequence, it couples to the test-mass motion.

In reality, the motion of the test mass and spacecraft will be coupled by the drag-free attitude control system (DFACS). The spacecraft motion is expected to be suppressed in on-ground processing [20]. For our purposes, we simply assume that the spacecraft (and the associated optical benches) perfectly follows a geodesic.

The laser beam then picks up an additional noise term $N_{23}^{\delta}(\tau)$ due to any deviation in the motion of the test-mass from geodesic, caused by spurious forces (cf. Appendix B). This noise represents the movement of the test mass towards the measuring optical bench, such that a positive value corresponds to a decrease in path length (see Fig. 3), and thus a positive phase shift. The noise term enters with a factor of 2, since the beam travels to the test mass and back.

Therefore, at the TMI photodiode, the phase components of the local beam (carrier and sideband) read

$$ \phi_{\text{tmi}_{12} \leftarrow 12}^{o}(\tau) = \phi_{12}^{o}(\tau), \quad (21a) $$

$$ \phi_{\text{tmi}_{12} \leftarrow 12}^{\epsilon}(\tau) = \phi_{12}^{\epsilon}(\tau) + \frac{\nu_0}{c} \left( N_{\text{tmi}_{12} \leftarrow 12}^{\text{ob}}(\tau) + 2N_{12}^{\delta}(\tau) \right), \quad (21b) $$

while the frequency offsets and fluctuations read

$$ \nu_{\text{tmi}_{12} \leftarrow 12}^{o}(\tau) = \nu_{12}^{o}(\tau), \quad (22a) $$
$$ \nu_{\text{tmi}_{12} \leftarrow 12}^{\epsilon}(\tau) = \nu_{12}^{\epsilon}(\tau) + \frac{\nu_0}{c} \left( \dot{N}_{\text{tmi}_{12} \leftarrow 12}^{\text{ob}}(\tau) + 2\dot{N}_{12}^{\delta}(\tau) \right). \qquad (22b) $$

### D. Adjacent beams

In this section, we study the propagation of a modulated laser beam generated by laser source 13 (attached to the adjacent optical bench), which travels through the optical fiber to the local optical bench 12, to finally interfere on the TMI and RFI photodiodes (see Fig. 2). We call it the *adjacent beam*. We express all phase and frequency quantities as functions of the TPS $\tau_1$.

Similarly to local beams, we neglect the propagation time for the adjacent beams, and model fluctuations in the optical path length by a noise term $N^{\text{ob}}(\tau)$. We model any nonreciprocal noise terms related to the propagation through the optical fibers by the backlink noise term $N_{12 \leftarrow 13}^{\text{bl}}(\tau)$, expressed as an equivalent path length change.

Therefore, the phase drifts and fluctuations of adjacent beams (carrier and sideband) at the TMI and RFI photodiodes read

$$ \phi_{\text{tmi/rfi}_{12} \leftarrow 13}^{o}(\tau) = \phi_{13}^{o}(\tau), \qquad (23a) $$

$$ \phi_{\text{tmi/rfi}_{12} \leftarrow 13}^{\epsilon}(\tau) = \phi_{13}^{\epsilon}(\tau) + \frac{\nu_0}{c} N_{12 \leftarrow 13}^{\text{bl}}(\tau) + \frac{\nu_0}{c} N_{\text{tmi/rfi}_{12} \leftarrow 13}^{\text{ob}}, \qquad (23b) $$

where $\phi_{13}^{o}(\tau)$ and $\phi_{13}^{\epsilon}(\tau)$ are the phase drifts and fluctuations of the laser beam produced by laser source 13, respectively. The equivalent frequency quantities are

$$ \nu_{\text{tmi/rfi}_{12} \leftarrow 13}^{o}(\tau) = \nu_{13}^{o}(\tau), \qquad (24a) $$

$$ \nu_{\text{tmi/rfi}_{12} \leftarrow 13}^{\epsilon}(\tau) = \nu_{13}^{\epsilon}(\tau) + \frac{\nu_0}{c} \dot{N}_{12 \leftarrow 13}^{\text{bl}}(\tau) + \frac{\nu_0}{c} \dot{N}_{\text{tmi/rfi}_{12} \leftarrow 13}^{\text{ob}}. \qquad (24b) $$

### E. Distant beams

Finally, we study the propagation of a modulated laser beam generated by laser source 21 (attached to the distant optical bench), which travels roughly 2.5 million kilometers in free space before it reaches the local optical bench 12. This *distant beam* eventually interferes on the ISI photodiode, see Fig. 2.

#### 1. Interspacecraft propagation

As described in Sec. III B, modulated beams are represented as the superposition of simple beams, each treated independently. Consequently, the same propagation equations apply to both carrier and sideband beams.

We shall derive the expression of a simple laser beam's phase $\Phi_{12 \leftarrow 21}(\tau)$ and frequency $\nu_{12 \leftarrow 21}(\tau)$ measured on receiver optical bench 12 (expressed in comoving time coordinate $\tau_1$) as a function of the same beam's phase $\Phi_{21}(\tau)$ and frequency $\nu_{21}(\tau)$ measured on emitter optical bench 21 (expressed in comoving time coordinate $\tau_2$). We write

$$ \Phi_{12 \leftarrow 21}(\tau) = \Phi_{21}(\tau - d_{12}(\tau)), \qquad (25) $$

where $d_{12}(\tau)$ is the proper pseudorange (PPR), which includes not only the light time of flight, but also conversions between reference frames associated to $\tau_1$ and $\tau_2$.

Since we model small in-band and large out-of-band effects independently, we need to decompose the PPR in a similar manner. We define $d_{12}^{o}(\tau)$ as slowly varying PPR offsets (e.g., due to constant path lengths and variations in orbital motion, relativistic effects, and coordinate transformations), and $d_{12}^{\epsilon}(\tau)$ as small in-band PPR fluctuations.

In our simulation, we only consider the effect of gravitational waves and neglect any other small in-band fluctuations of the PPRs (such as spacecraft jitter motion or variations of the interplanetary medium optical index). Therefore, if $H_{12}(\tau)$ denotes the integrated fluctuations of the PPR due to gravitational waves measured on MOSA 12, we have $d_{12}^{\epsilon}(\tau) = H_{12}(\tau)$. The total PPR now reads

$$ d_{12}(\tau) = d_{12}^{o}(\tau) + H_{12}(\tau). \qquad (26) $$

Applying this decomposition to Eqs. (6) and (25), we have

$$ \Phi_{12 \leftarrow 21}(\tau) = \nu_0 \cdot (\tau - d_{12}^{o}(\tau) - H_{12}(\tau)) \\ + \phi_{21}^{o}(\tau - d_{12}^{o}(\tau) - H_{12}(\tau)) \\ + \phi_{21}^{\epsilon}(\tau - d_{12}^{o}(\tau) - H_{12}(\tau)). \qquad (27) $$

We expand the previous equation to first order in both the small fluctuations $H_{12}(\tau)$, and $\phi_{21}^{\epsilon}(\tau)$ and neglect any second-order cross terms,

$$ \Phi_{12 \leftarrow 21}(\tau) = \nu_0 \cdot (\tau - d_{12}^{o}(\tau) - H_{12}(\tau)) \\ + \phi_{21}^{o}(\tau - d_{12}^{o}(\tau)) \\ - \nu_{21}^{o}(\tau - d_{12}^{o}(\tau))H_{12}(\tau) \\ + \phi_{21}^{\epsilon}(\tau - d_{12}^{o}(\tau)). \qquad (28) $$

We can again write the previous quantity as the sum of large phase drifts and small phase fluctuations, $\Phi_{12 \leftarrow 21}(\tau) = \nu_0\tau + \phi_{12 \leftarrow 21}^{o}(\tau) + \phi_{12 \leftarrow 21}^{\epsilon}(\tau)$, with

$$ \phi_{12 \leftarrow 21}^{o}(\tau) = \phi_{21}^{o}(\tau - d_{12}^{o}(\tau)) - \nu_0 d_{12}^{o}(\tau), \qquad (29a) $$

$$ \phi_{12 \leftarrow 21}^{\epsilon}(\tau) = \phi_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau)) \\ - [\nu_0 + \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))] H_{12}(\tau). \qquad (29b) $$
We write the equivalent instantaneous frequency $\nu_{12 \leftarrow 21}(\tau) = \nu_0 + \nu_{12 \leftarrow 21}^o(\tau) + \nu_{12 \leftarrow 21}^e(\tau)$ as the sum of a large frequency offsets and small frequency fluctuations,

$$ \nu_{12 \leftarrow 21}^o(\tau) = \nu_{21}^o(\tau - d_{12}^o(\tau))(1 - \dot{d}_{12}^o(\tau)) - \nu_0 \dot{d}_{12}^o(\tau), \quad (30a) $$

$$ \nu_{12 \leftarrow 21}^{\epsilon}(\tau) = \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))(1 - \dot{d}_{12}^{\epsilon}(\tau)) \\ - [\nu_0 + \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))] \dot{H}_{12}(\tau). \quad (30b) $$

Here, we have neglected first-order terms in $\dot{\nu}_A^o H_{12}(\tau)$, so these equations are only valid if the laser frequency is evolving slowly. This is discussed in more detail in Appendix C.

## 2. Distant beams at the interspacecraft interferometer photodiode

The received distant beam propagates inside the optical bench to interfere with the local beam at the ISI photodiode. As for the other beams, we only add a generic optical path length noise term $N_{isi12 \leftarrow 21}^{\text{ob}}(\tau)$.

We write the phase drifts and fluctuations of the distant beam at the ISI photodiode (valid for both carrier and sideband) as

$$ \phi_{isi12 \leftarrow 21}^{o}(\tau) = \phi_{21}^{o}(\tau - d_{12}^{o}(\tau)) - \nu_{0} d_{12}^{o}(\tau), \quad (31a) $$

$$ \phi_{isi12 \leftarrow 21}^{\epsilon}(\tau) = \phi_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau)) \\ - [\nu_0 + \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))] H_{12}(\tau) \\ + \frac{\nu_0}{c} N_{isi12 \leftarrow 21}^{\text{ob}}(\tau). \quad (31b) $$

Equivalently, the frequency offsets and fluctuations read

$$ \nu_{isi12 \leftarrow 21}^{o}(\tau) = \nu_{21}^{o}(\tau - d_{12}^{o}(\tau))(1 - \dot{d}_{12}^{o}(\tau)) - \nu_{0} \dot{d}_{12}^{o}(\tau), \quad (32a) $$

$$ \nu_{isi12 \leftarrow 21}^{\epsilon}(\tau) = \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))(1 - \dot{d}_{12}^{\epsilon}(\tau)) \\ - [\nu_0 + \nu_{21}^{\epsilon}(\tau - d_{12}^{\epsilon}(\tau))] \dot{H}_{12}(\tau) \\ + \frac{\nu_0}{c} \dot{N}_{isi12 \leftarrow 21}^{\text{ob}}(\tau). \quad (32b) $$

## F. Interferometers

### 1. Beatnote for simple beams

Using definitions given in Eq. (4), let us write the complex amplitude for two simple beams 1 and 2 interfering at a photodiode,

$$ E_1(\tau) = E_{1,0}(\tau) e^{i2\pi\Phi_1(\tau)}, \quad (33a) $$

$$ E_2(\tau) = E_{2,0}(\tau) e^{i2\pi\Phi_2(\tau)}. \quad (33b) $$

We ignore any effects due to spatial dimensions of the beam or the photodiode, and assume that such effects will be modeled as either an equivalent phase error in the readout signal, or as an independent quantity.<fnref n="3" />

The power of the total electromagnetic field measured near the photodiode is

$$ P(\tau) \propto |E_1(\tau) + E_2(\tau)|^2. \quad (34) $$

Substituting the expressions of the two beams yields

$$ P(\tau) \propto |E_{1,0}(\tau)|^2 + |E_{2,0}(\tau)|^2 \\ + 2E_{1,0}(\tau)E_{2,0}(\tau) \cos(2\pi(\Phi_1(\tau) - \Phi_2(\tau))). \quad (35) $$

The power near the photodiode has an oscillating component with a total phase of $\Phi_{\text{PD}}(\tau) = \Phi_1(\tau) - \Phi_2(\tau)$. We call this signal the *beatnote*.

Let us use the two-variable representation described in Eq. (6),

$$ \Phi_1(\tau) = \nu_0\tau + \phi_1^o(\tau) + \phi_1^e(\tau), \quad (36a) $$

$$ \Phi_2(\tau) = \nu_0\tau + \phi_2^o(\tau) + \phi_2^e(\tau), \quad (36b) $$

to express the total phase of the *beatnote* as the sum of large phase drifts and small phase fluctuations,

$$ \Phi_{\text{PD}}(\tau) = \phi_{\text{PD}}^{o}(\tau) + \phi_{\text{PD}}^{e}(\tau), \quad (37) $$

with

$$ \phi_{\text{PD}}^{o}(\tau) = \phi_{1}^{o}(\tau) - \phi_{2}^{o}(\tau), \quad (38a) $$

$$ \phi_{\text{PD}}^{e}(\tau) = \phi_{1}^{e}(\tau) - \phi_{2}^{e}(\tau). \quad (38b) $$

We simulate the equivalent instantaneous frequency defined as $\nu_{\text{PD}}(\tau) = \dot{\Phi}_{\text{PD}}(\tau)$. It can be written as

$$ \nu_{\text{PD}}(\tau) = \nu_{\text{PD}}^{o}(\tau) + \nu_{\text{PD}}^{e}(\tau), \quad (39) $$

where the *beatnote* frequency offsets $\nu_{\text{PD}}^o(\tau)$ and the *beatnote* frequency fluctuations $\nu_{\text{PD}}^e(\tau)$ are defined by

$$ \nu_{\text{PD}}^{o}(\tau) = \nu_{1}^{o}(\tau) - \nu_{2}^{o}(\tau), \quad (40a) $$

$$ \nu_{\text{PD}}^{e}(\tau) = \nu_{1}^{e}(\tau) - \nu_{2}^{e}(\tau). \quad (40b) $$

### 2. Beatnote polarity

A closer look at Eq. (35) shows that we do not have direct access to the total phase of the *beatnote* $\Phi_{\text{PD}}(\tau)$, but

<footnote n="3">³For example, DWS could be modeled as a direct measurement of beam tilt angles, with all beam angles represented by independent variables.</footnote>
only measure its cosine value. Therefore, the total phase can only be known up to a sign and a multiple of $2\pi$.

Physically, this sign ambiguity corresponds to the fact that the electrical signal does not contain any information about which of the two interfering laser beams is of higher frequency. In practice, however, the beatnote polarity can be determined at all times by applying a known frequency offset on the local laser beam and observing the resulting change in the beatnote frequency. In addition, once all lasers are locked, the beatnote polarities can simply be read from the frequency plan, as described in Sec. VI.

Therefore, we do not include the beatnote polarity ambiguity in our optical models, and we will instead assume that it is solved directly by the phasemeter, or in a first processing step on ground.

### 3. Beatnotes for modulated beams

We now study the electromagnetic field of two interfering modulated beams, labeled $k = 1, 2$. As derived in Sec. III B, we write both modulated beams as the sum of three independent simple beams, namely the carriers and the upper and lower sidebands,

$$E_k(\tau) = E_{k,c}(\tau) + E_{k,sb^+}(\tau) + E_{k,sb^-}(\tau), \quad (41)$$

with total phases

$$\Phi_{k,c}(\tau) = \nu_0\tau + \Phi_{k,c}^o(\tau) + \Phi_{k,c}^e(\tau), \quad (42a)$$

$$\Phi_{k,sb^+}(\tau) = \nu_0\tau + \Phi_{k,sb^+}^o(\tau) + \Phi_{k,sb^+}^e(\tau), \quad (42b)$$

$$\Phi_{k,sb^-}(\tau) = \nu_0\tau + \Phi_{k,sb^-}^o(\tau) + \Phi_{k,sb^-}^e(\tau). \quad (42c)$$

or the equivalent instantaneous frequencies

$$\nu_{k,c}(\tau) = \nu_0 + \nu_{k,c}^o(\tau) + \nu_{k,c}^e(\tau), \quad (43a)$$

$$\nu_{k,sb^+}(\tau) = \nu_0 + \nu_{k,sb^+}^o(\tau) + \nu_{k,sb^+}^e(\tau), \quad (43b)$$

$$\nu_{k,sb^-}(\tau) = \nu_0 + \nu_{k,sb^-}^o(\tau) + \nu_{k,sb^-}^e(\tau). \quad (43c)$$

The total power at the photodiode reads

$$|E_1(\tau) + E_2(\tau)|^2 = |E_{1,c}(\tau) + E_{1,sb^+}(\tau) + E_{1,sb^-}(\tau) + E_{2,c}(\tau) + E_{2,sb^+}(\tau) + E_{2,sb^-}(\tau)|^2. \quad (44)$$

Expanding this expression yields cross terms between all 6 terms, which correspond to beatnotes at their difference frequencies.

Because the sidebands are modulated at a frequency of about 2.4 GHz, most of these beatnote frequencies lie far outside of the phasemeters measurement bandwidth (approximately 5 MHz to 25 MHz).

Only three beatnotes lie inside this region:

(i) The carrier-carrier beatnote,

$$\Phi_{\text{PD},c}(\tau) = \Phi_{1,c}(\tau) - \Phi_{2,c}(\tau), \quad (45a)$$

$$\nu_{\text{PD},c}(\tau) = \nu_{1,c}(\tau) - \nu_{2,c}(\tau). \quad (45b)$$

(ii) The upper sideband-upper sideband beatnote,

$$\Phi_{\text{PD},sb^+}(\tau) = \Phi_{1,sb^+}(\tau) - \Phi_{2,sb^+}(\tau), \quad (46a)$$

$$\nu_{\text{PD},sb^+}(\tau) = \nu_{1,sb^+}(\tau) - \nu_{2,sb^+}(\tau). \quad (46b)$$

(iii) The lower sideband-lower sideband beatnote,

$$\Phi_{\text{PD},sb^-}(\tau) = \Phi_{1,sb^-}(\tau) - \Phi_{2,sb^-}(\tau), \quad (47a)$$

$$\nu_{\text{PD},sb^-}(\tau) = \nu_{1,sb^-}(\tau) - \nu_{2,sb^-}(\tau). \quad (47b)$$

Because the sidebands of the lasers on MOSAs 12, 23, and 31 (respectively 13, 32, and 21) are offset by 2.4 GHz (respectively, 2.401 GHz), and because we always interfere beams from both types of MOSA, these three beatnotes will always be offset by 1 MHz. Therefore, they can be tracked individually by the phasemeter.

Each of these beatnote frequencies can be decomposed again as a sum of large frequency offsets and small fluctuations, and we recover equations similar to Eq. (40). Therefore, the carrier and sideband parts of a modulated laser beam can be implemented as three distinct beams in the simulation, from which we form three beatnotes.

As described in the previous sections, we only include the carrier and upper-sideband laser beams in our model; as a consequence, we only compute the carrier-carrier and the upper sideband-upper sideband beatnotes.

### 4. Interspacecraft, test-mass, and reference interferometer beatnotes

To obtain the beatnote phases (or frequencies) measured by the ISI, TMI, and RFI, we can substitute in the previous equations the phases (or frequencies) of the interfering beams.

As discussed above, the beatnote polarities are arbitrary. As a convention, we will always write the beatnote phase (and frequencies) as the difference of the distant or adjacent beam phase (or frequency) and the local beam phase (or frequency),

$$\phi(\tau) = \phi_{\text{distant/adjacent}}(\tau) - \phi_{\text{local}}(\tau), \quad (48a)$$

$$\nu(\tau) = \nu_{\text{distant/adjacent}}(\tau) - \nu_{\text{local}}(\tau). \quad (48b)$$

Following the optical-bench design of Fig. 2, we have the following beatnote phase offsets and fluctuations, for both carriers and sidebands,
$$ \phi_{isi_{12}}(\tau) = \phi_{isi_{12} \leftarrow 21}(\tau) - \phi_{isi_{12} \leftarrow 12}(\tau), \quad (49a) $$

$$ \phi_{tmi_{12}}(\tau) = \phi_{tmi_{12} \leftarrow 13}(\tau) - \phi_{tmi_{12} \leftarrow 12}(\tau), \quad (49b) $$

$$ \phi_{rfi_{12}}(\tau) = \phi_{rfi_{12} \leftarrow 13}(\tau) - \phi_{rfi_{12} \leftarrow 12}(\tau), \quad (49c) $$

and similarly for beatnote frequencies.

## IV. PHASE READOUT, FREQUENCY DISTRIBUTION, AND CLOCK ERROR MODELING

We show in Fig. 4 an overview of the LISA phase readout chain, adapted from [27] (where technical details on the phase readout and frequency distribution system can be found). The optical beatnotes are converted to electrical signals by photoreceivers, which are then digitized by an analog-to-digital converter (ADC). The phase of these digital signals are then tracked by digital phase-locked loops (DPLLs).

The phasemeter is driven by an 80 MHz clock signal. Inside the phasemeter, the ADC samples the electrical beatnotes at the same rate, with an additional timing jitter intrinsic to the ADC. This ADC jitter results in a phase error in the measured beatnotes, which is expected above the requirements. To correct for the ADC jitter, an additional periodic pilot tone signal at 75 MHz is derived from the on-board clock and superimposed on each electrical signal fed to the phasemeter. The pilot tone phase is tracked alongside the main beatnotes in dedicated DPLL channels. By comparing the measured pilot tone phase against its nominal 75 MHz value, the ADC jitter can be corrected in the main beatnotes, such that the pilot tone becomes the effective reference clock signal for the phase measurements.

### A. Readout noise

We directly simulate the optical beatnote frequencies. Our simulated electrical signals are therefore the same quantities, with the addition of a readout noise term $N^{ro}(\tau)$. This readout noise accounts for both shot noise and any

<figure><img src="image_7.png" /><figcaption>FIG. 4. Overview of the phase readout chain of LISA, adapted from [27]. The photoreceiver converts the optical signal to electric signals, which are then digitized by the ADC. The intrinsic timing jitter of the ADC is corrected by a 75 MHz pilot tone superimposed to the photoreceiver signal. The phasemeter DPLLs then tracks the beatnote phases.</figcaption></figure>

errors due to front-end electronic in the photoreceivers; refer to Appendix B for more details.

### B. Phasemeter and pilot tone

Because the photoreceiver signals are already simulated as discrete beatnote frequency samples, we do not directly simulate the digitization process of the ADC nor the phase tracking by the DPLLs.

Furthermore, we do not simulate the pilot tone correction, but assume that it perfectly removes the ADC jitter. We do account for timing errors in the pilot tone itself, which are also expected above the requirements. These pilot tone errors will be corrected using the sidebands introduced in Sec. III B. Refer to the next Sec. IV for how we model the pilot tone and sideband signals.

### C. Frequency distribution and clock signals

#### 1. Frequency-distribution scheme

Most subsystems on board LISA are driven by timing signals derived from the USO. In our simulation model, we focus on processes for which timing is performance critical, which are summarized in Fig. 5.

Following the current mission design, each LISA spacecraft uses one dedicated clock (realized by an USO), from which all timing signals are derived. As described above, the timing reference for all phasemeter measurements is the pilot tone, which is derived from the USO by first up-converting its nominal frequency<fnref n="4" /> to 2.4 GHz, and then converting that signal to the desired $\nu_{PT} = 75$ MHz using frequency dividers. This conversion chain allows for a very stable phase relationship between the electrical pilot tone and the 2.4 GHz optical sideband [27], which are used in post-processing to reduce the timing errors of the pilot tone itself [22,28–32].

The 2.401 GHz sidebands used on right MOSAs, on the other hand, are less stable with respect to the pilot tone. This is acceptable, as additional clock noise in this signal can be corrected for using the sideband beatnotes in the RFI [22,32].

Lastly, any errors in the 80 MHz phasemeter clock are also corrected by the pilot tone correction, such that it is not performance critical and could either be directly synthesized from the USO or from the 2.4 GHz signal. This choice is currently irrelevant for our simulation since we directly simulate the pilot tone as reference clock for all measurements.

#### 2. Clock-signal model

We model the pilot tone signal as a periodic signal of the form

<footnote n="4">⁴No final decision has been made on the precise USO frequency that will be used for LISA.</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 5. Overview over the USO frequency distribution on one optical bench, based on [27]. Not depicted is a possible electrical comparison between the 2.4 GHz and 2.401 GHz signals. Note that the phasemeter clock is not performance critical, and could be synthesized from either the 2.4 GHz signal or directly from the USO, indicated by two possible connections in dotted lines.</figcaption></figure>

$$V_{\text{PT}}(\tau) = \cos(2\pi\nu_{\text{PT}}[\tau + q_1(\tau)]). \quad (50)$$

Here, $q_1(\tau)$ describes the timing deviations of the pilot tone generated on spacecraft 1 with respect to the TPS $\tau_1$, expressed in the latter. Note that the dominant noise source in the pilot tone generation is the USO itself [27], such that we assume the statistical properties of the pilot tone noise to be identical to those of the USO noise.

We further decompose $q_1(\tau)$ using two time series,

$$q_1(\tau) = q_1^o(\tau) + q_1^e(\tau), \quad (51)$$

to model large deterministic effects (such as clock frequency offsets and drifts) and small in-band stochastic fluctuations. As before, we do not simulate the 75 MHz signal itself, but only $q_1^o(\tau)$ and $q_1^e(\tau)$ [or rather $\dot{q}_1^o(\tau)$ and $\dot{q}_1^e(\tau)$ as the pilot tones fractional frequency fluctuations].

The clock signal is used to create the sidebands, as described in Sec. III C. The total phase of the sideband modulation signals is modeled as

$$\nu_{12}^m \cdot (\tau + q_1(\tau) + M_{12}(\tau)). \quad (52)$$

Here, $\nu_{12}^m$ is the constant nominal frequency<fnref n="5" /> of the modulating signal on optical bench 12. Imperfections in the frequency conversion between the pilot tone and the sidebands are modeled by an additional modulation noise term $M_{12}(\tau)$.

### D. Timer model

In order to model timestamping and pseudoranging (cf. Sec. VII), we not only need the frequency fluctuations of the local clock, but also the time shown by each

spacecraft timer. These times must be tracked down to at least ns precision while reaching values of around $10^8$ s at the end of the 10 years of extended mission. The use of double-precision floating-point numbers is not compatible with such a dynamic range. Therefore, we simulate offsets of that timer relative to the associated TPS $\delta\hat{\tau}_1^\tau(\tau) \equiv \delta\hat{\tau}_1(\tau)$, called timer deviations, which evolve slowly with time. The total clock time<fnref n="6" /> $\hat{\tau}_1^\tau(\tau)$ as a function of the TPS can then be computed by

$$\hat{\tau}_1^\tau(\tau) = \tau + \delta\hat{\tau}_1(\tau). \quad (53)$$

Timer deviations are closely related to the clock timing jitter,

$$\delta\hat{\tau}_1(\tau) = q_1(\tau) + \delta\hat{\tau}_{1,0}. \quad (54)$$

In this equation, $\delta\hat{\tau}_{1,0}$ accounts for the fact that we don't know the true time $\tau_{1,0}$ at which we turn on the timer, i.e., we cannot directly relate the initial phase of the clock signal $q_1(\tau_{1,0})$ to any external time frame.

### E. Signal sampling

#### 1. Signal sampling in terms of phase

The photoreceiver signals recorded, say, on spacecraft 1, are generated according to the TPS $\tau_1$. The measurements that are eventually telemetered, however, are recorded and timestamped with clock time $\hat{\tau}_1$. As a consequence, we need to resample the photoreceiver signals from the TPS to the clock time frame.

If a photoreceiver signal $\Phi_{\text{PD}}$ is expressed in terms of phase, this can be achieved following Eq. (3),

$$\Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) = \Phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)). \quad (55)$$

Therefore, we need to compute the TPS $\tau_1^{\hat{\tau}_1}(\tau)$ given a given clock time $\tau$. This quantity can be computed by writing Eq. (53) evaluated at $\tau_1^{\hat{\tau}_1}(\tau)$,

$$\hat{\tau}_1^\tau(\tau_1^{\hat{\tau}_1}(\tau)) = \tau_1^{\hat{\tau}_1}(\tau) + \delta\hat{\tau}_1(\tau_1^{\hat{\tau}_1}(\tau)). \quad (56)$$

We use Eq. (3) to rewrite the left-hand side, which gives, after rearranging,

$$\tau_1^{\hat{\tau}_1}(\tau) = \tau - \delta\hat{\tau}_1(\tau_1^{\hat{\tau}_1}(\tau)). \quad (57)$$

We can solve this implicit equation for $\tau_1^{\hat{\tau}_1}(\tau)$ iteratively, by computing

<footnote n="5">⁵By definition, these frequencies are at their nominal values. The real modulation signals will have a frequency offset due to the terms $q_1$ and $M_{12}$ in Eq. (15).</footnote>

<footnote n="6">⁶This timescale will be realized in practice by the so-called spacecraft elapsed time (SCET), which is the only timescale directly available onboard the satellites.</footnote>
$$ \delta\hat{\tau}_1^{(0)}(\tau) = \delta\hat{\tau}_1(\tau), \qquad (58a) $$

$$ \delta\hat{\tau}_1^{(n+1)}(\tau) = \delta\hat{\tau}_1(\tau - \delta\hat{\tau}_1^{(n)}(\tau)), \qquad (58b) $$

such that

$$ \lim_{n \to \infty} \delta\hat{\tau}_1^{(n)}(\tau) = \delta\hat{\tau}_1(\tau_1^{(n)}(\tau)). \qquad (59) $$

Since the timer deviations are evolving slowly, the iteration converges quickly. In our simulations, we stop after two iterations, such that

$$ \tau_1^{(n)} \approx \tau - \delta\hat{\tau}_1^{(n)}(\tau). \qquad (60) $$

We can then plug the previous equation in Eq. (55) to write all frame-independent measurements as a functions of the correct recording times $\Phi_{\text{PD}}^{\hat{\tau}_1}(\tau)$, given the same quantities expressed in the TPS. We find

$$ \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) \approx \Phi_{\text{PD}}^{\tau_1}(\tau - \delta\hat{\tau}_1^{(2)}(\tau)). \qquad (61) $$

This operation can be implemented with a time-varying fractional-delay filter (interpolation).

We introduce the timestamping operator $\mathbf{T}_1$, which shifts a signal $s(\tau)$ from the TPS to the clock time of spacecraft 1. Formally, its action is given by

$$ \mathbf{T}_1 s(\tau) = s(\tau - \delta\hat{\tau}_1^{(2)}(\tau)). \qquad (62) $$

Using this shorthand notation, Eq. (61) now reads

$$ \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) = \Phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)) \approx \mathbf{T}_1 \Phi_{\text{PD}}^{\tau_1}(\tau). \qquad (63) $$

Note that this is only valid for measurements expressed in phase, as frequencies are not frame-independent quantities.

## 2. Sampling errors in terms of frequency

The effect of sampling can also be expressed in terms of total frequency, where it manifests itself as a Doppler-like frequency shift.

In the following paragraph, we compute frequencies by taking the derivative of phase with respect to the clock time, since this is the time reference that the phasemeter will use to measure the signal frequency. From Eq. (55), and denoting function composition as $(\Phi_{\text{PD}}^{\hat{\tau}_1} \circ \tau_1^{\hat{\tau}_1})(\tau) = \Phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau))$, we have

$$ \nu_{\text{PD}}^{\hat{\tau}_1}(\tau) = \frac{d\Phi_{\text{PD}}^{\hat{\tau}_1}}{d\tau}(\tau) = \frac{d(\Phi_{\text{PD}}^{\tau_1} \circ \tau_1^{\hat{\tau}_1})}{d\tau}(\tau). \qquad (64) $$

Using the chain rule,

$$ \nu_{\text{PD}}^{\hat{\tau}_1}(\tau) = \nu_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)) \times \frac{d\tau_1^{\hat{\tau}_1}}{d\tau}(\tau). \qquad (65) $$

To compute the derivative of $\tau_1^{\hat{\tau}_1}(\tau)$, let us differentiate the defining implicit Eq. (57),

$$ \begin{aligned} \frac{d\tau_1^{\hat{\tau}_1}}{d\tau}(\tau) &= 1 - \frac{d(\delta\hat{\tau}_1 \circ \tau_1^{\hat{\tau}_1})}{d\tau}(\tau) \\ &= 1 - \frac{d\delta\hat{\tau}_1}{d\tau}(\tau_1^{\hat{\tau}_1}(\tau)) \times \frac{d\tau_1^{\hat{\tau}_1}}{d\tau}(\tau). \end{aligned} \qquad (66) $$

Using Eq. (54), we find $d\delta\hat{\tau}_1/d\tau = \dot{q}_1(\tau)$. Inserting this identity, we can rearrange the previous equation to get

$$ \frac{d\tau_1^{\hat{\tau}_1}}{d\tau}(\tau) = \frac{1}{1 + \dot{q}_1(\tau_1^{\hat{\tau}_1}(\tau))}, \qquad (67) $$

which finally yields for the total frequency,

$$ \nu_{\text{PD}}^{\hat{\tau}_1}(\tau) = \frac{\nu_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau))}{1 + \dot{q}_1(\tau_1^{\hat{\tau}_1}(\tau))} \approx \mathbf{T}_1 \left[ \frac{\nu_{\text{PD}}^{\tau_1}(\tau)}{1 + \dot{q}_1(\tau)} \right]. \qquad (68) $$

## 3. Sampling in two-variable decomposition

We now want to describe the effect of timing errors in the framework of two-variable decomposition. This will allow us to split the sampling errors derived previously into large deterministic offsets in the measurement timestamps, and small stochastic fluctuations that enter as an additional noise term. The latter represent what is often referred to as *clock noise*, see, e.g., [20].

However, we want to make it clear once more that this decomposition is entirely artificial. Both slow drifts and in-band clock noise describe the same physical process, namely the instability of the USO, on different time scales.

The sampling process applies to the total phase of each photoreceiver signal, given by Eq. (37) as

$$ \begin{aligned} \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) &= \Phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)) \\ &= \phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)) + \phi_{\text{PD}}^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)). \end{aligned} \qquad (69) $$

Since $\phi_{\text{PD}}^{\tau_1}(\tau)$ is very quickly evolving, small (first-order) timing fluctuations in $\tau_1^{\hat{\tau}_1}(\tau)$ must appear in the measurement described by $\phi_{\text{PD}}^{\tau_1}(\tau)$. Thus, we must account for the cross coupling between $\phi_{\text{PD}}^{\tau_1}(\tau)$ and $\phi_{\text{PD}}^{\tau_1}(\tau)$, and we cannot simply time shift both components individually.

We can insert Eqs. (51) and (54) into Eq. (57) to get

$$ \tau_1^{\hat{\tau}_1}(\tau) = \tau - \delta\hat{\tau}_{1,0} - q_1^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)) - q_1^{\tau_1}(\tau_1^{\hat{\tau}_1}(\tau)). \qquad (70) $$

We model clock-noise fluctuations $\dot{q}_1^\epsilon$ as band-limited noise, such that they remain small and we can expand the $\phi^\epsilon$ term in Eq. (69) to first order in $\dot{q}_1^\epsilon$,
$$ \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) = \phi_{\text{PD}}^{o}(\tau - \delta\hat{\tau}_{1,0} - q_1^o(\tau^{\hat{\tau}_1}(\tau))) + \phi_{\text{PD}}^{e}(\tau^{\hat{\tau}_1}(\tau)) \\ - \nu_{\text{PD}}^{o}(\tau - \delta\hat{\tau}_{1,0} - q_1^o(\tau^{\hat{\tau}_1}(\tau)))q_1^e(\tau^{\hat{\tau}_1}(\tau)). \quad (71) $$

Finally, we obtain the two variable-decomposition for the resampled photoreceiver phase.

$$ \phi_{\text{PD}}^{\hat{\tau}_{1,o}}(\tau) \approx \phi_{\text{PD}}^{o}(\tau - \delta\hat{\tau}_{1,o} - \mathbf{T}_1 q_1^o(\tau)), \quad (72a) $$

$$ \phi_{\text{PD}}^{\hat{\tau}_{1,e}}(\tau) \approx \mathbf{T}_1 \phi_{\text{PD}}^{e}(\tau) - \nu_{\text{PD}}^{o}(\tau - \delta\hat{\tau}_{1,o} - \mathbf{T}_1 q_1^o(\tau)) \mathbf{T}_1 q_1^e(\tau). \quad (72b) $$

For frequency data, we start with Eq. (68), and decompose again clock noise $\dot{q}_1$ into two variables, as explained in Sec. IV. We then expand it to first order in $\dot{q}_1^e$ to get

$$ \frac{1}{1 + \dot{q}_1^o(\tau^{\hat{\tau}_1}(\tau)) + \dot{q}_1^e(\tau^{\hat{\tau}_1}(\tau))} \\ \approx \frac{1}{1 + \dot{q}_1^o(\tau^{\hat{\tau}_1}(\tau))} - \frac{\dot{q}_1^e(\tau^{\hat{\tau}_1}(\tau))}{[1 + \dot{q}_1^o(\tau^{\hat{\tau}_1}(\tau))]^2}. \quad (73) $$

So in total, we have

$$ \nu_{\text{PD}}^{\hat{\tau}_1}(\tau) \approx \nu_{\text{PD}}^{\tau_1}(\tau^{\hat{\tau}_1}(\tau)) \left[ \frac{1}{1 + \dot{q}_1^o(\tau^{\hat{\tau}_1}(\tau))} - \frac{\dot{q}_1^e(\tau^{\hat{\tau}_1}(\tau))}{[1 + \dot{q}_1^o(\tau^{\hat{\tau}_1}(\tau))]^2} \right]. \quad (74) $$

We now expand $\nu_{\text{PD}}^{\tau_1}(\tau) = \nu_{\text{PD}}^{\tau_{1,o}}(\tau) + \nu_{\text{PD}}^{\tau_{1,e}}(\tau)$, and neglect the small coupling of $q_1^e(\tau)$ to the already small fluctuations $\nu_{\text{PD}}^{\tau_{1,e}}(\tau)$. We collect the terms to express the photodiode signal offsets $\nu_{\text{PD}}^{\hat{\tau}_{1,o}}(\tau)$ and fluctuations $\nu_{\text{PD}}^{\hat{\tau}_{1,e}}(\tau)$ after shifting to the clock time, using Eq. (61),

$$ \nu_{\text{PD}}^{\hat{\tau}_{1,o}}(\tau) \approx \frac{\mathbf{T}_1 \nu_{\text{PD}}^{\tau_{1,o}}(\tau)}{1 + \mathbf{T}_1 \dot{q}_1^o(\tau)}, \quad (75a) $$

$$ \nu_{\text{PD}}^{\hat{\tau}_{1,e}}(\tau) \approx \frac{\mathbf{T}_1 \nu_{\text{PD}}^{\tau_{1,o}}(\tau)}{1 + \mathbf{T}_1 \dot{q}_1^o(\tau)} - \frac{\mathbf{T}_1 \nu_{\text{PD}}^{\tau_{1,o}}(\tau) \mathbf{T}_1 \dot{q}_1^e(\tau)}{[1 + \mathbf{T}_1 \dot{q}_1^o(\tau)]^2}. \quad (75b) $$

To simplify further our equations, we define the frequency timestamping operator, which includes the rescaling by $1 + \dot{q}_1^o$. It is formally defined by its action on a signal $s(\tau)$,

$$ \dot{\mathbf{T}}_1 s(\tau) = \mathbf{T}_1 \left[ \frac{s(\tau)}{1 + \dot{q}_1^o(\tau)} \right] = \frac{\mathbf{T}_1 s(\tau)}{1 + \mathbf{T}_1 \dot{q}_1^o(\tau)}. \quad (76) $$

Now, photoreceiver frequency signals in the clock-time frame of spacecraft 1 read

$$ \nu_{\text{PD}}^{\hat{\tau}_{1,o}}(\tau) \approx \dot{\mathbf{T}}_1 \nu_{\text{PD}}^{\tau_{1,o}}(\tau), \quad (77a) $$

$$ \nu_{\text{PD}}^{\hat{\tau}_{1,e}}(\tau) \approx \dot{\mathbf{T}}_1 \left[ \nu_{\text{PD}}^{\tau_{1,e}}(\tau) - \frac{\nu_{\text{PD}}^{\tau_{1,o}}(\tau) \dot{q}_1^e(\tau)}{1 + \dot{q}_1^o(\tau)} \right]. \quad (77b) $$

## V. ONBOARD PROCESSING

In this section, we describe the processing steps the readout signals undergo on board the spacecraft, and in particular the filtering and downsampling steps. The sampling rates used in our simulation are shown schematically in Fig. 6. We then give the expression of the main measurement signals, which are the main outputs of the simulation.

### A. Filtering and downsampling

#### 1. Physics sampling rate

As described in Sec. IV and following the current mission design, the onboard phasemeters track the phase (or, equivalently, the instantaneous frequency) of sampled and digitized versions of the MHz beatnotes using DPLLs running at 80 MHz.

For performance reasons, we cannot simulate continuous analog signals nor DPLL signals at their real sampling rate. Instead, we use a discretized representation and rely on high-level models to capture the most significant effects. In our simulations, continuous quantities, as well as photoreceiver signals and beatnote measurements, are simulated at the physics rate

$$ f_s^{\text{phy}} = 16 \text{ Hz.} \quad (78) $$

Note that this physics rate matches the penultimate downsampling step of the real onboard decimation chain (described below), which is used by the DFACS.

<figure><img src="image_9.png" /><figcaption>FIG. 6. Overview of the real signal sampling rates, from continuous optical and electrical signals, to the cascaded downsampled signals down to the telemetered 4 Hz data. We also indicate the sampling rates used in the simulation (in gray boxes) to represent these signals; continuous and high-frequency signals are represented by discrete 16 Hz simulated quantities, while the telemetry data is simulated at their true 4 Hz rate.</figcaption></figure>
## 2. Antiaiasing filters

The current mission design suggests that the raw 80 MHz phasemeter beatnote signals are then filtered and down-sampled to various lower sampling rates, and ultimately to the final measurement rate of 4 Hz. This last measurement sampling rate is in line with the mission instrument design, and compatible with the limited telemetry budget and the required bandwidth for on-ground processing. The 4 Hz data are then telemetered down to Earth.

High-order digital low-pass finite impulse response (FIR) filters, as well as cascading filters, are expected to be used to prevent noise aliasing in the frequency band relevant for LISA data analysis, between $10^{-4}$ Hz and 1 Hz [33]. These filters must strongly attenuate the signals above the Nyquist frequency, while maintaining a gain close to unity and low phase distortion below 1 Hz. Their precise implementation is still under development.

In the simulation, we use a single digital symmetrical FIR filter to go from $f_s^{phy}$ to the final measurement sampling rate of

$$f_s^{\text{meas}} = 4 \text{ Hz.} \qquad (79)$$

The default implementation of the antiaiasing filter is described in Appendix D.

## 3. Decimation

Once the beatnote frequency measurements are filtered, we use a four-fold decimation (we select one sample out of four) to produce the final 4 Hz telemetry data. They are the main output of the simulation.

Analytically, we model the filtering and downsampling step with the continuous, linear filter operator F, which is applied to the beatnote frequency measurements.

## B. Telemetered beatnote measurements

We summarize here the downsampled, filtered beatnote measurements output by the phasemeter, i.e., the interspacecraft, test-mass, reference carrier and sideband beatnote frequencies. They are ultimately telemetered down to Earth.<fnref n="7" />

### 1. Beatnote measurement notation

For these beatnote measurements, we introduce a clear notation that uses the name of the associated interferometer and its index, complemented by the type of beam (carrier or sideband). The real phasemeter will only produce the total frequency or the total phase of the signal. For our studies, however, it is often useful to also have access to the underlying offsets and fluctuations in two separate variables, which is why we give here the signals in this form.

The simulation will provide an additional output for the total frequency, given as the sum of the two components.

For readability's sake, we drop all time arguments. We use delay operators to account for time shifts that appear when propagating signals. We denote $\mathbf{D}_{12}$ the delay operator associated with the PPR $d_{12}^o(\tau)$ defined in Sec. III E, such that for any signal $s(\tau)$,

$$\mathbf{D}_{12}s(\tau) = s(\tau - d_{12}^o(\tau)). \qquad (80)$$

Furthermore, we introduce the Doppler-delay operator, which is defined as

$$\dot{\mathbf{D}}_{12}s(\tau) = (1 - \dot{d}_{12}^o(\tau))s(\tau - d_{12}^o(\tau)). \qquad (81)$$

We also make use of the timestamping operators $\dot{T}_i$ introduced in Sec. IV E, and the downsampling and filtering operator F.

We will also use a shorthand notation for the beatnote frequency offsets in the TPS, which we define by

$$a_{12}^c \equiv \nu_{\text{isi}_{12},c}^o = \dot{\mathbf{D}}_{12}O_{21} - \nu_0 \dot{d}_{12}^o - O_{12}, \qquad (82a)$$

$$a_{12}^{\text{sb}} \equiv \nu_{\text{isi}_{12},\text{sb}}^{o} = a_{12}^{c} + \dot{\mathbf{D}}_{12}[\nu_{21}^{m}(1 + \dot{q}_{2}^{o})] - \nu_{12}^{m}(1 + \dot{q}_{1}^{o}), \qquad (82b)$$

$$b_{12}^c \equiv \nu_{\text{rfi}_{12},c}^o = O_{13} - O_{12}, \qquad (82c)$$

$$b_{12}^{\text{sb}} \equiv \nu_{\text{rfi}_{12},\text{sb}}^{o} = b_{12}^{c} + (\nu_{13}^{m} - \nu_{12}^{m})(1 + \dot{q}_{1}^{o}). \qquad (82d)$$

In addition, most of the laser-related terms $p_{12}$, $O_{12}$ will be determined by the laser-locking scheme, as described in Sec. VI.

### 2. Interspacecraft interferometer beatnote frequencies

The carrier-carrier beatnote frequency measurement in the ISI's contains the delayed distant and local laser frequency fluctuations $\dot{p}_{21}$ and $\dot{p}_{12}$, as well as the delayed distant and local optical-bench path length noises appearing as Doppler shifts $\dot{N}_{\text{isi}_{12}\leftarrow 21}^{\text{ob}}$ and $\dot{N}_{\text{isi}_{12}\leftarrow 12}^{\text{ob}}$. The effect of the gravitational-wave signal $\dot{H}_{12}$ also appears as an extra Doppler shifts on the distant beam. Lastly, the readout $\dot{N}_{\text{isi}_{12},c}^{\text{ro}}$ and clock noise $a_{12}^c \dot{q}_1^e / (1 + q_1^o)$ terms are added.

$$\text{isi}_{12,c}^{o} = \mathbf{F}\dot{\mathbf{T}}_{1}a_{12}^{c}, \qquad (83a)$$

$$\text{isi}_{12,c}^{o} = \mathbf{F}\dot{\mathbf{T}}_{1} \left\{ \begin{aligned} & \dot{\mathbf{D}}_{12}\dot{p}_{21} - (\nu_{0} + \mathbf{D}_{12}O_{21})\dot{H}_{12} \\ & + \frac{\nu_{0}}{c}\dot{N}_{\text{isi}_{12}\leftarrow 21}^{\text{ob}} - (\dot{p}_{12} + \frac{\nu_{0}}{c}\dot{N}_{\text{isi}_{12}\leftarrow 12}^{\text{ob}}) \\ & + \dot{N}_{\text{isi}_{12},c}^{\text{ro}} - \frac{a_{12}^{c}\dot{q}_{1}^{e}}{1 + \dot{q}_{1}^{o}} \end{aligned} \right\}, \qquad (83b)$$

<footnote n="7">⁷As mentioned before, there are other data streams, such as the angular readouts provided by DWS, which we do not model here. The measured pseudorange (MPR) measurements are described in Sec. VII.</footnote>
$$ \text{isi}_{12,c} = \text{isi}^{o}_{12,c} + \text{isi}^{e}_{12,c}. \qquad (83c) $$

The sideband-sideband beatnote frequency measurement is similar with two main differences. First, the distant and local laser frequency fluctuations are affected by the coupling of the modulation frequency with clock jitter and modulation noise $\nu_{21}^m(\dot{q}_2^e + \dot{M}_{21})$ and $\nu_{12}^m(\dot{q}_1^e + \dot{M}_{12})$. Secondly, the distant sideband beatnote frequency offsets are shifted by the modulation frequency affected by out-of-band clock errors $\nu_{21}^m(1 + \dot{q}_2^o)$. Overall, we get

$$ \text{isi}^{o}_{12,sb} = \mathbf{F}\dot{\mathbf{T}}_1 a^{sb}_{12}, \qquad (84a) $$

$$ \text{isi}^{e}_{12,sb} = \mathbf{F}\dot{\mathbf{T}}_1 \left\{ \dot{\mathbf{D}}_{12}(\dot{p}_{21} + \nu_{21}^m(\dot{q}_2^e + \dot{M}_{21})) - (\nu_0 + \mathbf{D}_{12}[O_{21} + \nu_{21}^m(1 + \dot{q}_2^o)])\dot{H}_{12} + \frac{\nu_0}{c}\dot{N}_{\text{isi}_{12}\leftarrow 21}^{\text{ob}} \right. \\ \left. - \left(\dot{p}_{12} + \nu_{12}^m(\dot{q}_1^e + \dot{M}_{12}) + \frac{\nu_0}{c}\dot{N}_{\text{isi}_{12}\leftarrow 21}^{\text{ob}}\right) + \dot{N}_{\text{isi}_{12},\text{sb}}^{\text{ro}} - \frac{a_{12}^{\text{sb}}\dot{q}_1^e}{1 + \dot{q}_1^o} \right\}, \qquad (84b) $$

$$ \text{isi}_{12,sb} = \text{isi}^{o}_{12,sb} + \text{isi}^{e}_{12,sb}. \qquad (84c) $$

### 3. Reference interferometer beatnote frequencies

The carrier-carrier beatnote frequency measurement in the RFIs contains the frequency fluctuations of the adjacent and local laser beams $\dot{p}_{13}$ and $\dot{p}_{12}$, as well as the associated optical-bench path length noises $\dot{N}_{\text{rfi}_{12}\leftarrow 13}^{\text{ob}}$ and $\dot{N}_{\text{rfi}_{12}\leftarrow 12}^{\text{ob}}$. The adjacent beam that travels through the optical fiber picks up the backlink noise $\dot{N}_{12}^{\text{bl}}$. The readout $\dot{N}_{\text{rfi}_{12},c}^{\text{ro}}$ and clock noise $b_{12}^c \dot{q}_1^e / (1 + q_1^o)$ terms are then added,

$$ \text{rfi}^{o}_{12,c} = \mathbf{F}\dot{\mathbf{T}}_1 b^{c}_{12}, \qquad (85a) $$

$$ \text{rfi}^{e}_{12,c} = \mathbf{F}\dot{\mathbf{T}}_1 \left\{ \dot{p}_{13} + \frac{\nu_0}{c}(\dot{N}_{\text{rfi}_{12}\leftarrow 13}^{\text{ob}} + \dot{N}_{12}^{\text{bl}}) \right. \\ \left. - \left(\dot{p}_{12} + \frac{\nu_0}{c}\dot{N}_{\text{rfi}_{12}\leftarrow 12}^{\text{ob}}\right) + \dot{N}_{\text{rfi}_{12},c}^{\text{ro}} - \frac{b_{12}^{c} \dot{q}_1^e}{1 + \dot{q}_1^o} \right\}, \qquad (85b) $$

$$ \text{rfi}_{12,c} = \text{rfi}^{o}_{12,c} + \text{rfi}^{e}_{12,c}. \qquad (85c) $$

The expression for the sideband-sideband beatnote frequency measurement follows the same logic, with the adjacent and local laser frequency fluctuations affected by the in-band clock and modulation noises $\nu_{13}^m(\dot{q}_1 + \dot{M}_{13})$ and $\nu_{12}^m(\dot{q}_1 + \dot{M}_{12})$,

$$ \text{rfi}^{o}_{12,sb} = \mathbf{F}\dot{\mathbf{T}}_1 b^{sb}_{12}, \qquad (86a) $$

$$ \text{rfi}^{e}_{12,sb} = \mathbf{F}\dot{\mathbf{T}}_1 \left\{ \dot{p}_{13} + \nu_{13}^m(\dot{q}_1 + \dot{M}_{13}) + \frac{\nu_0}{c}(\dot{N}_{\text{rfi}_{12}\leftarrow 13}^{\text{ob}} + \dot{N}_{12}^{\text{bl}}) \right. \\ \left. - \left( \dot{p}_{12} + \nu_{12}^m(\dot{q}_1 + \dot{M}_{12}) + \frac{\nu_0}{c}\dot{N}_{\text{rfi}_{12}\leftarrow 12}^{\text{ob}} \right) \right. \\ \left. + \dot{N}_{\text{rfi}_{12},\text{sb}}^{\text{ro}} - \frac{b_{12}^{\text{sb}}\dot{q}_1^e}{1 + \dot{q}_1^o} \right\}, \qquad (86b) $$

$$ \text{rfi}^{o}_{12,sb} = \text{rfi}^{o}_{12,sb} + \text{rfi}^{e}_{12,sb}. \qquad (86c) $$

### 4. Test-mass interferometer beatnote frequencies

The carrier-carrier beatnote frequency measurements in the TMI have the same form as for the RFI, with the exception of the additional local test-mass noise term $\dot{N}_{12}^{\delta}$,

$$ \text{tmi}^{o}_{12,c} = \mathbf{F}\dot{\mathbf{T}}_1 b^{c}_{12}, \qquad (87a) $$

$$ \text{tmi}^{e}_{12,c} = \mathbf{F}\dot{\mathbf{T}}_1 \left\{ \dot{p}_{13} + \frac{\nu_0}{c}(\dot{N}_{\text{tmi}_{12}\leftarrow 13}^{\text{ob}} + \dot{N}_{12}^{\text{bl}}) \right. \\ \left. - \left( \dot{p}_{12} + \frac{\nu_0}{c}(\dot{N}_{\text{tmi}_{12}\leftarrow 12}^{\text{ob}} + 2\dot{N}_{12}^{\delta}) \right) \right. \\ \left. + \dot{N}_{\text{tmi}_{12},c}^{\text{ro}} - \frac{b_{12}^{c} \dot{q}_1^e}{1 + \dot{q}_1^o} \right\}, \qquad (87b) $$

$$ \text{tmi}_{12,c} = \text{tmi}^{o}_{12,c} + \text{tmi}^{e}_{12,c}. \qquad (87c) $$

As mentioned previously, we do not model sideband-sideband beatnote measurements in the TMIs.

## VI. LASER LOCKING AND FREQUENCY PLANNING

As mentioned in Sec. III C, each laser source is either frequency-locked to a resonant cavity or phase locked to another laser source using a specific interferometric beatnote. In this section, we describe how we simulate these laser-locking control loops. We then list the various locking configurations available for LISA in its baseline configuration.

### A. Frequency planning

The beatnote frequencies that can be measured by the LISA phasemeters are limited to between 5 MHz and 25 MHz.<fnref n="8" /> As a consequence, all beatnote frequencies need to be controlled to fall in this range, which is achieved by introducing predetermined offset frequencies in the laser locking control loops. A set of these frequency offsets for all lasers over the whole mission duration is called a frequency plan.

<footnote n="8">⁸The exact frequency range remains to be defined. In addition, some margins are required for both the upper and lower bounds to account for the sideband beatnotes, which are offset by 1 MHz from the carrier beatnotes.</footnote>
The problem of finding such frequency plans has recently been studied ([34]), and exact solutions have been found. We will use these solutions as an input to the simulation.

## B. Locking condition

Laser locking is achieved by controlling the frequency of a locked laser, such that a given beatnote frequency $\nu_{\text{PD}}(\tau)$ remains equal to a preprogrammed reference value $\nu_{\text{plan}}(\tau)$ provided in the frequency plan.

We do not simulate the actual control loop, but instead directly compute the correct frequency offsets and fluctuations of the locked laser for this locking condition to be satisfied. In reality, the locking control loops will have finite gain and bandwidth, such that the locking beatnotes can still contain out-of-band glitches and noise residuals. Here, we consider the frequency lock to be perfect. This means that the locking beatnote offset is exactly equal to the desired value.

Locking control loops run according to their local clocks, such that the locking condition is fulfilled in the local clock time frame $\hat{\tau}_1$. In addition, the frequency-plan locking frequencies are interpreted as functions of the same local clock time frame. In terms of total phase, the result of this control is that the measured beatnote phase is controlled to be exactly equal to the frequency-plan phase,

$$ \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau) = \Phi_{\text{plan}}(\tau). \quad (88) $$

Note that the control loop operates on data delivered by the phasemeter at a high frequency of 80 MHz ([35]). As such, we simulate the locking before applying any filtering or downsampling.

The previous locking condition is expressed in the local time frame, but we really want to solve for it in the TPS. We can use Eqs. (3), (51), and (53) to relate the measured beatnote phase to its equivalent in the TPS,

$$ \Phi_{\text{PD}}^{\tau_1}(\tau) = \Phi_{\text{PD}}^{\hat{\tau}_1}(\tau + q_1^o(\tau) + q_1^e(\tau) + \delta\hat{\tau}_{1,0}). \quad (89) $$

Using this result, we can write the locking condition from Eq. (88) as

$$ \Phi_{\text{PD}}^{\tau_1}(\tau) = \Phi_{\text{plan}}(\tau + q_1^o(\tau) + q_1^e(\tau) + \delta\hat{\tau}_{1,0}). \quad (90) $$

We expand the previous equation to first order in $q_1^e(\tau)$,

$$ \Phi_{\text{PD}}^{\tau_1}(\tau) = \Phi_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}) \\ + \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})q_1^e(\tau). \quad (91) $$

The second-order term is proportional to the product $q_1^e(\tau)^2 \dot{\nu}_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})$ of the square of the clock fluctuations and the time derivative of the frequency-plan locking frequency. To evaluate the order of magnitude of

this term, we compute the average clock time deviation [36] after a time corresponding to its saturation frequency (described in Appendix B); we find a value of the order of $10^{-9}$ s. In addition, all currently available frequency plans verify $\dot{\nu}_{\text{plan}}(\tau) < 3$ Hz s$^{-1}$. Therefore, we neglect terms of the order of $3 \times 10^{-18}$ cycles, far below the $\mu$-cycle level of gravitational-wave signals.

From Eq. (91), one directly obtains the usual decomposition in phase drifts and fluctuations,

$$ \Phi_{\text{PD}}^{\tau_{1,o}}(\tau) = \Phi_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}), \quad (92a) $$

$$ \Phi_{\text{PD}}^{\tau_{1,e}}(\tau) = \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})q_1^e(\tau). \quad (92b) $$

Indeed, the current baseline is to use piecewise linear functions with daily inflexions as frequency-plan locking frequencies. As a consequence, the latter are slowly varying, i.e., only consist in large out-of-band frequency offsets, such that $\nu_{\text{plan}}(\tau) = \nu_{\text{plan}}^o(\tau)$. This also applies to the preprogrammed reference phase $\Phi_{\text{plan}}(\tau) = \phi_{\text{plan}}^o(\tau)$.

We denote the (local) locked laser phase drifts and fluctuations as $\phi_l^{\tau_{1,o}}(\tau)$ and $\phi_l^{\tau_{1,e}}(\tau)$, and the (distant or adjacent) reference laser phase drifts and fluctuations as $\phi_r^{\tau_{1,o}}(\tau)$ and $\phi_r^{\tau_{1,e}}(\tau)$. Using Eq. (37), we have

$$ \Phi_{\text{PD}}^{\tau_{1,o}}(\tau) = \phi_r^{\tau_{1,o}}(\tau) - \phi_l^{\tau_{1,o}}(\tau), \quad (93a) $$

$$ \Phi_{\text{PD}}^{\tau_{1,e}}(\tau) = \phi_r^{\tau_{1,e}}(\tau) - \phi_l^{\tau_{1,e}}(\tau) + N_{\text{PD}}^{\text{to}}(\tau). \quad (93b) $$

It is now straightforward to write the resulting locked-laser phase drifts and fluctuations,

$$ \phi_l^{\tau_{1,o}}(\tau) = \phi_r^{\tau_{1,o}}(\tau) - \Phi_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}), \quad (94a) $$

$$ \phi_l^{\tau_{1,e}}(\tau) = \phi_r^{\tau_{1,e}}(\tau) - \Phi_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})q_1^e(\tau) \\ + N_{\text{PD}}^{\text{to}}(\tau). \quad (94b) $$

For frequency, we start by taking the derivative of Eq. (90) and expand it once again in $q_1^e$,

$$ \begin{aligned} \nu_{\text{PD}}^{\tau_1}(\tau) = & (1 + \dot{q}_1^o(\tau))\nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}) \\ & + \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})\dot{q}_1^e(\tau) \\ & + (1 + \dot{q}_1^o(\tau) + \dot{q}_1^e(\tau)) \\ & \times \dot{\nu}_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})q_1^e(\tau). \end{aligned} \quad (95) $$

The last term of the previous equation has a similar form as the small correction neglected in Eq. (30b). We study such a term in Appendix C and find that it is several orders of magnitude below the main noise term. Therefore, we will neglect it in the rest of this derivation.

Using the same two-variable decomposition along with the frequency equivalent of Eq. (93),
$$ \nu_{\text{PD}}^{\tau_{1,o}}(\tau) = \nu_{r}^{\tau_{1,o}}(\tau) - \nu_{l}^{\tau_{1,o}}(\tau), \qquad (96a) $$

$$ \nu_{\text{PD}}^{\tau_{1,e}}(\tau) = \nu_{r}^{\tau_{1,e}}(\tau) - \nu_{l}^{\tau_{1,e}}(\tau) + \dot{N}_{\text{PD}}^{\text{ro}}(\tau), \qquad (96b) $$

we finally obtain the resulting locked-laser frequency offset and fluctuations,

$$ \nu_l^{\tau_{1,o}}(\tau) = \nu_r^{\tau_{1,o}}(\tau) - (1 + \dot{q}_1^o(\tau)) \times \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}), \qquad (97a) $$

$$ \nu_l^{\tau_{1,e}}(\tau) = \nu_r^{\tau_{1,e}}(\tau) - \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}) \times \dot{q}_1^e(\tau) + \dot{N}_{\text{PD}}^{\text{ro}}(\tau). \qquad (97b) $$

Note that these equations describe the locked laser at the photodiode. To properly simulate this effect, we need the locked lasers frequency at the laser source, which we denote here as $\bar{\nu}_l(\tau)$. In Sec. III C, we add to the local beam frequency fluctuations an optical path length noise term $\dot{N}_{\text{PD}}^{\text{ob}}(\tau)$ during its propagation from the laser source to the photodiode. As a consequence, we have

$$ \bar{\nu}_l^e(\tau) = \nu_r^e(\tau) - \nu_{\text{plan}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})\dot{q}_1^e(\tau) + \dot{N}_{\text{PD}}^{\text{ro}}(\tau) - \frac{\nu_0}{c}\dot{N}_{\text{PD}}^{\text{ob}}(\tau) \qquad (98) $$

for the local locked-laser fluctuations.

### C. Locking configurations

In total, 5 of the 6 lasers in the constellation are locked (directly or indirectly) to one primary laser. Each of the locked lasers is locked to either the adjacent laser, using the RFI, so that Eqs. (97a) and (98) read

$$ O_{12}(\tau) = \nu_{\text{rfi}_{12}}^{o}(-13)(\tau) - (1 + \dot{q}_1^o(\tau)) \times \nu_{\text{fplan,rfi}_{12}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}), \qquad (99a) $$

$$ \dot{p}_{12}(\tau) = \nu_{\text{rfi}_{12}}^{e}(-13)(\tau) - \nu_{\text{fplan,rfi}_{12}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})\dot{q}_1^e(\tau) + \dot{N}_{\text{rfi}_{12}}^{\text{ro}}(\tau) - \frac{\nu_0}{c}\dot{N}_{\text{rfi}_{12}}^{\text{ob}}(-12)(\tau), \qquad (99b) $$

or to the distant laser, using the ISI, such that we get

$$ O_{12}(\tau) = \nu_{\text{isi}_{12}}^{o}(-21)(\tau) - (1 + \dot{q}_1^o(\tau)) \times \nu_{\text{fplan,isi}_{12}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0}), \qquad (100a) $$

$$ \dot{p}_{12}(\tau) = \nu_{\text{isi}_{12}}^{e}(-21)(\tau) - \nu_{\text{fplan,isi}_{12}}(\tau + q_1^o(\tau) + \delta\hat{\tau}_{1,0})\dot{q}_1^e(\tau) + \dot{N}_{\text{isi}_{12}}^{\text{ro}}(\tau) - \frac{\nu_0}{c}\dot{N}_{\text{isi}_{12}}^{\text{ob}}(-12)(\tau). \qquad (100b) $$

These expressions can be substituted into the equations of Sec. V B to derive the telemetered beatnote measurements with locked lasers.

<figure><img src="image_10.png" /><figcaption>FIG. 7. Laser-locking configurations for laser 12 as primary laser.</figcaption></figure>
The LISA model presented here permits six distinct locking topologies. For each of them, we have the freedom to choose the primary laser, such that, in total, we have 36 possible locking configurations. We plot the six configurations with laser 12 as the primary laser in Fig. 7. The other 30 combinations can be deduced by applying permutations of the indices.

## VII. PSEUDORANGING

In addition to the GHz sideband modulation, each laser beam will also carry an additional modulation with a predetermined pseudorandom noise (PRN) code used for absolute ranging and timing synchronization. The basic measurement principle is to correlate the received PRN code in each ISI with a local copy generated on the receiving spacecraft. The result of this measurement is the measured pseudorange (MPR), which contains information on both the light travel time between the spacecraft and the clock desynchronization.

### A. Pseudoranging modulation

The PRN modulation is performed at a relatively high frequency of around 2 MHz, far outside our simulation bandwidth. We therefore do not model the actual phase modulation. This modulation also causes a small additional noise in our measurement band, at a level below 1 pm Hz⁻⁰.⁵ in units of displacement [37], which we do not model. In addition, we only model the PRN measurement in the ISI, and completely ignore the presence of the PRN codes in the other interferometers.

Instead, we model this measurement by directly propagating the time deviations of each spacecraft timer with respect to their TPSs, alongside the laser beams. The MPR is then computed as the difference between the received and local timer.

Similarly to the main interferometric measurements and as described in Sec. VA, pseudoranging simulation is performed at $f_s^{\text{phy}}$, while the MPRs are ultimately filtered and downsampled to a lower rate $f_s^{\text{meas}}$.

### B. Pseudoranging as a clock-time difference

We consider in the following paragraphs a beam received by optical bench 12 at the receiver TPS $\tau$, which was emitted from optical bench 21 at emitter TPS $\tau - d_{12}(\tau)$. Here, the PPR $d_{12}(\tau)$ contains the light time of flight, as well as the conversion between the two proper times.

Conceptually, the MPR measures the pseudorange, given as the difference between the time $\hat{\tau}_1^\tau(\tau)$ shown by the local clock of the receiving spacecraft at the event of reception of the beam, and the time $\hat{\tau}_2^\tau(\tau - d_{12}(\tau))$ shown by the local clock of the sending spacecraft at the event of emission of the beam. In reality, the MPR only measures the pseudorange up to the repetition period of the PRN code, which is around 1 ms. The full pseudorange is then recovered by

combining the MPR measurements with ground-based observations.

At the moment, we do not simulate this effect and assume that the MPR directly gives the pseudorange without ambiguity. In addition, we assume that the vacuum between the satellites is sufficiently good that we can neglect (or compensate for) any dispersion effects, such that the PRN code suffers exactly the same delay as the carrier and sidebands.

Thus, we can model the MPR as the difference

$$R_{12}(\tau) = \hat{\tau}_1(\tau) - \hat{\tau}_2(\tau - d_{12}(\tau)) + N_{12}^R(\tau), \quad (101)$$

where $N_{12}^R(\tau)$ is a ranging noise term modeling imperfections in the overall correlation scheme.

### C. Pseudoranging in terms of timer deviations

As explained in Sec. IV D, we do not simulate the total clock time $\hat{\tau}_1(\tau)$ for each spacecraft, but only deviations $\delta\hat{\tau}_1(\tau)$ from the associated TPS,

$$\hat{\tau}_1(\tau) = \tau + \delta\hat{\tau}_1(\tau) \quad \text{and} \quad \hat{\tau}_2(\tau) = \tau + \delta\hat{\tau}_2(\tau). \quad (102)$$

Inserting these definitions into Eq. (101) yields

$$R_{12}(\tau) = \delta\hat{\tau}_1(\tau) - [\delta\hat{\tau}_2(\tau - d_{12}(\tau)) - d_{12}(\tau)] + N_{12}^R(\tau). \quad (103)$$

Let us define the clock time of the sending spacecraft propagated to the photodiode of the distant interspacecraft interferometer as

$$\delta\hat{\tau}_{\text{isi}_{12}\leftarrow 2}(\tau) \approx \delta\hat{\tau}_3(\tau - d_{12}(\tau)) - d_{12}(\tau). \quad (104)$$

We can then express the MPR as the simple difference

$$R_{12}(\tau) \approx \delta\hat{\tau}_1(\tau) - \delta\hat{\tau}_{\text{isi}_{12}\leftarrow 2}(\tau) + N_{12}^R(\tau). \quad (105)$$

In our simulation, we make the additional assumption that $d_{12} \approx d_{12}^o$ for this measurement. This is valid since the terms contained in $d_{12}^e$ (in our simulation model, only $H_{12}$) create timing jitters much less than a nanosecond.

Notice that in Eq. (105), we compute the MPR as a function of the receiving TPSs, so that formally $R_{12} = R_{12}^\tau$. In reality, the MPR is measured according to the clock time of the receiving spacecraft, $R_{12}^\tau$. Similarly to all other measurements, we simulate this by first generating $R_{12}^\tau$ and then resampling the resulting time series to obtain $R_{12}^\tau$, as described in Sec. IV E.

## VIII. IMPLEMENTATION

The model presented in the previous sections has been implemented independently in two LISA Consortium simulators, namely LISA Instrument and LISANode.
<figure><img src="image_11.png" /><figcaption>FIG. 8. LISA Instrument simulation workflow. At each stage, entire time series are generated. In the end, the simulation products comprise the telemetered beatnotes and the MPRs.</figcaption></figure>

In this section, we briefly describe the structure of both simulators and highlight the key differences between them. Results obtained from these simulators are presented in Sec. IX.

A. LISA Instrument

LISA Instrument [38] is a Python-based implementation of the simulation model described in this paper. It is designed to facilitate fast exploratory studies, run quick or partial

simulations (instrumental effects and noises can easily be toggled on and off), and prototype new features.

LISA Instrument ships as a standalone Python package. As a consequence, it is easy to install, use, and integrate in traditional workflows, such as Jupyter Notebooks. LISA Instrument does not require a custom installation and can be used out-of-the-box on most computing clusters.

LISA Instrument relies strongly on traditional numerical libraries, such as NumPy and SciPy [39,40], and therefore benefits from fast optimized vectorized operations as it handles large arrays of data. Its runtime performance is studied and compared to that of LISA Node in Fig. 9.

LISA Instrument runs stage-by-stage simulations, where time series are generated for the entire simulation duration at each stage. The main stages of a simulation are represented in Fig. 8. First, time series are generated for all noises enabled in the simulation, following the prescription of Appendix B. LISA Instrument uses FIR and cascaded RC filters [41] to generate the noise time series. Then, local beam frequencies are computed (see Sec. III C). Local beams from locked-laser sources are obtained by substituting the results of locking condition equations found in Sec. VI C. These local beams are then propagated to obtain the adjacent and distant beam frequency time series (see Secs. III D and III E). Optical beatnotes and measured beatnote frequencies are obtained from the equations derived in Secs. III F and IV. MPRs are also computed according to the model described in Sec. VII. At this point, both beatnote frequencies and

<figure><img src="image_12.png" /><figcaption>FIG. 9. Runtime and memory performance of LISA Instrument (in blue) and LISA Node (in yellow, orange, and red for the various compiler optimization levels), for three instrumental configurations of increasing complexity (from left to right) and simulations durations. Dashed lines are used for extrapolated data.</figcaption></figure>
MPRs are expressed as functions of their respective TPSs; they are resampled at the next stage to their associated clock time frames, following the methodology given in Sec. IV E. Finally, all measurements are filtered and downsampled (cf. Sec. VA) to obtain the telemetered beatnote frequencies and MPR measurements described in Sec. V B.

A downside of this simple implementation is that memory usage increases drastically with the simulation length. Memory pressure can become limiting for long simulations (typically more than a few months, see Sec. VIII C) if many noises and instrumental effects are enabled. The alternative implementation of the same simulation model, described in the next section, is overall less flexible but is optimized for long simulations.

## B. LISANode

LISANode [42,43] is a simulation framework that allows the user to build modular simulation graphs out of atomic computational units, called nodes. This is realized using a mix of Python and C++, where the Python code is responsible for defining the graph structure and interconnecting the different nodes. The nodes themselves are implemented in C++, such that the final executable is a C++ command line program.

Using C++ offers the advantage that compiler optimizations produce a fast executable, and allows us to reuse legacy code from previous C++ LISA simulators, such as LISACode [21]. Naturally, the cost is reduced readability, usability, and slower development times. Another consequence is that LISANode needs to be compiled on each machine it runs. To work around this last difficulty, we offer containerization solutions, in the form of Docker and Singularity images that contain optimized compiled versions of LISANode along with the software environment necessary to run them. These images can be downloaded and used on local machines and on most computing infrastructures, and do not hinder the runtime performance of the simulations.

In LISA Instrument, data is generated for the whole simulation length. On the contrary, LISANode creates data one step at a time: a sample at time $t_n$ is computed for all quantities before repeating the same instructions for the next samples $t_{n+1}$. New samples are therefore simulated on the fly, keeping only in memory the data that is required for the current and future samples. This way, memory usage remains roughly constant regardless of the simulation length (see Sec. VIII C). This allows long simulations to run on memory-constrained machines.

## C. Runtime and memory performance

We evaluate the runtime and memory performance of LISA Instrument and LISANode for three instrumental configurations of increasing complexity. In the *simple configuration*, all noises but laser noise are neglected.

Most instrument effects are disabled, as we use a static constellation with constant arm lengths, do not lock the lasers, assume perfect clocks, and set $f_s^{\text{meas}} = f_s^{\text{phy}} = 4$ Hz (no filtering or downsampling). In the *intermediate configuration*, all noises and effects are activated except for clock errors (no resampling of the measurement to clock time frames). We use realistic orbits and frequency plan for the locking configuration N1-12. We filter all measurements and use the nominal sampling frequencies described in Sec. VA. Lastly, the *full configuration* includes the effects of imperfect clocks.

We run simulations of increasing durations, ranging from 1 h ($1.4 \times 10^4$ telemetered samples at 4 Hz for each channel) to 6 months ($6.2 \times 10^7$ samples). Missing points indicates that the simulation did not complete on our test machine (MacBook Pro M1, 2021, 64 GB of RAM) because of excessive memory pressure or runtime. We extrapolate the results to 1 year (dashed lines) using a linear (LISANode) or quadratic (LISA Instrument) fit to the existing data points.

We used the latest version of LISA Instrument and LISANode, and compiled three different LISANode executables to study the impact of optimizations; one with no compiler optimization (compiler flag -O0), one with some optimizations (-O1), and one with most optimizations (-O2) enabled. In each case, we measure the runtime and peak memory usage. Results are reported in Fig. 9. Note that we do not include compilation time for LISANode in these figures, which strongly depends on the chosen optimization level (around 5 s for -O0, 1 min for -O1, and 10 min for -O2).

In terms of runtime, as expected, LISA Instrument is significantly (up to several orders of magnitude) faster for simple simulations. For intermediary configurations, LISA Instrument remains faster than LISANode up to simulation length of a month. Considering full instrumental configurations, highly-optimized versions of LISANode are faster than LISA Instrument irrespective of simulation duration. This is especially true for simulations longer than a month, where LISANode runs roughly twice as fast. In addition, memory usage can become limiting for simulations of a day or longer on typical machines with a few GB of memory when using LISA Instrument. LISANode caps memory peak usage to low values of about 100 MB irrespective of the simulation length.

## D. Simulation parameters and simulation products

For reference, we give in Table I the list of options accepted by LISA Instrument to configure the simulations. They parametrize the instrumental configuration (type of orbits, choice of a laser locking configuration, design of the onboard filters, etc.), the various noise models (noise amplitudes and spectral shapes), and the length of the simulation. Note that similar options can be used with LISANode, with some slight variations in their names.
<table><caption>TABLE I. Simulation parameters available to configure LISA Instrument, alongside with their units and reference sections.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>size</td><td>Number of samples to simulate</td><td>...</td><td>...</td></tr><tr><td>dt</td><td>Measurement sampling period f<sub>s</sub><sup>meas</sup></td><td>s</td><td>Section VA</td></tr><tr><td>t0</td><td>Initial simulation time</td><td>s</td><td>...</td></tr><tr><td>physics_upsampling</td><td>Ratio f<sub>s</sub><sup>phy</sup>/f<sub>s</sub><sup>meas</sup> of physics and measurement sampling frequencies</td><td>...</td><td>Section VA</td></tr><tr><td>clockinv_tolerance</td><td>Convergence criterion for clock noise inversion</td><td>s</td><td>Section IVE</td></tr><tr><td>clockinv_maxiter</td><td>Maximum number of iterations for clock noise inversion</td><td>...</td><td>Section IVE</td></tr><tr><td>aafilter</td><td>Antialiasing filter design specifications</td><td>...</td><td>Section VA</td></tr><tr><td>orbits</td><td>Path to orbit file, or PPRs d<sub>ij</sub><sup>o</sup>(τ)</td><td>s</td><td>Section III E</td></tr><tr><td>gws</td><td>Path to gravitational-wave file, or link responses H<sub>ij</sub>(τ)</td><td>...</td><td>Section III E</td></tr><tr><td>lock</td><td>Laser locking configuration (e.g., N1-12)</td><td>...</td><td>Section VI</td></tr><tr><td>fplan</td><td>Path to frequency-plan file, or locking beatnote frequencies ν<sub>PD,r</sub>(τ)</td><td>Hz</td><td>Section VI</td></tr><tr><td>central freq</td><td>Laser central frequency ν<sub>0</sub></td><td>Hz</td><td>Section III B</td></tr><tr><td>laser_asds</td><td>Laser noise amplitude spectral density (ASD)</td><td>Hz Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>modulation_asds</td><td>Modulation noise ASD</td><td>s Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>modulation_freqs</td><td>Modulation frequencies ν<sub>ij</sub><sup>m</sup></td><td>Hz</td><td>Section III C</td></tr><tr><td>clock_asds</td><td>Clock noise ASD</td><td>Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>clock offsets</td><td>Clock offsets δτ<sub>i,0</sub> from TPS τ<sub>i</sub></td><td>s</td><td>Section IV D</td></tr><tr><td>clock_freqoffsets</td><td>Clock frequency offsets y<sub>i,0</sub></td><td>s<sup>-1</sup></td><td>Appendix B</td></tr><tr><td>clock_freqlindrifts</td><td>Clock frequency linear drifts y<sub>i,1</sub></td><td>s<sup>-2</sup></td><td>Appendix B</td></tr><tr><td>clock_freqquaddrifts</td><td>Clock frequency quadratic drifts y<sub>i,2</sub></td><td>s<sup>-3</sup></td><td>Appendix B</td></tr><tr><td>backlink_asds</td><td>Backlink noise ASD</td><td>m Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>backlink_fknees</td><td>Backlink noise knee frequency</td><td>Hz</td><td>Appendix B</td></tr><tr><td>testmass_asds</td><td>Test-mass noise ASD</td><td>m s<sup>-2</sup> Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>testmass_fknees</td><td>Test-mass noise knee frequency</td><td>Hz</td><td>Appendix B</td></tr><tr><td>oms_asds</td><td>Readout noise ASD</td><td>m Hz<sup>-0.5</sup></td><td>Appendix B</td></tr><tr><td>oms_fknees</td><td>Readout noise knee frequency</td><td>Hz</td><td>Appendix B</td></tr><tr><td>ranging biases</td><td>Ranging systematic bias N<sub>i</sub><sup>R,o</sup></td><td>s</td><td>Appendix B</td></tr><tr><td>ranging_asds</td><td>Ranging noise ASD</td><td>s Hz<sup>-0.5</sup></td><td>Appendix B</td></tr></tbody></table>

<table><caption>TABLE II. Simulation products, alongside their units, and reference equations. All quantities are output at the measurement sampling rate f<sub>s</sub><sup>meas</sup>.</caption><thead><tr><th>Dataset</th><th>Description</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>isi_carrier_offset</td><td>ISI carrier beatnote frequency offsets</td><td>Hz</td><td>Equation (83a)</td></tr><tr><td>isi_carrier_fluctuations</td><td>ISI carrier beatnote frequency fluctuations</td><td>Hz</td><td>Equation (83b)</td></tr><tr><td>isi_carriers</td><td>ISI carrier beatnote total frequency</td><td>Hz</td><td>Equation (83c)</td></tr><tr><td>isi_usb_offset</td><td>ISI upper-sideband beatnote frequency offsets</td><td>Hz</td><td>Equation (84a)</td></tr><tr><td>isi_usb_fluctuations</td><td>ISI upper-sideband beatnote frequency fluctuations</td><td>Hz</td><td>Equation (84b)</td></tr><tr><td>isi_usbs</td><td>ISI upper-sideband beatnote total frequency</td><td>Hz</td><td>Equation (84c)</td></tr><tr><td>rfi_carrier_offset</td><td>RFI carrier beatnote frequency offsets</td><td>Hz</td><td>Equation (85a)</td></tr><tr><td>rfi_carrier_fluctuations</td><td>RFI carrier beatnote frequency fluctuations</td><td>Hz</td><td>Equation (85b)</td></tr><tr><td>rfi_carriers</td><td>RFI carrier beatnote total frequency</td><td>Hz</td><td>Equation (85c)</td></tr><tr><td>rfi_usb_offset</td><td>RFI upper-sideband beatnote frequency offsets</td><td>Hz</td><td>Equation (86a)</td></tr><tr><td>rfi_usb_fluctuations</td><td>RFI upper-sideband beatnote frequency fluctuations</td><td>Hz</td><td>Equation (86b)</td></tr><tr><td>rfi_usbs</td><td>RFI upper-sideband beatnote total frequency</td><td>Hz</td><td>Equation (86c)</td></tr><tr><td>tmi_carrier_offset</td><td>TMI carrier beatnote frequency offsets</td><td>Hz</td><td>Equation (85a)</td></tr><tr><td>tmi_carrier_fluctuations</td><td>TMI carrier beatnote frequency fluctuations</td><td>Hz</td><td>Equation (87b)</td></tr><tr><td>tmi_carriers</td><td>TMI carrier beatnote total frequency</td><td>Hz</td><td>Equation (87c)</td></tr><tr><td>mprs</td><td>MPRs</td><td>s</td><td>Equation (105)</td></tr></tbody></table>
We also list in Table II the quantities output by both simulators, alongside their units, and reference equations.

# IX. RESULTS AND DISCUSSION

An example code snippet for simple simulation and on-ground processing is given in Appendix E for the current versions of LISA Instrument and PyTDI. In the following sections, we describe the results of more complete simulations.

## A. Telemetry measurements

We present here the results of numerical simulations performed with LISA Instrument. The simulations include all noises described in the previous sections, in addition to a gravitational-wave signal from the loudest verification binary listed in [44]. Its orbital period is about 569.4 s and its 4-year signal-to-noise ratio (SNR) is estimated at 113. We simulated three days (about $10^6$ s) of measurements at the final rate of $f_{\text{meas}} = 4$ Hz. We have scaled the amplitude of the gravitational-wave signal such that its 3-day SNR matches the expected 4-year SNR. The light travel times between the spacecraft are computed from orbits files provided by ESA [24]. They are treated as time varying but remain roughly constant over the simulation duration, with values for each link between 8.17 s and 8.32 s. Lasers are locked in the N1-12 configuration (cf. Sec. VI) with a frequency plan computed accordingly ([45]).

Fig. 10 shows the time evolution of the six ISI beatnotes in terms of total frequency. In the chosen locking configuration, ISI 31 and 21 beatnotes (in green and brown, respectively) are locking beatnotes, and therefore are piecewise linear functions entirely determined by the frequency plan (cf. Sec. VI). They do not contain any noise since we assume a perfect laser phase-lock loop at the frequencies we study. Conversely, the remaining ISI beatnotes (in blue, orange, red, and purple) are nonlocking. Therefore, they have large MHz trends driven by the frequency plan and the relative motion of the spacecraft, in addition to a number of noises, dominated by the ~30 Hz Hz$^{-0.5}$ laser noise.

As expected, the frequency plan ensures that the beatnotes remain in the valid range of the phasemeter, i.e., between 5 MHz and 25 MHz in absolute value.

The amplitude spectral densities of all carrier beatnote frequency fluctuations are presented in Fig. 11. We used a Python implementation the log-scale power spectral density (LPSD) method [46] developed by C. Vorndamme with Kaiser windows. We overlay the 10 pm noise reference curve (in black), which is a typical target noise level for metrology noise in a single LISA link [16]. Its power spectral density (PSD) in units of frequency reads

$$ \left( \frac{2\pi f}{1064 \text{ nm}} \frac{10 \text{ pm}}{\sqrt{\text{Hz}}} \right)^2 [1 + (2 \text{ mHz } f)^4]. \quad (106) $$

<figure><img src="image_13.png" /><figcaption>FIG. 10. Time series of the ISI carrier beatnote frequencies. Locking ISI carrier beatnote frequencies (green 31 and brown 21) are piecewise linear functions driven by the frequency plan, free of any small fluctuations (c.f. left focus on 31). Nonlocking ISI carrier beatnotes (blue 12, orange 23, red 13, and purple 32) have large trends driven by the frequency plan and Doppler effect, and small in-band fluctuations (cf. right focus on 32).</figcaption></figure>

We expect test-mass acceleration noise to remain above this reference curve at low frequencies.

The nonlocking ISI beatnotes (blue lines) are dominated by laser noise (at about $30 \text{ HzHz}^{-0.5}$), only modulated by the one or two-way transfer function. Further processing is required to reduce this laser noise to below the noise requirements, which will reveal the presence on the injected gravitational-wave signal, cf. Sec. IX B.

Nonlocking RFI beatnotes (orange lines) contain mostly readout noises, and therefore remain below the 10 pm noise reference curve (refer to Appendix B for the noise models used in the simulation). Nonlocking TMI beatnotes (purple lines) contain, in addition to the same readout noises, test-mass acceleration noises, which become dominant below ~5 MHz. At these frequencies, the test-mass acceleration noise is clearly above the noise reference curve, as expected. At higher frequencies, we see that the different TMIs have different noise levels. On optical benches where the RFI is used for locking most common noises in the beams cancel, and the TMI is dominated by its own readout noise only. On the adjacent optical benches, on the other hand, we see an increased noise level due to the fiber backlink noise added to the locked laser during propagation between the benches.

ISI and RFI locking beatnotes are represented as plain red lines. Since we assume perfect laser phase-lock loops, these beatnotes should be vanishing, and we measure here the numerical noise floor of our simulations at about $10^{-12} \text{ HzHz}^{-0.5}$, well below the expected gravitational-wave signals at about $\nu_0 \times 10^{-21} \approx 10^{-7}$ Hz. Such a low-numerical noise floor can be achieved despite the large
<figure><img src="image_14.png" /><figcaption>FIG. 11. Amplitude spectral densities of all carrier beatnote frequency fluctuations. Nonlocking ISI beatnotes (blue lines) are dominated by laser noise, while nonlocking RFI and TMI beatnotes (orange and purple lines) contain mostly readout and test-mass acceleration noises. The 10 pm reference curve is represented as a bold black line. Locking beatnotes (plain red lines) should be vanishing, but represent here the numerical noise floor at about $10^{-12}$ HzHz$^{-0.5}$. Nonlocking beatnote fluctuations computed from the total beatnote frequencies (dashed red lines) have a 1000x larger numerical noise floor and contain other numerical artifacts.</figcaption></figure>

dynamic range of the quantities in play thanks to the two-variable decomposition described in Sec. III B 3 (the precision of beatnote frequency fluctuations are only limited by the magnitude of the laser noise).

We compare these results to what can be obtained with a single-variable model. The same nonlocking beatnote fluctuations have been computed by linearly detrending the total beatnote frequencies (to remove large out-of-band trends), and are plotted as dashed red lines. We see a numerical noise floor between $10^{-8}$ and $10^{-9}$ HzHz$^{-0.5}$, leaving little margin with respect to the expected magnitude of the gravitational signals and secondary noises that we wish to simulate and study.

Next, Fig. 12 shows time series of the six MPRs, as described in Sec. VII. In addition to the expected light travel times of about 8 s, we can observe that they also include the differential initial timer offsets (of a few seconds) and clock drifts (a few tens of milliseconds per day).

Finally, we show in Fig. 13 the PSDs of the different clock-noise-related measurements we simulate, all converted to units of s Hz$^{-0.5}$. The MPR (blue) is dominated by a white noise down to the lowest frequencies. Around 0.1 mHz to 0.3 mHz, the noise level coincides with that of the clock noise measured by the ISI sideband measurements (orange). Here, we plot a signal combination rejecting

common mode noise between carrier and sideband, following Eq. (B1) of Ref. [32], such that the plotted curve is dominated by the actual clock noise at most of the frequency band. Similarly, using Eq. (B9) of Ref. [32], we can combine the RFI sideband beatnotes to give a measurement of the larger of the two modulation noise terms, labeled $\Delta M_1$ (green). We see that the modulation noise is orders of magnitude smaller than the clock noise in most of the band.

<figure><img src="image_15.png" /><figcaption>FIG. 12. Time series of the MPRs. We can observe large deterministic errors with respect to the expected light travel time of ≈8 s due to initial timer offsets and clock drifts.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 13. PSDs of clock-noise related measurements. We plot the MPR (blue line) alongside measurements derived from the ISI sidebands (orange line) and the RFI sidebands (green line). See main text for details.</figcaption></figure>

## B. Processed measurements

We have seen that the raw telemetered beatnotes can be grouped in three categories. Locking beatnotes do not contain any gravitational-wave signal or noise (assuming perfect laser locking) and are dominated by numerical noises in our simulations. Nonlocking RFI beatnotes are also signal free and are dominated by secondary (readout and test-mass) noises. Only nonlocking ISI beatnotes carry useful gravitational-wave information, but contain laser noise at many orders of magnitude above the expected signals, alongside other noise sources.

In order to detect and analyze the gravitational-wave signals, we must therefore reduce these sources of noise to reasonable levels. This is achieved by a processing technique called TDI, in which multiple interferometric readouts are time shifted and combined to cancel the main noise sources.

To demonstrate that this kind of processing is possible using our simulated data, we apply the algorithm described in [32], using `PyTDI` [47], to reduce the limiting noise

<figure><img src="image_17.png" /><figcaption>FIG. 14. Amplitude spectral densities of second-generation TDI combinations $X_2$, $Y_2$, and $Z_2$. Note that laser noise, overwhelming in the telemetered beatnote frequencies, is suppressed by many orders of magnitude, such that the injected verification binary (SNR of 113) is now clearly visible around 3.5 mHz.</figcaption></figure>

sources included in our simulation (laser and clock noise) to below the required level. Fig. 14 shows the spectra of the second-generation TDI combinations $X_2$, $Y_2$, and $Z_2$, in which the laser and clock noises have been suppressed. The gravitational-wave signal is clearly visible at the expected frequency of 3.5 mHz, with an SNR of about 100.

Following the conclusions of [32], the biggest contributors to the residual noise are the ISI readout noise at high frequencies and the test-mass acceleration noise at low frequencies; other nonsuppressed noise sources have a smaller but non-negligible contribution.

## X. CONCLUSION

In this paper, we proposed a model of the LISA measurement chain, which includes the propagation of optical signals (modulated laser beams, each containing a carrier and an upper sideband) across the constellation and on the optical benches; the phase readout of the different interferometers; as well as the on-board processing of the beatnote signals. We also included a high-level model for the MPR auxiliary measurements, which are used to estimate the interspacecraft distances, necessary for the on-ground processing (such as TDI). This model accounts for laser-locking control loops, and properly treats different time frames and clock errors.

We presented two implementations of the model, along with a comparison of their runtime and memory performance that highlights their respective advantages and drawbacks. LISA Instrument is a Python implementation that is easy to use, and very efficient for short simulations (a few months or less), while we recommend `LISANode` for longer and more complex simulations.

Some results obtained with LISA Instrument are presented to demonstrate the correctness of the implementation. In particular, we check that the beatnote measurements exhibit the expected behavior, in terms of total frequency or frequency fluctuations. For the latter, we show that our model keeps numerical noise to acceptable levels for the study of instrumental noises and gravitational-wave signals. We also check that laser and clock noise reduction by TDI performs as expected by computing the second-generation Michelson combinations. We confirm that the residual noises in these channels matches their expected levels, and that a typical gravitational signal hidden in the raw telemetry data becomes clearly apparent at the expected frequency and magnitude.

The injection of gravitational-wave signals is possible through the multiple interfaces of LISA Instrument and `LISANode` with other simulation tools. These tools include in particular `LISA Orbits` [24] for realistic spacecraft orbits, `LISA Glitch` [48] for injection of instrumental artifacts, `LISA GW Response` [49] for injection of gravitational signals, and `PyTDI` [47] for further on-ground processing.

At the time this paper is written, some important instrumental effects are still under development. While some of
them are already implemented in the simulation, we do not include them in this description of the model but refer to future dedicated publications. In particular, we do not include any dynamical effects or dynamical control loops to actuate on the spacecraft and test masses [50]. We also do not include tilt-to-length effects, which occur as an apparent path length change due to any misalignment of optical elements and that will be partially mitigated on ground using DWS measurements [51]. We do not simulate any of the ground-based observations that will be used to determine the spacecraft positions and velocities, as well as the offsets of the onboard clocks with respect to a global timescale. These are required inputs to the further processing and data analysis steps, and therefore will be included in a future version of the simulation.

Finally, we currently produce phasemeter measurements expressed as total frequencies in Hz. There are ongoing discussions to choose the best representation for telemetry data; in order to capture any effects related to this choice, we plan to update the simulators to use the official data format once it is agreed upon. In particular, if phase data must be produced, the model should be updated to include an initial phase for each optical beam and clock signal, which do not affect the current frequency data.

## ACKNOWLEDGMENTS

The authors thank the LISA Simulation Working Group and the LISA Simulation Expert Group for the lively discussions on all simulation-related activities. They would like to personally thank A. Petiteau, G. Heinzel, V. Müller, A. Hees, M. Staab, and J. N. Reinhardt for their insightful feedbacks. J. B. B. gratefully acknowledges support from UK Space Agency (Grant No. ST/X002136/1), the Centre National d'Études Spatiales (CNES), the Central National pour la Recherche Scientifique (CNRS) and the Université Paris-Diderot. J. B. B. has been supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, California Institute of Technology, administered by Universities Space Research Association under contract with NASA. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). O. H. gratefully acknowledges support by Centre National d'Études Spatiales (CNES) and by the Deutsches Zentrum für Luft- und Raumfahrt (DLR, German Space Agency) with funding from the Federal Ministry for Economic Affairs and Energy based on a resolution of the German Bundestag (Projects Ref. No. 50OQ1601 and No. 50OQ1801). This work was supported by the Programme National GRAM of CNRS/INSU with INP and IN2P3 cofunded by CNES.

## APPENDIX A: SYMBOL GLOSSARY

This paper defines a large number of quantities and uses or introduces many symbols to represent them. To facilitate

the reading, we have listed the main conventions in Table III and the main quantities in Table IV.

## APPENDIX B: NOISE MODELS

We describe here the different noise sources that we include in our simulations. For each noise, we give a short description and its mathematical expression. That includes spectral shapes in the form of their PSD for stochastic terms, as well as any deterministic effects.

Noise models are derived from allocations or current best estimates (CBEs) given in the performance model ([52]), where applicable.<fnref n="9" />

We give here a continuous description of these noise models; however, they are actually implemented as discrete noise sources at $f_s^{\text{phy}} = 16$ Hz.

### 1. Laser noise

Laser noise describes the optical phase fluctuations in the electromagnetic field of a free-running laser stabilized to a cavity (fluctuations in the field amplitude are not included here). It is given in the performance model by the allocation for the laser frequency stability in units of frequency,

$$S_{\hat{N}^p}(f) = (30 \text{ HzHz}^{-0.5})^2 \left[ 1 + \left( \frac{2 \times 10^{-3} \text{ Hz}}{f} \right)^4 \right]. \quad (B1)$$

### 2. Modulation noise

Modulation noise describes any mismatch in the phase of the modulation sidebands (transmitted to the distance optical bench) and the pilot tone (used as a local timing reference) [27].

Both the 2.4 GHz and the 2.401 GHz sideband signals used for modulation are generated from the local USO. The pilot tone is derived from the electrical 2.4 GHz signal using a series of low-noise frequency dividers. The optical modulation is performed using an electro-optical modulator (EOM) followed by a fiber amplifier.

The fiber amplifier is the dominating part for the 2.4 GHz signal [27]. We fit a rough model to the blue curve in Fig. (5.13) of Ref. [27] to obtain a timing jitter power spectral density of

$$S_M(f) = (10^{-14} \text{ s Hz}^{-0.5})^2 \left[ 1 + \left( \frac{1.5 \times 10^{-2} \text{ Hz}}{f} \right)^2 \right]. \quad (B2)$$

Note that the blue curve used for the fit corresponds to a 1 W fiber amplifier, which allows for lower noise levels than the more recent measurements cited in [32] for a 2 W

<footnote n="9">⁹Some noise sources do not match the CBE given in latest performance model, and we plan to update them as soon as possible to reflect the currently expected LISA performance.</footnote>
<table><caption>TABLE III. Summary of the conventional notations used in this paper.</caption><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>□<sup>t</sup>(τ)</td><td>Quantity expressed in the barycentric coordinate time (TCB) (global time frame)</td></tr><tr><td>□<sup>τ<sub>i</sub></sup>(τ)</td><td>Quantity expressed in the spacecraft proper time (TPS) i (related to t by relativistic corrections)</td></tr><tr><td>□<sup>τ<sub>z</sub></sup>(τ)</td><td>Quantity expressed in the clock time on board spacecraft i (related to τ<sub>i</sub> by instrumental imperfections)</td></tr><tr><td>□ or d/dr</td><td>Time derivative of quantity (in the specified time frame)</td></tr><tr><td>□<sub>i</sub></td><td>Quantity related to spacecraft i</td></tr><tr><td>□<sub>ij</sub></td><td>Quantity related to MOSA ij</td></tr><tr><td>□<sub>A→B</sub></td><td>Quantity measured on A propagated from B</td></tr><tr><td>□<sub>PD</sub></td><td>Quantity related to a generic photodetector (no indices)</td></tr><tr><td>□<sub>c</sub></td><td>Quantity related to the main carrier</td></tr><tr><td>□<sub>m</sub> or □<sub>m</sub></td><td>Quantity related to the phase-modulation signal</td></tr><tr><td>□<sub>sb</sub> or □<sub>sb+</sub></td><td>Quantity related to the upper sideband</td></tr><tr><td>□<sub>sb<sup>-</sup></sub></td><td>Quantity related to the lower sideband</td></tr><tr><td>□<sup>o</sup></td><td>Large out-of-band component of a quantity (offsets or drifts)</td></tr><tr><td>□<sup>e</sup></td><td>Small in-band component of a quantity (fluctuations)</td></tr></tbody></table>

amplifier. This noise is therefore underestimated in the current version of the simulation.

For the 2.401 GHz signal, we expect a higher noise level due to the electrical conversion chain, which can no longer be realized by simple frequency dividers. Following [32], we model this by increasing the modulation noise in the right-hand side optical benches by a factor 10.

### 3. Test-mass acceleration noise

Test-mass acceleration noise describes the optical path length variations due to the test-mass motion with respect to its nominal position inside its housing.

It is given in the performance model by the allocation value for the single test-mass acceleration noise in acceleration units. We include an extra factor 2 in order to account for the beam reflection onto the test mass, and neglect the high-frequency component because it is smaller than the optical metrology system (OMS) displacement noise (see below). Moreover, we whiten the noise at below $10^{-4}$ Hz to prevent numerical overflow. We get

$$S_{N^{\delta}}(f) = (2 \times 2.4 \times 10^{-15} \text{ m s}^{-2} \text{ Hz}^{-0.5})^2 \\ \times \left[ 1 + \left( \frac{0.4 \times 10^{-3} \text{ Hz}}{f} \right)^2 \right]. \quad (B3)$$

Note that this is an out-of-loop value, ignoring the coupling of test mass to spacecraft motion introduced by DFACS.

### 4. Backlink noise

Beams are transmitted between adjacent optical benches using optical fibres. During this transmission, the beams

can pick up an additional phase noise term. We model only the nonreciprocal noise terms, i.e., the difference between the phase shift of a beam propagating from optical bench ij to ik vs that of the beam propagating from ik to ij.

Backlink noise is given in the performance model by the allocation for the reference backlink in displacement,

$$S_{N^{\text{bl}}}(f) = (3 \times 10^{-12} \text{ m Hz}^{-0.5})^2 \\ \times \left[ 1 + \left( \frac{2 \times 10^{-3} \text{ Hz}}{f} \right)^4 \right]. \quad (B4)$$

We use the same value for the TMIs and RFIs.

### 5. Readout noise

We summarize as readout noise the equivalent positional readout error due to technical noise sources, such as shot noise.

The OMS displacement noise is given in the performance model by the allocation value for the overall displacement long-arm, test-mass, and reference noise entries; in terms of displacement,<fnref n="10" />

$$S_{N^{\text{ro}}}(f) = A^2 \left[ 1 + \left( \frac{2 \times 10^{-3} \text{ Hz}}{f} \right)^4 \right], \quad (B5)$$

where $A = 6.35 \times 10^{-12} \text{ m Hz}^{-0.5}$ for the ISIs, $1.42 \times 10^{-12} \text{ m Hz}^{-0.5}$ for the TMIs, and $3.32 \times 10^{-12} \text{ m Hz}^{-0.5}$ for the RFIs.

The performance model does not give values for the sideband beatnotes. We approximate them using $\epsilon = 0.15$.

<footnote n="10">¹⁰The overall displacement noise in the performance model summarizes multiple noise sources, some of which are already accounted for independently in this model. The values for the OMS displacement noise are therefore overestimates.</footnote>
<table><caption>TABLE IV. Symbols for quantities used in this model alongside the units used in the simulation.</caption><thead><tr><th>Symbol</th><th>Unit</th><th>Description</th></tr></thead><tbody><tr><td>c</td><td>m s<sup>-1</sup></td><td>Speed of light in a vacuum</td></tr><tr><td>ν<sub>0</sub></td><td>Hz</td><td>Optical frequency of the lasers (281.6 THz)</td></tr><tr><td>Φ<sub>ij</sub></td><td>Cycles</td><td>Total phase of local beam at laser source ij</td></tr><tr><td>ν<sub>ij</sub></td><td>Hz</td><td>Total frequency of local beam at laser source ij</td></tr><tr><td>O<sub>ij</sub></td><td>Hz</td><td>Carrier frequency offset of local beam at laser source ij</td></tr><tr><td>p<sub>ij</sub></td><td>Cycles</td><td>Carrier phase fluctuations of local beam at laser source ij</td></tr><tr><td>Φ<sub>ij,m</sub></td><td>Cycles</td><td>Total phase of modulating signal on MOSA ij</td></tr><tr><td>ν<sub>ij,m</sub></td><td>Cycles</td><td>Total frequency of modulating signal on MOSA ij</td></tr><tr><td>ν<sup>nl</sup><sub>ij</sub></td><td>Hz</td><td>Nominal frequency of the modulating signal on MOSA ij</td></tr><tr><td>Φ<sub>ij-ji</sub></td><td>Cycles</td><td>Total phase of beam ij propagated to MOSA ji</td></tr><tr><td>ν<sub>ij-ji</sub></td><td>Hz</td><td>Total frequency of beam ij propagated to MOSA ji</td></tr><tr><td>Φ<sub>ifo<sub>ij</sub>-kl</sub></td><td>Cycles</td><td>Total phase of the beam kl propagated to photodetector ifo (isi, rfi, or tmi)</td></tr><tr><td>ν<sub>ifo<sub>ij</sub>-kl</sub></td><td>Hz</td><td>Total frequency of the beam kl propagated to photodetector ifo (isi, rfi, or tmi)</td></tr><tr><td>Φ<sub>ifo<sub>ij</sub></sub></td><td>Cycles</td><td>Optical beatnote total phase of ifo (isi, rfi, or tmi)</td></tr><tr><td>ν<sub>ifo<sub>ij</sub></sub></td><td>Hz</td><td>Optical beatnote frequency of ifo (isi, rfi, or tmi)</td></tr><tr><td>τ̃<sup>τ<sub>i</sub></sup></td><td>s</td><td>Instrumental clock time on board spacecraft i as a function of the TPS</td></tr><tr><td>δτ̃<sup>δ</sup></td><td>s</td><td>Deviations of instrumental clock time on board spacecraft i from TPS i</td></tr><tr><td>τ̃<sup>τ<sub>i</sub></sup></td><td>s</td><td>TPS of spacecraft i as a function of instrumental clock time on board spacecraft i</td></tr><tr><td>F</td><td>...</td><td>Filter operator, modeling the filtering and decimation stages from 16 Hz to 4 Hz</td></tr><tr><td>T<sub>i</sub></td><td>...</td><td>Phase timestamping operator, transforming a phase quantity from TPS i to instrumental clock time i</td></tr><tr><td>dot T<sub>i</sub></td><td>...</td><td>Frequency timestamping operator, transforming a frequency quantity from TPS i to instrumental clock time i</td></tr><tr><td>D<sub>ij</sub></td><td>...</td><td>Delay operator, time-shifting a phase quantity by d<sub>ij</sub></td></tr><tr><td>dot D<sub>ij</sub></td><td>...</td><td>Doppler-delay operator, time-shifting a frequency quantity by d<sub>ij</sub> (including Doppler corrections)</td></tr><tr><td>a<sub>ij</sub></td><td>Hz</td><td>Shorthand notation for the optical beatnote frequency offsets in the ISI ij</td></tr><tr><td>b<sub>ij</sub></td><td>Hz</td><td>Shorthand notation for the optical beatnote frequency offsets in the RFI and TMI ij</td></tr><tr><td>ifo<sub>ij</sub></td><td>Hz</td><td>Frequency readout of the ifo (isi, rfi, or tmi)</td></tr><tr><td>q<sub>i</sub></td><td>s</td><td>Noise of instrumental clock on board spacecraft i with respect to TPS</td></tr><tr><td>M<sub>ij</sub></td><td>s</td><td>Modulation noise on MOSA ij</td></tr><tr><td>N<sup>ob</sup><sub>ifo<sub>ij</sub>-kl</sub></td><td>s</td><td>Noise of beam kl propagated to ifo (isi, rfi, or tmi) due to optical path length variations on the optical bench</td></tr><tr><td>N<sup>ro</sup><sub>ifo<sub>ij</sub></sub></td><td>s</td><td>Readout noise in ifo (isi, rfi, or tmi)</td></tr><tr><td>N<sup>bl</sup><sub>ifo<sub>ij</sub></sub></td><td>s</td><td>Backlink noise in ifo (rfi or tmi)</td></tr><tr><td>N<sup>δ</sup><sub>ij</sub></td><td>s</td><td>Test-mass displacement noise</td></tr><tr><td>N<sup>R</sup><sub>ij</sub></td><td>s</td><td>Ranging noise</td></tr><tr><td>d<sub>ij</sub></td><td>s</td><td>Proper pseudo-range (PPR) as the difference of TPSs i (at reception) and j (at emission)</td></tr><tr><td>R<sub>ij</sub></td><td>s</td><td>Measured pseudo-range (MPR) ij as the difference of clock times i (at reception) and j (at emission)</td></tr><tr><td>H<sub>ij</sub></td><td>s</td><td>Integrated fluctuations of the PPR ij due to gravitational waves</td></tr></tbody></table>

instead of 0.85 in the shot noise formula to account for the lower power level. This yields $A = 1.25 \times 10^{-11} \text{ m Hz}^{-0.5}$ for the ISIs, and $7.90 \times 10^{-12} \text{ m Hz}^{-0.5}$ for the RFIs.

Optical bench path length noise in terms of displacement is given by the performance model as

$$S_{N^{\text{ob}}}(f) = A^2, \qquad (\text{B6})$$

## 6. Optical bench path-length noise

Optical bench path-length noise summarizes different optical path length noises due to, for example, jitters of optical components in the path of the different beams.

where $A = 4.24 \times 10^{-12} \text{ m Hz}^{-0.5}$ for local beams in TMIs and $2 \times 10^{-12} \text{ m Hz}^{-0.5}$ for local beams in RFIs.
## 7. Ranging noise

Pseudoranging is performed by correlating local and distant PRN signals, cf. Sec. VII. Ranging noise describes the imperfection of the overall ranging measurement scheme in a single link due to technical noise sources. Note that this does not include any noise appearing in the pseudorange itself, such as clock noise or changes in the optical path length between the spacecraft.

Pseudoranging is given by an ad hoc model, combining a systematic bias $N_i^{R,o}$ and a zero-mean stochastic Gaussian white noise $N_i^{R,\epsilon}(\tau_i)$,

$$N_i^R(\tau_i) = N_i^{R,o} + N_i^{R,\epsilon}(\tau_i), \quad (B7)$$

with default values of $S_{N^{R,\epsilon}} = 0.9 \text{ m Hz}^{-0.5}$ and $N_i^{R,o} = 0.8$.

## 8. Clock noise

USOs on each spacecraft act as central time references for all onboard systems. As described in Sec. IV, we actually use the pilot tone as the timing reference for all phasemeter measurements. Clock noise here models any deviations of these pilot tones from the corresponding spacecraft proper time (TPS).

Clock noise is given by the model described in [32], expressed in terms of fractional frequency deviations as the sum of a random jitter, and constant deterministic frequency offset, linear drift, and quadratic drift,

$$\dot{q}_i(\tau) = \dot{N}_i^q(\tau) + y_{0,i} + y_{1,i}\tau + y_{2,i}\tau^2, \quad (B8)$$

$\dot{N}_i^q(\tau)$ is a random jitter, generated as a flicker noise with a PSD between $10^{-5}$ Hz and $f_s^{\text{phy}}/2 = 8$ Hz given by

$$S_{\dot{N}_i^q}(f) = (6.32 \times 10^{-14})^2 f^{-1}, \quad (B9)$$

The deterministic coefficients are

$$y_0 \approx 5 \times 10^{-7} \text{ s s}^{-1}, \quad (B10)$$

$$y_1 \approx 1.6 \times 10^{-14} \text{ s s}^{-2}, \quad (B11)$$

$$y_2 \approx 9 \times 10^{-23} \text{ s s}^{-3}. \quad (B12)$$

These values should be seen as orders of magnitude, and will be different for all three USOs.

## APPENDIX C: MAGNITUDE OF $\dot{\nu}$

To estimate the order of magnitude of the term $\dot{\nu}_A^o H_{12}(\tau)$ we neglected compared to the term $\nu_A^o \dot{H}_{12}(\tau)$ we included, we can observe the rate of change in the example frequency plan presented in Fig. 15. This is plotted in Fig. 15. As we can see, we have $\dot{\nu}_A^o < 3 \text{ Hz s}^{-1}$ for the whole 4-year duration. On the other hand, $\nu_A^o$ is of the order of

<figure><img src="image_18.png" /><figcaption>FIG. 15. Magnitude of beatnote frequency derivatives for all nine beatnotes in the example frequency plan. Data provided by G. Heinzel.</figcaption></figure>

10 MHz. We consider both $\nu_A^o$ and $\dot{\nu}_A^o$ as constant scaling factors for this estimate.

Note that $H_{12}(\tau)$ and $\dot{H}_{12}(\tau)$ are noise terms that we can evaluate in the frequency domain. We have

$$\mathcal{F}[\dot{H}_{12}](f) = 2\pi f \times \mathcal{F}[H_{12}](f). \quad (C1)$$

The usual LISA measurement band extends from $10^{-4}$ Hz to 1 Hz, such that even at the lower limit of $10^{-4}$ Hz, we have

$$\begin{aligned} \nu_A^o \mathcal{F}[\dot{H}_{12}](f) &\approx 10^4 \text{ Hz s}^{-1} \times \mathcal{F}[H_{12}](f) \\ &\gg \dot{\nu}_A^o \mathcal{F}[H_{12}](f) \approx 10^4 \text{ Hz s}^{-1} \times \mathcal{F}[H_{12}](f). \end{aligned} \quad (C2)$$

Note that the term $\nu_A^o \dot{H}_{12}(\tau)$ is already a very small correction to the dominant term $\nu_0 \dot{H}_{12}(\tau)$, such that we can safely neglect these additional terms.

## APPENDIX D: DEFAULT IMPLEMENTATION OF THE ANTIALIASING FILTER

By default, the antaliasing filter is designed from a Kaiser windowing function, with a transition band that extends from 1.1 Hz to 2.9 Hz, and minimum attenuation above 2.9 Hz of 240 dB. Note that the filter transition band extends above the Nyquist frequency, such that there will be a significant amount of aliasing during downsampling, as depicted in Fig. 16. However, since aliasing happens by reflection across the Nyquist frequency, any noise in the band $[f_s^{\text{meas}}/2, f_s^{\text{meas}} - 1 \text{ Hz}]$ will be aliased into the band $[1 \text{ Hz}, f_s^{\text{meas}}/2]$, such that it stays outside our measurement band of $[10^{-4} \text{ Hz}, 1 \text{ Hz}]$.

Analytically, we model this digital filter with the continuous, linear filter operator $\mathcal{F}$, which is applied to the beatnote frequency measurements. In the frequency domain, this is equivalent to multiplying our signals by the filter transfer function $\mathcal{F}[\mathcal{F}]$, pictured in Fig. 16 and given by
<figure><img src="image_19.png" /><figcaption>FIG. 16. Antialiasing filter transfer function magnitude. The transition band (gray) is chosen to avoid aliasing into the measurement band (blue). Unlike what is often presented, the frequency axis uses a linear scale here.</figcaption></figure>

$$ \mathcal{F}[\mathbf{F}](\omega) = \sum_{k=0}^{N} c_k e^{-i\omega k / f_s^{\text{phy}}}, \quad (\text{D1}) $$

where $c_k$ are the filter coefficients.

## APPENDIX E: EXAMPLE SIMULATION CODE SNIPPET

We give here an example simulation and processing code snippet for a simple instrumental setup and minimal processing. The simulation duration is fixed to 10 days. We inject the gravitational-wave signal of one of the strongest verification binaries expected in the LISA [44], as described in Sec. IX A. Note that the signal amplitude is rescaled such that the cumulated SNR over the simulation duration matches the full expected 4-year SNR. The relative frequency shifts induced by the gravitational-wave strain on the six LISA links are computed using LISA GW Response 2.1.2 [49].

```python
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import h5py
import scipy

from lisagwresponse import VerificationBinary
from lisainstrument
import Instrument
from pytdi import Data
from pytdi.michelson import X2, Y2, Z2

# Use a
standard set of Keplerian orbits
# These are provided by LISA Orbits
```

(Table continued)

```
(Continued)

orbit_file =
'keplerian-orbits.h5'
with h5py.File(orbit_file) as f:
    orbits_t0 = f_attrs['t0']

#
Simulation runs for 3 days at 4 Hz,
# and starts 10 s after the orbit file
dt = 0.25 # s
fs = 1
/ dt # Hz
duration = 60 * 60 * 24 * 3 # s
size = duration * fs # samples
t0 = orbits_t0 + 10 #
s

# Compute link responses to signal
source = VerificationBinary(
    period=569.4,
    distance=2089,
    masses=(0.8, 0.117),
    glong=57.7281,
    glat=6.4006,
    iota=60 *
    (np.pi / 180),
    orbits=orbit_file,
    size=size,
    dt=dt,
    t0=t0,
)

# Plot and write
the link responses to
disk
source.plot(source.t[:8000])
source.write('verification-binary.h5')

# Define the
instrumental setup, simulate
# and write the measurements to disk
instru = Instrument(
    orbits=orbit_file,
    gws='verification-binary.h5',
    laserShape='white',
    size=size,
    dt=dt,
    t0=t0,
)
instru.disable_all_noises(but='laser')
instru.write('measurements.h5')

#
Compute TDI Michelson channels
data = Data.from_instrument('measurements.h5')
```

(Table continued)
(Continued)

X2 =
X2.build(**data.args)(data.measurements)
Y2 = Y2.build(**data.args)(data.measurements)
Z2 =
Z2.build(**data.args)(data.measurements)

# Plot the TDI Michelson channels
psd = lambda
tseries: scipy.signal.welch(
    tseries,
    fs=fs,
    nperseg=2**16,
    window=('kaiser',
        30),
    detrend=None
)
freq, X2_psd = psd(X2)
freq, Y2_psd = psd(Y2)
freq, Z2_psd =

(Continued)

psd(Z2)
plt.loglog kHz, np��(10^2),
    label='X2')

plt.loglog kHz, np.sq(10^2),
    label='Y2')

(Continued)

# Plot the TDI Michelson channels
psd = lambda
tseries: scipy.signal.welch(
    tseries,
    fs=fs,
    nperseg=2**16,
    window=('kaiser',
        30),
    detrend=None
)
freq, X2_psd = psd(X2)
freq, Y2_psd = psd(Y2)
freq, Z2_psd =

(Continued)

The instrument and the measurements are computed using LISA Instrument 1.1.1 [38]. The spacecraft follow a set of Keplerian orbits with an average constellation arm length of 2.5 × 10⁹ m. The locking configuration N1-12 is used. Only a white laser frequency noise and the injected signal are present in the data, as all other effects are disabled. We compute the second-generation Michelson TDI channels X₂, Y₂, and Z₂ using PyTDI 1.2.1 [47], and plot them.

[1] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. **116**, 061102 (2016).

[2] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence, Phys. Rev. Lett. **116**, 241103 (2016).

[3] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2, Phys. Rev. Lett. **118**, 221101 (2017); **121**, 129901(E) (2018).

[4] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170608: Observation of a 19-solar-mass binary black hole coalescence, Astrophys. J. Lett. **851**, L35 (2017).

[5] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence, Phys. Rev. Lett. **119**, 141101 (2017).

[6] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. **119**, 161101 (2017).

[7] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X **9**, 031040 (2019).

[8] B. P. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW190425: Observation of a compact binary coalescence with total mass ~3.4M☉, Astrophys. J. Lett. **892**, L3 (2020).

[9] R. Abbott (LIGO Scientific and Virgo Collaborations), GW190412: Observation of a binary-black-hole coalescence with asymmetric masses, Phys. Rev. D **102**, 043015 (2020).

[10] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW190814: Gravitational waves from the coalescence of a 23 solar mass black hole with a 2.6 solar mass compact object, Astrophys. J. Lett. **896**, L44 (2020).

[11] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GW190521: A Binary Black Hole Merger with a Total Mass of 150M☉, Phys. Rev. Lett. **125**, 101102 (2020).

[12] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X **11**, 021053 (2021).

[13] R. Abbott et al. (LIGO Scientific and Virgo Collaborations), GWTC-2.1: Deep extended catalog of compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run, arXiv:2108.01045.

[14] R. Abbott et al. (LIGO Scientific, VIRGO, and KAGRA Collaborations), Observation of gravitational waves from two neutron star-black hole coalescences, Astrophys. J. Lett. **915**, L5 (2021).

[15] R. Abbott et al. (LIGO Scientific, VIRGO, and KAGRA Collaborations), GWTC-3: Compact binary coalescences observed by LIGO and Virgo during the second part of the third observing run, arXiv:2111.03606.

[16] P. Amaro-Seoane et al. (LISA Collaboration), Laser interferometer space antenna, arXiv:1702.00786.

[17] N. J. Cornish and R. W. Hellings, The effects of orbital motion on LISA time delay interferometry, Classical Quantum Gravity **20**, 4851 (2003).
[18] L. J. Rubbo, N. J. Cornish, and O. Poujade, Forward modeling of space borne gravitational wave detectors, Phys. Rev. D **69**, 082003 (2004).

[19] M. Vallisneri, Synthetic LISA: Simulating time delay interferometry in a model LISA, Phys. Rev. D **71**, 022001 (2005).

[20] M. Otto, Time-delay interferometry simulations for the laser interferometer space antenna, Ph.D. thesis, Leibniz University, Hannover, 2015.

[21] A. Petiteau, G. Auger, H. Halloin, O. Jeannin, E. Plagnol, S. Pireaux, T. Regimbau, and J.-Y. Vinet, LISACode: A scientific simulator of LISA, Phys. Rev. D **77**, 023002 (2008).

[22] O. Hartwig and J.-B. Bayle, Clock-Jitter reduction in LISA time-delay interferometry combinations, Phys. Rev. D **103**, 123027 (2021).

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

[24] J.-B. Bayle, A. Hees, M. Lilley, and C. Le Poncin-Lafitte, LISA Orbits (2022).

[25] A. Wolski, Theory of electromagnetic fields, in CERN Accelerator School: Course on RF for Accelerators (2011).

[26] L. Wissel, A. Wittchen, T. S. Schwarze, M. Hewitson, G. Heinzel, and H. Halloin, Relative-Intensity-Noise Coupling in Heterodyne Interferometers, Phys. Rev. Appl. **17**, 024025 (2022).

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

[28] R. Hellings, G. Giampieri, L. Maleki, M. Tinto, K. Danzmann, J. Hough, and D. Robertson, Heterodyne laser tracking at high Doppler rates, Opt. Commun. **124**, 313 (1996).

[29] R. W. Hellings, Elimination of Clock Jitter noise in space born laser interferometers, Phys. Rev. D **64**, 022002 (2001).

[30] M. Tinto, F. B. Estabrook, and J. W. Armstrong, Time delay interferometry for LISA, Phys. Rev. D **65**, 082003 (2002).

[31] M. Tinto and O. Hartwig, Time-delay interferometry and clock-noise calibration, Phys. Rev. D **98**, 042003 (2018).

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

[33] O. Gerberding, Phase readout for satellite interferometry, Ph.D. thesis, Leibniz University, Hannover, 2014.

[34] G. Heinzel, LISA Consortium internal technical note (2018).

[35] K. Yamamoto (private communication).

[36] W. Riley, P. L. U. Time, and F. Division, *Handbook of Frequency Stability Analysis*, NIST Special Publication (U.S. Department of Commerce, National Institute of Standards and Technology, Washington, 2008).

[37] J. J. Esteban Delgado, Laser ranging and data communication for the laser interferometer space antenna, Ph.D. thesis, Leibniz University, Hannover, 2012.

[38] J.-B. Bayle, O. Hartwig, and M. Staab, LISA Instrument (2022).

[39] C. R. Harris et al., Array programming with NumPy, Nature (London) **585**, 357 (2020).

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

[41] S. Plaszczyński (Planck HFI), Generating long streams of $1/f^a$ noise, Fluct. Noise Lett. **07**, R1 (2007).

[42] J.-B. Bayle, Simulation and data analysis for LISA: Instrumental modeling, time-delay interferometry, noise-reduction performance study, and discrimination of transient gravitational signals, Ph.D. thesis, Paris University VII, APC, 2019.

[43] J.-B. Bayle, O. Hartwig, A. Petiteau, and M. Lilley, LISA Node (2022).

[44] T. Kupfer, V. Korol, S. Shah, G. Nelemans, T. R. Marsh, G. Ramsay, P. J. Groot, D. T. H. Steeghs, and E. M. Rossi, LISA verification binaries with updated distances from Gaia Data Release 2, Mon. Not. R. Astron. Soc. **480**, 302 (2018).

[45] G. Heinzel (private communication).

[46] M. Tröbs and G. Heinzel, Improved spectrum estimation from digitized time series on a logarithmic frequency axis, Measurement **39**, 120 (2006).

[47] M. Staab, J.-B. Bayle, and O. Hartwig, PYTDI (2022).

[48] J.-B. Bayle, E. Castelli, and N. Korsakova, LISA Glitch (2022).

[49] J.-B. Bayle, Q. Baghi, A. Renzini, and M. Le Jeune, LISA GW Response (2022).

[50] H. Inchauspé, M. Hewitson, O. Sauter, and P. Wass, New LISA dynamics feedback control scheme: Common-mode isolation of test mass control and probes of test-mass acceleration, Phys. Rev. D **106**, 022006 (2022).

[51] 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).

[52] LISA Consortium internal technical note (2019).