# Authenticated time for detecting GNSS attacks

Marco Spanghero  Networked Systems Security Group
KTH Royal Institute of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>Kewei Zhang  Networked Systems Security Group
KTH Royal Institute of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>Panagiotis Papadimitratos  Networked Systems Security Group
KTH Royal Institute of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>
###### Abstract

Information cross-validation can be a powerful tool to detect manipulated,
dubious GNSS data. A promising approach is to leverage time obtained over
networks a mobile device can connect to, and detect discrepancies between the
GNSS-provided time and the network time. The challenge lies in having reliably
both accurate and trustworthy network time as the basis for the GNSS attack
detection. Here, we provide a concrete proposal that leverages, together with
the network time servers, the nearly ubiquitous IEEE 802.11 (Wi-Fi)
infrastructure. Our framework supports application-layer, secure and robust
real time broadcasting by Wi-Fi Access Points (APs), based on hash chains and
infrequent digital signatures verification to minimize computational and
communication overhead, allowing mobile nodes to efficiently obtain
authenticated and rich time information as they roam. We pair this method with
Network Time Security (NTS), for enhanced resilience through multiple sources,
available, ideally, simultaneously. We analyze the performance of our scheme
in a dedicated setup, gauging the overhead for authenticated time data (Wi-Fi
timestamped beacons and NTS). The results show that it is possible to provide
security for the external to GNSS time sources, with minimal overhead for
authentication and integrity, even when the GNSS-equipped nodes are mobile,
and thus have short interactions with the Wi-Fi infrastructure and possibly
intermittent Internet connectivity, as well as limited resources.

##  I Introduction

Cost effective, compact, high precision Global Navigation Satellite Systems
(GNSS) receivers enable the integration of precise location and timing
services in diverse, mobile platforms, ranging from cars to smart-phones to
Internet of Things devices. While bringing extensive benefits to users, the
intrinsic lack of security in GNSS [1] exposes new attack surfaces. Civilian-
grade GPS and GLONASS do not guarantee the authenticity of the received
signals. Even though this problem is being addressed by GALILEO’s Navigation
Message Authentication (NMA) ([2, 3, 4]), it will not preclude all the attacks
(e.g., Distance Decreasing [5, 6]). At the same time, easy accessibility to
software-defined radios and open source software simulators, reduce the cost
of platforms capable of mounting advanced attacks ([7, 8, 9, 10]). This is a
significant threat, successfully demonstrated in simulated attack scenarios
and it can disrupt the availability, accuracy or correctness of the GNSS
information.

This is why defensive methods have been investigated. Beyond modifying GNSS,
e.g. with NMA, various defense mechanisms can be implemented on the receiver
side ([11, 12, 13, 14, 15]). Countermeasures based on Inertial Measurement
Units (IMUs) ([16, 17]), signal processing such as Automatic Gain Control
(AGC) monitoring and phase measurements across multiple antennas, are few
examples of defenses based on integrating external sensors or exploiting
physical characteristics of the signal ([18, 19, 20, 21, 22]).

Most of the developed defense approaches do not leverage the likely network
connectivity of the GNSS equipped platforms. In fact, modern devices can have
access to several sources of information for validating GNSS obtained location
and timing. Such a capability can be valuable as it could be readily deployed
to safeguard GNSS-based position (and timing). The recent results in [23]
exploit the co-existence of multiple time sources to cross-validate time
information obtained from the GNSS receiver, demonstrating how a combination
of different time sources can be used to detect deviations in GNSS time and
detect GNSS attacks that caused such time deviations. In the context of
protecting mobile phones from spoofing, side-information based approaches,
such as time sources, inertial sensors and A-GPS (Assisted GPS) information
[24] are also investigated.

Even though existing approaches, such as [23], raise the bar, it is paramount
to have trustworthy information from the leveraged alternate sources. If not,
in the presence of an adversary that can attack both the GNSS side (e.g.,
spoofing) and the Internet side (e.g., impersonating or manipulating the
alternate source of time information), the detection would fail. Time sources
such as the Network Time Protocol (NTP), 3G/4G Radio Base Stations (RBSs) and
WiFi Access Points (APs) can provide accurate timestamps. But they do not
immediately offer protection of the time information, allowing powerful
attackers to feed the victim receiver with tampered time information. The lack
of authentication in NTP frames and WiFI beacons can be exploited by a
malicious actor that forges messages to arbitrarily skew the time of
requesting clients; it can block legitimate replies and substitute them with
bogus ones [25], aiming to overtake the attack detection system.

Beyond the aforementioned vulnerability, using NTP and Internet-based time
services assumes the GNSS-equipped device is able to frequently connect to the
remote infrastructure. In high mobility scenarios, the interaction between
different entities is limited in time by their relative speed; establishing a
link with an AP, ideally a trusted one, might be hard for highly mobile nodes.
Cellular internet access can be more effective, yet it can be challenging in
dense network areas or in rural environments. Cellular RBS physical layer can
provide accurate time information (Primary Synchronization Signals (PSS) [26,
27]), but exploiting this information comes at a significantly higher
implementation complexity for the receiver. Furthermore, as the cellular
infrastructure relies on GNSS to maintain synchronization between RBSs, an RBS
may be under the influence of an attacker we want to defend the mobile nodes
from. Hence, a connection-less approach, based on opportunistic information
could be less complex to implement and applicable, particularly in an urban
environment, with dense APS deployment.

It is, of course, critical to authenticate the timing information from the
time servers and the APs. To address this, we investigate the use of Network
Time Security (RFC 7380 [28]) and we propose an authentication scheme for WiFi
beacons, studying different configurations for authenticated time broadcast.
The result is a framework that supports authenticated time provision to cross-
validate GNSS-obtained time, leveraging mainstream mobile systems. Building on
the approach of [23] and [24], our scheme supports different options to
safeguard time information. In particular, we investigate how to achieve this
with low computation and communication cost, and reasonable credential
management complexity. We are interested in a solution that, while being easy
to deploy, remains effective even with intermittent Internet connectivity.
Existing numerous timing servers, ubiquitous, redundant WiFi APs for
resilience to misbehavior, and low cost cryptographic protection are the
ingredients of our proposal.

The paper is organized as follows: Section II presents the system and attacker
model, Section III discusses how to safeguard external time information
towards detecting GNSS attacks. Section IV presents implementation details and
the computational overhead of our proposal. Section V discusses future
developments and concludes.

##  II System and Attack Model

We consider a mobile client that interacts with a fixed terrestrial network
and obtains GNSS information through a GNSS receiver. Clients can access
Internet over WiFi APs, possibly authenticated. Clients receive broadcast
messages from APs within range, without necessarily being associated with
those APs. The APs advertise WiFi networks with 802.11 protocol compliant
beacons and in a single hop network, where all the receivers obtain a wireless
broadcast message (beacon) practically simultaneously.

Typically, GNSS attackers in the literature influence victim GNSS receivers.
We extend this model and consider an adversary that can also attack the
terrestrial network(s) the GNSS-enabled devices connect to. The attacker can
attempt to manipulate the alternative sources of information (WiFi beacons and
NTP frames) in a coordinated effort towards deceiving the GNSS receiver,
notably towards defeating the GNSS cross-validation approach.

We do not limit the capacity and objectives of the attacker in terms the
deviations it can cause to the GNSS Position Velocity and Timing (PVT)
solution. In addition, we do not presume which specific methods can be used by
the adversary.

Let $\rho$ be a n-vector of pseudo-ranges, $H$ be a n x 4 matrix of receiver
state observations and $x$ be the PVT solution for the corresponding
observation, nominally $\bm{x}=[x,y,z,T_{GNSS}]$, $\bm{v}$, n x 1 vector, be a
measurement noise. Let $\bm{f}$ be an n x 1 vector that corresponds to the
adversarial action; then the following models a generic GNSS receiver under
adversarial influence:

$\bm{\rho}=\bm{H}\cdot\bm{x}+\bm{f}+\bm{v}$

Adversaries able to craft $\bm{f}$ to selectively tamper only with the
location, or only with the velocity, part of the PVT solution, while
intentionally not modifying the time information, cannot be thwarted by the
method considered here or by any other time-based spoofing detection method
[29].

The attacker, provided with WiFi and, possibly, other radios, and Internet
access can attempt forging 802.11 beacon frames. The content of the beacon
frames can be either malicious (containing misleading timing information) or
meaningless, to consume resources of the receiver. The attacker can spoof or
meacon single beacons or sequences of beacons and NTP frames to mask the
misbehavior in the GNSS information.

The attacker and the victim can be subject to the same propagation delay when
receiving messages from an AP. Furthermore, the adversary may have higher
computational power than the clients and the APs, but she cannot obtain system
credentials and cryptographic keys.

##  III Safeguarding Time-based GNSS verification

Augmentation information from externally provided time references can be
useful only under the assumption that the terrestrial networking
infrastructure is secure. Consider Figure 1, where an adversary controls part
of the infrastructure and can emulate/forge bogus timing information, in a
coordinated effort in parallel to the GNSS attack. Within the attacker’s
influence area, mobile clients are subject to the GNSS spoofer actions and
they receive a combination of legitimate and malicious time information,
without being able to distinguish the former from the latter.

Figure 1: Mobile node transversing several AP coverage zones. One node can be
in range of multiple APs at the same time. The adversarial zone is influenced
by the combined effort of an attacker spoofing GNSS receiver while
disseminating bogus beacon information, or even introducing rogue APs.

Wi-Fi APs can offer secure links to mobile clients in range and to the
Internet, provided that the client has enough time to establish a connection
and that it has the necessary credentials (i.e., Wi-Fi password, or
certificate). Such a link can be leveraged to access Internet time services,
to cross-check the GNSS time with that provided by the Network Time Protocol
(NTP) infrastructure. If the attacker can influence only the GNSS receiver of
the receiver, this approach is valid and can detect adversaries that tamper
with the GNSS, if the modification of the GNSS time is larger than the
accuracy provided by NTP. An advanced attacker that hijacks the NTP server
connection could mask the adversarial action on the GNSS receiver by forging
bogus NTP replies. This is possible even if the connection between the AP and
the mobile client is secured, assuming the attacker knows which NTP servers
are used by the client and the interaction between the client and the NTP
server is unauthenticated. In case an attacker is able to tamper with the NTP
information, Network Time Security protocol (NTS) can be used to provide
authenticated NTP frames.

Although this raises the bar for a potential attacker, there are limitations
due to the structure of the Wi-Fi protocol. As an example, mobile clients are
able to connect to only one AP at the time. There are instances when the
mobile device is not able to connect to any remote time provider. Clearly, the
mobile node can connect to the APs it has valid credentials for. Furthermore,
Wi-Fi associations are slow and the authentication process introduces latency.
Depending on the mobility of the client, it might be infeasible for it to
establish a secure connection.

In such situations, with limited (secure) access to the Internet, the Wi-Fi
infrastructure can still be useful, in an opportunistic manner. The client
does not connect to the APs but it exploits the management and advertisement
information distributed by the APs, more specifically Wi-Fi beacons, which
provide a rich source of timing information. Defined in Section 9.3.3 of the
IEEE Wireless LAN MAC Specification [30] as management frames, they are used
to advertise networks and AP capabilities. Additionally, they contain a
64-bit, microsecond resolution timestamp that can be used to synchronize
devices within the range of the AP; providing a precise source of relative
time information.

On the other hand, as beacons are not authenticated, based on the current
standard and deployment status they provide a source of time of low
trustworthiness. It is relatively simple for an attacker to spoof/meacon
existing AP beacons or emulate new APs with timestamps that can support the
adversarial action on the GNSS side. This would still require a sophisticated
adversary to orchestrate a coordinated attack on the GNSS receiver and the
terrestrial infrastructure. Interestingly, this might be comparably easier
than hijacking the connection to an NTP server. To leverage unauthenticated
Wi-Fi beacons in an attempt to validate the GNSS-provided time, the mobile
client might collect information from a diverse set of APs. This comes at the
cost of elevated overhead and would require assumptions on the maximum number
of APs the adversary is capable of controlling/emulating. Overall, a scheme
difficult to be effective.

Authenticated timestamps would provide an additional layer of security for the
GNSS attack detection. Several options can be considered. If the mobile node
can associate itself with the AP, it could use a shared symmetric key to
authenticate beacons. Establishing such a symmetric key can be done in several
standard ways, but it is subject to the latency of establishing a secure
connection with the AP. Moreover, even so, it would allow any of the mobile
nodes to transmit such beacons.

The APs can authenticate broadcast beacons. A naive solution would be to use
asymmetric/public key cryptography to digitally sign every message/beacon. The
creation and verification of a digital signature for each beacon would
introduce significant overhead due to the relatively high beacon rate
($\approx$ $10\text{\,}\mathrm{beacon}\text{\,}{\mathrm{s}}^{-1}$). This would
be costly for the APs and not sustainable for constrained mobile devices,
enabling clogging or resource depletion Denial of Service (DoS) attacks [31].

###  III-A Efficient beacon authentication

The approach used here is to predominantly apply symmetric key cryptography
for broadcast authentication with the necessary and infrequent use of
asymmetric key digital signatures. An approach such as [32] provides efficient
broadcast authentication. To achieve lightweight authentication of beacons
(and timestamps) here, the hash chain elements are used as symmetric keys.
Each beacon is protected by an authentication code ($MAC^{n}$) calculated with
the hash chain element that will be disclosed at the following transmission
cycle, the subsequent beacon. Upon receipt of that latter beacon, the mobile
nodes in range can authenticate the previous, sufficiently early received,
beacon. Beacons very close to (and trivially after) the key disclosure are
rejected, given the single hop broadcast.

Fig. 2 shows the generation of the hash chain and the order of the use of its
elements. The last element, the anchor, is digitally signed by the AP, to bind
a hash chain to the specific AP. This essentially allows any client that moves
within the AP range to obtain a digitally signed anchor for the hash chain and
validate any subsequent beacon authentication code. The public key that
corresponds to the private key, used for digital signing, is certified by a
trusted third party, a Certification Authority or Public Key Infrastructure
(PKI). For simplicity, here, each AP identity is bound to its public key. The
resultant AP certificate can either be made available over the Internet or it
can be periodically broadcast by the AP, attached to a subset of the
transmitted beacons (leveraging an idea from [33, 34]).

Figure 2: Chain structure and generated payloads. Infrequent digital
signatures are applied to the hash chain elements

The hash chain is pre-computed by the AP and it can cover long periods of
time. Due to the high frequency of the beacons and given that each hash chain
element is used as a one-time key material (thus requiring one element per
each beacon timestamp), the length of the hash chain can be significant. As
mobile clients can join the coverage area of an AP at any point in time, if
the client is provided only with the initial anchor to the hash chain, it
might be forced to perform a large number of hash calculations before
confirming the validity of any hash chain element and message. The AP can
digitally sign intermediate anchors, to be disclosed by the AP as the usage of
the hash chain progresses. By verifying the digital signature of the
intermediate anchor, the mobile client calculates a much lower number of
hashes initially. Different strategies are possible to distribute anchors,
similarly to the distribution of the AP certificates.

We can exploit the regularity of the beacon emission. The hash chain is
generated in advance, alongside with the signed anchors. MACs are computed and
packets are assembled and scheduled for transmission at precise times,
allowing the AP to execute all the cryptographic operations outside the
transmission critical path.

Finally, a reminder that the scheme does not operate in isolation but together
with possibly less frequent access to internet time servers. This provides a
time reference to mobile nodes, useful for this broadcast authentication as
well as the ability to detect replays of sequences of authenticated beacons.
To be more effective against such attacks, additional information on the AP
and the network context can be used.

##  IV Test environment Setup and results

The method described in Section III-A is implemented in a laboratory test
setup. The authentication scheme for the Wi-Fi beacons allows mobile clients
to detect APs, obtain digitally signed anchors and verify the beacons payload.
Figure 3 presents the implemented system architecture and the interaction of
its components. The authentication scheme is implemented at the application
layer and we evaluate a deployment of the system with one client, one AP
(supporting our authentication scheme) and several APs without authentication.
We pre-install the certificate of the AP in the client. This does not affect
the generality of the setup, as certificates and anchors can be broadcasted.

The test system is composed of the following components:

  * •

Host machine: Intel I7 CPU running a Linux system whose kernel supports high
precision timing

  * •

Host Wi-Fi card: Intel Corporation Wireless-AC 9260

  * •

Host virtual Wi-Fi card: mac80211_hwsim kernel module

  * •

Wi-Fi beacons: from surrounding APs in the test area

  * •

Secure Wi-Fi beacons: dedicated beacon emulator in the test area

  * •

NTP servers: three servers in Sweden (provided by Netnod [35])

  * •

NTS servers: three servers in Sweden (provided by Netnod [35])

  * •

GPS PVT rate: $1\text{\,}\mathrm{H}\mathrm{z}$

  * •

Ublox EVK-6T evaluation kit [36]

To emulate a mobility scenario, the client interacts with the AP by receiving
the beacon broadcast information over Wi-Fi, in a non connected state.
Additionally, it communicates with the NTS/NTP infrastructure, over a wireless
link.

Figure 3: Framework for secure time broadcasting and Wi-Fi beacons
authentication.

We developed a setup that allows to emulate a beacon transmitter. We rely on
high accuracy software timers to emit beacons at a constant rate of 100TU. To
transmit custom Wi-Fi beacon frames, the wireless cards are configured in
”monitor” mode. This allows leveraging the packet injection capabilities of
the wireless interface to send raw frames without directly interacting with
the driver. To achieve the highest accuracy, it would be necessary to modify
the hardware driver for the Wireless interface. To the best of our knowledge,
in the current 802.11 implementation in Linux, beacons are configured in user
space as static templates. The information is received in the kernel module
where channel state information is added and ultimately it is scheduled for
transmission by the driver, where the final timestamp is added. This static
structure does not allow directly supporting our scheme, as information
contained in the Wi-Fi beacons needs to be updated at each transmission cycle.
Nevertheless, there is no limitation with the current standard to have dynamic
beacon payloads, provided that the compulsory information is included.

To test the concept, the mac80211_hwsim Wi-Fi simulation driver is used. This
kernel module allows setting a virtual Wi-Fi devices that can either be
configured as a Station (STA) or Access Point (AP). The virtual Wi-Fi device
supports packet injection and behave like ideal Wi-Fi devices. The major
limitation is that it does not take into account real propagation medium or
coexistence with other wireless networks. The generation, transmission and
detection of the beacons uses Scapy framework [37] and is Python3.7 based.

The specialized authentication payload is included in the Wi-Fi beacon
structure, at the Vendor Specific tag fields. This is allowed by the protocol,
as long as the assembled packet is smaller than the maximum size (2320 octets
[30]). Timestamps are placed in the dedicated mandatory field of the beacon
consisting of 8 octets, corresponding to a 64 bit microsecond timestamp. All
clients that do not support the Vendor Specific fields we define, simply
ignore them; our implementation can coexist with existing wireless networks
without disrupting their behavior. The implementation does not depend on a
specific signature method and it could be extended to support different
signatures.

The integration of Wi-Fi APs with authenticated NTS servers provides a solid
basis to detect time-shifting GNSS spoofing attacks and more generally
assesses the health of the GNSS information. We evaluate the potential
overhead introduced by the authentication of the NTP information in NTS and by
our time authentication scheme, to understand whether the same accuracy as the
non authenticated options is achievable. We expect that authentication causes
minimal, if any, degradation.

###  IV-A Evaluation results

No substantial performance degradation was observed while using NTS for
retrieving absolute timestamps from a remote server over 4G networks when
compared to NTP in the same conditions. Table IV-A shows the results of a
6-hour test against three NTS servers provided by NetNod in Sweden. Although
there is a small degradation in the offset of the synchronizing client when
using NTS, the standard deviation of RTD and offset of the two methods are
compatible, and so is their accuracy.

Source |  RTD |  Offset
---|---|---
nts.sth2.ntp.se 2a01:3f7:2:202::201 |  $\mu:2.825\mathrm{e}{-2}$ |  $\mu:1.228\mathrm{e}{-3}$
$\sigma:7.200\mathrm{e}{-3}$ |  $\sigma:4.331\mathrm{e}{-3}$
nts.sth1.ntp.se 2a01:3f7:2:52::10 |  $\mu:2.765\mathrm{e}{-2}$ |  $\mu:9.380\mathrm{e}{-4}$
$\sigma:7.151\mathrm{e}{-3}$ |  $\sigma:3.309\mathrm{e}{-3}$
nts.sth.ntp.se 2a01:3f7:2:62::10 |  $\mu:2.774\mathrm{e}{e-2}$ |  $\mu:1.042\mathrm{e}{-3}$
$\sigma:6.970\mathrm{e}{-3}$ |  $\sigma:3.131\mathrm{e}{-3}$

TABLE I: Performance of NTS over cellular 4G connection (measurements in
seconds ($\mathrm{s}$))

Source |  RTD |  Offset
---|---|---
sth1.ntp.se 2a01:3f7:2:1::1 |  $\mu:2.673\mathrm{e}{-2}$ |  $\mu:3.177\mathrm{e}{-4}$
$\sigma:6.512\mathrm{e}{-3}$ |  $\sigma:3.409\mathrm{e}{-3}$
sth2.ntp.se 2a01:3f7:2:2::1 |  $\mu:2.631\mathrm{e}{-2}$ |  $\mu:4.223\mathrm{e}{-4}$
$\sigma:5.526\mathrm{e}{-3}$ |  $\sigma:2.975\mathrm{e}{-3}$
npt3.sptime.se 2001:6b0:42:3::123 |  $\mu:3.931\mathrm{e}{-2}$ |  $\mu:-5.802\mathrm{e}{-3}$
$\sigma:1.375\mathrm{e}{-3}$ |  $\sigma:6.985\mathrm{e}{-3}$

TABLE II: Performance of NTP over cellular 4G connection (measurements in
($\mathrm{s}$))

(a) Histogram of NTS Round Trip Delay distribution for three NTS servers in
Sweden. Measurement details are provided in Table IV-A

(b) Histogram of NTS Offset distribution for three NTS servers in Sweden.
Measurement details are provided in Table IV-A

Figure 4: NTS performance plots

(a) Histogram of NTP Round Trip Delay distribution for three NTP servers in
Sweden. Measurement details are provided in Table IV-A

(b) Histogram of NTP Offset distribution for three NTP servers in Sweden.
Measurement details are provided in Table IV-A

Figure 5: NTP performance plots

The results presented for NTS and NTP pertain to clients (mobile nodes)
connected to the Internet for extended periods of time. If only short lived
connections are available, the achievable accuracy depends on the local clock
stability. Intuitively, when a remote time source is not available, the mobile
client can still use its local clock. If a remote time source is available at
any time, it can be used to continuously correct the local clock to provide a
stable reference. If the remote time source is intermittently available, the
local clock will drift depending on its quality.

The security improvements obtained from NTS come at an overhead due to the
additional operations required by the NTS protocol. The setup of the
cryptographic keys introduces an initial overhead, even if it is a one-time
operation. The Message Authentication Codes computation for each NTS
transaction introduces limited overhead. According to the analysis in [38],
based on a previous version of the NTS standard, NTS generates approximately
six times more computational load compared to NTP, for the same number of
transmitted and received packages. On the other hand, only a fraction of this
overhead is caused by the additional authentication (approximately 3%). In the
context of this work, the computational overhead is one of the considered
aspects. The accuracy of NTS compared to NTP and the increased level of
robustness still make NTS a preferable solution for GNSS spoofing detection.
Additionally, the polling cycle can be regulated to limit the overhead caused
by the NTS protocol.

Considering the proposed extension for Wi-Fi, the scheduling of the beacons
for transmission preserves the time information accuracy. To achieve accuracy
within the standard’s tolerance, the beacon scheduling system must be
implemented at the driver level; while the pre-computation and assembly of the
payloads can be performed at the application level.

According to the 802.11 physical layer specification, each beacon is
transmitter at a specific Target Beacon Transmission Time (TBTT). If the
medium is available at the requested TBTT (in out case 100TU), the scheduled
beacon is transmitted. If the medium is not available at the requested TBTT,
it will be accessed as per the 802.11 physical layer specification (section
11.1.2.1). In this case the beacon may be delayed because of the Carrier Sense
Multiple Access (CSMA) deferral. All the subsequent beacon transmissions will
be scheduled at the nominal TBTT.

The standard (section 11.1.2.4) mandates a minimum accuracy of the Time
Synchronization Function (TSF) timer to be 100ppm. This limits the accuracy of
the synchronization between STA and AP. Specifically, for out application, it
limits the accuracy of the timestamps in the beacons. The synchronization
accuracy of $\pm 0.01\%$ over one TBTT of 100TU is within
$25\text{\,}\mathrm{\SIUnitSymbolMicro s}$, accounting for other delays due to
the transmission in the physical layer.

The current setup is based on a software simulator able to provide a realistic
beacon distribution, with emission interval (TBTT and TSF synchronization
error) standard deviation comparable to an hardware transmitter. This effect
is currently simulated by adding a 100ppm error to the beacon timestamps.
Nevertheless, the current implementation shows the possibility of achieving
authentication with low overhead, given the rate of the beacons. Preliminary
measurements show that the results in terms of accuracy (and in terms of
effectiveness in detecting the GNSS attacks) will be very close to the ones
presented in [23], where the implementation had no security features.

Similarly to the measurements performed in [23], the GNSS time is modified by
progressively increasing time offsets until a total bias on the real GNSS time
of $120\text{\,}\mathrm{ms}$ is obtained at time $t_{0}$. After $t_{0}$, the
time bias is kept constant. A detection threshold based on the beacon
timestamp accuracy is defined ($\epsilon_{Wi-
Fi}=$25.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, obtained from the simulated
beacon transmitter), given by the standard deviation of the received beacon
timestamps. Different observation windows are identified, corresponding to
increasing beacon intervals ($T_{Window}={1,3,5}$\mathrm{s}$$). Figure 6 shows
the GNSS spoofing detection capability of the present implementation of the
setup, with different observation windows. Only beacons that pass the
authentication test are used to detect potential attackers.

(a) Observation window of $1\text{\,}\mathrm{s}$

(b) Observation window of $3\text{\,}\mathrm{s}$

(c) Observation window of $5\text{\,}\mathrm{s}$

Figure 6: GNSS spoofing detection with authenticated beacon timestamps.
Results based on a simulated transmitter. The accuracy of the detection system
improves with longer observation windows

We briefly discuss the computational cost of the proposed Wi-Fi augmentation
scheme. The cost of an HMAC operation is low in modern platforms, considerably
lower than verifying a digital signature. For a standard Intel i7 laptop CPU,
there is one order of magnitude difference in the delay between calculating
one HMAC and verifying one signature. The amortization of the few signature
verifications over a large number of beacons reduces the cost of verification,
by reducing the number of the most computationally expensive operation. Table
III presents the typical delay of each cryptographic operation, when performed
on an Intel i7-8750H mobile CPU.

As cost metric, we consider the time spent by the client to authenticate the
beacons time information. We define $C_{SIG}$, cost of validating one digital
signature, $C_{HMAC}$, cost of calculating one HMAC and $C_{HASH}$, cost of
computing one hash element. If the receiver does not have any prior
information on the state of the AP, the cost of the first authentication is
calculated as $C^{\prime}=C_{SIG}+C_{HMAC}+k\cdot C_{HASH}$; plus the cost to
validate the AP certificate if needed. The cost of first verification
increases linearly with the distance $k$ of the current hash element element
from the most recent digitally signed anchor, but it is bounded by the
frequency of the signed anchors. In the worst case scenario, the client is not
able to retrieve any digitally signed anchor and has to compute the entirety
of the $H_{C}$. All the subsequent authentications have a cost
$C^{\prime\prime}=C_{HMAC}+C_{HASH}$; or multiple $C_{HASH}$ in case of loss
of some beacons. Clients that frequently transition between coverage areas of
different APs are subject to a first verification every time they enter a new
AP area.

Operation |  Op # |  Avg. Delay |  Total
---|---|---|---
Signature_verify (RSA1024) |  1 |  $52\text{\,}\mathrm{\SIUnitSymbolMicro s}$ |  $52\text{\,}\mathrm{\SIUnitSymbolMicro s}$
HMAC (SHA256) |  |  $7.67\text{\,}\mathrm{\SIUnitSymbolMicro s}$ |  $7.67\text{\,}\mathrm{\SIUnitSymbolMicro s}$
HASH (SHA256) |  $k$ |  $1.64\text{\,}\mathrm{\SIUnitSymbolMicro s}$ |  $k\cdot$$1.64\text{\,}\mathrm{\SIUnitSymbolMicro s}$
TABLE III: Cost of cryptographic operations for authenticating Wi-Fi beacons
at the receiver. Computations are performed on an Intel i7 cpu.

For example, consider an urban scenario, with the client moving at an average
speed of $15\text{\,}\mathrm{km}\text{\,}{\mathrm{h}}^{-1}$, and a trajectory
of $100\text{\,}\mathrm{m}$ within the coverage area of an AP. The total
transit time for the AP coverage is 24.03$\mathrm{s}$ and during this period
the total number of beacons a client can receive by one AP is $\approx$
234$\mathrm{beacon}$. For intermediate signed anchor periods longer than the
travel time of the client within the AP coverage area (e.g., the signed anchor
is updated every $60\text{\,}\mathrm{s}$, corresponding roughly to
$600\text{\,}\mathrm{beacon}$), the client needs to verify on average only one
anchor for every new AP, with a worst case of $k=599$. In this case, the time
spent for authenticating the beacons is $\approx$3.21\text{\,}\mathrm{ms}$$,
corresponding to $0.013\%$ of the transit time.

For APs with hash chains valid for long periods of time (i.e., one hash chain
valid for $24\text{\,}\mathrm{h}$), the usage of intermediate anchors
effectively reduce the cost of joining the coverage area of an AP. Without
intermediate anchors, a client joining the coverage area close to the hash
chain expiration time would be penalized. The distribution of the intermediate
anchors can be performed directly by the AP, at regular intervals, or through
the Internet. The former solution requires larger beacon payloads, while the
latter supposes (at least intermittent) network connectivity.

##  V Conclusions

Alternative sources of time are beneficial in detecting time-manipulating GNSS
spoofers. Opportunistic time providers, Wi-Fi beacons, and dedicated timing
infrastructures, NTP, provide largely available, precise sources of time
information. While proving effective in detecting misbehaving GNSS receivers,
i.e., when the client is able to obtain enough information, unauthenticated
time information can be relied upon.

Different options to secure augmentation time information are investigated,
notably Network Time Security (NTS) and modified Wi-Fi beacons to support
authentication. This scheme imposes limited overhead, does not disrupt the
normal operation of the Wi-Fi APs and can be deployed readily.

Furthermore, considering NTP and NTS, there is no significant performance loss
due to the authentication of the NTP frames, making NTS a valid solution to
provide GNSS-based receivers secure Internet time.

##  VI Acknowledgment

This work supported by the Swedish Foundation for Strategic Research (SSF)
SURPRISE project and the KAW Academy Fellowship Trustworthy IoT project. In
addition, we would like to thank Dr. Hongyu Jin for the constructive feedback
and discussion during the development of this work.

## References

  * [1] J. Volpe, “Vulnerability assessment of the transportation infrastructure relying on the global positioning system,” 2001.
  * [2] The European Union, “Galileo Open Service Signal-In-Space Interface Control Document (OS SIS ICD), Issue 1.2,” 2015.
  * [3] I. Fernández-Hernández, V. Rijmen, G. Seco-Granados, J. Simon, I. Rodríguez, and J. D. Calle, “A navigation message authentication proposal for the Galileo open service,” _NAVIGATION_ , vol. 63, no. 1, pp. 85–102, 2016.
  * [4] P. Walker, V. Rijmen, I. Fernandez-Hernandez, G. Seco-Granados, J. Simón, J. Calle, and O. Pozzobon, “Galileo open service authentication: a complete service design and provision analysis,” in _Proceedings of the 28th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2015)_ , 2015, pp. 3383–3396.
  * [5] K. Zhang and P. Papadimitratos, “GNSS receiver tracking performance analysis under distance-decreasing attacks,” in _Proceedings of the 2015 International Conference on Location and GNSS (ICL-GNSS 2015)_ , Gothenburg, Sweden, 2015, pp. 1–6.
  * [6] K. Zhang and P. Papadimitratos, “On the Effects of Distance-decreasing Attacks on Cryptographically Protected GNSS Signals,” in _Proceedings of the 2019 International Technical Meeting of The Institute of Navigation (ION ITM 2019)_ , Reston, Virginia, 2019, pp. 363–372.
  * [7] K. C. Zeng, S. Liu, Y. Shu, D. Wang, H. Li, Y. Dou, G. Wang, and Y. Yang, “All your GPS are belong to us: Towards stealthy manipulation of road navigation systems,” in _27th USENIX Security Symposium_ , 2018, pp. 1527–1544.
  * [8] D. Schmidt, K. Radke, S. Camtepe, E. Foo, and M. Ren, “A survey and analysis of the GNSS spoofing threat and countermeasures,” 2016.
  * [9] K. Zeng, Y. Shu, S. Liu, Y. Dou, and Y. Yang, “A practical GPS location spoofing attack in road navigation scenario,” in _HotMobile 2017 - Proceedings of the 18th International Workshop on Mobile Computing Systems and Applications_. New York, New York, USA: Association for Computing Machinery, Inc, feb 2017, pp. 85–90. [Online]. Available: http://dl.acm.org/citation.cfm?doid=3032970.3032983
  * [10] L. Huang and Q. Yang, “Low-Cost GPS Simulator GPS Spoofing by SDR,” _DEFCON’15_ , 2015. [Online]. Available: https://www.researchgate.net/publication/286330869_Low-cost_GPS_simulator_-_GPS_spoofing_by_SDR
  * [11] P. Papadimitratos and A. Jovanovic, “Method to secure GNSS based locations in a device having GNSS receiver,” Apr. 17 2012, US Patent 8,159,391.
  * [12] M. L. Psiaki and T. E. Humphreys, “GNSS Spoofing and Detection,” _Proceedings of the IEEE_ , vol. 104, no. 6, pp. 1258–1270, 2016.
  * [13] P. Papadimitratos and A. Jovanovic, “Protection and Fundamental Vulnerability of GNSS,” in _IEEE International Workshop on Satellite and Space Communications (IEEE IWSSC)_ , Toulouse, France, October 2008, pp. 167–171.
  * [14] K. Zhang and P. Papadimitratos, “Secure Multi-constellation GNSS Receivers with Clustering-based Solution Separation Algorithm,” in _IEEE Aerospace Conference_ , Big Sky, MT, USA, March 2019, pp. 1–9.
  * [15] ——, “Safeguarding NMA Enhanced Galileo OS Signals from Distance-Decreasing Attacks,” in _International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+)_ , Miami, FL, USA, September 2019, pp. 4041–4052.
  * [16] J. T. Curran and A. Broumendan, “On the use of low-cost IMUs for GNSS spoofing detection in vehicular applications,” in _Proceedings of the ITSNT_ , 2017.
  * [17] S. Khanafseh, N. Roshan, S. Langel, F.-C. Chan, M. Joerger, and B. Pervan, “GPS spoofing detection using RAIM with INS coupling,” in _2014 IEEE/ION Position, Location and Navigation Symposium-PLANS 2014_. IEEE, 2014, pp. 1232–1239.
  * [18] P. Y. Montgomery, T. E. Humphreys, and B. M. Ledvina, “Receiver-autonomous spoofing detection: Experimental results of a multi-antenna receiver defense against a portable civil GPS spoofer,” in _Proceedings of the ION International Technical Meeting_ , 2009, pp. 124–130.
  * [19] M. L. Psiaki, S. P. Powell, and B. W. O’hanlon, “GNSS spoofing detection using high-frequency antenna motion and carrier-phase data,” in _Proceedings of the 2013 International Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2013)_ , 2013, pp. 2949–2991.
  * [20] F. Bastide, D. Akos, C. Macabiau, and B. Roturier, “Automatic gain control (AGC) as an interference assessment tool,” in _ION GPS/GNSS 2003, 16th International Technical Meeting of the Satellite Division of The Institute of Navigation_ , 2003, pp. pp–2042.
  * [21] H. Borowski, O. Isoz, S. Lo, F. Eklöf, and D. Akos, “Detecting false signals with automatic gain control,” _GPS World_ , vol. 23, no. 4, pp. 38–43, 2012.
  * [22] M. L. Psiaki, B. W. O’hanlon, S. P. Powell, J. A. Bhatti, K. D. Wesson, and T. E. Humphreys, “GNSS spoofing detection using two-antenna differential carrier phase,” 2014.
  * [23] K. Zhang, M. Spanghero, and P. Papadimitratos, “Protecting gnss-based services using time offset validation,” in _2020 IEEE/ION Position, Location and Navigation Symposium (PLANS)_ , 2020, pp. 575–583.
  * [24] S. Ceccato, F. Formaggio, G. Caparra, N. Laurenti, and S. Tomasin, “Exploiting side-information for resilient gnss positioning in mobile phones,” in _2018 IEEE/ION Position, Location and Navigation Symposium (PLANS)_.
  * [25] A. Malhotra, I. E. Cohen, E. Brakke, and S. Goldberg, “Attacking the Network Time Protocol,” 2017. [Online]. Available: http://dx.doi.org/10.14722/ndss.2016.23090
  * [26] 3GPP, “Radio subsystem synchronization TS 05.10,” https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=265, 2003\.
  * [27] ——, “GSM/EDGE Radio subsystem synchronization TS 45.010,” https://portal.3gpp.org/desktopmodules/Specifications/SpecificationDetails.aspx?specificationId=2712, 2019\.
  * [28] D. F. Franke, D. Sibold, K. Teichel, M. Dansarie, and R. Sundblad, “Network Time Security for the Network Time Protocol,” Internet Engineering Task Force, Internet-Draft draft-ietf-ntp-using-nts-for-ntp-28, Mar. 2020, work in Progress. [Online]. Available: https://datatracker.ietf.org/doc/html/draft-ietf-ntp-using-nts-for-ntp-28
  * [29] T. Humphreys, J. Bhatti, D. Shepard, and K. Wesson, “The Texas Spoofing Test Battery: Toward a Standard for Evaluating GPS Signal Authentication Techniques,” Tech. Rep.
  * [30] IEEE, “Ieee standard for information technology–telecommunications and information exchange between systems local and metropolitan area networks–specific requirements part 11: Wireless lan medium access control (mac) and physical layer (phy) specifications,” _IEEE Std 802.11-2012 (Revision of IEEE Std 802.11-2007)_ , pp. 1–2793, 2012.
  * [31] H. Jin and P. Papadimitratos, “Dos-resilient cooperative beacon verification for vehicular communication systems,” _Ad Hoc Networks_ , vol. 90, p. 101775, November 2019, recent advances on security and privacy in Intelligent Transportation Systems. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S1570870518307108
  * [32] A. Perrig, R. Canetti, D. X. Song, and J. D. Tygar, “Efficient and Secure Source Authentication for Multicast,” in _Proceedings of the Network and Distributed System Security Symposium,{NDSS} 2001, San Diego, California, {USA}_. The Internet Society, 2001. [Online]. Available: https://www.ndss-symposium.org/ndss2001/efficient-and-secure-source-authentication-multicast/
  * [33] G. Calandriello, P. Papadimitratos, J.-P. Hubaux, and A. Lioy, “On the Performance of Secure Vehicular Communication Systems,” _IEEE Transactions on Dependable and Secure Computing (IEEE TDSC)_ , vol. 8, no. 6, pp. 898–912, November 2011.
  * [34] H. Jin and P. Papadimitratos, “Proactive Certificate Validation for VANETs,” in _IEEE Vehicular Networking Conference (IEEE VNC)_ , Columbus, OH, USA, December 2016.
  * [35] “Netnod’s NTP service. Keeping Sweden in sync. Guaranteed.” [Online]. Available: https://www.netnod.se/ntp/network-time-protocol-%28ntp%29-services-?gclid=EAIaIQobChMIjM61o5ep6gIVwYSyCh1qyQK_EAAYASACEgIL7vD_BwE
  * [36] ublox, “EVK-6 u-blox Evaluation Kits User Guide,” https://www.u-blox.com/sites/default/files/products/documents/EVK-6_UserGuide_(28GPS.G6-EK-10040).pdf, Accessed 2020-02-01.
  * [37] Philippe Biondi and the Scapy community, “Scapy,” 2020. [Online]. Available: https://scapy.net/
  * [38] M. Langer, K. Teichel, D. Sibold, and R. Bermbach, “Time synchronization performance using the network time security protocol,” in _2018 European Frequency and Time Forum (EFTF)_ , 2018, pp. 138–144.


# Doppler effect in TianQin time-delay interferometry

Lu Zheng  Shutao Yang  Xuefeng Zhang<EMAIL_ADDRESS>MOE Key Laboratory
of TianQin Mission, TianQin Research Center for Gravitational Physics $\&$
School of Physics and Astronomy, Frontiers Science Center for TianQin,
Gravitational Wave Research Center of CNSA, Sun Yat-sen University (Zhuhai
Campus), Zhuhai 519082, China

###### Abstract

The current design of space-based gravitational wave detectors utilizes
heterodyne laser interferometry in inter-satellite science measurements.
Frequency variations of the heterodyne beatnotes are predominantly caused by
the Doppler effect from relative satellite motion along lines of sight.
Generally considered to be outside the measurement band, this Doppler
frequency shift appears to have been overlooked in numerical simulations of
time-delay interferometry (TDI). However, the potential impact on the
implementation of TDI should be assessed. The issue is particularly relevant
to TianQin that features geocentric orbits, because of strong gravity
disturbances from the Earth-Moon system at frequencies $<1\times 10^{-4}$ Hz.
In this proof-of-principle study, based on high-precision orbital data
obtained from detailed gravity field modeling, we incorporate the Doppler
shift in the generation of TianQin’s beatnote phase signals. To remove the
large-scale Doppler phase drift at frequencies $<1\times 10^{-4}$ Hz, we
develop a high-performance high-pass filter and consider two possible
processing sequences, i.e., applying the filter before or after TDI
combinations. Our simulation results favor the former and demonstrate
successful removal of the low-frequency gravity disturbances for TianQin
without degrading the TDI performance, assuming 10 m pseudo-ranging
uncertainty. The filtering scheme can be used in developing the initial noise-
reduction pipeline for TianQin.

††preprint: APS/123-QED

##  I  Introduction

Gravitational waves (GWs) are excellent information carriers and are expected
to provide a unique and new method for observing and studying the universe.
With the development of technology, GWs have been detected by the LIGO
detectors for the first time in 2015 utilizing laser interferometry [1]. LIGO
have an armlength of 4 km, and their most sensitive frequency band is about
$10-10^{3}$ Hz [2]. Limited by available armlengths and seismic noise, in
order to detect GWs in the mHz frequency band, it is necessary to send GW
detectors to space. For this reason, space-based GW detector missions such as
LISA and TianQin have been proposed [3, 4].

The TianQin mission is expected to deploy three identical satellites in a
circular orbit with the radius about $10^{5}$ km of the Earth, which will form
an approximately equilateral triangle constellation with an armlength of about
$1.7\times 10^{5}$ km [5, 6]. This interference armlength enables the TianQin
observatory to have a sensitive detection capability in the 0.1 mHz - 1 Hz
frequency band. There are a wide variety of important astrophysical and
cosmological sources in this frequency band, including massive binary black
holes, extreme mass ratio inspirals (EMRIs), compact galactic binaries and
stochastic GW backgrounds [7, 8, 9, 10].

Like other space-based interferometric GW detectors, TianQin suffers from
overwhelming laser phase noise due to the coupling with unequal armlengths
caused by the relative motion of the satellites. After many theoretical
studies and experimental verifications, time delay interferometry (TDI) is
considered to be an effective method for deducting laser phase noise [11, 12,
13, 14, 15, 16, 17]. As a data post-processing technique, TDI can construct
virtual equal-armlength interferometers to reduce laser phase noise. The LISA
team has developed multiple sets of TDI simulation software [18, 19, 20, 21,
22, 23], such as LISA Simulator, Synthetic LISA, LISACode, and LISANode, and
the capabilities have been improved over time. As possibly the first
simulator, LISA Simulator was based on an old interferometry layout [18].
Synthetic LISA and LISACode are more in line with the real physical scenario
[19, 20]. LISANode is a new prototype simulator which includes an up-to-date
instrumental configuration [24], various noise sources, and improved TDI
algorithms [23, 25, 26]. Important lessons can be learned from LISA.
Nevertheless, to handle unique features of TianQin, we have felt the
necessities to develop a TDI simulation program specifically for the mission.
Some earlier efforts were made, for instances, in studying the science cases
of EMRIs [9] and in evaluating residual unequal armlengths after TDI [27],
and, more recently, in assessing the space plasma noise propagation in TDI
[28]. The new MATLAB-based program is named TQTDI and has been developed in
tandem with our orbit simulator TQPOP (TianQin Quadruple Precision Orbit
Propagator [29]).

TDI is currently a necessary method for deducting laser phase noise that is
several orders of magnitude higher than GW signals in space GW detection, but
there are other important effects to be reckoned with in TDI data processing.
Because of phase locking offsets and the Doppler shift of the laser frequency
due to relative motion between satellites, the beatnote frequencies can vary
by $\sim 10$ MHz (The bandwidth for phase measurement is about 5 to 25 MHz
[30, 31]). In the raw output of phasemeters, small GW-induced phase
fluctuations are buried under large phase ramps, but most of previous TDI
simulations in literature seem to have disregarded the Doppler shift in
generating the science beatnote signals. Though the Doppler effect is
considered out-of-band especially for LISA, the potential impact on the
effectiveness of TDI deserves careful studies. Additionally, the initial noise
reduction pipeline should take the issue into account in order to handle
realistic data, since it only has access to total phases or frequencies [25,
32].

Proximity of the TianQin satellites to the Earth has its benefits, such as
easier data communication, availability of the Global Navigation Satellite
Systems, etc. However, the complex gravity environment of the Earth-Moon
system also sets a lower bound at $1\times 10^{-4}$ Hz to the detection
frequency band of TianQin [29, 33], which poses a risk of interfering with TDI
processing.

In order to explore whether or not the Doppler effect affects TianQin’s TDI
processing, we need to compare the results of TDI before and after deducting
the Doppler shift. Due to the huge magnitude difference between GW signals and
the Doppler effect when expressed in phase, double-precision (16 digits) can
no longer meet the precision requirement. As a remedy, our simulator has been
developed using quadruple-precision arithmetic (34 digits) at necessary
places, which is a technique we have previously applied for TQPOP [29].
Because the orbital Doppler effect is very close to the detection frequency
band of TianQin [29, 33], the method of deducting it also needs to be
carefully designed and tested to ensure that it does not affect TDI
performance and GW detection.

Ground-based GW interferometers already in operation, such as advanced LIGO,
also suffer from various environmental noises, among which seismic and
Newtonian noises in the raw data is rather high at low frequencies below 10 Hz
[34, 35, 36]. In order to facilitate subsequent scientific data analysis, one
approach adopted by advanced LIGO is using high-pass filter in gstlal
calibration pipeline to remove the huge low frequency noise [37, 38].
TianQin’s case is somewhat similar but with the complication of TDI. Hence we
may also take a similar approach in the pre-data processing for TianQin to
remove this orbital noise.

The paper is organized as follows. In Sec. II we introduce the cause of the
Doppler effect and deduce the TDI formulas with the effect. In Sec. III, first
we analyze the characteristics of the orbital Doppler effect of the TianQin
satellites, which shows the feasibility of using a high-pass filter to deduct
the effect. Then we briefly illustrate the design method of the filter used in
our research and the simulation tool TQTDI we developed is described. Section
IV shows the simulation results, including the data processing order and the
effect of ranging errors, and theoretical deduction of the influence of
Doppler effect on TDI processing is also included. Last, conclusions are
presented in Sec. V.

##  II  Theoretical description

The TianQin constellation is shown in Fig. 1, which includes three spacecrafts
that are labeled SC$i$ (For $i$ = 1, 2, 3, same below). The lengths of laser
links between satellites are denoted as $L_{i}$ and $L_{i^{\prime}}$,
respectively, according to the counter-clockwise and clockwise propagation
directions. Each of the three spacecrafts carries two optical benches (OB)
marked as $i$ and $i^{\prime}$. In the split measurement scheme [22], each OB
will output three sets of measurement data from the science interferometer,
test-mass (TM) interferometer and reference interferometer (The derivation in
this paper neglects clock sideband measurements). The measurement data are
denoted as $s_{i}$, $\epsilon_{i}$, and $\tau_{i}$, respectively, which
contain both GW signals and noises. TDI algorithm is to linearly combine the
measurements and their delayed terms to suppress the laser phase noise and
retain the GW signals. There are a variety of TDI combinations, such as
Michelson combinations X (1st generation), X2 (2nd generation), Sagnac
combinations $\alpha$, $\beta$, and so on [12, 39]. The detailed composition
of $s_{i}$, $\epsilon_{i}$, and $\tau_{i}$ has been deduced in the previous
literature [22], and for simplicity, here we only include a few crucial
physical quantities: GW signals, laser phase noise and Doppler frequency shift
in the unit of phase for the derivation of the TDI formula.

Figure 1:  Labeling conventions used in TianQin split-interferometer.

Without loss of generality, we take OB1 as an example to derive the TDI
combinations with the Doppler effect. The total phase of a photon originating
from the laser source of OB1 at the time $t$ is given as

$\Phi_{1}(t)=\omega_{1}t+p_{1}(t),$ (1)

where $\omega_{1}$ is the center laser frequency and $p_{1}$ is the phase
fluctuation of the pre-stabilized laser. The total phase of the laser emitted
from SC2 when it reaches SC1 at the time $t$ can be expressed as

$\Phi_{2}^{\prime}(t)=\omega_{2^{\prime}}^{\prime}t+p_{2^{\prime}}(t-\tau_{3})+H_{1}(t),$
(2)

where $\omega_{2^{\prime}}^{\prime}$ is the center frequency of the laser on
OB2’ with the Doppler shift, $p_{2^{\prime}}\left(t-\tau_{3}\right)$ is the
phase fluctuation of the laser on OB2’ with the delay along the arm $L_{3}$,
and $H_{1}(t)$ is the phase changes caused by GWs. Hence, the phase difference
formed by the interference of the two laser beams on OB1 can be expressed as

$\Delta\Phi_{1}(t)=(\omega_{2^{\prime}}^{\prime}-\omega_{1})t+p_{2^{\prime}}(t-\tau_{3})-p_{1}(t)+H_{1}(t).$
(3)

In Eq. (3), we regard the term $(\omega_{2^{\prime}}^{\prime}-\omega_{1})t$ as
being solely caused by the Doppler effect, and assume that the fixed frequency
difference introduced during phase locking has been removed. In reality, the
beatnote signal measured by the phasemeter contains this Doppler term
($(\omega_{2^{\prime}}^{\prime}-\omega_{1})t$), and it might affect the GW
detection due to the complex gravitational field environment of the geocentric
orbit. Therefore, we will add this term in the following derivation, and focus
on examining it in this work.

By introducing the delay operator $D_{i}f(t)=f(t-\tau_{i})$, the science
interferometer signals of the two optical benches in SC1 can be written as

$\begin{split}s_{1}&=\phi_{1}^{d}+D_{3}p_{2^{\prime}}-p_{1}+H_{1},\\\
s_{1^{\prime}}&=\phi_{1^{\prime}}^{d}+D_{2^{\prime}}p_{3}-p_{1^{\prime}}+H_{1^{\prime}},\end{split}$
(4)

where $\phi_{1}^{d}=\int(\omega_{2^{\prime}}^{\prime}-\omega_{1})dt$ and
$\phi_{1^{\prime}}^{d}=\int(\omega_{3}^{\prime}-\omega_{1^{\prime}})dt$
represent accumulated Doppler phases. The multiplication is converted to an
integral here because the frequency of the received laser beam sent by the
remote satellite changes due to the relative movement between the satellites.
Note that in this short derivation alone we ignore OB displacement noise,
clock noise, readout noise, and other noises here [22], and only consider
laser phase noise, the Doppler effect, and GW signals. In this case, the TM
and reference interferometer signals become trivial. All other interferometer
signals from SC2 and SC3 can be constructed by cyclic permutation of the
unprimed and primed indices: 1 → 2, 2 → 3 and 3 → 1 as well as 1’ → 2’, 2’ →
3’ and 3’ → 1’.

Then, we follow the standard procedure to construct the intermediate TDI
variables $\eta_{1}$ and $\eta_{1^{\prime}}$ to remove the frequency
fluctuations of the primed lasers [22, 21]:

$\begin{split}\eta_{1}&=H_{1}+\phi_{1}^{d}+D_{3}p_{2}-p_{1},\\\
\eta_{1^{\prime}}&=H_{1^{\prime}}+\phi_{1^{\prime}}^{d}+D_{2}p_{3}-p_{1}.\end{split}$
(5)

Similarly, the other four intermediate variables can be obtained by cyclic
permutation of the indices. With these variables, we can construct various TDI
combinations to remove the laser noises, and the following is the expression
for the Michelson-X combination [17]:

$\begin{split}X=&[D_{2^{\prime}}D_{2}-I](\eta_{1}+D_{3}\eta_{2^{\prime}})\\\
-&[D_{3}D_{3^{\prime}}-I](\eta_{1^{\prime}}+D_{2^{\prime}}\eta_{3}),\end{split}$
(6)

where $I$ is the identity operator. The above equation, like the derivation in
previous literature [17], will deduct most of the laser noise, leaving only
the residual term:

$\delta
X^{p}=D_{3^{\prime}}D_{3}D_{2}D_{2^{\prime}}p_{1}-D_{2}D_{2^{\prime}}D_{3^{\prime}}D_{3}p_{1}.$
(7)

One can see from Eq. (4) to Eq. (6) that the TDI combination works in the same
way for both $H_{i}$ and $\phi_{i}^{d}$, and does not normally have an
additional suppression on the latter. Therefore, we need to analyze and
evaluate TianQin’s Doppler effect jointly with TDI to determine whether or not
it affects GW detection, and if so, how to deduct it.

##  III  Simulation setup

Figure 2 shows a preliminary TianQin data transmission chain. In the
satellites, the phasemeter outputs are downsampled to about 10 Hz, and
transmitted to the Earth. On the ground, the raw data are preprocessed to
reduce the instrument noises through clock synchronization, TDI, and other
means [32]. Then the resulting data products are analysed to extract GW
signals for science application. Our TQTDI program simulates the steps after
the downsampling and before the scientific data processing.

Figure 2:  Preliminary TianQin data transmission chain.

###  III.1 Orbital Doppler effect

The Doppler frequency shift is a main component of the phasemeter measurement
data. TianQin’s floating test masses are subject to complex gravitational
perturbations from the Earth, Moon, and other celestial bodies in the solar
system. The perturbations are the primary cause of the orbital Doppler effect.
In previous work, our team have determined the amplitudes and frequencies of
the gravity disturbances for TianQin’s orbit, and evaluated its impact on the
range acceleration noise [29, 33]. The result shows that the Earth-Moon’s
gravity field dominates at low frequencies and is out of the preliminary
detection band, which makes it easier to remove the Doppler effect. But the
evaluation has also shown the insufficiency of double-precision arithmetic for
the numerical simulation because of the enormous dynamic range of the Doppler
effect that is about 20 orders of magnitude higher than the noise requirement.
To overcome this difficulty, we have developed TDI simulation program
supporting quadruple-precision arithmetic in necessary modules, which allows
for 34 significant digits in representing simulated data.

TQPOP is a quadruple precision orbit simulation program for TianQin developed
by our team [29], which provides a relative truncation error of $<10^{-20}$
and sub-pm/Hz1/2 precision. Various sources of gravitational perturbations are
included in the simulator, such as the Earth and Moon’s static gravity, the
Earth’s solid and ocean tides, the Sun’s point mass and oblateness $J_{2}$ and
so on. Using orbits of the TianQin satellites calculated by TQPOP, the Doppler
effect in terms of phase variations can be obtained, and the results are shown
in Fig. 3. One can see that the amplitude spectral density (ASD) is mainly
distributed below $1\times 10^{-4}$ Hz. It suggests that a straightforward way
to remove the Doppler effect is by high-pass filtering.

Figure 3:  Time and frequency domain diagrams of the orbital Doppler effect of
the TianQin satellites (90 days). Above: The phase variations associated with
the Doppler effect generated for the six laser links, cf. Fig. 1 for
notations. Bottom: The amplitude spectral density (ASD) of the Doppler effect
and the noise requirement on a single link.

###  III.2 High-pass filter design

A high level of linearity of the phase response [16] is needed to meet the
stringent requirement for data fidelity in GW science data processing. Thereby
a FIR (Finite impulse response) filter that can meet the linear phase
requirement is adopted in this work. In order to completely filter out the
Doppler effect without affecting GW signals, the FIR filter needs to have the
characteristics of large stopband attenuation, narrow transition width, and
small peak-to-peak ripple in the passband, which makes the design more
difficult. Based on the Parks-McClellan optimal equiripple algorithm [40, 41],
we develop a quadruple-precision high-pass filter program with MATLAB.
Applying the filter multiple times in series can remove the Doppler effect by
more than ten orders of magnitude while preserving GW signals.

###  III.3 Program modules

The simulation flow of TQTDI is shown in Fig. 4, it can be roughly divided
into two parts: raw data generation and preprocessing including TDI and the
Doppler effect removal. Two possible preprocessing sequences, i.e., applying
the filter before or after TDI combinations, are considered. The raw data,
namely the beatnote signals, are dominated by the Doppler effect, but also
contains GW signals and various noises. For data generation, one needs to know
the satellite orbits and GW sources. The delay time required in the TDI
processing can be estimated from the orbits [21].

Figure 4:  Simulation modules of TQTDI. Noises: Generating laser phase noises,
TM displacement noises, readout noises, and more. GW sources: Providing GW
strain based on formulas or file reading, and source information. Orbit:
Generating satellite orbit data, calculated by Kepler’s formula or importing
quadruple-precision orbit data from TQPOP. GW response: Calculating TianQin’s
responses from satellite orbits and GW sources. Doppler: Calculating the
orbital Doppler effect. Ranging estimation: Calculating the laser travel time
between two satellites based on the orbits (general relativistic effects
ignored), with ranging errors if needed. Beatnote signals: Synthesizing
science interferometer signals, TM interferometer signals, and reference
interferometer signals of six OBs. TDI: Performing TDI processing on the
beatnote signals according to the estimated delay time to cancel the laser
phase noise. Filter: Filtering out orbital Doppler effects by a high-pass
filter. The order of filtering and TDI may swap.

##  IV  Simulation results

Using TQTDI, we simulate quadruple-precision raw data with a low sampling rate
of 0.02 Hz and a duration of about 90 days. This is chosen partly to save
computation time. Then the laser phase noise and Doppler effect are removed by
TDI processing and high-pass filtering, respectively. To demonstrate the
effectiveness of the preprocessing, a double white dwarf GW signal, with a
strain $h=10^{-18}$ and frequency $f_{0}$ = 1 mHz, is selected as a test
signal in the simulation. Noise parameters are acquired from [4], e.g., 10
Hz/Hz1/2 for pre-stabilized lasers. Phase-locking, which induces a linear
growth trend in the output phase data, has not been added to our simulation,
so the detrending is omitted. The results are discussed below.

###  IV.1 TDI and filtering without ranging error

The simulation of this subsection is the ideal situation, that is, one has the
perfect information of the absolute ranging of each two satellites. Figure 5
shows the results of TDI processing first and then high-pass filtering. The
science interferometer signal (blue curve) shows that the laser phase noise
dominates at higher frequencies while the Doppler effect dominates at lower
frequencies. Apparently, the double white dwarf GW test signal is completely
overwhelmed. After TDI processing, the laser noise is suppressed by about ten
orders of magnitude, and the test GW signal appears, as shown by the red curve
in the figure. The residual noise at higher frequencies is consistent with the
theoretical value of the residual secondary noise after TDI-X [42], indicating
that the TDI-X combination can suppress the laser phase noise below the
secondary noise at least in the frequency band of 0.1–1 mHz [27]. Note that
the expected violation of the noise requirement above 1 mHz [27] is not
clearly reflected in the 90-day ASD due to non-stationarity of the laser phase
noise that is proportional to time-varying residual unequal armlengths [11].

Figure 5:  Phase ASD of the simulation results of 90 days. The blue curve is
the science interferometer signal of OB1. The red and yellow curves are TDI-X
and TDI-X after high-pass filtering, respectively. The black dashed line is
the TDI-X residual secondary noise requirement of TianQin. The insert is a
zoomed-in plot around $10^{-3}$ Hz and the pick is the test GW signal.

Because of the huge dynamic range of the Doppler effect, quadruple-precision
data are used throughout the above simulation, which consume more memory and
time than double-precision. In addition, the magnitude of the Doppler effect
is much larger than GW signals, which may affect subsequent scientific
analysis if not removed. From the first and middle graphs in Fig. 6, one can
only see the characteristics of the orbital motion in the time domain, but no
obvious traces of various noises or GW signals. So the output data of TDI
processing is high-pass filtered. From the yellow curve in Fig. 5, it can be
seen that the GW signal becomes apparent and the dynamic range of the data is
reduced after filtering. Now the data can be converted into double-precision
without losing accuracy. Meanwhile, one can also see the sinusoidal GW signal
in the time domain, as shown in the bottom graph of Fig. 6.

Figure 6:  Time domain results of the simulation. The blue, red, and yellow
curves show the corresponding time series of the ASD curves in Fig. 5. The
insert of the bottom graph displays an expanded waveform.  Figure 7:  Phase
ASD of the simulation results for filtering before TDI. The blue and red lines
are the original science interferometer signal and the filtered signal,
respectively. TDI-X from the filtered data is plotted by the yellow line. The
residual secondary noise requirement is indicated by black dashed line.

The simulation results for exchanging TDI processing and filtering sequence
are shown in Fig. 7. The laser phase noise is also suppressed below the
secondary noise like the previous case in Fig. 5. Therefore, it can be
concluded that when there is perfect knowledge of the armlengths, the sequence
of TDI processing and filtering out the Doppler effect does not affect the
results.

We also simulate several different TDI combinations, and the results of the
combinations $\rm X,\ X_{2}$ and $\alpha$ are shown in Fig. 8. By comparing
the simulation and theoretical results, one can see that the second generation
Michelson combinations can suppress the laser phase noise below the secondary
noise in the tested frequency band, but the first generation Sagnac
combination cannot.

Figure 8:  Phase ASDs of different TDI combinations, namely, 1st (yellow) and
2nd (green) generation Michelson combinations and 1st generation Sagnac
$\alpha$ combination (blue). The dashed lines represent the corresponding
theoretical residual secondary noises.

###  IV.2 TDI and filtering with ranging error

Next, we consider the case with ranging error, and use imperfect knowledge of
the armlengths in the simulation. Pseudo-random codes ranging technique [24]
is normally used to determine the absolute range between satellites. Here we
assume that the ranging errors of the six links are uncorrelated zero-mean
Gaussian white noises [45, 32, 31], and simulate several scenarios with
various noise levels. By adding white noise series with different RMS (Root-
Mean-Square) to the laser travel times in TDI, multiple cases with different
ranging errors can be simulated.

Figure 9 shows the simulation results of TDI processing and then filtering,
and the ranging error of $10^{-1}\ \rm m,\ 10^{-3}\ m,\ 10^{-5}\ m,\ 10^{-7}\
m$ (RMS) are chosen. It can be seen that the green curve coincides with the
black dashed line at $>10^{-4}$ Hz, which means the ranging error needs to
reach $10^{-7}$ m to make the TDI processing achieve the expected performance.
However, in the scenarios without the Doppler effect, namely swapping the
order of TDI processing and filtering, the former still shows good suppression
even when the ranging error reaches 10 m as displayed in Fig. 10. The
discrepancy in the above simulation results implies that an increase in the
residual noise may be caused by the Doppler effect not removed in the data. To
confirm this, we have the following theoretical derivation.

Figure 9:  The simulation results of TDI processing without removing the
Doppler effect in the science interferometer signals, and with ranging errors
$10^{-1}\ \rm m,10^{-3}\ m,10^{-5}\ m,10^{-7}\ m$, respectively. The solid
blue curve is the science interferometer signal with the Doppler effect, and
other colored solid curves are the results of TDI-X with different ranging
errors. The residual secondary noise requirement is shown by the black dashed
line.

When there is delay error, $D\phi(t)$ will be $\phi(t-\tau-\delta\tau)$, where
$\delta\tau$ represents the time error associated with the ranging error. Then
performing a Taylor expansion with respect to $\delta\tau$, one can get the
following formula:

$\phi(t-\tau-\delta\tau)\approx\phi(t-\tau)-\dot{\phi}(t-\tau)\delta\tau.$ (8)

In the scenario where $\dot{\phi}(t-\tau)$ changes very slowly relative to
$\delta\tau$ as in the case of the Doppler effect, we can assume that
$\dot{\phi}(t-\tau)$ is a constant, denoted by $\dot{\phi}$. Applying the
Fourier transform, we obtain the error term:

$\tilde{\phi}_{err}={\dot{\phi}}\,\delta\tilde{\tau}.$ (9)

The above equation shows that the error emerges from the coupling of the
derivative of the measured phase $\dot{\phi}$ and the time error
$\delta\tilde{\tau}$. The main component of the interferometer signal comes
from the Doppler effect before the high-pass filtering, and from laser phase
noise after the filtering. In Fig. 7, the maximum value of the blue curve is
nearly ten orders of magnitude higher than that of the red one. That is to say
that $\dot{\phi}$ is reduced by nearly ten orders of magnitude after the
filtering, which can explain why the ranging error can be relaxed when the
Doppler effect is not present.

Using a similar method, one can deduce that for the Doppler effect in TDI-X
combination, the error term introduced due to the coupling with the ranging
error is

$\delta
X_{err}^{d}\approx-2\dot{\phi}_{1}^{d}(t)(\delta\tau_{2}+\delta\tau_{2^{\prime}})+2\dot{\phi}_{1^{\prime}}^{d}(t)(\delta\tau_{3}+\delta\tau_{3^{\prime}}).$
(10)

In Eq. (10), we assume the delay time $\tau_{i}=\tau$ and
$\dot{\phi}_{i}^{d}(t-N\tau)=\dot{\phi}_{i}^{d}(t)$ for $N=1,\ 2,\ 3$. Because
$\dot{\phi}_{i}^{d}(t)$ changes quite slowly that
$\dot{\phi}_{i}^{d}(t-N\tau)-\dot{\phi}_{i}^{d}(t)$ is five orders of
magnitude lower than $\dot{\phi}_{i}^{d}(t)$. Similarly, considering
$\dot{\phi}_{i}^{d}(t)$ as a constant, and assuming that all ranging noises
are uncorrelated and have the same frequency characteristics, then in the
Fourier domain we obtain

$\delta\tilde{X}_{err}^{d}\approx 4\dot{\phi}^{d}\delta\tilde{\tau}.$ (11)

For TianQin, $\dot{\phi}_{i}^{d}\sim 10^{7}$ Hz, and by Eq. (11), we can
estimate that the residual ranging error compatible with the Doppler effect
needs to reach about $10^{-7}$ m for TDI-X combination to be below the
residual secondary noise over the test frequency band. However, achieving such
precision is currently unfeasible with pseudo-ranging techniques, which can
achieve accuracy at the centimeter level [43, 44, 45].

Therefore, under the current technology, we recommend removing the Doppler
phase drift before performing TDI combinations in the data preprocessing in
order to alleviate the requirement on the ranging error.

Figure 10:  The simulation results of TDI with the Doppler effect filtered out
and with the ranging error is $10^{2}\ \rm m,10^{1}\ m,10^{0}\ m,10^{-1}\ m$.
The high-pass filtered science interferometer signal is plotted by the solid
blue curve, and other colored solid curves are the results of TDI-X with
different ranging errors. The black dashed line is the residual secondary
noise requirement. The insert shows an expanded view.

##  V  Conclusion

In this proof-of-principle study, based on the sub-pm/Hz1/2 precision orbits
for the TianQin satellites, we develop a quadruple-precision TDI simulation
program to account for the large-scale Doppler effect in the science beatnote
phase signals. A high-performance high-pass filter compatible with standard
TDI processing is designed to remove the Doppler phase drift at low
frequencies so that GW signals may emerge in the time domain. Our simulations
further demonstrate that TDI processing without removing the Doppler effect
(i.e., the orbital noise) beforehand would impose a much higher requirement on
the pseudo ranging error, and the explanations are given (see Sec. IV.2 and
Appendix). Hence, we recommend filtering out the Doppler phase drift before
performing TDI combinations in TianQin’s initial noise-reduction pipeline
under the existing pseudo ranging capability. With this scheme, the Earth-
Moon’s gravity disturbance at frequencies $<1\times 10^{-4}$ Hz is expected to
have little impact on the GW detection at frequencies $>1\times 10^{-4}$ Hz.

As an added note, the study has assumed the two-way pseudo-ranging for
implementing TDI as a baseline for TianQin. However, this is not the only
option, and may be complemented by, e.g., TDI ranging [46, 47, 48]. Thereby it
would be interesting to consider this latter scenario in future work.
Moreover, one might also consider finding possible usages for the low-
frequency gravity-field information discarded by the filtering.

###### Acknowledgements.

The authors thank Yi-De Jing, Chengjian Luo, Bobing Ye, Jun Luo and anonymous
referee for helpful discussions and comments. X. Z. is supported by the
National Key R&D Program of China (Grant No. 2020YFC2201202 and
2022YFC2204600).

##  Appendix A Derivation of Eqs. (10) and (11)

The result of the Doppler term after being processed by TDI-X can be obtained
by replacing the $\eta_{i}$ in Eq. (6) with $\phi^{d}_{i}$. When there are
ranging errors, it can be written as

$\begin{split}X^{d}&=\phi_{1}^{d}(t-\tau_{2}-\tau_{2^{\prime}}-\delta\tau_{2}-\delta\tau_{2^{\prime}})-\phi_{1}^{d}(t)\\\
&+\phi_{2^{\prime}}^{d}(t-\tau_{3}-\tau_{2}-\tau_{2^{\prime}}-\delta\tau_{3}-\delta\tau_{2}-\delta\tau_{2^{\prime}})\\\
&-\phi_{2^{\prime}}^{d}(t-\tau_{3}-\delta\tau_{3})\\\
&-\phi_{1^{\prime}}^{d}(t-\tau_{3^{\prime}}-\tau_{3}-\delta\tau_{3}-\delta\tau_{3^{\prime}})+\phi_{1^{\prime}}^{d}(t)\\\
&-\phi_{3}^{d}(t-\tau_{2^{\prime}}-\tau_{3^{\prime}}-\tau_{3}-\delta\tau_{2^{\prime}}-\delta\tau_{3}-\delta\tau_{3^{\prime}})\\\
&+\phi_{3}^{d}(t-\phi_{2^{\prime}}-\delta\tau_{2^{\prime}})\end{split}$ (12)

Using Taylor expansion, we get

$\begin{split}X^{d}&\sim\phi_{1}^{d}(t-\tau_{2}-\tau_{2^{\prime}})-\dot{\phi_{1}^{d}}(t-\tau_{2}-\tau_{2^{\prime}})(\delta\tau_{2}+\delta\tau_{2^{\prime}})\\\
&-\phi_{1}^{d}(t)+\phi_{2^{\prime}}^{d}(t-\tau_{3}-\tau_{2}-\tau_{2^{\prime}})\\\
&-\dot{\phi_{2^{\prime}}^{d}}(t-\tau_{3}-\tau_{2}-\tau_{2^{\prime}})(\delta\tau_{3}+\delta\tau_{2}+\delta\tau_{2^{\prime}})\\\
&-\phi_{2^{\prime}}^{d}(t-\tau_{3})+\dot{\phi_{2^{\prime}}^{d}}(t-\tau_{3})(\delta\tau_{3})\\\
&-\phi_{1^{\prime}}^{d}(t-\tau_{3^{\prime}}-\tau_{3})+\dot{\phi_{1^{\prime}}^{d}}(t-\tau_{3^{\prime}}-\tau_{3})(\delta\tau_{3}+\delta\tau_{3^{\prime}})\\\
&+\phi_{1^{\prime}}^{d}(t)-\phi_{3}^{d}(t-\tau_{2^{\prime}}-\tau_{3^{\prime}}-\tau_{3})\\\
&+\dot{\phi_{3}^{d}}(t-\tau_{2^{\prime}}-\tau_{3^{\prime}}-\tau_{3})(\delta\tau_{2^{\prime}}+\delta\tau_{3}+\delta\tau_{3^{\prime}})\\\
&+\phi_{3}^{d}(t-\tau_{2^{\prime}})-\dot{\phi_{3}^{d}}(t-\tau_{2^{\prime}})(\delta\tau_{2^{\prime}})\end{split}$
(13)

Then the error term introduced by the coupling of Doppler effect and ranging
error is:

$\begin{split}\delta
X^{d}_{err}&=-\dot{\phi_{1}^{d}}(t-\tau_{2}-\tau_{2^{\prime}})(\delta\tau_{2}+\delta\tau_{2^{\prime}})\\\
&-\dot{\phi_{2^{\prime}}^{d}}(t-\tau_{3}-\tau_{2}-\tau_{2^{\prime}})(\delta\tau_{3}+\delta\tau_{2}+\delta\tau_{2^{\prime}})\\\
&+\dot{\phi_{2^{\prime}}^{d}}(t-\tau_{3})(\delta\tau_{3})+\dot{\phi_{1^{\prime}}^{d}}(t-\tau_{3^{\prime}}-\tau_{3})(\delta\tau_{3}+\delta\tau_{3^{\prime}})\\\
&+\dot{\phi_{3}^{d}}(t-\tau_{2^{\prime}}-\tau_{3^{\prime}}-\tau_{3})(\delta\tau_{2^{\prime}}+\delta\tau_{3}+\delta\tau_{3^{\prime}})\\\
&-\dot{\phi_{3}^{d}}(t-\tau_{2^{\prime}})(\delta\tau_{2^{\prime}})\end{split}$
(14)

First we assume the delay time $\tau_{i}=\tau$ and
$\dot{\phi}_{i}^{d}(t-N\tau)=\dot{\phi}_{i}^{d}(t)$ for $N=1,\ 2,\ 3$.
Furthermore, because the Doppler effect on the same arm is approximately
equal, we have

$\begin{split}\dot{\phi}_{1}^{d}(t)\approx\dot{\phi}_{2^{\prime}}^{d}(t)\\\
\dot{\phi}_{1^{\prime}}^{d}(t)\approx\dot{\phi}_{3}^{d}(t)\end{split}$ (15)

In this way, Eq. (14) can be approximated as Eq. (10). By assuming
$\delta\tilde{\tau}_{2}=\delta\tilde{\tau}_{2^{\prime}}=\delta\tilde{\tau}_{3}=\delta\tilde{\tau}_{3^{\prime}}=\delta\tilde{\tau}$,
the Fourier transform of Eq. (10) yields Eq. (11).

##  Appendix B Verification of Eqs. (10) and (11)

Here we verify Eqs. (10) and (11) with numerical simulations. The two solid
lines in Fig. 11 are the simulation results, and the two dashed lines are the
results calculated by Eq. (10) and Eq. (11), respectively. As can be seen from
Eq. (9) and Eq. (11), the error introduced by TDI-X is 4 times that of a
single delay operator. So we simulate $4D_{2}\phi^{d}_{1}(t)$ and TDI-X
combination with a ranging error of $10^{-3}\ \rm m$, of which the ASDs are
shown by the blue and red solid lines in Fig. 11. It can be seen that their
amplitudes are the same at $>1\times 10^{-4}$ Hz, and are consistent with the
theoretical calculation of the dashed lines. The zoomed-in image below shows
that the red solid line and the yellow dashed line coincide, and therefore it
is shown that the approximations in the deriving Eqs. (10) and (11) are
reasonable.

Figure 11:  Comparison between theoretical calculation from Eqs. (10, 11) and
numerical simulation. The blue and red solid lines simulated by TQTDI are the
ASD of 4 times the Doppler term $\phi_{1}^{d}(t)$ with the delay $\tau_{2}$
and TDI-X combination, both with a ranging error of $10^{-3}$ m, respectively.
Yellow dashed line is the ASD calculated from Eq. (10). Taking the RMS of the
derivative of each Doppler term as the value of $\dot{\phi}^{d}$, and
substituting it into Eq. (11) give the green dashed line.

## References

  * Abbott  _et al._ [2016] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Observation of gravitational waves from a binary black hole merger, Physical Review Letters 116, 061102 (2016).
  * Abbott  _et al._ [2009] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. Amin, S. Anderson, W. Anderson, M. Arain, _et al._ , LIGO: the laser interferometer gravitational-wave observatory, Reports on Progress in Physics 72, 076901 (2009).
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, _et al._ , Laser interferometer space antenna, arXiv preprint arXiv:1702.00786 (2017).
  * Luo  _et al._ [2016] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, _et al._ , TianQin: a space-borne gravitational wave detector, Classical and Quantum Gravity 33, 035010 (2016).
  * Ye  _et al._ [2019] B.-B. Ye, X. Zhang, M.-Y. Zhou, Y. Wang, H.-M. Yuan, D. Gu, Y. Ding, J. Zhang, J. Mei, and J. Luo, Optimizing orbits for TianQin, International Journal of Modern Physics D 28, 1950121 (2019).
  * Tan  _et al._ [2020] Z. Tan, B. Ye, and X. Zhang, Impact of orbital orientations and radii on TianQin constellation stability, International Journal of Modern Physics D 29, 2050056 (2020).
  * Wang  _et al._ [2019] H.-T. Wang, Z. Jiang, A. Sesana, E. Barausse, S.-J. Huang, Y.-F. Wang, W.-F. Feng, Y. Wang, Y.-M. Hu, J. Mei, _et al._ , Science with the TianQin observatory: Preliminary results on massive black hole binaries, Physical Review D 100, 043003 (2019).
  * Huang  _et al._ [2020] S.-J. Huang, Y.-M. Hu, V. Korol, P.-C. Li, Z.-C. Liang, Y. Lu, H.-T. Wang, S. Yu, J. Mei, _et al._ , Science with the TianQin observatory: Preliminary results on galactic double white dwarf binaries, Physical Review D 102, 063021 (2020).
  * Fan  _et al._ [2020] H.-M. Fan, Y.-M. Hu, E. Barausse, A. Sesana, J.-d. Zhang, X. Zhang, T.-G. Zi, J. Mei, _et al._ , Science with the TianQin observatory: Preliminary result on extreme-mass-ratio inspirals, Physical Review D 102, 063016 (2020).
  * Mei  _et al._ [2021] J. Mei, Y.-Z. Bai, J. Bao, E. Barausse, L. Cai, E. Canuto, B. Cao, W.-M. Chen, Y. Chen, Y.-W. Ding, _et al._ , The TianQin project: current progress on science and technology, Progress of Theoretical and Experimental Physics 2021, 05A107 (2021).
  * Tinto and Armstrong [1999] M. Tinto and J. W. Armstrong, Cancellation of laser noise in an unequal-arm interferometer detector of gravitational radiation, Physical Review D 59, 102003 (1999).
  * Armstrong  _et al._ [1999] J. Armstrong, F. Estabrook, and M. Tinto, Time-delay interferometry for space-based gravitational wave searches, The Astrophysical Journal 527, 814 (1999).
  * Estabrook  _et al._ [2000a] F. Estabrook, M. Tinto, and J. Armstrong, Time-delay analysis of LISA gravitational wave data: Elimination of spacecraft motion effects, Physical Review D 62, 042002 (2000a).
  * Tinto  _et al._ [2003] M. Tinto, D. A. Shaddock, J. Sylvestre, and J. Armstrong, Implementation of time-delay interferometry for LISA, Physical Review D 67, 122003 (2003).
  * De Vine  _et al._ [2010] G. De Vine, B. Ware, K. McKenzie, R. E. Spero, W. M. Klipstein, and D. A. Shaddock, Experimental demonstration of time-delay interferometry for the laser interferometer space antenna, Physical Review Letters 104, 211103 (2010).
  * Schwarze  _et al._ [2019] T. S. Schwarze, G. F. Barranco, D. Penkert, M. Kaufer, O. Gerberding, and G. Heinzel, Picometer-stable hexagonal optical bench to verify LISA phase extraction linearity and precision, Physical review letters 122, 081104 (2019).
  * Tinto and Dhurandhar [2021] M. Tinto and S. V. Dhurandhar, Time-delay interferometry, Living Reviews in Relativity 24, 1 (2021).
  * Rubbo  _et al._ [2004] L. J. Rubbo, N. J. Cornish, and O. Poujade, Forward modeling of space-borne gravitational wave detectors, Physical Review D 69, 082003 (2004).
  * Vallisneri [2005] M. Vallisneri, Synthetic LISA: Simulating time delay interferometry in a model lisa, Physical Review D 71, 022001 (2005).
  * Petiteau  _et al._ [2008] A. Petiteau, G. Auger, H. Halloin, O. Jeannin, E. Plagnol, S. Pireaux, T. Regimbau, and J.-Y. Vinet, LISAcode: A scientific simulator of LISA, Physical Review D 77, 023002 (2008).
  * Otto [2015] M. Otto,  _Time-delay interferometry simulations for the laser interferometer space antenna_ , Ph.D. thesis, Hannover: Gottfried Wilhelm Leibniz Universität Hannover (2015).
  * Otto  _et al._ [2012] M. Otto, G. Heinzel, and K. Danzmann, TDI and clock noise removal for the split interferometry configuration of LISA, Classical and Quantum Gravity 29, 205003 (2012).
  * Bayle [2019] 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, Université de Paris; Université Paris Diderot; Laboratoire Astroparticules … (2019).
  * Heinzel  _et al._ [2011] G. Heinzel, J. J. Esteban, S. Barke, M. Otto, Y. Wang, A. F. Garcia, and K. Danzmann, Auxiliary functions of the LISA laser link: ranging, clock noise transfer and data communication, Classical and Quantum Gravity 28, 094008 (2011).
  * Bayle  _et al._ [2021] J. B. Bayle, O. Hartwig, and M. Staab, Adapting time delay interferometry for LISA data in frequency, Physical Review D 104, 023006 (2021).
  * Hartwig and Bayle [2021] O. Hartwig and J.-B. Bayle, Clock-jitter reduction in LISA time-delay interferometry combinations, Physical Review D 103, 123027 (2021).
  * Zhou  _et al._ [2021] M.-Y. Zhou, X.-C. Hu, B. Ye, S. Hu, D.-D. Zhu, X. Zhang, W. Su, and Y. Wang, Orbital effects on time delay interferometry for TianQin, Physical Review D 103, 103026 (2021).
  * Jing  _et al._ [2022] Y.-D. Jing, L. Zheng, S. Yang, X. Zhang, L. Lu, B. Tang, and W. Su, Plasma noise in TianQin time-delay interferometry, Physical Review D 106, 082006 (2022).
  * Zhang  _et al._ [2021] X. Zhang, C. Luo, L. Jiao, B. Ye, H. Yuan, L. Cai, D. Gu, J. Mei, and J. Luo, Effect of Earth-Moon’s gravity on TianQin’s range acceleration noise, Physical Review D 103, 062001 (2021).
  * Barke [2015] S. Barke,  _Inter-spacecraft frequency distribution for future gravitational wave observatories_ , Ph.D. thesis, Gottfried Wilhelm Leibniz Universität Hannover (2015).
  * Bayle and Hartwig [2023] J.-B. Bayle and O. Hartwig, Unified model for the lisa measurements and instrument simulations, Physical Review D 107, 083019 (2023).
  * Hartwig [2021] O. Hartwig,  _Instrumental modelling and noise reduction algorithms for the Laser Interferometer Space Antenna_ , Ph.D. thesis (2021).
  * Luo and Zhang [2022] C. Luo and X. Zhang, Effect of Earth-Moon’s gravity on TianQin’s range acceleration noise. ii. impact of orbit selection, Physical Review D 105, 102007 (2022).
  * Aasi  _et al._ [2015] J. Aasi, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Advanced LIGO, Classical and quantum gravity 32, 074001 (2015).
  * Fiori  _et al._ [2021] I. Fiori, A. Effler, P. Nguyen, F. Paoletti, R. M. Schofield, and M. C. Tringali, Environmental noise in gravitational-wave interferometers, in  _Handbook of Gravitational Wave Astronomy_ (Springer, 2021) pp. 1–72.
  * Nguyen  _et al._ [2021] P. Nguyen, R. Schofield, A. Effler, C. Austin, V. Adya, M. Ball, S. Banagiri, K. Banowetz, C. Billman, C. Blair, _et al._ , Environmental noise in advanced LIGO detectors, Classical and Quantum Gravity 38, 145001 (2021).
  * Viets  _et al._ [2018] A. Viets, M. Wade, A. Urban, S. Kandhasamy, J. Betzwieser, D. A. Brown, J. Burguet-Castell, C. Cahillane, E. Goetz, K. Izumi, _et al._ , Reconstructing the calibrated strain signal in the advanced LIGO detectors, Classical and Quantum Gravity 35, 095015 (2018).
  * Cannon  _et al._ [2021] K. Cannon, S. Caudill, C. Chan, B. Cousins, J. D. Creighton, B. Ewing, H. Fong, P. Godwin, C. Hanna, S. Hooper, _et al._ , GstLAL: A software framework for gravitational wave discovery, SoftwareX 14, 100680 (2021).
  * Shaddock  _et al._ [2003] D. A. Shaddock, M. Tinto, F. B. Estabrook, and J. Armstrong, Data combinations accounting for lisa spacecraft motion, Physical Review D 68, 061303 (2003).
  * McClellan  _et al._ [1973] J. McClellan, T. Parks, and L. Rabiner, A computer program for designing optimum FIR linear phase digital filters, IEEE Transactions on Audio and Electroacoustics 21, 506 (1973).
  * Parks and Burrus [1987] T. W. Parks and C. S. Burrus,  _Digital filter design_ (Wiley-Interscience, 1987).
  * Estabrook  _et al._ [2000b] F. Estabrook, M. Tinto, and J. Armstrong, Time-delay analysis of LISA gravitational wave data: Elimination of spacecraft motion effects, Physical Review D 62, 042002 (2000b).
  * Esteban  _et al._ [2011] J. J. Esteban, A. F. García, S. Barke, A. M. Peinado, F. G. Cervantes, I. Bykov, G. Heinzel, and K. Danzmann, Experimental demonstration of weak-light laser ranging and data communication for lisa, Optics express 19, 15937 (2011).
  * Wang  _et al._ [2015] Y. Wang, G. Heinzel, and K. Danzmann, First stage of lisa data processing. ii. alternative filtering dynamic models for lisa, Physical Review D 92, 044037 (2015).
  * Xie  _et al._ [2020] X. Xie, T. Geng, Q. Zhao, Y. Lv, H. Cai, and J. Liu, Orbit and clock analysis of bds-3 satellites using inter-satellite link observations, Journal of geodesy 94, 1 (2020).
  * Tinto  _et al._ [2005] M. Tinto, M. Vallisneri, and J. Armstrong, Time-delay interferometric ranging for space-borne gravitational-wave detectors, Physical Review D 71, 041101 (2005).
  * Francis  _et al._ [2015] S. P. Francis, D. A. Shaddock, A. J. Sutton, G. De Vine, B. Ware, R. E. Spero, W. M. Klipstein, and K. McKenzie, Tone-assisted time delay interferometry on GRACE Follow-On, Physical Review D 92, 012005 (2015).
  * Page and Littenberg [2021] J. Page and T. B. Littenberg, Bayesian time delay interferometry, Physical Review D 104, 084037 (2021).