# Prediction of recurrence free survival of head and neck cancer using PET/CT
radiomics and clinical information

###### Abstract

The 5-year survival rate of Head and Neck Cancer (HNC) has not improved over
the past decade and one common cause of treatment failure is recurrence. In
this paper, we built Cox proportional hazard (CoxPH) models that predict the
recurrence free survival (RFS) of oropharyngeal HNC patients. Our models
utilise both clinical information and multimodal radiomics features extracted
from tumour regions in Computed Tomography (CT) and Positron Emission
Tomography (PET). Furthermore, we were one of the first studies to explore the
impact of segmentation accuracy on the predictive power of the extracted
radiomics features, through under- and over-segmentation study. Our models
were trained using the HEad and neCK TumOR (HECKTOR) challenge data, and the
best performing model achieved a concordance index (C-index) of 0.74 for the
model utilising clinical information and multimodal CT and PET radiomics
features, which compares favourably with the model that only used clinical
information (C-index of 0.67). Our under- and over-segmentation study confirms
that segmentation accuracy affects radiomics extraction, however, it affects
PET and CT differently.

Index Terms—  Head and neck cancer, Medical imaging, PET/CT, Multimodal,
Radiomics

##  1 Introduction

###  1.1 Head and neck cancer

Head and neck cancer (HNC) is the seventh most prevalent cancer in the world.
Approximately 660,000 people are diagnosed with the disease annually [1].
Common causes of treatment failure and death in HNC is distance metastases and
second primary cancers. Survival and treatment outcome of patients can
potentially be improved by identifying prognostic factors that indicate tumour
aggressiveness or risk of recurrence during diagnosis. In doing so, patients
can be stratified into different risk groups and, this can be taken into
account when chemotherapy regimens and radiation doses are being planned for
patients [2].

###  1.2 Medical imaging

Medical imaging such as 2-deoxy-2-[$\prescript{18}{}{\mathbf{F}}$] fluoro-D-
glucose (FDG) Positron Emission Tomography (PET) and Computed Tomography (CT)
contain large amounts of quantitative mineable information including clinical
information, such as prognostic features. Radiomics is the process of
extracting quantitative information from images, and using them with an aim to
aid clinical decision making [2, 3, 4].

###  1.3 Radiomics

Radiomics features are extracted from the segmented region of interest (e.g.
tumour) of the image [5]. These features are hypothesised to have prognostic
power and, analysis by [4] found that some radiomics features have
associations with gene-expression pattern in tumours. Furthermore, [6] have
also observed that CT and PET based radiomics have capabilities to predict the
control rate of the local tumour in HNC. These examples demonstrate a wide
range of predictive power radiomics hold. There are different types of
radiomics features, such as shape, first-order and second-order statistics.
Shape features represent information such as the volume and sphericity of the
tumour. First-order statistics describe histogram-based properties of the
image such as the skewness. Lastly, second-order statistic features provide
textural information such as tumour heterogeneity [2, 3].

###  1.4 Aims

This study aims to create multimodal model that can predict the recurrence
free survival (RFS) of patients diagnosed with HNC (located in the oropharynx
region). Our contributions are as follows. First, a novel model was built
using multimodal radiomics extracted from FDG-PET/CT and non-imaging data
(clinical information) [7]. Second, this study also explored the effects of
segmentation accuracy on the quality of the radiomics features extracted for
RFS prediction to understand the importance of segmentation accuracy in the
multimodal RFS workflow [3, 7]. We were one of the first studies to explore
the impact of segmentation accuracy on the predictive power of the extracted
multimodal radiomics features for RFS prediction.

##  2 Data and Methods

###  2.1 Data

The HEad and neCK TumOR (HECKTOR) challenge [7] provided FDG-PET and CT images
of patients diagnosed with oropharyngeal HNC. Furthermore, for each patient,
the PET images were registered with CT images. In addition to image data,
clinical information - centre ID, gender, age, weight, tobacco, alcohol
consumption, performance status, HPV status, surgery and/or chemotherapy were
provided for majority of the patients. The data included in the challenge was
built from 7 different centres with 524 subjects.

Moreover, the training dataset images were provided with the ground truth
labels for primary Gross Tumour Volumes (GTVp) and metastatic lymph nodes
(GTVn). A subset of the cases ($n=488$) contained RFS information. As this
study aims to build RFS prediction models, this subset of dataset was used for
the study. This training dataset was further split into training (85 %),
validation (7.5 %) and testing (7.5 %) sets. Further information about the
challenge dataset can be found in [7].

###  2.2 Preprocessing

Fig. 1: An exemplar (a) PET (b) CT and (c) Segmentation of a patient from
HECKTOR challenge [7]

The clinical dataset was pre-processed for survival analysis as some
participants had missing values. For binary data, (tobacco, alcohol, HPV
status and surgery), missing values were set as zero, negative status, denoted
as 0, was replaced with -1 and positive status, denoted as 1, was kept as it
is. The missing data for weight was set as 75 kg as instructed by HECKTOR 2022
and performance status and centre ID were dropped from the dataset [7, 8].

PET, CT, and their annotations (segmentation labels) were resampled so that
they have the same isotropic voxel spacing of 2 mm $\times$ 2 mm $\times$ 2
mm. During this process, Bspline and nearest neighbour interpolator were used
for PET/CT and annotations respectively. Furthermore, the field of view of
PET, CT and their annotations were different for some patients (e.g. Fig. 1).
To ensure uniformity in size, PET, CT and labels were cropped using bounding
boxes. The bounding boxes were determined by finding the field of view that
covers both PET and CT for each patient using code from [7]. Lastly, the
annotation of GTVp and GTVn were merged as a single label, so the radiomics
features from GTVp and GTVn could be extracted together.

Fig. 2: Radiomics workflow (Image adapted from [5] and method adapted from
[9])

###  2.3 Radiomics workflow

The radiomics workflow (shown in Fig. 2) typically involves, image
acquisition, image segmentation, extraction of radiomics features, feature
selection, and model fitting using the extracted features [5]. The pyradiomics
package [10] was used to extract features from the tumour regions in PET and
CT (107 features each). After feature extraction, the radiomics features of
PET and CT were compared, so that features with different values could be
relabelled as PET- and CT-specific features. Features that had the same value
(e.g. volume) were considered common features. In the next step, features were
selected using the least absolute shrinkage and selection operator (LASSO)[11]
with 8-fold cross-validation. The selected features were then fed into the Cox
proportional hazard (CoxPH) model [9, 12].

####  2.3.1 Cox proportional hazard model

Cox proportional hazard (CoxPH) models are used to model the effects of an
explanatory variable on survival. The mathematical expression of the CoxPH
model is shown in Eq. (1), and it expresses the hazard at time $t$ for an
individual that have the set of explanatory variables
$\\{\mathbf{X}|\mathbf{X}=(X_{1},X_{2},...X_{p})\\}$. $h_{0}(t)$ is the
baseline hazard function and, the coefficient ($\beta_{i}$) of the explanatory
variables are estimated by maximising the likelihood function [12]. The CoxPH
model was used in this study to model the RFS, which is the time without any
recurrence measured in days after the end of treatment [7].

Here, the dataset was used to evaluate the predictive power of the model,
which was evaluated by calculating the concordance index (C-index). The
C-index takes a value between 0 and 1, where 1 indicates a perfect concordance
[13]. Lastly, different combinations of imaging and non-imaging features were
fed into the LASSO and CoxPH analysis to find the optimal combinations of
features for the prediction model [9].

$h(t,\mathbf{X})=h_{0}(t)\exp{\left(\sum_{i=1}^{p}{\beta_{i}X_{i}}\right)}$
(1)

###  2.4 Impact of segmentation quality

As shown in Fig. 2, segmentation is one of the crucial processes in the
radiomics workflow as features will be extracted from the region of interest
(e.g. tumour area) [5]. Therefore, to study how the accuracy of segmentation
affects the quality of the radiomics features, morphological operations of the
ground truth annotations were performed to study over-segmentation (using
dilation), and under-segmentation (using erosion). By changing the extend of
dilation/erosion (defined by $r$ as the radius of structuring element in voxel
units), the resulting segmentations were used to extract the radiomics
features and fed into the LASSO and CoxPH analysis. In this study, we consider
the segmentation accuracy to be represented by the deviation
(erosion/dilation) annotations from the ground truth, while the quality of the
radiomics features is represented by the C-index value on RFS.

##  3 Results

The CoxPH regression model for predicting RFS was built using different
combinations of radiomics features and clinical information and the results
are shown in Tab. 1. The best performing model achieved a C-index of 0.74, and
this was built using both PET and CT features and clinical information. Of all
the ground truth models, PET+clinical model had the lowest performance,
C-index of 0.65. This was even lower than the clinical only model.

One example of radiomics features extracted and selected for the CoxPH
regression model was sphericity, which describes the roundness of the tumour
compared to a sphere. Another example was gray level size zone matrix small
area emphasis. This can be used to understand how much fine textures there are
within the tumour [10].

The result from the under-segmentation and over-segmentation study is shown in
Tab. 1 from the third row onward. As shown, the eroded CT+clinical model had a
comparable performance to the the ground truth CT model, which suggests that
the CT radiomics features were not affected by under-segmentation. However,
CT+clinical model was more sensitive to over-segmentation. On the other hand,
PET+clinical model performance was affected by under-segmentation and had the
best performance when annotations were over-segmented by $r$=2 (with C-index
of 0.68).

Lastly, the volume of some tumours was too small when the ground truth
annotations were eroded by $r$=2 or more, so radiomics features (both for PET
and CT) could not be extracted, therefore no results are presented.

CoxPH Regression Model | C-index (Train) | C-index (Test)
---|---|---
Clinical | 0.63 | 0.67
Ground truth (CT $+$ Clinical) | 0.77 | 0.73
Ground truth (PET $+$ Clinical) | 0.74 | 0.65
Ground truth (CT$+$PET$+$Clinical) | 0.78 | 0.74
Eroded r=1 (CT $+$ Clinical) | 0.76 | 0.74
Eroded r=1(PET $+$ Clinical) | 0.73 | 0.61
Eroded r=1(CT$+$PET$+$Clinical) | 0.77 | 0.69
Dilated r=1 (CT$+$ Clinical) | 0.75 | 0.68
Dilated r=1 (PET $+$ Clinical) | 0.72 | 0.62
Dilated r=1 (CT$+$PET$+$Clinical) | 0.76 | 0.64
Dilated r=2 (CT$+$Clinical) | 0.75 | 0.68
Dilated r=2 (PET$+$Clinical) | 0.73 | 0.68
Dilated r=2 (CT$+$PET$+$Clinical) | 0.77 | 0.67

Table 1: RFS prediction using different combination of imaging (derived from
CT and/or PET) and non-imaging (clinical) features. Ground truth, eroded
(under-segmentation) and dilated (over-segmentation) means that the radiomics
features were extracted using the ground truth annotations, eroded, and
dilated annotations, respectively.

##  4 Discussion

CoxPH model built using PET/CT and clinical features had the best performance
(C-index of 0.74). This performance was better than the clinical only model
(C-index of 0.67), which shows that additional information from radiomics
features improved the prediction accuracy of RFS. However, it is important to
note that PET+clinical model performed the worst out of all the ground truth
models. This suggests that PET radiomics features may not be contributing
substantially to the survival predictions. This can potentially be explained
by a couple of reasons.

The first potential reason is that the bin width used to discretise the images
to extract radiomics. Image discretisation is done before feature extraction
and the default bin width equal to 25 provided by pyradiomics [10] was used in
this study. However, it is important to remember that PET and CT have a
different range of measurement, and PET has a smaller range than CT, which
means that bin width could have been too large for PET. If this was the case,
we may have lost some information when extracting PET radiomics features [5].
This will be further explored and fine-tuned for future extraction process.

Another potential reason is that many of the ground truth annotations were
initially drawn on the CT of the PET/CT scan. This means that the ground truth
annotations may not accurately cover the tumours in PET [7]. Interestingly,
this may be true considering how PET radiomics model performed better when
annotations were dilated by r=2 (C-index of 0.68). This result may suggest
that the area around the tumour boundary may contain useful information for
survival prediction. This can be further supported by the fact that out of all
the PET models, the eroded model had the lowest performance (C-index of 0.61).
Normally, the area around the tumour boundary is quite blurry in PET due to
its low resolution and many existing segmentation methods often ignore the
blurred edges [14]. This result suggests that the information from PET may
contain valuable information and should be considered when developing a
segmentation algorithm for radiomics analysis.

Furthermore, it is interesting to see that unlike PET features, CT features
were not affected by under-segmentation. However, their performance dropped
for over-segmentation. These preliminary findings suggest that in order to
build a successful multimodal prediction model using radiomics features, it
may be better to have a separate segmentation for PET and CT images [15].
Further investigation should be done to understand why over- and under-
segmentation affects CT and PET features differently. One potential way to
explore this, is to find selected features that change in value due to under-
and over- segmentation. This can also give us a better insight into
understanding which features are more sensitive to segmentation accuracy [5].
Overall, the over- and under-segmentation study showed that quality of the
radiomics features are affected by the segmentation accuracy.

##  5 Conclusion

This study highlights two promising key results. Firstly, it has shown that
the inclusion of multimodal radiomics features improves the efficacy of the
CoxPH model in predicting RFS. Secondly, it has shown that segmentation
accuracy does affect the quality of the radiomics features extracted. For the
first result, it will be interesting to see if the addition of radiomics
features - handcrafted (e.g. number of tumour masses) [16] and deep learning
features [9] will help improve the model further. Moreover, it will be
intriguing to build a model using only PET or CT radiomics, so that it can be
used understand how well image data is able to predict RFS. For the second
result, the study has shown that the impact of over- and under-segmentation
differs depending on whether the features are extracted from PET or CT.
Further study should be done to understand if having a separate segmentation
algorithm for PET and CT is better than having one segmentation for both. This
has important consequences for multi-modal imaging, as it highlights the
importance of tailoring segmentation algorithms depending on the type of
medical images. Lastly, further study needs to be done to improve our
understanding on how exactly segmentation accuracy affects the radiomics
features. Introducing errors that reflect the inter- and intra-observer
segmentation variability seen in real settings will help us give a better
insight into the implications of segmentation accuracy in radiomics studies
[5]. This can be implemented by introducing wide range of segmentation errors,
such as inclusion of reader’s bias and omission errors [17].

##  6 Compliance with ethical standards

This research study was conducted retrospectively using human subject data
made available in open access by (HECKTOR2022) [7]. Ethical approval was not
required as confirmed by the license attached with the open access data.

##  7 Acknowledgments

This work was supported by EPSRC grant number EP/S0240-93/1, GE HealthCare and
Centre for Doctoral Training in Sustainable Approaches to Biomedical Science:
Responsible and Reproducible Research (SABS : $R^{3}$) Doctoral Training
Centre, University of Oxford.

## References

  * [1] M. Gormley, G. Creaney, A. Schache, K. Ingarfield, , and D. I. Conway,  “Reviewing the epidemiology of head and neck cancer: definitions, trends and risk factors,”  Br. Dent. J., vol. 233, no. 9, pp. 780–786, 2022.
  * [2] M. Vallieres et al.,  “Radiomics strategies for risk assessment of tumour failure in head-and-neck cancer,”  Sci. Rep., vol. 7, 2017.
  * [3] S. Rizzo, F. Botta, S. Raimondi, D. Origgi, C. Fanciullo, A. G. Morganti, and M. Bellomi,  “Radiomics: the facts and the challenges of image analysis,”  Eur. Radiol. Exp., vol. 2, pp. 1–8, 2018.
  * [4] Hugo. J.W.L. Aerts et al.,  “Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach,”  Nat. Commun, vol. 5, 2014.
  * [5] J.E. van Timmeren, D. Cester, S. Tanadini-Lang, H. Alkadhi, and B. Baessler,  “Radiomics in medical imaging-”how-to” guide and critical reflection,”  Insights into imaging, vol. 11, pp. 1–16, 2020.
  * [6] M. Bogowicz, O. Riesterer, L. S. Stark, G. Studer, J. Unkelbach, M. Guckenberger, and S. Tanadini-Lang,  “Comparison of pet and ct radiomics for prediction of local tumor control in head and neck squamous cell carcinoma,”  Acta oncologica, vol. 56, no. 11, pp. 1531–1536, 2017.
  * [7] V. Andrearczyk et al.,  “Overview of the hecktor challenge at miccai 2022: Automatic head and neck tumor segmentation and outcome prediction in pet/ct,”  in Head and Neck Tumor Segmentation and Outcome Prediction in third challenge, HECKTOR 2022, V. Andrearczyk, V. Oreiller, M. Hatt, and A. Depeursinge, Eds. 2023, LNCS, pp. 1–30, Springer Nature Switzerland.
  * [8] K. Wang et al.,  “Recurrence-free survival prediction under the guidance of automatic gross tumor volume segmentation for head and neck cancers,”  in Head and Neck Tumor Segmentation and Outcome Prediction, third challenge, HECKTOR 2022, V. Andrearczyk et al., Eds. 2023, LNCS, pp. 144–153, Springer Nature Switzerland.
  * [9] E. Bourigault, D. R. McGowan, A. Mehranian, and B. W. Papież,  “Multimodal pet/ct tumour segmentation and prediction of progression-free survival using a full-scale unet with attention,”  in Head and Neck Tumor Segmentation and outcome prediction, second challenge, HECKTOR 2021, pp. 189–201. Springer, 2022.
  * [10] J. J. M. van Griethuysen et al.,  “Computational radiomics system to decode the radiographic phenotype,”  Cancer Res., vol. 77, no. 21, pp. e104–e107, 2017.
  * [11] R. Tibshirani,  “Regression shrinkage and selection via the lasso,”  J. R. Stat. Soc. B (Methodol), vol. 58, no. 1, pp. 267–288, 1996.
  * [12] D. G. Kleinbaum and M. Klein,  Survival analysis a self-learning text,  Springer, third edition, 1996.
  * [13] Cameron Davidson-Pilon,  “lifelines: survival analysis in python,”  Journal of Open Source Software, vol. 4, no. 40, pp. 1317, 2019.
  * [14] L. Li, W. Lu, Y. Tan, and S. Tan,  “Variational pet/ct tumor co-segmentation integrated with pet restoration,”  IEEE trans radiat plasma med sci, vol. 4, pp. 37–49, 2019.
  * [15] S. Steyaert, M. Pizurica, D. Nagaraj, P. Khandelwal, T. Hernandez-Boussard, A. J. Gentles, and O. Gevaert,  “Multimodal data fusion for cancer biomarker discovery with deep learning,”  Nat Mach Intell, vol. 5, pp. 351–362, 2023.
  * [16] L. Rebaud, T. Escobar, F. Khalid, K. Girum, and I. Buvat,  “Simplicity is all you need: out-of-the-box nnunet followed by binary-weighted radiomic model for segmentation and outcome prediction in head and neck pet/ct,”  in Head and neck tumor segmentation and outcome prediction in third challenge, HECKTOR 2022, pp. 121–134. Springer, 2023.
  * [17] Ş. Vădineanu, D. M. Pelt, O. Dzyubachyk, and K. J. Batenburg,  “An analysis of the impact of annotation errors on the accuracy of deep learning for cell segmentation,”  in Medical Imaging with Deep Learning, 2022.


resizeenv

# A quantum-assisted master clock in the sky: global synchronization from
satellites at sub-nanosecond precision

Sage Ducoing<EMAIL_ADDRESS>Department of Physics and Astronomy, Louisiana
State University, Baton Rouge, LA 70803, USA  Ivan Agullo  Department of
Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803, USA
James E. Troupe  Xairos Systems, Inc., Denver, CO 80124, USA  Stav Haldar
<EMAIL_ADDRESS>Department of Physics and Astronomy, Louisiana State
University, Baton Rouge, LA 70803, USA

###### Abstract

This article develops a protocol to synchronize clocks on board a network of
satellites equipped with quantum resources. We show that, in such a
constellation, satellites reinforce each other’s sync capabilities, forming a
common clock that is more stable and precise than its constituents. We
envision the resulting network as a master clock able to distribute time
across the globe, providing the basis for a future quantum global navigation
satellite system or a space-based quantum network. As an example of its
capabilities, we show that a constellation of 50 satellites equipped with
modest quantum resources, and distributed amongst 5 orbits at an altitude of
500 km, allows synchronization of clocks spread across the globe at sub-
nanosecond precision.

##  I Introduction

The limits of science and technology depend to a large extent on our ability
to measure, hold and distribute time at high precision. Space-based platforms
equipped with quantum technologies have been envisaged as promising candidates
to revolutionize the scale, precision and security of sync networks Khatri
_et al._ (2021); Sidhu  _et al._ (2021); California-Polytechnic-State
University (Last accessed: Dec. 26 2021).

One concrete protocol for time distribution from satellites using quantum
resources was recently introduced in Ref. Haldar  _et al._ (2023a), based on
previous work Ho  _et al._ (2009); Lamas-Linares and Troupe (2018); Lee  _et
al._ (2019). In this protocol, a satellite acts as an intermediary to
synchronize clocks located at two ground stations on Earth’s surface. This
protocol utilizes the interchange of pairs of entangled photons, and takes
advantage of the tight time-of-birth correlations ($\lesssim 10$ femtoseconds)
associated with photon pairs produced via spontaneous parametric down
conversion. We refer to this as a Quantum Clock Synchronization (QCS)
protocol. This QCS protocol between stationary clocks was experimentally
demonstrated in Ref. Lee  _et al._ (2019) first and more recently in Ref.
Spiess  _et al._ (2021).

This article proposes the use of a network of satellites interconnected and
synchronized by means of an extension of this protocol, as a promising
infrastructure to distribute time globally across the Earth. We support the
proposal on an exhaustive study of the capabilities of such a network, based
on computer simulations able to recreate the dynamics of the network in real
time. The analysis shows that such a constellation of inter-communicated
satellites acquires collective advantages that enhanced the ability of
individual satellites to distribute time to clientele at the Earth surface.
The advantages concern both precision and area coverage.

The utilization of satellites as intermediaries benefits from the low loss
communication channel provided by free-space Khatri  _et al._ (2021); Sidhu
_et al._ (2021); Oi  _et al._ (2017); California-Polytechnic-State University
(Last accessed: Dec. 26 2021); Villar  _et al._ (2020); Yin  _et al._ (2020);
Dai  _et al._ (2020). This contrasts with the large losses of traditional
fiber-optic methods, which make them inefficient for communication beyond a
few 100 kms Sidhu  _et al._ (2021); Khatri  _et al._ (2021).

But, shifting to a satellite comes with its own difficulties. On the one hand,
satellites are in relative motion among each other and with ground stations,
and the relative velocities affect the achievable precision of the protocol.
On the other hand, it is impractical to put highly stable, large atomic clocks
onboard satellites, limiting the sync capabilities of each individual
satellite. Our study shows that the collective behaviour of multiple
satellites brings in advantages which compensate these limitations.

We conclude that any city on the Earth’s surface can sync continuously
throughout the day with the “master clock” that the network of satellites
defines, establishing synchronization among ground station at the sub-
nanosecond precision. Such a precision at global scales is currently beyond
the capabilities of classical protocols, including the GPS, which provides
$\gtrsim 40$ ns of precision (95% of the time) O’Connor  _et al._ (2019); Zhu
and Groten . Furthermore, the interchange of entangled photons can be used to
provide an added layer of security through Bell’s inequality violating
polarization correlations Lamas-Linares and Troupe (2018); Ho  _et al._
(2009).

Our quantitative analysis is done assuming off-the-shelf low size, weight and
power (SWaP) clocks onboard satellites, aided by already available, or near-
term quantum devices, such as single photon detectors and entanglement sources
(spontaneous parametric down conversion, SPDC, sources). On the other hand, we
introduce some approximations and simplifications in the analysis. Although
the loss model, described in Section III, provides a realistic approximation
for our purposes, it ignores the effects of cloud coverage, variable noise
levels during the day and night, pointing errors, and atmospheric turbulence.
Additionally, we approximate the satellite orbits as circular. Nonetheless,
our simulation techniques are flexible enough to include these effects in an
extended analyses, which we leave for future work. We do not include
relativistic corrections or corrections arising from frequency drifts between
clocks arising from other sources. We justify these simplifications by
noticing that the optimal duration of synchronization event turns of to be of
the order of $\lesssim 1\rm ms$, an interval that is short enough for the
accumulated effect from first order relativistic corrections and other
frequency drifts to be well below the sync precision.

For a broader perspective on clock synchronization and its role in
technological and scientific applications we point the reader to Ref. Troupe
_et al._ (2022) by some of us. A brief review of the precision levels and
network scales achievable by other state-of-the art clock synchronization
techniques McGrew  _et al._ (2018); whi ; Delva  _et al._ (2012); Deschênes
_et al._ (2016); Bergeron  _et al._ (2019) is also provided in Refs. Troupe
_et al._ (2022); Haldar  _et al._ (2023a, b). We also point the readers to
other proposals for QCS, i.e., the use of quantum resources to synchronize
clocks in Refs. Yurtsever and Dowling (2002); Giovannetti  _et al._ (2001,
2004); Valencia  _et al._ (2004a); Jozsa  _et al._ (2000); Ilo-Okeke  _et al._
(2020); Lafler and Lanning (2022). Proposals of quantum network of clocks made
in Refs. Kómár  _et al._ (2014) and Nichol  _et al._ (2022) are particularly
relevant.

The rest of this paper is organized as follows. Section II summarizes the
basics of the concrete QCS protocol we consider and its application to moving
clocks described in Lamas-Linares and Troupe (2018); Ho  _et al._ (2009);
Haldar  _et al._ (2023a, b). These tools form the basis on which the rest of
this article will rest. In Section III we briefly review the simulation
methods we use. In Section IV we discuss the simplest case of a network with
inter-satellite sync capabilities, namely a ring of synchronized satellites,
and we analyze its advantages and quantify its capabilities through relevant
figures of merit. This simple configuration will form the backbone of our
proposed master clock. In Section V, we increase the complexity of the network
by adding satellites in different orbits, and re-evaluate figures of merit to
quantify the advantage of this new functionality. Section VI is devoted to
provide a test case for synchronizing a global network of ground stations, by
choosing six cities across the globe. Section VII contains a summary and a
discussion of our results, identifying some key benefits of adding inter-
satellite QCS, both at the level of the single orbit and between orbits
(master clock).

##  II The QCS protocol and the concept of satellite precision shadow

###  II.1 The protocol

We start by briefly describing the QCS protocol for two stationary clocks, as
first introduced in Ho  _et al._ (2009), and show how their offset can be
determined using photon time-stamp cross-correlation functions. We then review
how to extend the protocol to moving clocks, and how to evaluate sync
precision and network scale.

We consider two clocks, Alice’s (A) and Bob’s (B), which we initially assume
tick at the the same frequency. The task of the protocol is to evaluate the
constant time offset between A and B. We must thus assume the clocks to be
stable enough to maintain their frequencies up to a required precision while
the protocol is executed. Suppose the constant offset between A and B is equal
to $\delta$, with B ahead of A. Both Alice and Bob are equipped with SPDC
sources that generate pairs of entangled photons, and two single photon
detectors to detect photons generated locally and by the other party,
respectively. $t_{i}^{(A)}$ and $t_{j}^{(B)}$ denote times at which photons
detections events $i$ and $j$ occur locally at Alice’s and Bob’s clocks
respectively. Let $\Delta t_{AB}$ and $\Delta t_{BA}$ be the travel time for
light between A and B, and vice-versa, respectively. The round trip time is
therefore $\Delta t=\Delta t_{AB}+\Delta t_{BA}$. The protocol described here
rests on the assumption of reciprocity, namely the equality between $\Delta
t_{AB}$ and $\Delta t_{BA}$. In classical protocols using radio frequency
signals, and for the case where Alice and Bob are a ground station and a
satellite in a LEO (low Earth orbit), respectively, the fluctuations in signal
multi-path propagation can create channel non-reciprocity of the order a
nanosecond. In contrast, for optical photons this figure goes down to the
order of 10’s of femtoseconds Belmonte  _et al._ (2017); Taylor  _et al._
(2020), which justifies our assumption. Now, to obtain the absolute time
difference between clocks, $\delta$, consider first a photon pair produced at
Alice’s location. One of the photons of the pair is detected locally at
Alice’s detector and the other member of the pair travels to Bob, accumulating
a travel time $\Delta t_{AB}$. For any particular pair event produced at
Alice’s site, the difference between the time labels recorded by Alice and Bob
will be $t^{(A)}_{i}-t^{(B)}_{i}=\Delta t_{AB}+\delta$. Similarly, for any
pair produced at Bob’s site $t^{(B)}_{j}-t^{(A)}_{j}=\Delta t_{BA}-\delta$.
The value of $\delta$ can be obtained from the knowledge of these two time
differences. In practice, however, one finds the added difficulty that Alice
and Bob receive many photons from ambient noise, and need a way to identify
which photons belong to an entangled pair (several photon pairs lose one
partner, therefore, also act like ambient noise). Correlation functions are a
useful tool to disentangle the signal from the noise. Consider, first, pairs
produced at Alice’s location. The detection times can be translated into
distributions by

$\displaystyle a(t)=\sum_{i}\delta\left(t-t^{(A)}_{i}\right)\,,\ \ \ \
b(t)=\sum_{j}\delta\left(t-t^{(B)}_{j}\right),$ (1)

where $\delta(x)$ is Dirac’s delta, and $i$ and $j$ index arbitrary detection
events, which can arise either from photons belonging to entangled pairs or
from other detector triggers such as stray light, dark counts, etc. We now
compute cross-correlation as:

$c_{AB}(T)=(a\star b)(T)=\int_{0}^{t_{\rm Aqc}}a(t)\,b(t+T)\,\mathrm{d}t\,,$
(2)

where $t_{\rm Aqc}$ is the duration of the protocol. The function $c_{AB}(T)$
counts the total number of detection events at Alice and Bob that are
separated precisely by an interval $T$. For a sufficiently large number of
detected entangled pairs (large signal-to-noise), this correlation function
would show a peak at $T=T_{AB}=\Delta t_{AB}+\delta\rm{.}$ Likewise, if we
consider pairs created on Bob’s site, we can construct another cross-
correlation

$c_{BA}(T)=(b\star a)(T)=\int_{0}^{t_{\rm Aqc}}b(t)\,a(t+T)\,\mathrm{d}t,$ (3)

which will show a peak at $T=T_{BA}=\Delta t_{BA}-\delta\rm{.}$ From these, we
can extract both the round trip time $\Delta t$ and the absolute time
difference between clocks without any knowledge of the length of the path
between Alice and Bob:

$\displaystyle\Delta t=T_{AB}+T_{BA}\,,\ \ \ \ \ \
\delta=\frac{1}{2}\left(T_{AB}-T_{BA}\right),$ (4)

where in the last equation we have used channel reciprocity, $\Delta
t_{AB}=\Delta t_{BA}$.

The quantum entanglement between in photon pairs created at Alice’s and/or
Bob’s locations does not play an important role in determining $\delta$,
except for the fact that the production of entangled pairs by SPDC ensures
that the photons in each pair are generated within a time window typically of
a few $100$ fs—several orders of magnitude smaller than the synchronization
precision we aim to obtain Valencia  _et al._ (2004b). In contrast, the
entanglement between photon pairs plays an important role in increasing the
security of this protocol to target malicious attacks by utilising well
established quantum cryptography techniques, such as Bell tests.

Additional difficulties appear when Alice and Bob are in relative motion. If
the distance between the two parties change in time, detection events
corresponding to different pairs of entangled photons will not accumulate in a
single peak of the correlation function. The peak will spread. If the spread
happens fast enough, the peak will be buried under the noise. On the contrary,
if synchronization can be achieved before noisy events accumulate in the
correlation function, the protocol would be successful. There is, therefore, a
competition between the duration of the protocol, $t_{Acq}$, and the noise
level. Note that $t_{Acq}$ cannot be made arbitrarily small, since enough
photons must be collected to get statistically significant peaks for the
cross-correlation functions. This effect is very similar to adding a relative
frequency difference between the clocks at A and B, which has been studied in
detail in Ref. Ho  _et al._ (2009).

The optimal precision this protocol can achieve when applied to moving clocks
can be obtained in terms of the following parameters: the ebit rate of the
SPDC sources, $R$; the total loss in the channel between the clocks A and B
(including detector losses), $\eta$; the radial relative velocity,
$v_{rel}^{rad}$; and the minimum (average) number of photons exchanged between
A and B needed for the cross-correlation peaks to be statistically
significant, $N_{\rm min}$.

Denoting the optimal sync precision by $t_{\rm bin}$, it is given by Haldar
_et al._ (2023b)

$t_{\rm bin}=\frac{N_{\rm min}\,v_{rel}^{rad}}{c\,R\,\eta},$ (5)

where $c$ is the speed of light in vacuum. Since $t_{\rm bin}$ gives the
smallest detectable time offset Haldar  _et al._ (2023b), a larger $t_{\rm
bin}$ means a worse sync precision.

Note that link transmissivity $\eta$ and the relative radial velocity
$v_{rel}^{rad}$ are functions of the distance between A and B and,
consequently, they change continuously in time. Throughout this paper, we will
use Eqn. (5) to determine the precision at which two moving clocks can
synchronize, applied to either satellite-satellite or satellite-ground
station. Distances between the parties will be computed by numerically
simulating the dynamics of a constellation of satellites distributed in
different orbits.

###  II.2 Precision shadow and sync traces

We introduce next a visual tool that becomes handy in determining the limits
of performance and scale of the QCS network. We use the term “satellite
precision shadow” for this tool. Given a satellite, the level of precision at
which it can sync with ground stations in Earth’s surface is determined by the
distance between the two parties, as we just explained above. Therefore, if we
fix a value of the desired sync precision, $t_{\rm bin}$, Eqn. (5) can be used
to determine the patch on Earth’s surface under the satellite at which
synchronization at such precision can be established. This is patch defines
what we call precision shadow of the satellite. Any two ground stations
located within the precision shadow at an instant can sync at optimal
precision $t_{\rm bin}$ using the satellite as a mediator.

But the previous definition of precision shadow does not take into account the
precision of the clock on board the satellite. Improving the precision of that
clock implies that a satellite can “hold time” for longer and, consequently,
be able to sync cities that are further apart. In other words, the size of the
precision shadow increases, elongating along the direction of the motion of
the satellite. The precision of the satellite’s clock can be quantified by
means of the “holdover time” $\tau$, defined as the average interval it takes
for the clock to accumulate a deviation equal to the desired sync precision
$t_{\rm bin}$.

As an example, Fig. 1 shows the precision shadow of a satellite in a 500 km
LEO polar orbit, when the satellite is just above the intersection of the
Prime Meridian and the Equator. We have chosen $t_{\rm bin}=1\,ns$ as sync
precision target. The plot shows the effect of increasing the satellite’s
clock holdover times $\tau$. Details of the simulations are spelled out in
more detail in the next section.

To illustrate the application of the sync protocol under consideration via our
numerical simulations, we choose two cities, New York and Salt Lake City in
this example, and compute the maximum precision at which they can sync by
means of a single LEO satellite in a 500 km polar orbit, as a function of
time. We call this type of plot _sync trace_ Haldar  _et al._ (2023a), since
it monitors the sync precision between two ground stations in the course of
time. Fig. 2 shows the sync trace between New York and Salt Lake City. The
maximum sync precision at a given instant $t$ is determined as the minimum
value of the two precisions at which the ground stations can individually sync
with the satellite, within an interval $\tau$ around $t$ (recall, $\tau$ is
the satellite’s clock holdover). Obviously, the maximum sync precision
achievable is given by the precision of the satellite’s clock—assuming that
the ground stations use better clocks. This absolute maximum precision is
indicated by a horizontal green line in Fig. 2, which we take to be 1 ns in
this example. This maximum can be improved just equipping satellites with
better clocks. As an aid to intuition, we mention here that off-the-shelf
Rubidium atomic clocks can hold time to 1 ns precision for roughly $10$
minutes DAR (2019). This will be a common feature of all sync traces shown in
this paper.111The precision of the clock on board satellite can be quantified
by the Allan deviation curve Oelker  _et al._ (2019), which is the standard
measure of a clock’s stability is often obtained empirically. Here, for our
simulations it suffices to take a simplistic, worst case scenario approach,
corresponding to precision $C$ for an interval $\tau$, and zero beyond that
interval.

Figure 1: 1 ns precision shadow of a single LEO satellite in a 500 km altitude
polar orbit when it is located just above the intersection of the Prime
Meridian and the Equator. All points within the grey shadow can sync with the
satellite at a precision better than 1 ns at this instant of time. As the
holdover increases from $\tau=0$ min (top) to $\tau=5$ min (center) and
finally $\tau=10$ min (bottom), the shadow elongates along the satellite’s
direction of motion, increasing the network scale.

The improvement in sync outcomes as the holdover time $\tau$ is increased for
the New York-Salt Lake City ground station pair ($\approx$ 3500 km apart) can
be seen in Fig. 2. For $\tau=0$ there was no connection due to visibility
constraints. As holdover increases, and thus also the extent of the satellite
shadow, some low-precision sync is possible twice per day. For $\tau=10$
minutes we finally see the capability of sub-nanosecond sync between the two
cities.

Having established the capabilities of a single satellite, we will now
investigate the effects of adding satellites, both with and without the
ability to synchronize amongst themselves, on ground-level sync outcomes.

Figure 2: Sync traces between New York and Salt Lake City, using a single LEO
satellite in a 500 km polar orbit, both in terms of precision and the number
of connections every day. As the holdover time $\tau$ is increased from (top)
$\tau=0$ min to (center) $\tau=5$ min, and finally to (bottom) $\tau=10$ min,
we start seeing sync events, first at worse than nanosecond precision, but
then, for 10 min holdover, we see sub-nanosecond level sync precision. The
green (horizontal) line indicates the cut-off introduced by the satellite
clock’s precision within the holdover time, which has been chosen to be 1 ns.

##  III Simulation methods

Our simulation consists mainly of two parts: (1) simulating the dynamics of
the satellite constellation and ground the stations and; (2) simulating the
lossy quantum communication channel between the two parties attempting
synchronization.

Regarding dynamics, we simulate the trajectory of a constellation of
satellites following circular orbits (extension to non-circular is possible),
and compute at every time step of the simulation the relative positions and
distances of each satellite and ground station. As described below, we use
these distances to compute the efficiency of the free space communication
channel.

Furthermore, we check the visibility of satellites from ground stations and
satellites among themselves, in the following manner. For a ground station
satellite pair, we check that the angle between the line joining the ground
stations to the satellite and the ground’s station zenith is less than
$\pi/2$. Whereas for satellite-satellite visibility, the distance from the
center of the Earth to the center of the line joining the two satellites, must
be greater than $R_{E}+\Delta_{\rm th}$, where $\Delta_{\rm th}$ is the
atmospheric thickness. We assume $\Delta_{\rm th}=10$ km. This condition makes
the satellite-satellite links impervious to effects of atmospheric scattering.

From the position of each satellite and ground station, we can evaluated the
relative radial velocity between any parties at every time step in the
simulation. The relative radial velocity between two parties (A,B) is given
by:

$v_{\rm rel}^{\rm
rad}=(\mathbf{r_{A}}\times\mathbf{\omega_{A}}-\mathbf{r_{B}}\times\mathbf{\omega_{B}}).\frac{(\mathbf{r_{A}}-\mathbf{r_{B}})}{|(\mathbf{r_{A}}-\mathbf{r_{B}})|}\,,$

where, $\mathbf{r_{A}},\mathbf{r_{B}}$ are their position vector with respect
to the center of the Earth, and $\mathbf{\omega_{A}},\mathbf{\omega_{B}}$ are
their angular velocities.

Regarding the communication channel, we take into account the various losses
by calculating effective efficiency factors. For the ground station-satellite
links, we must distinguish between the uplink (ground station to satellite)
and downlink (satellite to ground station) $\eta_{(up)}$ and $\eta_{(dwn)}$
respectively, since we assume different telescope radii for the two
(atmospheric turbulence and scattering effects would also be slightly
asymmetric, but our model does not include these details). $\eta_{(up)}$ is
the probability that a photon pair generated out of the SPDC source at the
ground station will lead to a double detection event, i.e., one partner will
be detected at A locally and the other at B after travelling through free
space and the atmosphere. Similarly for $\eta_{(down)}$. For satellite-
satellite links $\eta_{(up)}=\eta_{(down)}$.

We model the satellite-ground quantum communication channel as follows. For
concreteness, let us focus on a downlink channel. The following will hold
similarly for an uplink channel. We assume clear skies and approximate the
downlink channel as only lossy (background noise is accounted for at the
detectors). That is, photons are either transmitted through the channel or
lost in transmission. The dominant sources of loss are (1) beam spreading
(free-space diffraction loss), (2) atmospheric absorption/scattering, and (3)
non-ideal photodetectors on the satellite and on the ground. We characterize
these loss mechanisms by their transmittance values, which is the fraction of
the received optical power to the transmitted power (which is also equal to
the probability to transmit/detect a single photon). Let these transmittance
values be, respectively,

$\eta^{(dwn)}_{fs}(L),\quad\eta^{(dwn)}_{atm}(L),\quad\kappa_{sat},\quad\kappa_{gs}\,,$

where the superscripts refer to the downlink, $fs$ refers to free-space
diffraction loss, $atm$ to atmospheric loss, $L$ is the link distance
(physical distance) between satellite and receiver (which in turn depends on
the satellite altitude, $h$, position of the satellite in its orbit, and
position of the ground station on Earth’s surface), $\kappa$ denotes non-ideal
detection efficiencies for the onboard satellite detector array ($sat$) as
well as the detector array at the ground station ($gs$), and all transmittance
values are less than or equal to $1$. Simple analytic formulae are used to
estimate the free-space and atmospheric transmittance values in accordance
with Khatri  _et al._ (2021). Therefore, the overall efficiencies are given
as:

$\displaystyle\eta_{(dwn)}$ $\displaystyle=$
$\displaystyle\eta^{(dwn)}_{fs}\eta^{(dwn)}_{atm}\kappa_{sat}\kappa_{gs}\,,$
(6) $\displaystyle\eta_{(up)}$ $\displaystyle=$
$\displaystyle\eta^{(up)}_{fs}\eta^{(up)}_{atm}\kappa_{sat}\kappa_{gs}\,.$ (7)

The satellite-satellite link efficiencies can also be similarly evaluated. One
key difference, given our assumptions, is that these channels are free of
atmospheric scattering, i.e. $\eta^{(up)}_{atm}=\eta^{(dwn)}_{atm}=1$.

For convenience, we list in Table 1 other operational parameter entering in
our simulations along this paper, and the values for them we have used to
illustrate our techniques.

With this information, we use Eqn.(5) to evaluate the time sync precision
between the two parties at any time step of the simulation.

Altitude of orbits | h = 500 km
---|---
Operational wavelength |  $\lambda$ = 810 nm
Radii of telescopes | ($r_{sat}$, $r_{gs}$) = (10 cm, 60 cm)
Detector efficiencies | ($\kappa_{sat}$,$\kappa_{gs}$) = (0.5, 0.5)
Source rate |  $R$ = $10^{7}$ entangled-pairs/s
Table 1: Various operational parameters for the simulations

##  IV Single orbit master clock

In this section, we begin by discussing the simplest case of inter-satellite
synchronization: synchronizing a ring of satellites following the same orbit.
Inter-satellite synchronization is illustrated in Fig. 3. It will form the
backbone of our proposed master clock.

Figure 3: Schematic for a 3-link time sync involving two satellite-ground
station links and one satellite-satellite link. Each link has sync precision
$P_{i}=-{\rm log_{10}}(t_{bin,i})$, and the effective sync precision between
the two ground station clocks (GS1/GS2) is restricted to the minimum precision
over all links. These satellite-satellite links can arise through either
intra-orbit or inter-orbit synchronization, which are explored in detail in
Sections IV and V.

###  IV.1 Intra-orbit synchronization

Consider a ring of satellites in circular orbit around the Earth. When the
satellites do not sync with one another, for two cities to synchronize they
must simultaneously lie within the precision shadow of the same satellite.
This requirement greatly constrains the possibilities of synchronization,
given the limited size of the shadows. It is challenging to sync cities
separated by more than $\approx$ 1000 km, unless satellites are equipped with
high quality clocks (see Figs. 1). Another option is allowing the (poor/medium
quality) satellite clocks to synchronize with each other. With this
functionality, for two cities to synchronize, they do not need to connect to a
common satellite, but must sync to any satellite in the constellation within a
time window $\tau$.

In this subsection, we analyze the necessary conditions for adding this
functionality —intra-orbit synchronization— as well as the consequent
improvement in precision and network scale of the QCS network.

An obvious condition for intra-orbit sync is satellite to satellite
visibility. If we place the first satellite at the point (0, 0, $R_{E}+h$),
with respect to the Earth’s center, and where $R_{E}$ is the Earth’s radius,
the angle the second satellite makes with the positive $z$-axis must be less
than $\pi-2\theta_{s}$. Here $\theta_{s}$, which we call the shadow angle, is

$\theta_{s}=\arcsin(\frac{R_{E}}{R_{E}+h}).$ (8)

If we consider polar orbits at $h=500$ km altitude, this corresponds to
$\theta_{s}\approx 68^{\circ}$. For an orbit of equally spaced satellites at
this altitude, nine satellites suffice for the visibility condition to be
satisfied at any instant, and each satellite is always in view from the
nearest neighbor. Increasing the orbital height decreases the required number
of satellites for continuous visibility.

For two satellites A and B in the same (circular) orbit, the orbital
velocities are equal, i.e., $\mathbf{\omega_{A}}=\mathbf{\omega_{B}}$, so
$v_{\rm rel}^{\rm rad}=0$. Therefore, according to Eqn. (5), the sync
precision has no lower bound and is only limited in practice by the time-
stamping precision. Further, once nearest-neighbor satellites sync with each
other, all satellites become synchronized transitively. Thus, to form a master
clock comprised of a ring of polar satellites sharing a common time, the
visibility condition becomes the only requirement which, as mentioned above,
is solved by adding enough satellites to the orbit.

###  IV.2 Synchronization precision and network size

Once the visibility condition between satellites has been established, the
task of syncing two ground stations GS1 and GS2 via two satellites SAT1 and
SAT2 proceeds as follows (see Fig. 3 for an illustration).

We first find the minimum precision that GS1 can get at time $t$ from all
satellites in the orbit, and then the minimum precision that GS2 can get from
the orbit within a period $(t-\tau/2,t+\tau/2)$, where $\tau$ is the holdover
time of clocks on board satellites (assuming all identical). The larger of
these two values is the sync precision for the city pair. This protocol uses
that satellites can sync among themselves at better precision than satellites
with ground stations, because they do not move relative to each other.

We consider now an example of polar orbit with 10 satellites at 500 km
altitude. In Figs. 4 and 5, we compare the regions of achievable nanosecond
sync precision for a disconnected orbit — no inter-satellite sync — to an
orbit with intra-orbit sync functionality. We show the result for different
holdover times $\tau=0,5,$ and $10$ min. The differently colored precision
shadows in Fig. 4 indicate that they are disconnected and independent shadows,
hence sync between two cities occurs only when they fall within one single
shadow patch (one single color). Once the intra-orbit synchronization
functionality is added (Fig. 5), all shadow patches merge together to form a
single bigger shadow. With this functionality established, it is clear that
the network scale (geographic extent) is enhanced substantially. Without any
holdover capabilities, each shadow is $\approx 8^{\circ}$ in width, and
increases by about $3.5^{\circ}$ per added minute of holdover. We see in Fig.
5 that in order to have continuous latitudinal connection (bottom plot) with
this configuration of ten satellites, we require 10 minutes of holdover time.

Before discussing detailed examples, we further stress here the advantage
achieved via intra-orbit connectivity. Equipped with modest clocks with 5-10
mins of holdover (off-the-shelf Rubidium atomic clocks) we can provide high
precision synchronization between cities on the opposite ends of the globe
(latitude-wise). For example, 1 ns sync precision is achievable between New
York City in the US and Puerto Montt in Chile with 5 mins of holdover time
(see Fig. 6 (middle)), and between New York City and Beijing, China, with 10
mins of holdover time (see Fig. 10 (top)), with the simple 1 orbit
configuration discussed above. This corresponds to a sync network range of
roughly 10000 km (considering geodesic distances between the cities). In
contrast, for disconnected satellites this scale could not surpass 4000 km
with similar holdover times. The intra-orbit functionality changes the
constellation capabilities from continental to global scale.

Figure 4: 1 ns precision shadows for a constellation of 10 satellites in a 500
km polar orbit along the Prime Meridian. In this plot, the satellites are not
in sync with each other, and consequently the different patches are
independent. Each plot shows a different value of the holdover times: $\tau=0$
min (top), $\tau=5$ min (middle), and $\tau=10$ min (bottom).

Figure 5: 1 ns sync precision shadow for an orbit of 10 satellites in a 500 km
altitude polar orbit along the Prime Meridian. The satellites are synchronized
with each other, and consequently the the individual patches are part of a
unique shadow. Different holdover times are shown: $\tau=0$ min (top),
$\tau=5$ min (middle), and $\tau=10$ min (bottom).

###  IV.3 Sync Examples

We will now present some concrete examples of QCS using city pairs along
similar longitudes, and compare the sync capabilities with and without intra-
orbit connectivity. We choose cities along the same longitude because adding
satellites to the same orbit only increases the coverage area of our QCS
network along the direction of the satellites’ motion (latitude). As above, we
simulate an orbit of 10 equally-spaced satellites in a polar orbit along the
Prime Meridian, and choose other parameters as in Table 1. The precision is
recorded over a period of one day, with a simulation time step of $1\rm s$.
The following “figures of merit” are evaluated: average precision of sync
events at precision better than 1ns, the largest time gap between sync events
with precision better than 1ns, and the percentage of the day that ground
stations are synchronize with precision better than 1ns. For cases in which
either or both of the cities are not visible to any satellite, we use the term
“invis.” in the data. Similarly, when sync events are possible but are not
better than 1ns in precision, we report the average precision above 10 $\mu$s
and neglect the remaining figures. Table 2 contains these figures of merit for
sync events between New York/Montreal (New York City-Montreal), New York/Port
Au Prince (New York City-Port Au Prince), and New York/Puerto Montt (New York
City-Puerto Montt), for holdover times of 0, 2, 5, and 10 minutes. Fig. 6
shows precisions over one day for holdover time $\tau$ = 5 min for the same
city pairs. Note that the recorded and plotted precision values are now
reported as $-\rm log_{10}(t_{bin})$, so an increase in this value indicates
an ability to sync clocks with better precision. For example, for nanosecond
and microsecond precision this quantity is equal to 9 and 6, respectively.

New York City-Montreal Figures of Merit
---
| Disconnected | Connected
Holdover | 0 min | 2 min | 5 min | 10 min | 0 min | 2 min | 5 min | 10 min
Average Precision | invis | 9.99 | 10.67 | 11.07 | invis | 9.99 | 10.67 | 11.17
Largest Gap (min) | invis | 494 | 491 | 486 | invis | 494 | 491 | 486
Connected Fraction (%) | invis | 8 | 18 | 33 | invis | 8 | 18 | 33

New York City-Port Au Prince Figures of Merit
---
| Disconnected | Connected
Holdover | 0 min | 2 min | 5 min | 10 min | 0 min | 2 min | 5 min | 10 min
Average Precision | invis | 8.30 | 9.13 | 10.72 | invis | 9.02 | 10.04 | 11.18
Largest Gap (min) | invis | n/a | 551 | 527 | invis | 560 | 538 | 527
Connected Fraction (%) | invis | n/a | 2 | 16 | invis | 0.20 | 10 | 26

New York City-Puerto Montt Figures of Merit
---
| Disconnected | Connected
Holdover | 0 min | 2 min | 5 min | 10 min | 0 min | 2 min | 5 min | 10 min
Average Precision | invis | invis | invis | invis | invis | 9.04 | 10.02 | 11.19
Largest Gap (min) | invis | invis | invis | invis | invis | 577 | 495 | 490
Connected Fraction (%) | invis | invis | invis | invis | invis | 2 | 13 | 32

Table 2: Figures of merit for the clock synchronization of city pairs along
similar longitudes for various holdover times. The precision numbers are
represented as absolute values in the log scale (higher is better). We use
“invis.” to mean that one or both cities is invisible to any satellite in the
orbit, thus synchronization is not possible. When the cities are visible but
the average achievable sync precision is below the threshold of 1ns, we change
the threshold to 10 $\mu$s and do not report the largest gap or connected
fraction. This is done to indicate that, while the cities are visible to a
common satellite, they are unable to achieve a precision greater than 1ns. The
three regimes of intra-orbit connectivity advantages are clear from these
results: New York City-Montreal does not significantly benefit from
connectivity as the cities always lie in the shadow of a single satellite; New
York City-Port Au Prince has some increase in all figures; the most
significant benefit occurs with New York City-Puerto Montt as neither city is
visible to the shadow of a single satellite simultaneously, thus intra-orbit
connection results in an exponential benefit in average sync precision.

Figure 6: Sync precision versus time for various city pairs with 5 minutes of
holdover time. The simulation spans one day. The green (horizontal) line
represents the 1 ns cut-off, introduced by the satellite clock’s precision
within the holdover time. Figures on the left represent the sync precision
when the satellites in an orbit are unsynchronized, whereas on the right we
consider a constellation with intra-orbit sync capabilities. Clearly New York
City-Montreal(top) pair retains the same precision regardless of intra-orbit
connectivity as explained above. New York City-Port Au Prince (center)
significantly increases in average precision as well as portion of the day
where 1ns sync is possible. New York City-Puerto Montt (bottom) can sync only
when the satellites are sync with each other as they will never lie in the
same shadow segment.

These results show that in regimes where pairs of ground stations fall outside
of the same shadow, as in Fig. 6(e), or on the edges of the same shadow, as in
Fig. 6(c), equipping the orbit with satellite to satellite sync capability
greatly improves the achievable precision and network scale. On the other
hand, pairs of cities which are close enough to fall in the shadow of a single
satellite do to get any extra benefit from this new functionality, as
expected. These results show that benefits that could be achieved by simply
equipping satellites with better clocks can also be achieved, and even
surpassed, by allowing satellites to sync amongst themselves.

##  V Multi-orbit master clock

In this section, we introduce another level of functionality to the satellite
constellation, namely synchronization between satellites in different orbits—
inter-orbit synchronization. This is motivated by the requirement to increase
the scale of the network longitude-wise, since adding holdover and allowing
intra-orbit sync only allows the extension of the network scale along the
satellite’s direction of motion, which for polar orbits is approximately along
fixed longitudes. Thus, to provide global coverage, inter-orbit
synchronization is a necessary feature. At the same time, adding this
functionality is more challenging than intra-orbit connectivity, because
satellites in different orbits move at high velocities relative to each other,
unlike ones in the same orbit which are at relative rest. Therefore, the high
precision sync achieved between satellites in the same orbit is only possible
for nearby orbits. We will explore and quantify this difficulty, as well as it
advantages.

Figure 7: Possible links between satellites. The solid red line represents an
intra-orbit link between satellites 1 and 2 in the same orbit, allowing a
given orbit to maintain a constant sync precision between all members. The
yellow dashed line is an inter-orbit link between a conjugate satellite pair
(i.e. two satellites with the same number in a neighboring orbit) and the
green dashed line is an inter-orbit link between two non-conjugate satellites.
The conjugate pairs are the best candidates for inter-orbit sync due to their
closer proximity (and thus less loss) than other pairs. Note that the
simulations introduce an orbital offset of 0.25∘ (too small to be shown in
this figure) to prevent satellite collision at the poles.

###  V.1 Inter-orbit synchronization

Let us begin by describing the relative dynamics of satellites in neighbouring
orbits, since, analogous to the intra-orbit case, if nearest-neighbour orbits
are synced all orbits are synced by transitivity. Let us consider two polar
orbits both at the same altitude. Since such orbits intersect at the poles, a
small offset must be added in the initial polar angles of the satellites to
avoid catastrophe. Next we must ensure visibility between the satellites. This
is again analogous to the single orbit case. Since 9 satellites were needed to
ensure visibility amongst satellites in a ring, translating to a polar angle
separation of $40^{\circ}$ along the orbit, the two orbits must also be no
more than $40^{\circ}$ apart azimuthally for satellites in different orbits to
see each other near the Equator. Since the Equator is the point of maximum
separation between the orbits, this condition ensures visibility throughout
the satellites’ journey in their respective orbits. Next, we come to the
question of sync precision. We know from Eqn. (5), that high precision can be
achieved if the two clocks are very close ($\eta$, the channel transmissivity,
is high) and if they move slowly relative to each other. Intuitively, the
pairs of satellites that are conjugate to each other, i.e., say the pair of
satellite labeled with the 3 in each orbit (see Fig.7 ), are the ones that
remain closest neighbours throughout the orbital motion. Two points in the
orbit are of special interest: (1) near the Poles, where the distances are
minimum, but relative radial velocity is maximum, and (2) near the Equator,
where conjugate satellites are at maximum separation but their velocities
become almost parallel and perpendicular to the line joining them, barring the
effect of the small offset. Thus, conjugate pairs of satellites are good
candidates for ensuring sync between orbits. Nonetheless, this trade-off
between separation and relative radial velocity makes their sync dynamics non-
trivial.

Further, if any conjugate pair between orbits can sync, the entire orbits are
synced since we have already defined the conditions needed for intra-orbit
sync and determined that it is continuously available at a high level of
precision. Thus, the sync precision of the master clock is limited by the
inter-orbit sync capabilities, which we shall discuss now with a concrete
example.

Let us now imagine a QCS constellation with 5 orbits (this ensures constant
visibility between nearest-neighbour orbits via conjugate satellites) and 10
satellites in each orbit (this ensures nearest-neighbour intra-orbit
visibility). All orbits are polar and at a constant altitude of 500 km. The
polar angle offset between satellites in neighbouring orbits is chosen to be
$0.25^{\circ}$. Our aim is to determine the precision of this master clock as
a function of time.

To this end, let us first look at the sync between conjugate satellites for a
fixed pair of next-to-each-other orbits, say orbit 1 and 2. We see in Fig. 8
that high precision is obtained roughly once in 4.5 minutes. This can be
understood in the following way. The members of each conjugate pair can sync
with high precision four times in its orbit – twice at the Poles and twice at
the Equator. The period of the orbit is roughly 90 min, and there are 10
satellites per orbit. This gives us a gap of around 2.25 minutes between
peaks. But, the behaviour of the conjugate pairs providing service to the
Eastern Hemisphere is identical to the ones on the other side, and thus peaks
in sync traces arising from satellites numbered 1-5 overlap with those arising
from satellites numbered 6-10. Therefore, we see roughly a 4-5 minutes gap
between successive peaks. Therefore, if the satellite clocks had a holdover of
roughly 5 minutes the satellites in the constellation be in continuous sync at
sub-nanosecond level precision.

Figure 8: Sync precision traces between orbit numbers 1 and 2 in a QCS
constellation. The constellation consists of 10 satellites in 5 polar orbits
each at 500 km altitude. The different colors represent precision between
different conjugate pairs of satellites (E.g., precision between the pair of
satellites labelled as 1 of the orbits 1 and 2 is the brown colored curve, the
pair of satellites labelled as 2 is pink, etc.) All the peaks in the plot
provide sub-nanosecond precision between orbit numbers 1 and 2. The gaps
between these peaks are roughly 4-5 minutes, thus, a holdover time of 5
minutes for satellite clocks is enough to maintain sub-nanosecond precision
continuously.

On the other hand, there is also the effect of offset between orbits. Let us
look at the precision between the 1st satellite of orbit pairs 1-2, 2-3, 3-4
and 4-5. From Fig. 9 (top) it seems that all satellite pairs achieve high
precision at the same time. But zooming into one of these peaks, in Fig. 9
(bottom) we find that these peaks are all slightly shifted w.r.t. each other
due to the initial offset. Thus, again, to achieve high precision between
different orbits, we need some holdover for the satellite clocks. In this
case, a holdover of around 30 seconds is sufficient, which is much smaller
than the holdover requirement discussed above and shown in Fig. (8).

Figure 9: Effect of polar angle offset between conjugate satellites in
neighbouring orbits on the sync precision. (Top) All orbit pairs sync at high
precision and apparently also at the same time. B ut when we zoom into the
peaks (Bottom) the effect of the polar offset becomes clear. For a small
offset of $0.25^{\circ}$, the peaks for precision between different orbit
pairs are displaced w.r.t. each other by not more than 30 seconds. Hence a
holdover of 5 minutes required for covering the gaps between sync of different
conjugate pairs is enough to also compensate this effect.

In summary, we conclude that a master clock with sync precision at least at
the sub-nanosecond level can be built with constellation of 50 LEO satellites
equipped with modest clocks with a small holdover time of $\approx 5$ min.
Here, we must re-stress that the working principle of the master clock is to
reinforce low-stability clocks forming a network through quantum links. For
example, with the constellation just mentioned, clocks that can hold time at
sub-nanosecond level for only a few minutes are combined to build a master
clock that holds a common time indefinitely. As a consequence, the coverage
area —network scale— on the ground also extends many-fold. Such a master clock
is a primitive to a quantum network and imbibes many of the latter’s essential
features.

###  V.2 Sync Examples

To illustrate the enhancement in sync outcomes provided by inter-orbit sync
functionality, we present an example with three city pairs along similar
latitudes, namely New York-Salt Lake City, New York-Madrid, and New York-
Beijing. The constellation is the same as discussed in the previous
subsection, and we assume a modest holdover time of 10 minutes. This choice is
also motivated by the fact that it leads to a continuous 1 ns precision shadow
for an orbit (see Fig. 5 (bottom)).

We are now interested in the fraction of the day for which sub-nanosecond sync
is possible. Similar to the examples presented in the single-orbit case, we
choose city pairs representing three distinct regimes of improvement over the
case where no inter-orbit sync is assumed (but intra-orbit synchronization).
We refer to the former as _connected_ and the latter as _disconnected_ case in
this context. In our example, although New York and Beijing have the largest
separation distance of the pairs presented, they always lie within the
precision shadow of the same orbit (on the opposite sides of the Earth). Thus,
their ability to synchronize is not affected by whether or not the orbits can
synchronize amongst themselves. This is evident in Table 3, as both the
disconnected and connected cases have 100$\%$ daily sync coverage at the 0.1
ns level (see Fig. 10 (top)). This example also provides an interesting use
case for a single orbit master clock.

The next case showing an intermediate benefit from inter-orbit sync is seen
through the New York-Salt Lake City pair. The cities are close enough to be
visible to the same orbit but achieve sub-nanosecond sync precision only about
$60\%$ of the day without inter-orbit connectivity. We see this visually in
Fig. 10 (middle) as the precision peaks are discontinuous in the disconnected
case (left) and become fully continuous once the orbits can sync (right). So
by allowing inter-orbit sync, we gain an extra $40\%$ of the day in which sub-
nanosecond sync is possible. Finally, a drastic improvement is achieved for
the the New York-Madrid case. It is not possible to sync these cities without
inter-orbit synchronization, as the two cities are never simultaneously
visible to the same orbit. However, since either city is always visible to
some orbit, when the orbits are continuously synced, sub-nanosecond sync
between these cities can be achieved throughout the day (see Fig. 10
(bottom)).

| New York-Salt Lake | New York-Madrid | New York-Beijing
---|---|---|---
| Disconn. | Conn. | Disconn. | Conn. | Disconn. | Conn
Largest Gap (min) | 65 | 0 | invis | 0 | 0 | 0
Conn. Frac. (%) | 60 | 100 | invis | 100 | 100 | 100

Table 3: Figures of merit for the clock synchronization of city pairs along
similar latitudes for the “master clock” configuration.

Figure 10: Comparison of sync precision values without (left) and with (right)
inter-orbit synchronization for city pairs along similar latitudes for a
constellation made of 5 orbits of 10 satellites each with $\tau$= 10 minutes.
The pairs are New York-Beijing (top), New York-Salt Lake City (middle), and
New York-Madrid (bottom). The horizontal green (solid) line represents the 1
ns cut-off, introduced by the satellite clock’s precision within the holdover
time. The plots show that, synchronization at 1 ns precision is continuously
available between the selected cities. A modest improvement is seen in the
sync connectivity between New York and Salt Lake City, and a drastic
improvement between New York and Madrid.

##  VI Global Network

Now that we have discussed the advantages of inter-satellite synchronization
and proposed the design of a master clock constellation, it remains to be
shown that this clock can indeed support sub-nanosecond sync precision between
any two cities on Earth. Rather than carefully selected cities along constant
latitudes/longitudes as previously done, we now calculate the sync precision
for pairs of cities randomly spread across 6 continents: Seattle, USA; London,
United Kingdom; Rio Grande, South America; New Delhi, India; Cape Town, South
Africa; and Sydney, Australia. The constellation is still the same previously
discussed.

Our simulations show that sub-nanosecond sync can be achieved for 100% of the
day for all city pairs for such a configuration, thus providing a continuously
connected global sync network. Fig. 11 affirms this for 3 city pairs, where
$<1$ ns sync precision is achievable over the entire day. However, as
previously discussed, the constellation was carefully selected to have a
limiting precision of 1 ns. Thus, we can conservatively conclude that a global
timing network at 1 ns sync precision and continuous coverage is achievable.

Nevertheless, these sync traces show that $10^{-10}-10^{-11}\rm s$ level sync
precision could be achieved while maintaining continuous global connectivity,
if the master-clock induced nanosecond cut-off could be loosened. One way to
do this is by adding more satellites to each orbit, which will reduce the
holdover requirement for global coverage and also for inter-orbit QCS (see
Fig. 8).

Figure 11: Sync precision over one day for a global network: Seattle-Cape Town
(top;16,415 km separation), Sydney-New Delhi (middle; 10,422 km), Rio Grande-
London (bottom; 7,801 km). This example provides a further proof of concept
that given the previously defined master clock configuration, sub-nanosecond
sync precision is achievable 100% of the day, regardless of city separation
distance.

##  VII Conclusions

Enhancing the precision of time synchronization not only pushes the boundaries
of fundamental science but also has important technological use cases Troupe
_et al._ (2022), both in the classical and quantum realms. Classical
synchronization techniques such as those using optical frequency combs have
been demonstrated to obtain femtosecond level precision, but they currently
lack scalability Deschênes  _et al._ (2016). The utilization of quantum
resources presents an alternative avenue, accompanied by complementary
advantages, and is thus deserving of exploration. For instance, QCS comes with
an added layer of quantum security which is not available in any classical
protocol. Thus, in this paper, by proposing a method to extend QCS to a global
scale, we provide a forward looking solution that would form the basis of
future secure communications, navigation and ranging. We show that our scheme
provides continuous global synchronization at the sub-nanosecond level using a
quantum-assisted master clock in the sky. This is done by envisioning
synchronization between satellites in the constellation. To this end, we
consider inter-satellite links, show their feasibility, and quantitatively
prove their ability to provide continuous and global time synchronization
coverage.

Our analyses has been done with some reasonable simplifying assumptions, as
mentioned in Section I. We are aware that, in a real experiment, the prevalent
levels of noise and jitter of channels, detectors, and time-stampers may be
higher than our estimates. Current free-space QCS demonstrations require
$\approx 1\;\rm s$ of acquisition time Lafler and Lanning (2022); Spiess  _et
al._ (2021), which is much more than the $\approx 1\;\rm ms$ we use in this
work, and which helps bypass the effects of relative velocities. Nonetheless,
our analyses provides useful trends both qualitative and quantitative that
will inform design parameters for an implementation. Although we underestimate
noise, we also use very conservative parameters, representing off-the-shelf
SPDC sources and detectors. Entanglement sources with rates of $\approx
10^{9}-10^{10}\;s^{-1}$ are now available and we anticipate that they would
soon be space ready Anwar  _et al._ (2022); Lohrmann  _et al._ (2020).
Similarly, SNSPDs (superconducting nanowire single-photon detectors) can
provide efficiencies much beyond 90%, far outperforming simple semiconductor
based detectors we have assumed Natarajan  _et al._ (2012). These improvements
will allow the acquisition time to stay within the millisecond range, even
with increased channel noise/loss, justifying our results in this work.
Finally, here we have only focused on the use of bipartite entanglement for
clock synchronization. Generation of multipartite entanglement has seen
considerable technological progress recently Deng  _et al._ (2007); Pan  _et
al._ (2012); Krenn  _et al._ (2017); Bergamasco  _et al._ (2017) and it would
be useful to explore further how such states could be used as a tool to
synchronise more than two parties at the same time (also see Kómár  _et al._
(2014)).

The main conclusions of our analysis can be summarised as follows:

  * •

QCS links between satellites in the same orbit increase the latitudinal extent
of the network. Using this minimal functionality, a high precision, single
orbit master clock can be set up, since there are no relative velocities
between satellites in the same orbit. Therefore, this is a functionality that
can be achieved without the requirement of any holdover for the satellite
clocks. Thus, surprisingly, very low stability clocks can reinforce each other
within the same orbit, and then provide high precision synchronization on the
ground.

  * •

By allowing different orbits to synchronize amongst themselves, the extent of
the network increases longitudinally. It should be noted that all other
constellation enhancements, such as increasing satellite clock holdover,
number of satellites, etc., could only increase the network scale along the
direction of the satellite’s motion. Hence, to achieve global coverage, inter-
orbit synchronization is essential.

  * •

Finally, as an example we have shown that a master clock comprising of 50
satellites distributed uniformly in five 500 km polar orbits, and equipped
with both intra and inter orbit QCS links is able to provide sub-nanosecond
precision globally, uninterrupted. This ability is beyond what is currently
achievable via the GPS.

In conclusion, the above analysis makes the proposed QCS constellation a
promising precursor to a quantum version of GPS. Our proposal can also be
considered as a “quantum-assisted network of clocks”. Individual nodes have
low timing stability, but together they form an highly precise master clock
network.

###### Acknowledgements.

We have benefited from discussions with Rachel McDonald. I.A., S.D. and S.H.,
are supported by the NSF grants PHY-2110273, by the RCS program of Louisiana
Boards of Regents through the grant LEQSF(2023-25)-RD-A-04, by funds from
Xairos Inc., and by the Hearne Institute for Theoretical Physics.

## References

  * Khatri  _et al._ (2021) Sumeet Khatri, Anthony J. Brady, Renée A. Desporte, Manon P. Bart, and Jonathan P. Dowling, “Spooky action at a global distance: analysis of space-based entanglement distribution for the quantum internet,” NPJ Quantum Information 7, 4 (2021).
  * Sidhu  _et al._ (2021) Jasminder S. Sidhu, Siddarth K. Joshi, Mustafa Gündogan, Thomas Brougham, David Lowndes, Luca Mazzarella, Markus Krutzik, Sonali Mohapatra1, Daniele Dequal, Giuseppe Vallone, Paolo Villoresi, Alexander Ling, Thomas Jennewein, Makan Mohageg, John Rarity, Ivette Fuentes amd Stefano Pirandola, and Daniel K. L. Oi, “Advances in space quantum communications,” IET Quantum Communication 1, 1–11 (2021).
  * California-Polytechnic-State University (Last accessed: Dec. 26 2021) California-Polytechnic-State University, “CubeSat Design Specification Revision 14,” https://www.cubesat.org/cds-announcement (Last accessed: Dec. 26 2021).
  * Haldar  _et al._ (2023a) Stav Haldar, Ivan Agullo, Anthony J. Brady, Antía Lamas-Linares, W. Cyrus Proctor, and James E. Troupe, “Towards global time distribution via satellite-based sources of entangled photons,” Phys. Rev. A 107, 022615 (2023a).
  * Ho  _et al._ (2009) Caleb Ho, Antía Lamas-Linares, and Christian Kurtsiefer, “Clock synchronization by remote detection of correlated photon pairs,” New Journal of Physics 11, 045011 (2009).
  * Lamas-Linares and Troupe (2018) Antía Lamas-Linares and James Troupe, “Secure quantum clock synchronization,” SPIE 10547, Advances in Photonics of Quantum Computing, Memory, and Communication XI 10547 (2018).
  * Lee  _et al._ (2019) Jianwei Lee, Lijiong Shen, Alessandro Cerè, James Troupe, Antia Lamas-Linares, and Christian Kurtsiefer, “Symmetrical clock synchronization with time-correlated photon pairs,” Applied Physics Letters 14, 101102 (2019).
  * Spiess  _et al._ (2021) Christopher Spiess, Sebastian Töpfer, Sakshi Sharma, Andrej Kržič, Meritxell Cabrejo Ponce, Uday Chandrashekara, Nico Lennart Döll, Daniel Rieländer, and Fabian Steinlechner, “Clock synchronization with correlated photons,”  (2021).
  * Oi  _et al._ (2017) Daniel KL Oi, Alex Ling, Giuseppe Vallone, Paolo Villoresi, Steve Greenland, Emma Kerr, Malcolm Macdonald, Harald Weinfurter, Hans Kuiper, Edoardo Charbon, and Rupert Ursin, “CubeSat quantum communications mission,” EPJ Quantum Technology 4 (2017), 10.1140/epjqt/s40507-017-0060-1.
  * Villar  _et al._ (2020) Aitor Villar, Alexander Lohrmann, Xueliang Bai, Tom Vergoossen, Robert Bedington, Chithrabhanu Perumangatt, Huai Ying Lim, Tanvirul Islam, Ayesha Reezwana, Zhongkan Tang, Rakhitha Chandrasekara, Subash Sachidananda, Kadir Durak, Christoph F. Wildfeuer, Douglas Griffin, Daniel K. L. Oi, and Alexander Ling, “Entanglement demonstration on board a nano-satellite,” Optica 7, 734 (2020).
  * Yin  _et al._ (2020) Juan Yin, Yu-Huai Li, Sheng-Kai Liao, Meng Yang, Yuan Cao, Liang Zhang, Ji-Gang Ren, Wen-Qi Cai, Wei-Yue Liu, Shuang-Lin Li, Rong Shu, Yong-Mei Huang, Lei Deng, Li Li, Qiang Zhang, Nai-Le Liu, Yu-Ao Chen, Chao-Yang Lu, Xiang-Bin Wang, Feihu Xu, Jian-Yu Wang, Cheng-Zhi Peng, Artur K. Ekert, and Jian-Wei Pan, “Entanglement-based secure quantum cryptography over 1,120 kilometres,” Nature 582, 501–505 (2020).
  * Dai  _et al._ (2020) Hui Dai, Qi Shen, Chao-Ze Wang, Shuang-Lin Li, Wei-Yue Liu, Wen-Qi Cai, Sheng-Kai Liao, Ji-Gang Ren, Juan Yin, Yu-Ao Chen, Qiang Zhang, Feihu Xu, Cheng-Zhi Peng, and Jian-Wei Pan, “Towards satellite-based quantum-secure time transfer,” Nature Physics 16, 848–852 (2020).
  * O’Connor  _et al._ (2019) A.C. O’Connor, M.P. Gallaher, K. Clark-Sutton, D. Lapidus, Z.T. Oliver, T.J. Scott, D.W. Wood, M.A. Gonzalez, E.G. Brown, and J. Fletcher, “Economic benefits of the global positioning system (GPS),” Report sponsored by the National Institute of Standards and Technology (2019).
  * (14) S. Y. Zhu and E. Groten, “Relativistic effects in GPS,” in  _GPS-Techniques Applied to Geodesy and Surveying_ (Springer-Verlag) pp. 41–46.
  * Troupe  _et al._ (2022) James Troupe, Stav Haldar, Ivan Agullo, and Paul Kwiat, “Quantum clock synchronization for future nasa deep space quantum links and fundamental science,” (2022), arXiv:2209.15122 .
  * McGrew  _et al._ (2018) W. F. McGrew, X. Zhang, R. J. Fasano, S. A. Schäffer, K. Beloy, D. Nicolodi, R. C. Brown, N. Hinkley, G. Milani, M. Schioppo, T. H. Yoon, and A. D. Ludlow, “Atomic clock performance enabling geodesy below the centimetre level,” Nature 564, 87–90 (2018).
  * (17) “Wrin1588 · wiki · projects / white rabbit standardization,” https://ohwr.org/projects/wr-std/wiki/wrin1588.
  * Delva  _et al._ (2012) P. Delva, F. Meynadier, C. Le Poncin-Lafitte, P. Laurent, and P. Wolf, “Time and frequency transfer with a microwave link in the ACES/PHARAO mission,” 2012 European Frequency and Time Forum (2012).
  * Deschênes  _et al._ (2016) Jean-Daniel Deschênes, Laura C. Sinclair, Fabrizio R. Giorgetta, William C. Swann, Esther Baumann, Hugo Bergeron, Michael Cermak, Ian Coddington, and Nathan R. Newbury, “Synchronization of distant optical clocks at the femtosecond level,” Physical Review X 6 (2016).
  * Bergeron  _et al._ (2019) Hugo Bergeron, Laura C. Sinclair, William C. Swann, Isaac Khader, Kevin C. Cossel, Michael Cermak, Jean-Daniel Deschênes, and Nathan R. Newbury, “Femtosecond time synchronization of optical clocks off of a flying quadcopter,” Nature Communications 10 (2019).
  * Haldar  _et al._ (2023b) Stav Haldar, Ivan Agullo, and James E. Troupe, “Synchronizing clocks via satellites using entangled photons: Effect of relative velocity on precision,”  (2023b), arXiv:2209.15122 .
  * Yurtsever and Dowling (2002) Ulvi Yurtsever and Jonathan P. Dowling, “Lorentz-invariant look at quantum clock-synchronization protocols based on distributed entanglement,” Phys. Rev. A 65, 052317 (2002).
  * Giovannetti  _et al._ (2001) Vittorio Giovannetti, Seth Lloyd, and Lorenzo Maccone, “Quantum-enhanced positioning and clock synchronization,” Nature 412, 417 (2001).
  * Giovannetti  _et al._ (2004) Vittorio Giovannetti, Seth Lloyd, Lorenzo Maccone, Jeffrey H. Shapiro, and Franco N. C. Wong, “Conveyor-belt clock synchronization,” Physical Review A 70, 043808 (2004).
  * Valencia  _et al._ (2004a) Alejandra Valencia, Giuliano Scarcelli, and Yanhua Shih, “Distant clock synchronization using entangled photon pairs,” Applied Physics Letters 85, 2655 (2004a).
  * Jozsa  _et al._ (2000) Richard Jozsa, Daniel S. Abrams, Jonathan P. Dowling, and Colin P. Williams, “Quantum clock synchronization based on shared prior entanglement,” Phys. Rev. Lett. 85, 2010–2013 (2000).
  * Ilo-Okeke  _et al._ (2020) Ebubechukwu O. Ilo-Okeke, Louis Tessler, Jonathan P. Dowling, and Tim Byrnes, “Entanglement-based quantum clock synchronization,” in  _FIFTH INTERNATIONAL CONFERENCE ON QUANTUM TECHNOLOGIES (ICQT-2019)_ (AIP Publishing, 2020).
  * Lafler and Lanning (2022) Randy Lafler and R. Nicholas Lanning, “Quantum time transfer for freespace quantum networking,”  (2022).
  * Kómár  _et al._ (2014) P. Kómár, E. M. Kessler, M. Bishof, L. Jiang, A. S. Sørensen, J. Ye, and M. D. Lukin, “A quantum network of clocks,” Nature Physics 10, 582–587 (2014).
  * Nichol  _et al._ (2022) B. C. Nichol, R. Srinivas, D. P. Nadlinger, P. Drmota, D. Main, G. Araneda, C. J. Ballance, and D. M. Lucas, “An elementary quantum network of entangled optical atomic clocks,” Nature 609, 689–694 (2022).
  * Belmonte  _et al._ (2017) A. Belmonte, M.T. Taylor, L. Holberg, and J. Kahn, “Effect of atmospheric anisoplanatism on earth-to-satellite time transfer over laser communication links,” Optics Express 25, 15676 (2017).
  * Taylor  _et al._ (2020) M.T. Taylor, A. Belmonte, L. Holberg, and J. Kahn, “Effect of atmospheric turbulence on timing instability for partially reciprocal two-way optical time transfer links,” Phys. Rev. A 101, 033843 (2020).
  * Valencia  _et al._ (2004b) Alejandra Valencia, Giuliano Scarcelli, and Yanhua Shih, “Distant clock synchronization using entangled photon pairs,” Applied Physics Letters 85, 2655 (2004b).
  * DAR (2019) “gps.gov,” https://www.gps.gov/governance/advisory/meetings/2019-06/burke.pdf (2019), [Accessed 18-11-2023].
  * Oelker  _et al._ (2019) E. Oelker, R.B. Hutson, C.J. Kennedy, L. Sonderhouse, T. Bothwell, A. Goban, D. Kedar, C. Sanner, J.M. Robinson, G.E. Marti, D.G. Matei, T. Legero, M. Giunta, R. Holzwarth, F. Riehle, U. Sterr, and J. Ye, “Demonstration of $4.8\times 10^{-17}$ stability at 1 s for two independent optical clocks,” Nature Photonics 12, 714–419 (2019).
  * Anwar  _et al._ (2022) Ali Anwar, Chithrabhanu Perumangatt, Aitor Villar, Alexander Lohrmann, and Alexander Ling, “Development of compact entangled photon-pair sources for satellites,” Applied Physics Letters 121 (2022), 10.1063/5.0109702.
  * Lohrmann  _et al._ (2020) Alexander Lohrmann, Chithrabhanu Perumangatt, Aitor Villar, and Alexander Ling, “Broadband pumped polarization entangled photon-pair source in a linear beam displacement interferometer,” Applied Physics Letters 116 (2020), 10.1063/1.5124416.
  * Natarajan  _et al._ (2012) Chandra M. Natarajan, Michael G. Tanner, and Robert H. Hadfield, “Superconducting nanowire single-photon detectors: physics and applications,”  (2012), 10.1088/0953-2048/25/6/063001, arXiv:1204.5560 .
  * Deng  _et al._ (2007) Z. J. Deng, M. Feng, and K. L. Gao, “Preparation of entangled states of four remote atomic qubits in decoherence-free subspace,” Phys. Rev. A 75, 024302 (2007).
  * Pan  _et al._ (2012) Jian-Wei Pan, Zeng-Bing Chen, Chao-Yang Lu, Harald Weinfurter, Anton Zeilinger, and Marek Żukowski, “Multiphoton entanglement and interferometry,” Rev. Mod. Phys. 84, 777–838 (2012).
  * Krenn  _et al._ (2017) Mario Krenn, Armin Hochrainer, Mayukh Lahiri, and Anton Zeilinger, “Entanglement by path identity,” Phys. Rev. Lett. 118, 080401 (2017).
  * Bergamasco  _et al._ (2017) N. Bergamasco, M. Menotti, J. E. Sipe, and M. Liscidini, “Generation of path-encoded greenberger-horne-zeilinger states,” Phys. Rev. Appl. 8, 054014 (2017).