%% APPENDIX %%

\section{Pseudothermal source and experimental setup}
\label{appendix:setup:detailed}

\begin{figure*}
  \centering
  \includegraphics[width=0.8\linewidth]{figures/setup/20240627_source_multiemitter.eps}
  \caption{\label{setup:multiemitter}
  Detailed experimental setup using a pseudothermal bunched light source. BS: 50:50 beamsplitter; Att: attenuator; HWP: half-wave plate; QWP: quarter-wave plate; APD: avalanche photodiode.
  }
\end{figure*}

% \subsection{Optical setup}

The primary source used in the frequency tracking measurement is a pseudo-thermal light source, where light from a laser running above the lasing threshold is split into two separate paths, one of which is delayed beyond the coherence time of the laser. If the two paths are perfectly indistinguishable (i.e., in polarization, spatial mode, and brightness), they act as independent light emitters and exhibit bunching effects.

From the Wiener-Khintchine theorem~\cite{khintchine1934correlations}, the Lorentzian spectral profile from a coherent laser light source corresponds to a $g^{(1)}$ with a Laplacian timing profile of the form,

\begin{equation*}
g^{(1)}(\tau) := \frac{\langle{}{E^*(t)E(t+\tau)}\rangle{}}{\langle{}E^*(t)E(t)\rangle{}}
= \exp\left(-\frac{|\tau|}{\tau_c}\right)
\end{equation*}

\noindent{}\\
as a function of coherence time $\tau_c$.

Thermal light can be modelled as a large collection of independent emitters of light~\cite{loudon2000quantum}. In the case of $\nu$ emitters, we have
\begin{align*}
\langle{}E^*(t)E(t+\tau)\rangle = \nu{}\langle{}E_i^*(t)E_i(t+\tau)\rangle
\end{align*}
and
\begin{align*}
\langle{}E^*&(t)E^*(t+\tau)E(t+\tau)E(t)\rangle \\
&= \nu{}\langle{}E_i^*(t)E_i^*(t+\tau)E_i(t+\tau)E_i(t)\rangle \\
&+ \nu{}(\nu{}-1)\left[\langle{}E_i^*(t)E_i(t)\rangle + \left|\langle{}E_i^*(t)E_i(t+\tau)\rangle\right|\right],
\end{align*}
% \noindent{}\\
and we can recover the Siegert relation~\cite{siegert1943fluctuations} and determine the $g^{(2)}$ of the light source

\begin{align*}
g^{(2)}(\tau) &:= \frac{\langle{}E^*(t)E^*(t+\tau)E(t+\tau)E(t)\rangle}{\langle{}E^*(t)E(t)\rangle^2} \\
&= 1 + \frac{\nu-1}{\nu}\left|g^{(1)}(\tau)\right|^2 \\
&= 1 + \frac{1}{2}\exp\left(-\frac{2|\tau|}{\tau_c}\right),\quad{} \nu=2,
\end{align*}

\noindent{}\\
which saturates at a value of $g^{(2)}(\tau=0)=1.5$.

Our experimental setup, per Fig.\,\ref{setup:multiemitter}, uses light from a distributed-feedback (DFB) laser diode of wavelength 780\,nm, coupled into 780-HP fiber to project into a single optical mode. The coherence length of the laser is approximately 200\,m; a delay fiber of length 400\,m was used in the delay arm, corresponding to a 2\,\textmu{}s delay. Polarization rotation in the delay arm is compensated using a free-space link with a HWP and QWP, while attenuation is achieved by slight decoupling of the collection mode, before recombination using a fused-coupler beamsplitter (BS). The $g^{(2)}(\tau=0) = 1.44(3)$ achieved in the setup is attributed to the fact that the laser is not fully coherent~\cite{yeo2024crossg2} and other experimental imperfections. % TODO: Verify the coherence length of the laser.

A second BS is used to further split the light into two symmetrical channels for downstream detection using fiber-pigtailed active-quenched Si avalanche photodiodes (Excelitas SPCM-800-10-FC). Attenuation in each channel is achieved by cascading fiber beamsplitters as well as fiber decoupling at the mating sleeves, and is measured to be relatively stable over multiple days. Fiber spools for delay was not used to avoid potential complications from timing delays attributed to fiber length changes.

Time tagging is subsequently performed by timestamp devices (S-Fifteen Instruments TDC2) with 4\,ps nominal timing resolution and 20\,ps 1-$\sigma$ timing jitter.
The clock to each timestamp device is supplied by external 10\,MHz clock distribution boards with quartz crystal oscillators and without any onboard temperature stabilization.

% The powermeter may have underestimated the power supplied.

% \subsection{Software setup}
% \label{appendix:setup:software}

% The timestamps are fed via a USB 2.0 uplink to two separate computers, which supplements coarse timing accuracy using the Network Time Protocol (NTP) on the order of milliseconds. The two computers rely on different master NTP servers for timing synchronization to emulate more realistic network conditions.
% Peak finding and tracking is subsequently performed to obtain a time difference with nanosecond resolution. The bidirectional communication for coincidence peak agreement is performed over direct TCP/IP sockets, using \texttt{qcrypto} with the \texttt{fpfind}~\cite{fpfind} add-on for inline software-based active frequency compensation. Each pairwise event corresponds to a timing difference with nanosecond resolution.

% The software responsible for peak tracking in \texttt{qcrypto} is \texttt{costream}, which was supplied with the parameters \texttt{-w 1024 -u 1024 -Q 50000 -a 4096} that accounted for the relatively broad coincidence peak. \texttt{costream} performs a two-pointer search for pairwise coincidence events, within a coincidence tracking window of width 256\,ns to account for frequency drifts.

% Talk abit about the encoding of the timestamp events and the Ethernet link as well.
% Timestamps are 64-bits blah.
% Initial time synchronization accuracy of 10~ms facilitated by NTP, before further refinement by the peak finding program.
% We make an additional small note that the accuracy of the initial coarse timing synchronization is also important, and constraints the choice of $N\delta{}t > 2\tau$ in order for the measured timing difference after the circular convolution to be accurate.


% \subsection{Concurrent frequency measurement}

The actual frequency offset between the two clocks was measured concurrently with the clock synchronization experiment by mixing separate copies of the $10$\,MHz signals and passing through a low pass filter with $4$\,GHz cut-off frequency. This mixed signal was sampled by an oscilloscope at a rate of $2.5$\,ksamples/s over $10$\,s, before performing a $2^{14}$-bin FFT (nominal resolution of approximately $0.15$\,Hz) with a Hann window.

The frequency offset error is then calculated by measuring the difference between the servoed frequency offset (from the frequency estimation step in the clock synchronization) and the measured frequency offset with linear interpolation. The histogram of offset errors is fitted using a Gaussian probability density function, obtaining a standard deviation of $3.2$\,ppb. % 3.23506471

% This is not an aircon transition.
% The spike and

% This standard deviation is dependent on presence of the laboratory air-conditioning, increasing to $3.6$~ppb between the $10$~hour and $22$~hour mark when the air-conditioning is switched off, and decreasing to $2.6$~ppb when the air-conditioning is switched on. This can be ascribed to the larger  temperature fluctuations in the environment resulting in a stronger frequency drifts in the quartz oscillators.

% %
% % 38 seconds
% This measurement was started on Monday, 26 Aug 2024 11.33am, with a mode hop observed at the 11h 38min mark. The transition was at 11.11pm.
% The second one occured at 23h 10m 28s, so at time 10.43m % ()
%  8 hour mark = 7.30pm
% 22 hour mark = 9.30am






























% \newpage
% \section{Using FFT for cross-correlation}
% \label{appendix:fft}

% Identification of the timing delay between two parties involves histogramming of the time difference between all timestamps measured by both parties, and then identifying the cluster of time difference. After binning in the histogram, this appears as a peak in the histogram. For a set of $m$ and $n$ timestamps,

% \begin{align}
% (f \star g)[n] = \sum_m f[m]g[n+m]
% \end{align}

% \noindent\\
% This has a time complexity of $\mathcal{O}(mn)$.

% A more computationally efficient method to calculate the time lag relies on the circular convolution theorem to compute the discrete cross-correlation using the Fourier Transform, i.e.

% \begin{align}
% (f \star g)[n] = \mathcal{F}^{-1}\left\{
%   \overline{\mathcal{F}\{f\}} \cdot{} \mathcal{F}\{g\}
% \right\}[n]
% \end{align}

% In discrete space, the discrete Fourier Transform can be efficiently implemented as a Fast Fourier Transform (FFT), which determines the final time complexity of the cross-correlation of $\mathcal{O}(n \log{n})$. This is more efficient for larger number of elements $n$, at the expense of $\mathcal{O}(n)$ memory.

% We note that the histogramming method works well when the number of elements to be computed is small, by constaining the search to a limited coincidence window such that $m$ is as small as possible.

% We additionally note that algorithms faster than $\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.
% which are unavailable in time randomized strings... but maybe not that true.



























% \newpage
\section{Peak finding probabilities using Poisson statistics vs Gaussian assumption}
\label{appendix:poissonprobability}

\begin{figure*}

  % Top left
  \begin{minipage}[l]{0.45\linewidth}
    \includegraphics[width=\textwidth]{figures/peakfinding/20240919_peakfinding_old_2poisson.png}
  \end{minipage}
  \hspace{0em}
  % Right
  \begin{minipage}[!r]{0.45\textwidth}
    \includegraphics[width=\textwidth]{figures/peakfinding/20240919_peakfinding_old_1gaussian.png}
  \end{minipage}
  \caption{\label{finding:gaussian}
    \textbf{Left:} Probabilities to find the correct peak position by solving Eqn.\,\ref{eqn:coincnofreq:binoverlap}, given the singles detection rate $s_1=s_2=100$\,kcounts/s, coincidence rate $c=650$\,counts/s, and bin overlap of $\nu=0.5$. Since there is no frequency offset, the peak finding probabilities reaches 100\% with increasing number of FFT bins.\\
    \textbf{Right:} Probabilities to find the correct peak position using the normal approximation for $X_i$ in Eqn.\,\ref{eqn:coincnofreq:binoverlap}, with the same parameters.
  }
\end{figure*}

The corresponding probabilities for Poisson model and the corresponding Gaussian approximation model are presented in Fig.\,\ref{finding:gaussian}. We find that the previous models assuming Gaussian-distributed bin values incorrectly predict the success probability of peak finding: in particular, it generally underestimates the probability for low time bin widths, and overestimates for small number of FFT bins.

FFT can be efficiently implemented if the number of bins $N$ are either in powers of 2, or has factorizations with small primes. Here we use $N=2^q, q\in\mathbb{Z}^{+}$ for its logarithmic scale.




























% \newpage
\section{Derivation of max-order distribution}
\label{appendix:derivation}

Given that the FFT is performed with N bins, the accidental coincidences in each of these bins follow the Poisson distribution with mean $\lambda{}$. The maximum value across all these bins increases with an increasing number of bins. The distribution of the max-order statistic (i.e., the maximum of all bins) with respect to the number of coincidences $k$ and number of bins $N$ of mean value $\lambda{}$ is denoted $f_{(N)}(k|\lambda)$.

Since each bin is assumed to be independent and identically distributed, the corresponding max-order cumulative distribution function (CDF) can be expressed as the product of the CDF of individual bins,
\begin{align*}
F_{(N)}(k|\lambda{}) = \underbrace{F(k|\lambda{})\times\ldots{}\times{}F(k|\lambda{})}_{N\ \text{times}} = \left[F(k|\lambda{})\right]^N.
\end{align*}

We therefore express the max-order probability mass function (PMF) in terms of the Poisson PMF and CDF of a single bin,
\begin{align*}
  f_{(N)}(k|\lambda)
  &= F_{(N)}(k|\lambda{}) - F_{(N)}(k-1|\lambda{}) \\
  &= \left[F(k|\lambda{})\right]^N - \left[F(k-1|\lambda{})\right]^N \\
  &= \left[F(k|\lambda{})\right]^N - \left[F(k|\lambda{}) - f(k|\lambda{})\right]^N.
\end{align*}

Notably, this form is already tenable for direct computation without further simplification, even though computing the difference of $n^{\text{th}}$-powers generally causes catastrophic cancellation due to floating-point rounding errors.
The max-order Poisson PMF distribution width heuristically scales roughly with $\sqrt{\lambda{}}$, which requires the PMF to be accurate to at least $1/\sqrt{\lambda{}}$, e.g., $\lambda{}\le 10^{4}$ needs at least $\sim{}10^{-2}$ accuracy.
Since $F(k|\lambda{}) \in [0,1]$ and the fact that exponentiation can be easily applied for large $N$ using numerical techniques (such as exponentiation-by-squaring), the rounding errors can be minimized to near the floating-point precision, e.g., $\sim{}10^{-16}$ for 64-bit floats.

For larger $\lambda{} > 10^{4}$ (i.e., high accidental coincidence rates per time bin), the normal approximation for each bin remains appropriate,
\begin{align*}
f(k) \sim{} \text{P}(k|\lambda{}) \approx{} \mathcal{N}(x=k|\mu=\lambda{},\sigma=\sqrt{\lambda{}}),
\end{align*}
and the corresponding max-order probability distribution function (PDF) follows a more tractable form for numerical computation,
\begin{align*}
  f_{(N)}(x|\mu,\sigma)
  &= \frac{d}{dx} F_{(N)}(x|\mu,\sigma) \\
  &= \frac{d}{dx} \left[F(x|\mu,\sigma)\right]^N \\
  &= N f(x|\mu,\sigma) \left[F(x|\mu,\sigma)\right]^{N-1}.
\end{align*}
This in fact corresponds to the first order term for the discrete case after a binomial expansion, noting that $f(k|\lambda{}) \le F(k|\lambda{})$, with
\begin{align*}
  f_{(N)}&(k|\lambda) \\
  &= \left[F(k|\lambda{})\right]^N - \left[F(k|\lambda{}) - f(k|\lambda{})\right]^N \\
  &= \left[F(k|\lambda{})\right]^N - \left[\left[F(k|\lambda{})\right]^N \right. \\
  &\quad \left. - N f(k|\lambda{})\left[F(k|\lambda{})\right]^{N-1} + \ldots\right] \\
  &= N f(k|\lambda{})\left[F(k|\lambda{})\right]^{N-1} - \mathcal{O}\left(f^2F^{N-2}\right).
\end{align*}

% Not too relevant...
% We further note that, under the statistical significance framework~\cite{ho2009clock} of normalizing the bin distributions (into the standard normal distribution), the associated peak finding probability can be directly computed as an integral, given the desired statistical significance threshold $S_\text{th}$. This is given by
% \begin{align*}
% \text{Prob}\left[\text{Peak}>S_\text{th}\right] = \int_0^\infty{} dz \int_{-\infty{}}^\infty{} dx\; Nf(x) F(x)^{N-1}\cdot{} \phi(x - z - S_\text{th}).
% \end{align*}






























% \newpage
\section{Peak tracking derivation}
\label{appendix:math:tracking}

\begin{figure*}
  \centering
  \includegraphics[width=0.65\linewidth]{figures/peaktracking/20240627_tracking_50h.png}
  \caption{\label{result:longertracking}
    Tracking measurement over 50 hours, similar to that of Fig.\,\ref{fig:result:main}.
  }
\end{figure*}

We perform active frequency compensation by estimating the clock frequency offset from the set of timestamps.
The timing difference $\tau_i := t_i - t'_i$, between a timestamp pair $\{t_i,t'_i\}$ from both parties, is continuously served by searching for photon pair detection events within a prescribed coincidence tracking window.

The frequency offset $\Delta{}u$ is given by the ratio between measured elapsed time $\Delta{}t_i := t_i - t_{i-1}$ with respect to some reference elapsed time, in this case the elapsed time measured by the other peer $\Delta{}t'_i$,

\begin{align*}
\frac{\Delta{}t_i}{\Delta{}t'_i} = 1 + \Delta{}u.
\end{align*}

Rewriting in terms of the measured successive timing difference, we can estimate the frequency offset by measuring the rate of change in the timing difference, i.e.,

\begin{align*}
\Delta{}u_i = \frac{\tau_i - \tau_{i-1}}{\Delta{}t'_i}.
\end{align*}

Active frequency compensation is therefore achieved by performing a timing correction for each timestamp $t_i$ using the estimated frequency offset,

\begin{align}
t_i &\rightarrow{} \Delta{}t_i ( 1 + \Delta{}u) + t_{i-1} = t_i + \Delta{}t_i \Delta{}u,
\end{align}

\noindent{}\\
where the overall frequency offset $\Delta{}u$ accumulated is also given by
$\Delta{}u = \prod_0^i{\left(1 + \Delta{}u_i\right)} - 1$.

The peak tracking is intrinsically noisy due to the signal being dominated by accidental coincidences.
We use an exponential moving average filter given by,
\begin{align*}
\tau{}'_i = \alpha\tau_{i} + (1 - \alpha)\tau'_{i-1},\quad{} \tau'_0 = \tau_0,
\end{align*}
which behaves like a low pass filter to smooth the signal. The cofficient $\alpha$ represents the relative weight of timing offset $\tau_i$ with respect to the accumulated average $\tau'_{i-1}$, and can be adjusted depending on the observed signal rate.

Given a unit step impulse, the time it takes to reach $1 - 1/e$ of the signal is associated with a time constant $\beta$ related to the coefficient $\alpha$,
\begin{align*}
\alpha = 1 - \exp{\left(-\frac{\overline{\Delta{}t}}{\beta}\right)} \approx{} \frac{\overline{\Delta{}t}}{\beta},\qquad{} \beta \gg \overline{\Delta{}t},
\end{align*}
after Taylor expansion of the exponential, where $\overline{\Delta{}t}$ is the average separation time between consecutive timestamp events. The singles rate $s$ in our experiment is 200\,kcounts/s, so $\overline{\Delta{}t} = 1/s = 5$\,\textmu{}s. We set the time constant $\beta = 50$\,ms, which is 4 orders of magnitude longer than $\overline{\Delta{}t}$ for averaging.






% \graphic[]{figures/tracking/20240801_simulation.png}{
%   Simulation of peak tracking.
% }



% \section{Peak searching}
% \label{appendix:searching}

% This cross-correlation is efficienctly performed using Fast Fourier Transform (FFT), relying on the circular convolution theorem. More details in Appendix~\ref{appendix:fft}. The choice of parameters dominates whether the cross-correlation yields a peak, specifically number of time bins $N = 2^q$, and initial width of each bin $\delta{}t$.
% An exhaustive 2D search of this parameter space will thus derive whether a source can be used for clock synchronization.







% \section{Peak tracking}
% \label{appendix:tracking}



% \section{Alternative setup}

% Alternative setups were used.

% \graphic[label=setup:all]{figures/setup/20240627_sources.pdf}{
%   Three different light sources used to test the model. \textbf{Top:} Multi-emitter source, created by an interferometer delay exceeding the coherence length of the laser. \textbf{Middle:} Sub-threshold laser thermal source, with the spectrum narrowed down using a grating monochromator and etalon. \textbf{Bottom:} Spontaneous parametric downconversion pair source, which generates photon pairs. BS: Beamsplitter, DM: Dichroic mirror, ppKTP: Periodically-poled potassium titanyl phosphate crystal.
% }

% \textbf{Subthreshold laser source:} Our setup is relatively simple \ref{setup:all} with the thermal light generated using a subthreshold distributed-feedback (DFB) laser diode at 780nm. It is subsequently filtered with a combination of a grating monochromator and a temperature-tuned etalon to increase its coherence time. The spectral profile is determined by the etalon which uses a fused silica (Suprasil 311) substrate with a tuning parameter of 4GHz/K. The optical coating has a reflectivity of 97\% on both surfaces giving a finesse of 103. The etalon has a thickness of 1mm resulting in a free spectral range of 103GHz and a bandwidth of 1GHz FWHM corresponding to a coherence time of around 1ns~\cite{tan2014measuring}. A polarizer is also used to clean up the laser mode. We use a 50:50 fiber beam splitter is then used to split and send the light to two parties each with their own single-photon detectors and FPGA-based timestamp devices. Both detectors are actively-quenched avalanche photodiodes with an efficiency of about 15\% at 780nm with a timing jitter of around 40ps. The timestamp devices have a resolution of 4ps with a timing jitter of around 20ps with timings derived from a common quartz oscillator. Cross correlation is then done using the recorded time stamps to generate the unnormalized $\text{g}^{(2)}(\tau)$.

% \textbf{Multi-emitter coherent laser source:} \note{Description goes here...}

% \textbf{Spontaneous parametric down-conversion (SPDC) pair source:} \note{Description goes here...}

% \graphic[label=setup:thermal:g2]{figures/setup/g2_thermal_nanowire.pdf}{
%   Second order correlation measurement of the thermal source using two superconducting nanowire detectors. Clock synchronization is done by tracking the change in peak position.
% }



% \graphic[label=setup:multiemitter2]{figures/setup/thermalsetup-eps-converted-to.pdf}{
%   The thermal light source is a sub-threshold laser diode, filtered by a narrow bandpass filter and two temperature-controlled etalons of different thickness.
% }


% \graphic[label=setup:pairs]{figures/setup/20231123_source_qkd_horizontal_boxed.pdf}{
%   \note{Placeholder, to update with better diagram}. Correlations using pair source.
% }


































% \vspace{2em}
% \newpage
\section{Longer time scale measurement}

Fig.\,\ref{result:longertracking} shows a longer experimental measurement spanning over 50\,hours, with a similar drift behaviour observed in the reconstructed frequency offset. The frequency offsets were not measured concurrently, so no tracking accuracy estimation is available for this dataset.

% Why this wasn't part of the main figures
Measurement was performed using a different pair of Si APDs instead (Micro Photon Devices SPAD), which had an order of magnitude lower system efficiency and resulted in a lower $-90$\,dB attenuation channel possible at the same detection rate parameters.
The measurement was also performed using suboptimal parameters for peak tracking that dropped off too many coincidences.

The timestamps were disciplined with onboard voltage-controlled temperature-compensated quartz crystal oscillators, so the frequency offset between the two timestamps generally remained the same, as well as nearly an order of magnitude better frequency stability.

% The measured frequency difference is low on the order of 10\,ppb due to the fact that the two timestamps shared the same air conditioned environment, but this technique will still work for frequency drifts larger than 100\,ppb.
%Peak tracking will start to fail at frequency drifts larger than X, we intend to resolve this issue.













































% \newpage
% \section{A comment on the number of coincidences reported by costream}
% \label{appendix:costream}

% The expected number of background accidentals, considering the total\dots
