% Journals to target:
% - PR Applied / PRX Quantum
% - APL
% - NJP
% - RSI

\documentclass[
    pra,
    superscriptaddress,
    twocolumn,
    aps,
    10pt,
]{revtex4-2}

% Imported packages
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{mathtools}
\usepackage{amsfonts}
% \usepackage[figure,figure*]{hypcap}  % let references link to figure instead of caption

\begin{document}

\title{Clock Synchronization without Frequency References under\\Weak Timing Correlations}

\author{Justin Yu Xiang Peh}
\author{Darren Ming Zhi Koh}
\affiliation{Centre for Quantum Technologies, 3 Science Drive 2, Singapore 117543}
\author{Zifang Xu}
\affiliation{Department of Physics, National University of Singapore, 2 Science Drive 3, Singapore, 117542}
\author{Xi Jie Yeo}
\author{Peng Kian Tan}
\affiliation{Centre for Quantum Technologies, 3 Science Drive 2, Singapore 117543}
\author{Christian Kurtsiefer}
\affiliation{Centre for Quantum Technologies, 3 Science Drive 2, Singapore 117543}
\affiliation{Department of Physics, National University of Singapore, 2 Science Drive 3, Singapore, 117542}

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

%\keywords{Suggested keywords}%Use showkeys class option if keyword
%display desired

\maketitle
\textbf{
Clock synchronization is necessary for communication and distributed computing tasks.
Previous schemes based on photon timing correlations use pulsed light or photon pairs for their strong timing correlations.
In this work, we demonstrate successful synchronization of quartz clocks using weakly time-correlated photons of 180\,ns coherence time from a pseudothermal bunched light source.
A synchronization timing jitter of 10\,ns is achieved over symmetric -102\,dB optical loss channels between two parties, over a span of 25\,hours.
We also derive a model to accurately estimate the coherence peak searching success probabilities.
}\\

\noindent

% Future directions:
% - PIN diode w/ thermal in same spectrum
%
% Justifications:
% - Pulse energy high


% Things to put inside:
% 1. Experimental setup picture
% 2. g(2) of 780 nm
% 3. Current vs g(2) height
% 4. Reconstruction of peak
% 5. Allan deviation to compare atomic clock with optical oscillator

% Add mini-page if necessary
% \begin{minipage}[c]{0.45\linewidth}
% \end{minipage}































\section{Introduction}

% Clock synchronization is necessary for communication and distributed computing tasks. Previous schemes based on measuring the timing correlation of light utilise strongly time-correlated photon pairs generated from spontaneous parametric down-conversion (SPDC). In this work, we demonstrate synchronization of quartz clocks using weakly time-correlated photons, after an initial coarse timing synchronization over Network Time Protocol (NTP). A synchronization accuracy of $10$~ns is achieved over symmetric optical loss channels of -$101.8(0.4)$~dB on each side, using photon bunching of $180$~ns coherence time.
% Remote clocks can be synchronized by sharing a timing correlated photon source between parties and measuring its second-order coherence $g^{(2)}(\tau)$, with the correlation peak timing corresponding to the timing delay $\tau_0$ between clocks. This has been demonstrated with highly time-correlated photon pairs. We show that light sources with weaker timing correlation, such as photon bunching in thermal light, can be used for clock synchronization as well, despite their low correlation peak of $g^{(2)}(\tau) \le{} 2$.

% Why do clock synchronization, and how is it currently done?
Clock synchronization is used in everyday tasks such as navigation and distributed computing. This is commonly implemented using the Network Time Protocol or global navigation satellite system (e.g., GPS) time synchronization, achieving precision of milliseconds or tens of nanoseconds, respectively~\cite{mills2006computer,renfro2022analysis}. Quantum communication protocols also require clock synchronization, but on the order of nano- to pico-seconds, and is typically achieved either using pulsed laser light~\cite{giorgetta2013optical} or with classical signals on dedicated optical or electronic channels~\cite{gerritsWhiteRabbitassistedQuantum2022,alshowkanSynchronizingQuantumLocal2022}.

Modern systems can use the resources typically in the protocol itself to perform the clock synchronization, such as photon pairs from spontaneous parametric down-conversion (SPDC)~\cite{lee2019symmetrical,ho2009clock}. This is possible with the use of frequency standards such as Rubidium (Rb) clocks which provide a long term frequency stability of \textless1\,ppb/day. In comparison, crystal oscillators without temperature stabilization have frequency stability of only 100\,ppb/day.

To remove the dependency on Rb frequency standards, it was proposed to find and track the timing difference between photon pairs due to their strong timing correlation on the order of picoseconds~\cite{ho2009clock,spiessClockSynchronizationCorrelated2023}, by identifying an initial coincidence peak with low timing resolution then iteratively apply frequency corrections to improve resolution. Weak coherent photon pulses were also proposed for clock frequency transfer~\cite{spiessClockSynchronizationPulsed2023}. In both cases, the cross-correlation peak is very strong, i.e., $g^{(2)}(\tau=0) \gg 2$, with resolution generally limited only by the timing jitter of the generation and detection optoelectronics.

% Our proposition: Clock synchronization with thermal light, and its benefits
Other sources of timing-correlated light are thermal light sources which has been used in application such as ghost imaging~\cite{bennink2002two} and range finding~\cite{tan2023practical}.
These timing correlations arise from temporal photon bunching, also known as the Hanbury-Brown-Twiss effect~\cite{brown1956correlation}.
In particular, the use of thermal light opens up the potential for distributed clock synchronization due to photon bunching being preserved across arbitrary partitioning.

However, unlike photon pairs, identification of this peak is more challenging due to the low signal (i.e., $g^{(2)}(\tau)\le{}2$) as well as fluctuations about the background of $g^{(2)} = 1$. There is an earlier proposal that suggests the use of thermal light for clock synchronization, relying on the use of low efficiency two-photon absorption in single photon detectors to resolve the bunching characteristic~\cite{zhu2013new,xu2020determination}; such a scheme has yet to be demonstrated.

% What sexy results did we achieve?
Here, we demonstrate clock synchronization using one such weak timing correlation light source, i.e., light from pseudothermal source with a $g^{(2)}(\tau)=1.44$ and coherence time $\tau_c = 180(6)$\,ns, over a \mbox{-102}\,dB transmission channel with an accuracy of 10\,ns, and continuously track the frequency drift over a period of 25\,hours to demonstrate its stability.
Our scheme only requires single photon detection and crystal oscillators as reference clocks.
We additionally derive the probability of achieving clock synchronization over transmission channels, which is applicable across different time-correlated photon sources, including SPDC pair sources and thermal sources.

% When performing quantum communication protocols, the existing quantum resource can be used as well for clock synchronization, can rely on quantum resources However, uses such as Quantum Communication with a high clock rate require a higher level of precision on the order of picoseconds\cite{li2023high}. This is typically done through the exchange of timing correlated light such as pulsed light or photon pairs from spontaneous parametric down-conversion (SPDC)\cite{}. However, these optical systems can be complex and costly to operate due to the many components involved.

































\section{Clock synchronization with bunched light}



% np.max(np.abs(np.diff(ys))) = 17.695
% np.std(ys) = 7.767015788677711
\begin{figure*}

  % Top left
  \begin{minipage}[l]{0.35\linewidth}
    \vspace{1em}
    \includegraphics[width=\textwidth]{figures/setup/20240627_source_multiemitter_simple.eps}
    \vspace{0.5em}

    \includegraphics[width=1\textwidth]{figures/setup/20240903_g2_dataset024a1_v1.png}
  \end{minipage}
  \hspace{0em}
  % Right
  \begin{minipage}[!r]{0.60\textwidth}
    \includegraphics[width=\textwidth]{figures/peaktracking/20240827_tracking_50h.png}
  \end{minipage}
  \caption{\label{fig:result:main}
    \textbf{Top left:} Simplified experimental setup for clock synchronization. Bunched (pseudothermal) light is sent to two different parties through separate -102\,dB channels, before detection by Silicon avalanche photodetectors (Si APDs) and timestamping by time taggers running on independent quartz crystal oscillators. Singles count rate $s_1$ and $s_2$ of approximately 200\,kcounts/s are recorded on each side. BS: beamsplitter; Att: attenuators.\\
    \textbf{Bottom left:}
    Plot of $g^{(2)}(\tau)$ of the pseudothermal light source, with the curve fit in grey corresponding to $g^{(2)}(0)=1.44(1)$ peak and coherence time $\tau_c=180(6)$\,ns. Error bars correspond to Poissonian errors from counting statistics. \\
    \textbf{Right:} Long-term trace of frequency offset, timing offset, and event rates during a $25$\,hour clock synchronization run between independent time taggers running on separate quartz crystal oscillators, after an initial frequency correction of $4.0$\,ppm. Correlation peak tracking is performed using a $256$\,ns coincidence window. The sharp spike and dip in event rates are attributed to laser mode hopping.
  }
\end{figure*}

% Main result achieved
Timing synchronization was achieved using bunched light between two remote parties over symmetrical loss channels, each with \mbox{-102}\,dB average transmission, for more than 25\,hours using only simple crystal oscillators, see Fig.\,\ref{fig:result:main}.

% Source description
The source of bunched light is a pseudothermal source~\cite{yeo2024twoemitter}, established using a laser of 780\,nm wavelength sent into an unbalanced Mach-Zehnder interferometer with optical delay longer than the laser coherence time $\tau_c$; phase fluctuations in the laser manifests the bunching effect, see Appendix~\ref{appendix:setup:detailed} for a detailed explanation.
This results in an intensity correlation that can be measured through the second-order coherence function $\text{g}^{(2)}(\tau)$, given by
\begin{align}
  \text{g}^{(2)}(\tau)=1+\frac{1}{2}\exp\left(-\frac{2|\tau|}{\tau_{\text{c}}}\right)
  \label{eqn:g2source}
\end{align}
due to the Lorentzian spectral profile of the laser, where the coherence time $\tau_{\text{c}}$ is inversely proportional to the frequency spectral width $\Delta f$ of the source~\cite{loudon2000quantum}.
The brightness of the source is 1.55(2)\,mW, with a measured coherence peak of $g^{(2)}(\tau=0)=1.44(1)$ and coherence time of $\tau_c = 180(6)$\,ns.

The light from the source is shared between $n=2$ parties with optical fiber channel tranmissions of $-101.9(4)$\,dB and $-102.2(4)$\,dB using optical attenuators.
Photon arrival events are detected using Silicon avalanche photodetectors (APDs) on each side, with count rates of 192\,kcounts/s and 182\,kcounts/s without correcting for dark counts and afterpulsing.

% servo -> report / served
The detection events are read by independent time taggers disciplined by different free-running 10\,MHz quartz crystal oscillators.
The generated timestamps are then exchanged for coincidence peak finding and tracking in real-time.
The peak tracking program continuously serves a timing offset between each party by pairing photon detection events between both parties within a $256$\,ns coincidence window, over the span of the 25\,hour measurement.
The average coincidence rate of $8,500$\,events/s is consistently higher than the accidental count rate of $7,000$\,events/s, indicative of successful frequency tracking.

Due to the presence of clock frequency drift, typically from temperature fluctuations and electronic noise, the frequency offset between the clocks changes over time as well. We are able to reconstruct this offset, by monitoring the drift in served timings using individual samples spanning 10.74\,s, with a resolution of 0.537\,s, shown in Fig.\,\ref{fig:result:main}.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{figures/peaktracking/20240923_freqtracking_measured.png}
  \caption{
    \label{result:freq:accuracy}
    A concurrent measurement of the actual frequency offset between the two external $10$\,MHz clocks, as well as the corresponding error in the reconstructed frequency offset. The accuracy of the frequency servoed by the peak tracking algorithm is $3.2$\,ppb on average.
  }
\end{figure}

The actual frequency offset was obtained via a heterodyne measurement using identical copies of the clock signals with integration time of 10\,s, shown in Fig.\,\ref{result:freq:accuracy}.
A maximum instantaneous frequency difference of 35.4\,ppb, corresponding to an average drift of 3.3\,ppb/s, was observed.%TODO: Get the actual value!

We estimate the error in the reconstructed signal by performing linear interpolation on the measurement and subsequent differencing.
The reconstructed frequency offset is found to be in good agreement with the measured frequency offset, with a root-mean-squared error of 3.2\,ppb averaged over the 25\,hour span.


































\section{Peak searching}

Clock synchronization can be decomposed into two parts: first identifying the initial clock frequency and time difference between two parties, then tracking of said peak to monitor the timing drift.
The former is performed by distributing photons with time-correlated statistics to each party, and finding the coincidence peak using cross-correlation to identify the time delay.

Efficient peak identification relies on the circular convolution theorem to compute the cross-correlation between two sets of detection timestamps $a[k]$ and $b[k]$,
\begin{align}
    g^{(2)}(\tau)&\sim(a \star b)[k] = \mathcal{F}^{-1}\!\left\{
      \overline{\mathcal{F}\{a\}} \cdot{} \mathcal{F}\{b\}
    \right\}\![k],
\end{align}
under the discretization $\tau = k \delta{}t$ with $k \in \mathbb{Z}$ and time resolution $\delta{}t$, using the Fourier Transform $\mathcal{F}$ and its inverse~\cite{ho2009clock}\footnote{
  We additionally note that algorithms faster than FFT's $\mathcal{O}(n \log{n})$ time complexity are also possible~\cite{calderaroFastSimpleQubitBased2020a,krauseClockOffsetRecovery2024}, but require additional encoding of time-synchronization strings within the light pulses.
}.
The success of peak identification is constrained by the singles rate and the rate of true coincidences, and additionally requires the choice of optimal parameters for the Fast Fourier Transform (FFT) parameters for computing the time delay: specifically, the number of time bins $N = 2^q$ ($q\in\mathbb{Z}^+$), and initial bin width $\delta{}t$.

\begin{figure}
  \centering
  \includegraphics[width=\linewidth]{figures/peakfinding/20240919_peakfinding_old_3poissonfreq2.png}
  \caption{
    \label{fig:finding:surfaceplot}
    Simulated probabilities to find the correct peak position by solving the required true coincidence equation found in Eqn.\,\ref{eqn:peakfinding}, given the singles detection rate $s_1=s_2=100$\,kcounts/s, coincidence rate $c=650$\,counts/s, bin overlap $\nu=0.5$, and frequency offset error of up to $\Delta{}u=100$\,ppb.
    These parameters correspond to the setup in Fig.\,\ref{fig:result:main} with an additional $3$\,dB attenuation per channel, and shows good agreement with measurements.
  }
\end{figure}

We model the peak finding probability --- in both cases, with and without a reliable frequency reference between both parties --- and perform an exhaustive parameter scan to identify appropriate FFT parameters in Fig.\,\ref{fig:finding:surfaceplot}. Our model provides better estimations of the peak finding probability compared to previous works, by avoiding the normal approximation to the noise in the FFT (see Appendix~\ref{appendix:poissonprobability}).


\subsection{No frequency offset}

% where $\mathcal{F}$ is the $\mathcal{O}(N \log{N})$ Fast Fourier Transform (FFT). The correlation noise floor is attributed to accidental coincidences $C_\text{acc}$ in each of the $N$ detection time bin of width $\delta{}t$, and follows the probability distribution $X_{i\in\{1,\ldots{},N\}} \sim \text{Poisson}(\lambda=C_\text{acc}\equiv{}s_1s_2T\delta{}t)$ after an integration time of $T=N\delta{}t$. The required coincidence rate above background in order to be identified as the correlation peak is therefore
% \begin{align}\label{eqn:peakfinding}
%   c > \frac{1}{2kT}\left(\max\{X_i\}-X_1\right),
% \end{align}

% % \vspace{0.1ex}
% \noindent


The minimum acquisition time required for the cross-correlation is $T = N\delta{}t$ to obtain a flat cross-correlation noise floor.
This noise floor arises from accidental coincidences --- attributed to noise sources, such as coincidences with uncorrelated detection events and dark counts --- with an expected value of $C_a = (s_1 s_2 \delta{}t) T$, where $s_1$ and $s_2$ are the detection event rates on each channel.
The detection events are well-approximated by a Poisson distribution after binning, so the accidental coincidences in each time bin also follows a Poisson distribution $\text{P}$, i.e., $X_{k\in\{1,\ldots{},N\}} \sim \text{P}(\lambda=C_a)$.
% ~\footnote{TODO: This argument flow may need more refinement, and may be a little math involved.}
% Need to substantiate this... but there are papers mentioning Poissonian statistics. Intuitively, conduct Bernoulli trials, then look for coincident events, which has probability p^2. This is effectively another Bernoulli trial -> binomial distribution for number of successes -> poisson distribution under infinitisimal time steps.

The maximum observed value across all time bins $X_{(N)} \equiv \max\{X_k\}$ is therefore given by the max-order Poisson distribution (derived in Appendix~\ref{appendix:derivation}) whose probability distribution is
\begin{equation*}
  f_{X_{(N)}}(x) = \left[F_X(k|\lambda{})\right]^N - \left[F_X(k|\lambda{}) - f_X(k|\lambda{})\right]^N,
\end{equation*}
where $f_X$ and $F_X$ correspond to the probability mass and cumulative distribution functions of a single bin.

The coincidence rate above background accidentals in a single time bin (denoted $c_e$) required to be identified as the highest peak in the cross-correlation is thus
\begin{equation}
c_e > \frac{1}{T}\left(X_{(N)} - X\right).
\label{eqn:coincnofreq}
\end{equation}

% We make a small note, that the success probability of finding the coincidence peak is determined solely by the average singles rate $\sqrt{s_1 s_2}$ and the excess coincidence rate above background $c_e$, given the same reference clock supplied to both parties. We can define an analogous heralding efficiency metric $\eta_h$ commonly used in the context of a photon pair source, whose equivalent is $\eta_h = c_e/\sqrt{s_1 s_2}$.

The coincidence rate per time bin can be maximized by setting the timing resolution $\delta{}t$ to be of the same scale as the coherence time $\tau_c$, so that most of the coincidence events fall within the same time bin, i.e., $\delta{}t \sim \tau_c$.
Some of these events may fall into an adjacent time bin instead due to off-centered bins (since the exact time offset is not known \textit{a priori}): this introduces a smudge factor into Eqn.\,\ref{eqn:coincnofreq} representing the degree of bin overlap $\nu \in [0.5, 1]$, yielding
\begin{equation}
c_e > \frac{1}{\nu{}T}\left(X_{(N)} - X\right).
\label{eqn:coincnofreq:binoverlap}
\end{equation}

% There is also an upper bound, but its not very useful information in realistic scenarios:
%while being limited by the mean separation between photon detections on a single detector, i.e.
% \begin{equation}
% \delta{}t \sim \tau_c,\quad{} \delta{}t \ll \frac{1}{\max\{s_1, s_2\}},
% \end{equation}

\subsection{Variable frequency offset}

In the case of two separate clocks with slightly different clock frequencies, the timing delay between photon arrivals between each party can drift by $\Delta{}\tau_{\text{drift}} \approx T\Delta{}u(t)$ due to the non-zero clock-frequency offset $\Delta{}u(t)$ after some elapsed time $T$.
Under a sufficiently small $T$, the frequency offset can be approximated as a constant, i.e., $\Delta{}u(t) = \Delta{}u$.

In order to maximize the number of coincidences in a single time bin, the time bin should ideally be as wide as the timing drift, i.e., $\delta{}t > \Delta{}\tau_\text{drift}$, or in other words, $N < 1/\Delta{}u$. We model this as an additional smudge factor $\mu \equiv \max{\{1, N\Delta{}u\}}$, and together with Eqn.\,\ref{eqn:coincnofreq:binoverlap} yields the minimum required true coincidence rate for successful peak finding given by
\begin{equation}
c_e > \frac{1}{\mu\nu{}T}\left(X_{(N)} - X\right).
\label{eqn:peakfinding}
\end{equation}
The surface of this equation is plotted in Fig.\,\ref{fig:finding:surfaceplot}, by performing Monte Carlo simulations for a specific set of parameters and across different $N$ and $\delta{}t$.

The reduction in coincidence counts due to the presence of frequency offset can be mitigated either by choosing a smaller $N$, or by performing a frequency precompensation~\cite{spiessClockSynchronizationCorrelated2023} on the set of timestamps $\{t_i\}$ to reduce the apparent $\Delta{}u$ between the two clocks and allow for larger $N$ values. The compensation is given by the mapping
\begin{align}
  t_i \,&\rightarrow{}\, \Delta{}t_i ( 1 + \Delta{}u) + t_{i-1} = t_i + \Delta{}t_i \Delta{}u,
  \label{eqn:timingcorrection}
\end{align}
where $\Delta{}t_i \equiv t_i - t_{i-1}$ is the separation between consecutive timestamps.

In practice, while the frequency offset between quartz clocks can be high ($\sim{}\!10$\,ppm), the short-term stability of the clocks themselves are much higher ($<\!10$\,ppb/s), so after an initial precompensation reference, a small frequency precompensation step size of $100$\,ppb is sufficient to identify most peaks within 3 peak searching attempts~\footnote{Using offsets of $\Delta{}u_0$ and $\Delta{}u_0 \pm 100$\,ppb.}.

Once a initial peak has been found, the frequency and timing resolution can be further improved by repeating the respective corrections with progressively smaller time bins, until the desired resolution has been reached~\cite{ho2009clock}.
Our optimized implementations of peak finding in Python, and frequency compensation in C, are available under GPLv2 license in GitHub~\cite{fpfind}.





























\section{Peak tracking}


\begin{figure}
  \includegraphics[width=0.99\linewidth]{figures/peaktracking/20240826_noisytracking4.png}
  \vspace{-2ex}

  \includegraphics[width=\linewidth]{figures/peaktracking/20240825_accuracy4.png}
  \caption{\label{fig:tracking}
    Characterization of peak tracking accuracy under a constant frequency offset $\Delta{}u=10$\,ppb, with a time constant $\beta=50$\,ms chosen for the exponential moving average filter used in the experiment.\\
    \textbf{Top:} Measurement of reported timing offsets (blue) against the actual offset (black) over a period of 30\,seconds. The inset highlights a $\approx{}\!30$\,ns timing error over a 0.5-second window, which suggests that a histogram fit does not improve accuracy. \\
    \textbf{Bottom:} Measurements of average offset error and jitter across different time constants $\beta$, over a measurement period of 10\,minutes. Timing jitter is equivalent to the root-mean-squared error, and corresponds to 10\,ns with $\beta=50$\,ms (marked in grey dashed line). Peak tracking fails with time constants less than 5\,ms.
  }
\end{figure}


% Where does clock drift come from
The use of quartz crystal oscillators as frequency references result in changes in frequency offset over time, due to their long-term $\Delta{}u(t)$ stability of up to 100\,ppb/day.
This causes a drift in the correlation peak position over time which needs to be tracked.

Tracking can be performed by looking for coincidences within a sufficiently wide coincidence window, such that drifts in the peak are captured.
% This process is similar to the pair finding process used for pair sources.
However, when using bunched light for peak tracking, coincidences events are effectively dominated by background accidentals due to the low $g^{(2)}(\tau)$ of the light source.
Directly returning each time difference $\tau_i$ found will result in the peak being quickly lost.

% Using EMA to track the peak successfully
We apply a smoothing operation using an exponential moving average filter, to introduce damping against events far from the current estimated peak position.
This effectively minimizes noise fluctuations and allows for the peak to be tracked.

% How we quantify accuracy, what are the results, and why
In order to quantify the accuracy of the peak tracking, we supply a constant $10$\,ppb frequency offset between the two timestamps using a function generator and derive the true timing offset.
Using a time constant $\beta=50$\,ms for the moving average filter, a timing jitter of 10\,ns can be obtained, as shown in Fig.\,\ref{fig:tracking}, which is an order of magnitude better than the $\tau_c=180(6)$\,ns coherence time of the source.
The jitter was tested to remain the same even under larger frequency offsets of 50\,ppb as well.

A tracking lag is also present and can be seen in Fig.\,\ref{fig:tracking}, where the measurement best-fit line is slightly displaced from the true timing offset. This is attributed to the frequency offset being non-zero and constant, and manifests as a mean timing error of around $6$\,ns at $\beta=50$\,ms. While this lag scales proportionally to the frequency offset, this lag is expected to not contribute significantly due to the random fluctuations of the frequency offset about its mean, in the case of free-running clocks.

We also make a small note regarding the use of timing offset histograms to improve timing resolution.
For photon pair and pulsed sources, higher timing resolution can be achieved by collecting a histogram within a small timing window and extracting the peak from a normal distribution curve fit, due to the narrow peak coherence signature of $\tau_c\ll{}1$\,ns.
However, for the bunched source with a significantly noisier coherence peak signature, this technique does not improve the timing jitter, as can be seen in the served timing offset in the inset of Fig.\,\ref{fig:tracking} being consistently far from the mean.

% Keeping the frequency offset low, and why
Active frequency compensation is also performed online to keep the frequency offset small, to avoid peak tracking loss since the correlation peak is less likely to drift out of the coincidence window.
We do this by estimating the frequency offset $\Delta{}u$ from the rate of change in timing drift (see Appendix~\ref{appendix:math:tracking}), then performing timing compensation for each timestamp using the estimated frequency offset, identical to that of Eqn.\,\ref{eqn:timingcorrection}.









































\section{Conclusion}

We demonstrated successful clock synchronization between two parties over a symmetrical -102\,dB transmission channel, using a pseudothermal bunched light source of coherence time $\tau_c = 180(6)$\,ns. Peak tracking was performed online over a span of 25\,hours with active frequency compensation and an exponential moving average filter of time constant $\beta=50$\,ns, achieving an overall timing jitter of 10\,ns.

Previous papers on clock synchronization performed using photon pair sources rely on a peak significance metric for peak finding~\cite{ho2009clock,spiessClockSynchronizationCorrelated2023}, used as a threshold for quantifying the probability of the cross-correlation peak being attributed to noise.
While it remains a useful metric for estimating peak location, this underestimates the peak finding probabilities under low signal conditions with small time bins.
We develop a model that accounts for the Poisson nature of the coincidences, and calculate instead the true probability of the cross-correlation peak being the signal, as well as derive optimal FFT parameters for a given true coincidence and singles rate.
This remains applicable even when using any other sources of timing correlations, including photon pairs and thermal light.

The clock synchronization scheme in this work can be directly applied to protocols that use the photon bunching as a resource, or indirectly by compensating the reported timings with the measured timing and frequency offsets.
Direct frequency compensation can also be performed by actively correcting the quartz oscillator frequencies, e.g., by means of a voltage error signal, so that the clock signals themselves can be utilized as part of a clock distribution network.

We can also take advantage of the fact that the second-order coherence is preserved across arbitrary partitioning of the light, to distribute the signal amongst multiple parties in a star topology.
Since splitting light into two separate channels incurs an additional insertion loss of $3$\,dB per channel, clock distribution to $2^n$ parties can be achieved with only $-3n$\,dB of additional loss per channel, e.g., a 128-party setup with this source incurs about $-80$\,dB.

This work paves the way towards clock synchronization using telecommunication O-band and C-band bunched light, which will be able to propagate with minimal chromatic dispersion over longer distances, e.g., 500\,km of G.652 telecommunication fiber with 0.2\,dB/km optical loss at 1550\,nm.
This additionally opens up the possibility of using erbium-doped fiber amplifiers to amplify the correlation signal.
While timing offset changes as a result of thermal expansion in optical fibers will need to be addressed, this effect is fundamentally indistinguishable from a clock frequency drift in this scheme and can hence be compensated for in the same manner.
%In particular, the use of the $g^{(2)}$ of thermal light in a symmetric configuration makes it invulnerable to chromatic dispersion effects.


\section{Acknowledgements}
The author would like to thank Ng Boon Long for helpful discussions in setting up the concurrent frequency measurement.
% This work is not supported by any grant.

\bibliographystyle{unsrt}
\bibliography{bibliography.bib}

% \onecolumngrid
% \newpage

\appendix
\input{appendix}

\end{document}
%
% ****** End of file apssamp.tex ******