1 Introduction↩︎

Hong-Ou-Mandel (HOM) interference [1] is an intriguing quantum optical phenomenon that has become a crucial tool in many applications in quantum metrology, telecommunications, and in the wider field of quantum technology [2]. In its best-known form, it describes the behavior of two indistinguishable photons that, coupled into two separate inputs of a 50:50 beam splitter (BS), bunch together and escape from the same output port. In fact, the measurement of the average coincidences between photons exiting from different BS ports as a function of a distinguishability parameter \(\eta\) shows a drop in coincidence events registered with \(\eta\) approaching 0. For example, if the only distinguishability parameter is time, the graph of coincidences versus differences in times of arrival \(\Delta t\) shows a dip at \(\Delta t=0\). The dependence of the coincidence counts on a suitable physical parameter (e.g. the delay \(\Delta t\)) can be exploited to measure the values difference itself (e.g. \(\Delta t\)) that, in turn, is the cause of the distinguishability.

Advances in single photon detection [3], [4] and light engineering [5]–[7] have strongly improved the performances of HOM-based experiments. A pioneering work [8] has recently reported attosecond precision (4.7 as) and accuracy (0.5 as). This work, as well as other metrological applications of the HOM effect, e.g. [9], [10], have exploited an optimization of the setup informed by quantum estimation theory. In the “ordinary” HOM experiments described above, the coincidence rate approaches a constant if the dimensionless distinguishability parameter is greater than one. For example, in the case of time delay, the relevant quantity is the ratio between the delay \(\Delta t\) and the biphoton coherence time \(\tau\). Consequently, traditional HOM-based interferometry is blind to quantum interference when the two photons are distinguishable at the detectors, i.e. \(\Delta t/\tau>1\). To overcome this limitation, and due to the quantum nature of HOM interference [11], [12], an initial distinguishability of the two photon pathways can be erased right before detection, fully retrieving the interference pattern (a setup historically known as "quantum eraser" [13]–[15]). In particular, given the parameter of interest (time of arrival in our case), resolving in the conjugate variable (frequency) restores the indistinguishability between the two photons also for values of the delay \(\Delta t\gg\tau\) as a direct consequence of Heisenberg uncertainty principle. Considerable improvements in single photon detectors and time tagging electronics capabilities have indeed made it possible to devise experimental configurations where the measurements are taken by resolving the photons’ degrees of freedom. The potential of resolved coincidence measurement has been shown as a proof of principle  [16]–[18], and its reach has been extended to quantum communication [19], quantum optical coherence tomography  [20], [21], imaging [22], [23], entanglement generation [24], [25], boson sampling [26]–[28], generation and detection of structured coalescence states [29] and precise measurement [22], [30]–[32].

For the goal of precision time measurements, an approach employing frequency entangled photons (discrete photon entanglement) through estimation of biphoton beat notes has also demonstrated promising results [33]. However, only non-classical states of light [30], [33] have been investigated in this context. In particular, the use of heralded photon pairs imposes strict limitations in terms of bandwidth and number of generated photon pairs. Here we follow a different route and experimentally demonstrate a scheme recently discussed as a theoretical proposal by Triggiani et al. [34], [35]. The theory of Refs. [34], [35] shows that, for sufficiently high resolution of the detectors, resolved HOM-based measurement can be employed to estimate arbitrarily large values of the delay or the transverse displacement between the photons.

To the best of our knowledge, our work is the first demonstration of a time-delay estimation scheme with frequency-resolved HOM interference of coherent light. In addition, and in contrast to previous works, we used, as frequency-resolving detectors, dispersive gratings and array of Single-Photon Avalanche Diodes (SPADs) on both outputs rather than filters or dispersive fibers. The use of coherent light implies better performance in terms of photons bandwidth and generated photon rate if compared with heralded photons for example generated by pumping spontaneous parametric downconversion crystals. Moreover the use of SPAD arrays offer better performance in term of acquisition time, since it is possible to manage a larger photon flux. Our theoretical and experimental results show that, despite finite resolution, the frequency-resolving technique allows us to estimate time delays much larger than the coherence time of the photons without loss of sensitivity. This result demonstrates the disruptive potential of frequency-resolved HOM metrology.

2 Material and Methods↩︎

2.1 Theory↩︎

We consider a two-photon interference setup equipped with single-photon frequency-resolving detectors, where two identical photons impinge on the two faces of a balanced beam splitter at a given temporal delay \(\Delta t\). The probability of detecting the frequencies \(\omega_{1},\omega_{2}\) for photons bunching (B) or antibunching (A) at the output ports of the beam splitter is given by [34] \[P_{A/B}(\omega_{1},\omega_{2})= \frac{1}{2}f(\omega_{1})f(\omega_{2})\left(1\mp\mathcal{V}\cos((\omega_{1}-\omega_{2})\Delta t\right), \label{eq:InfRes}\tag{1}\] where \(f\) is the spectral distribution of the single photon, \(\mathcal{V}\) is the visibility of the interference fringes, and \(\Delta t\) is the delay between the photons. For simplicity, we will consider \(f\) as a Gaussian spectral distribution with variance \(\sigma^2=\frac{1}{4\tau^2}\), with \(\tau\) coherence time of the photon. Noticeably, Eq. 1 corresponds to the scenario in which the frequency resolution of the detectors, quantified by the smallest difference \(\delta_\omega\) in the frequencies that can be distinguished, is large enough so that \(\delta_\omega \Delta t\ll 1\) and \(\delta_\omega \tau\ll 1\). If this condition is not verified, the continuous probability density distribution in Eq. 1 should be integrated over the frequency range associated to each pixel of the cameras, resulting in a discrete probability mass function \[P_{A/B}^{\omega_{01},\omega_{02}}=\int_{\omega_{01}-\frac{\delta_\omega}{2}}^{\omega_{01}+\frac{\delta_\omega}{2}}\int_{\omega_{02}-\frac{\delta_\omega}{2}}^{\omega_{02}+\frac{\delta_\omega}{2}}\mathrm{d}\omega_1\mathrm{d}\omega_2\;P_{A/B}(\omega_{1},\omega_{2}), \label{eq:Discrete}\tag{2}\] where \(\omega_{0i}\), \(i=1,2\), are the central frequencies associated with two given pixels. Eq. 2 exactly reproduces the more realistic scenario of finite resolution and can be computed analytically, but takes, in general, a complicated form. We can simplify it by assuming that, within the integration intervals of width \(\delta_\omega\), \(f(\omega_i)\) is essentially constant, which is equivalent to requiring that \(\delta_\omega/\sigma\ll 1\). We can thus approximate \[P_{A/B}^{\omega_{01},\omega_{02}}\simeq \frac{1}{2}f(\omega_{01})f(\omega_{02})\int_{\omega_{01}-\frac{\delta_\omega}{2}}^{\omega_{01}+\frac{\delta_\omega}{2}}\int_{\omega_{02}-\frac{\delta_\omega}{2}}^{\omega_{02}+\frac{\delta_\omega}{2}}\mathrm{d}\omega_1\mathrm{d}\omega_2\;\left(1\mp\mathcal{V}\cos((\omega_{1}-\omega_{2})\Delta t\right),\] that, after elementary steps, reduces to \[P_{A/B}^{\omega_{01},\omega_{02}}\simeq \frac{1}{2}f(\omega_{01})f(\omega_{02})\delta_\omega^2\left(1\mp\mathcal{V}\mathop{\mathrm{sinc}}^2\left(\frac{\Delta t \delta_\omega}{2}\right)\cos((\omega_{01}-\omega_{02})\Delta t)\right), \label{eq:ResProb}\tag{3}\] where \(\mathop{\mathrm{sinc}}(x)=\sin(x)/x\). In the next section, we show that this approximation agrees with the experimental data collected. On the other hand, by integrating \(P_{B/A}(\omega_{1},\omega_{2})\) over the whole range of frequencies \(\omega_1,\omega_2\), we can evaluate the usual probabilities of observing bunching or antibunching photons in a standard HOM interference experiment with bucket detectors, yielding the well-known HOM dip \[P_{A/B}^\mathrm{NR}=\frac{1}{2}(1\mp\mathcal{V}|\mathcal{F}(\Delta t)|^2), \label{eq:PNR}\tag{4}\] where \(\mathcal{F}\) is the Fourier transform of the spectral distribution \(f\).

Before proceeding with the analysis of the precision of our estimation setup, it is interesting to understand the predicted effect of a finite frequency resolution \(\delta_{\omega}\) on the interference pattern at a more fundamental level. For the delay values \(\Delta t\) so that \(\delta_\omega \Delta t\ll1\), Eq. 3 reduces to Eq. 1 . In this scenario, the bunching and antibunching probability distributions present a constant quantum beat, ultimately modulated by the spectral distribution of the photons. Remarkably, this is true even for \(\Delta t/\tau\gg 1\), that is, for values of the delay \(\Delta t\) for which a non-resolved approach would see no interference: despite the finite frequency resolution, the measurement completely erases the temporal distinguishability between the photons at the detection. For higher values of \(\Delta t\), the amplitude of the beats is reduced and the effective visibility \(\mathcal{V}\mathop{\mathrm{sinc}}^2(\Delta t \delta_\omega/2)\simeq \mathcal{V}(1-\Delta t^2 \delta_\omega^2/12)\) slowly decreases with \((\Delta t \delta_\omega)^2/12\), as the finite frequency resolution still partially erases the temporal distinguishability. Only for \(\Delta t\delta_\omega \gg 1\) the beats vanish, with the probability of bunching and antibunching reducing to Eq. 4 .

We now analyze the theoretical bounds on the precision for the estimation of the delay \(\Delta t\), first recovering some known results in recent literature. The ultimate precision \(\delta\Delta t\) achievable in the estimation of \(\Delta t\) for separable photons after \(N\) repeated measurements is given by the quantum Cramér-Rao bound (QCRB) \[\delta\Delta t_{Q}=\frac{1}{\sqrt{N\cdot H}}=\sqrt{\frac{2}{N}}\tau, \label{eq:QFI}\tag{5}\] where \(H\), the quantum Fisher information [36], does not depend on the time delay to be estimated. Remarkably, it is known that frequency-resolved two-photon interference detection can saturate such an ultimate precision for perfect interference visibility (\(\mathcal{V}=1\)) and high frequency resolution (\(\delta_\omega\ll\sigma\) and \(\delta_\omega\Delta t\ll1\)) [34]. In this case, the best precision achievable, given by the Cramér-Rao bound (CRB) \(\delta\Delta t_R=\frac{1}{\sqrt{N\cdot F_\mathrm{R}}}\), where \(F_{\mathrm{R}}\) is the Fisher information associated with the resolving approach [37], reduces to Eq. 5 independently of the value of the delay \(\Delta t\).

For imperfect visibilities \(\mathcal{V}<1\), although the optimality is lost, it is still possible to estimate delay arbitrarily larger than the coherence time of the photons with the frequency-resolving scheme [34]. On the other hand, if a non-resolved measurement scheme is adopted, the best uncertainty attainable, given by the CRB \(\delta\Delta t_\text{NR}=\frac{1}{\sqrt{N\cdot F_\text{NR}}}\) where \(F_\text{NR}\) is the Fisher information associated with non-resolved two-photon interference measurements, can be calculated for Gaussian wave packets as [34] \[F_\text{NR}=\frac{H}{2}\frac{\mathcal{V}^2}{\exp{\left(\frac{\Delta t^2}{2\tau^2}\right)-\mathcal{V}^2}}\frac{\Delta t^2}{\tau^2}. \label{eq:fnr}\tag{6}\] Noticeably, even for perfect visibility \(\mathcal{V}=1\), the Fisher information \(F_\mathrm{NR}\) depends on the value of the delay to be estimated, decreasing exponentially with \(|\Delta t/\tau|\). This is a well-known drawback of sensing schemes based on standard two-photon interference, since these techniques are blind to the interference between photons whose wave packets do not overlap in time. This drawback is overcome with the use of frequency-resolving cameras, but the effect of a finite resolution remained so far an open question.

We show the effect of a finite frequency resolution \(\delta_\omega\) on the precision of the estimation in Figure 1, where we compare the QFI \(H=\frac{1}{2\tau^2}\) and the FI \(F_\mathrm{R}\) and \(F_\mathrm{NR}\), with the FI \(F_{\delta_\omega}\) associated with a setup with finite frequency resolution \(\delta_\omega\), evaluated in section 2 of Supplement 1. We can see from Figure 1 that, despite a more realistic finite frequency resolution, the estimation uncertainty is smaller than the non-resolving setup for every value of the delay. This is especially true for \(\Delta t\geqslant\tau\), where \(F_{\delta_\omega}\) decreases with \(\Delta t\) at a much slower rate than \(F_{\mathrm{NR}}\), effectively allowing for the estimation of delays that are much larger in magnitude of the coherence time of the photons. This is our main theoretical result. We refer to section 2 of Supplement 1 for the derivation of the finite frequency-resolution CRB. In the following, we present the experimental setup that we have developed and employed to test our theoretical predictions.

2.2 Experimental Setup↩︎

The presented optical setup relies on a Mach-Zehnder interferometer (MZI), which employs a pulsed laser as its light source. A delay line is inserted in one arm of the MZI to introduce a time-imbalance between the two arms, whereas in the other arm a moving mirror removes first order interference by randomizing the relative phase between photons through changes in the optical path. The photons exiting the two ports of the second BS travel along optical paths whose difference is set to half the repetition time of the source, to be later recombined and sent to a frequency-resolving single photon detector array, as shown in Figure 2.

The photon source is a pulsed laser with a central wavelength of \(531.5\,\mathrm{nm}\), a line width of \(\sim 0.5\,\mathrm{nm}\) and a repetition rate of \(40\,\mathrm{MHz}\). After an attenuator and a half-wave plate, the laser is coupled to the single-mode fiber, part of a 50:50 fiber beam splitter (BS1). An output of this BS passes through an optical delay line (ODL), in which photons undergo a variable path in free space and are later coupled back in fiber. The other output is transmitted by a free-space polarizing beam splitter (PBS) and a quarter-wave plate (QWP) to be then reflected back by a mirror mounted on a piezoelectric actuator (PA) to randomize the relative phase between the photons. The reflected photons pass back through the QWP and are reflected by the PBS to be finally coupled back into an optical fiber. Prior to the second optical fiber beam splitter (BS2), the polarization controllers (PC) erase the residual distinguishability in polarization, thus maximizing the interference visibility.

The randomization of photons relative phase is obtained by driving the piezoelectric actuator with a \(1\,\mathrm{kHz}\) triangular wave that has peak to peak voltage \(V_\text{pp}=30\,\mathrm{V}\). This voltage corresponds to a change in the optical path of approximately \(32\cdot\lambda/2\), and ensures, together with the choice of waveform, a uniform distribution of phases. Furthermore, with these settings, the measurement duration \(t_\text{meas}\approx6.1\,\mathrm{s}\) is much longer than the \(31.3\,\mathrm{\mu s}\) it takes the piezoelectric actuator to vary the optical path by one \(\lambda/2\).

The photons exiting the BS2 travel on two optical fibers of different lengths, so that they reach a beam combiner (BC) at a time interval of about \(12.5\,\mathrm{ns}\), corresponding to half the repetition time of the laser. In this way, both output of BS2 go to the same identical spectrometer. We indicate the two BS outputs with “early” and “late”.

The spectrometer consists of a beam expander (BE) with a magnification factor of \(M=5.33\), followed by a ruled reflective diffraction grating with 1200 grooves/mm and a plano-convex lens with a focal length of \(750\,\mathrm{mm}\). At the focal plane of the lens, 8 single-photon avalanche diodes (SPADs) from a 32-SPAD array, developed at the Polytechnic University of Milan  [38], are positioned for detection. The SPAD array features a center-to-center distance of \(250\,\mathrm{\mu m}\) between adjacent devices, corresponding to an expected beat frequency of \(\Delta\omega\simeq1.8\,\mathrm{THz}\). The main source of uncertainty arises from the \(50\,\mathrm{\mu m}\) diameter of the sensitive area of the SPADs. The SPADs exhibit a timing jitter of \(45\,\mathrm{ps}\) (FWHM). Their output signals are processed by a custom timing module, also developed by the same research group, operating in time-tagging mode with a resolution of \(1.5\,\mathrm{ps}\). Synchronization is provided by a \(40\,\mathrm{MHz}\) reference clock signal from the laser head.

3 Results and discussion↩︎

The goal of our experiment is to test our theoretical prediction of the error \(\delta\Delta t_{\delta_\omega}\) in the estimation of the delay \(\Delta t\) with the finite-resolution frequency-resolving approach. To do so, we have estimated the experimental uncertainty \(\delta\Delta t_\text{exp}\) multiplied by the square root of the total number of coincidence events obtained in our setup \(\sqrt{N}\), and compared it to our theoretical prediction \(\delta\Delta t_{\delta_\omega}\) multiplied by \(\sqrt{N}\) (whose corresponding FI is shown in Figure 1).

We have at first performed a control experiment with bucket detectors, that is, without resolving the frequencies of the photon, to retrieve the already known result of the sensitivity for delay estimation in a standard HOM experiment. In the non-resolved case, after removing the grating, we have acquired the HOM dip shown in the upper panel in Figure 3 by recording the average coincidence count \(C(\Delta t)\) on the output ports of BS2 as a function of the relative time delay \(\Delta t\) between the photons impinging on BS2, controlled through an optical delay line (ODL). Each value of \(C(\Delta t)\), at a given delay \(\Delta t\), is calculated by the average of \(n_r=10\) repeated measurements \(c_i(\Delta t)\), while we have estimated the experimental uncertainty \(\delta C\) associated with the count \(C(\Delta t)\) as \(\frac{\delta c}{\sqrt{n_r}}\), where \(\delta c\) is the mean squared error of \(c_i(\Delta t)\). The upper panel of Figure 3 shows the acquired average coincidence count \(C(\Delta t)\) in a \(2\,\mathrm{ns}\) time window (black points) and a superimposed fitted curve \(C_\mathrm{fit}\left(\Delta t\right) = A \left( 1-\mathcal{V} \exp\left(-\frac{\Delta t^2}{4\tau^2}\right) \right)\), justified by Eq. 4 , where the amplitude \(A\), the visibility \(\mathcal{V}\) and the coherence time \(\tau\) are the fitting parameters. Since our purpose is to test the theoretical prediction, which assumes perfect knowledge of all the parameters of the model, we consider the fitting parameters exact and without error. In a real estimation procedure, this would be achieved with a one-time precise calibration of the setup prior to the estimations.

To calculate the experimental uncertainty \(\delta \Delta t_{\mathrm{exp}}^\mathrm{NR}(\Delta t)\) over \(N=n_pn_r\) pairs of observed photons, where \(n_p=C(\Delta t)\) is the number of observed photons pairs and \(n_r\) the number of repeated measurements as a function of the delay, similarly to the procedure used in [39], we employ the fitted curve \(C_\text{fit}(\Delta t)\) as a model to estimate the delay set by the ODL by inversion. We obtain the experimental uncertainty \(\delta \Delta t_\text{exp}^{NR}\) on the measure of delay \(\Delta t_\text{exp}^\mathrm{NR}\) with our non-resolving setup by error propagation from the uncertainty in the photon coincidence \(C_\text{fit}\), namely

\[\begin{align} \delta \Delta t_\text{exp}^\mathrm{NR} = \left| \frac{\partial C_\text{fit}}{\partial \Delta t} \right|^{-1} \delta C . \label{deltatexp} \end{align}\tag{7}\] The lower panel of Figure 3 shows the experimental uncertainty \(\sqrt{N}\delta \Delta t_\text{exp}^\mathrm{NR}\) (black points), the expected uncertainty \(\sqrt{N}\delta \Delta t_{NR}\) (red curve) given by the Fisher information \(F_\mathrm{NR}\) in Eq. 6 , and the ultimate quantum CRB \(\sqrt{N}\delta\Delta t_{Q}\) (blue curve) shown in Eq.@eq:eq:QFI , per pair of observed photons. It is clear that the uncertainty associated with the adopted non-resolved scheme rapidly diverges for large delays \(\Delta t\gg\tau\), with \(\tau\simeq 0.44~\mathrm{ps}\) given by the fit of \(C_\text{fit}\).

With the introduction of the diffraction grating to disperse radiation, as previously explained, we have performed an acquisition using eight pixels of the SPAD array by time tagging the photons arrival in each pixels. The parameters used in the acquisition are reported in the experimental setup section. Using a specifically developed python software (see section 1 of Supplement 1), we calculated the matrix of coincidences \(C^{B}_{i,j}\) between the pixel \(i\) and \(j\) around the same detectors time difference representing the same output port of the BS2 HOM beam splitter (bunching), and matrix \(C^{A}_{i,j}\) of coincidence around 12.5 ns time difference representing different output ports (anibunching). The coincidences are calculated using repeated measurement to estimate the experimental uncertainty as in the non resolved case. Each pixel \(i\) of the array, as explained in the experimental setup paragraph, corresponds to a given wavelength range \(\delta\lambda\) centered at \(\lambda_i\) with its associated optical frequency \(\omega_i\). The black points in Figure 4 (5) show the counted coincidences \(C^{A}_{i,j}\)(\(C^{B}_{i,j}\)) for six(four) different values of pixel pairs \(i,j\) corresponding to optical frequency pairs (\(\omega_i,\omega_j\)). The six (four) plots demonstrate the presence of quantum beating in the coincidence counts between pixels at different separations \(k=|i-j|\), corresponding to different frequency differences \(\Delta \omega_k=|\omega_i-\omega_j|\). We have fitted the data with \[C^{A/B}_{k,\text{fit}}\left(\Delta t\right) = \mathcal{N} \left( 1\mp\mathcal{V} \mathop{\mathrm{sinc}}^2\left(\frac{\Delta t \delta}{2}\right)\cos(\Delta\omega_k\Delta t)\right), \label{eq:ResFit}\tag{8}\] justified by our model in Eq. 3 . The fits produces beating frequency \(\Delta \omega_k\), that we report in the table, proportional to \(k\) as expected. From Figure 4 and Figure 5 it is clear that the oscillations are still visible outside the coherence time \(\tau \simeq 0.44~\mathrm{ps}\), which we show enables the estimation of time delays even for photons that do not overlap in time. Similarly to the analysis for the non-resolving setup, we consider the fitted curves in Eq. 8 as exact and without error, since we aim to compare the experimental uncertainties with the associated CRB, which assumes perfect knowledge on the parameters of the model.

We determined the experimental uncertainty \(\delta \Delta t_\text{exp}\) for the frequency-resolving case through a maximum-likelihood estimation approach. Employing the fitted curves in Eq. 8 as probability mass function generating the measured data \(C^X_k\), with \(X=A,B\) and \(k=|i-j|\), we evaluate the log-likelihood function \[\mathcal{L}(\Delta t)=\sum_{\substack{X=A,B\\k}}C_k^X\ln C^X_{k,\mathrm{fit}}(\Delta t) \label{eq:LogLike}\tag{9}\] and thus the uncertainty over \(N=n_r \sum\limits_{\substack{X=A,B\\k}}C_k^X\) observed coincidences associated with the value of \(\Delta t\) that maximizes \(\mathcal{L}(\Delta t)\) is given by \[\delta \Delta t_\mathrm{exp}=\sum\limits_{\substack{X=A,B\\k}}\delta C_k^X \left|\frac{\frac{\partial^2\mathcal{L}(\Delta t)}{\partial C_k^X\partial\Delta t}}{\frac{\partial^2\mathcal{L}}{\partial\Delta t^2}}\right|, \label{eq:UncRes}\tag{10}\] caused by the fluctuations \(\delta C_k^X\) of the experimental data \(C_k^X\) (more details in section 3 of Supplement 1). Figure 6 shows \(\sqrt{N}\delta \Delta t_\text{exp}\) (black points), quantum CRB (blue line), CRB for the non-resolved approach as obtained in Eq. 6 (red line), and the CRB for finite frequency resolution that we developed (green line), per pair of observed photons. The precision obtained in our frequency-resolving setup does not diverge as rapidly as the non-resolving approach confirming a good time resolution also for non-overlapping pulses as predicted by our theoretical model.

Some values of \(\delta\Delta t_\mathrm{exp}\) \(\sqrt{N}\) lie below the theoretical bounds of the CRB and the QCRB. This behavior can be easily explained with two observations. First, retrieving the fitted parameters in Eq. 8 and performing the estimation of the uncertainties from the same data renders our model overfitted and the estimated uncertainties too optimistic. Second, the expression in Eq. 10 has a more convoluted dependence on the expected coincidence counts compared to Eq. 7 . In conjunction with a more delicate measurement procedure, this caused larger statistical fluctuations. However, these complications are expected to affect equally the estimated value of \(\delta\Delta t_\mathrm{exp}\) independently of the true value of the delay \(\Delta t\). In other words, the overall dependence of \(\delta\Delta t_\mathrm{exp}\) on increasing values \(\Delta t\), that is the focus of our study and our main result, is not affected.

4 Conclusion↩︎

We have proposed an experimental setup capable to perform frequency-resolved Hong-Ou-Mandel interferometry and proved its reliability with classical light. We have validated pioneering theoretical works [34], [35] and showcased the metrological advantage achievable with the frequency-resolved configuration in realistic experimental conditions, i.e. with finite resolutions, while performing time-delay estimations.

The results reported, far from being limited to the demonstration of intriguing fundamental properties of light, have the potential to become a benchmark for innovations in a variety of applications beyond metrology, such as boson sampling, imaging and sensing.

Funding. European Union - Next Generation EU: NRRP Initiative, Mission 4, Component 2, Investment 1.3 - D.D. MUR n. 341 del 15.03.2022 – Next Generation EU (PE0000023 - "National Quantum Science and Technology Institute). European Union - Next Generation EU, Missione 4 Componente 1, PRIN 2022, project title "QUEXO", CUP: D53D23002850006. Italian Space Agency (Subdiffraction Quantum Imaging SQI 2023-13-HH.0).
Acknowledgments. We thank Vincenzo Tamma for useful discussions. We wish to acknowledge Vincenzo Buompane, Paolo Baiocco and Graziano Spinelli for technical support.

5 Software for the calculation of the coincidence↩︎

The code to detect coincidences is written in Python. The algorithm [alg:coinc] reports the pseudocode, where \(t_a\) and \(t_b\) are the vectors containing the temporal tags of the photons arriving on the two pixels, or groups of pixels, between which the coincidences are calculated. At each step of the while loop, a coincidence is detected only if \(|\Delta t|<\text{cw}\), namely if the absolute value of the difference between the arrival time of the two photons is smaller than the coincidence window. If this happens, the code considers two photons “used” and moves on to the next ones in both lists (lines 11\(-\)​12). If the coincidence condition is not met, then the script only discards the photon that arrived first (lines 14 or 16). For example, if the photon that arrived at time \(t_a[i]\) is the first of the two and there is a coincidence, in the next step the coincidence between \(t_a[i+1]\) and \(t_b[j]\) will be verified.

6 Derivation of the frequency-resolved CRB for finite resolutions↩︎

The Fisher information \(F\) associated with the estimation of a parameter \(\phi\) from the measurement of a random variable \(X\) taking values \(x\) from a probability function \(P(x|\phi)\) is given by [37] \[F = \mathbb{E}_X\left[\left(\frac{\partial}{\partial\phi}\ln P(X|\phi)\right)^2\right],\] where \(\mathbb{E}_X\) denote the expectation value over the random variable \(X\). Applied to our exact probability mass function \(P_{A/B}^{\omega_{01},\omega_{02}}\) found in Eq. (2) of the main document, we get \[F_{\delta_\omega}=\sum_{\substack{\omega_{01},\omega_{02}\\X=A,B}}\frac{1}{P_{A/B}^{\omega_{01},\omega_{02}}}\left(\frac{\partial P_{A/B}^{\omega_{01},\omega_{02}}}{\partial\Delta t}\right)^2, \label{eq:FIdeltaApp}\tag{11}\] where the summation over \(\omega_{01},\omega_{02}\) is performed over the central frequencies of all the pixels in the camera, which are assumed to be enough in number to cover the whole spectrum \(f(\omega)\) of the photons. Although possible, the analytical evaluation of \(F_{\delta_\omega}\) from Eq. 11 is cumbersome and does not yield an informative expression. The plot of the Fisher information \(F_{\delta_\omega}\) in fig. 1 of the main document has been numerically retrieved by setting the central frequencies equally distant \(\{\omega_{01},\omega_{02}\}=\{n\delta_\omega,m\delta_\omega\}\) and summing over \(n,m\in\{-5,-4,\dots,4,5\}\), which covers the whole spectrum \(f\) for \(\delta_\omega=1\)THz and \(\tau=\frac{1}{2\sigma}=0.44~\mathrm{ps}\).

7 Maximum-likelihood estimator and sensitivity in Eq. 11↩︎

The maximum-likelihood estimator \(\Delta t_\mathcal{L}(\{C_k^X\})\) is defined as the value of \(\Delta t\) that, given the set of coincidence counts \(\{C_k^X\}\) observed, maximizes the likelihood function \[\mathcal{L}(\Delta t)=\sum_{\substack{X=A,B\\k}}C_k^X\ln C^X_{k,\mathrm{fit}}(\Delta t), \label{eq:DerLogApp}\tag{12}\] already shown in Eq. (10) of the main document, with \(k\) running over the set of pixel pairs, where we assumed that \(C^X_{k,\mathrm{fit}}\) are the exact probabilities that generated the coincidence counts \(\{C_k^X\}\). Then, by definition, \[0=\frac{\partial \mathcal{L}(\Delta t)}{\partial \Delta t}\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}}=\sum_{\substack{X=A,B\\k}}C_k^X\frac{\partial}{\partial\Delta t}\ln C^X_{k,\mathrm{fit}}(\Delta t)\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}}.\] If we differentiate this expression with respect to the outcomes \(C_k^X\), remembering that \(\Delta t_\mathcal{L}\equiv\Delta t_\mathcal{L}(\{C_k^X\})\) depends on the outcomes, we have \[0=\frac{\mathrm{d}}{\mathrm{d}C_k^X} \frac{\partial \mathcal{L}(\Delta t)}{\partial \Delta t}\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}} =\frac{\partial^2 \mathcal{L}(\Delta t)}{\partial C_k^X\partial \Delta t}\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}} +\frac{\mathrm{d}\Delta t_\mathcal{L}}{\mathrm{d}C_k^X}\frac{\partial^2\mathcal{L}(\Delta t)}{\partial \Delta t^2}\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}}.\] We can thus evaluate \[\frac{\mathrm{d}\Delta t_\mathcal{L}}{\mathrm{d}C_k^X}=-\frac{\frac{\partial^2 \mathcal{L}(\Delta t)}{\partial C_k^X\partial \Delta t}}{\frac{\partial^2\mathcal{L}(\Delta t)}{\partial \Delta t^2}}\Bigg\vert_{\Delta t=\Delta t_\mathcal{L}},\] that can be easily calculated from Eq. 12 , and thus get \[\delta \Delta t_{\delta_\omega} = \sum_{\substack{X=A,B\\k}}\left|\frac{\mathrm{d}\Delta t_\mathcal{L}}{\mathrm{d}C_k^X}\right|\delta C_k^X,\] found in Eq.(11) of the main document.

References↩︎