1. Introduction

What is the ultimate limit of optical resolution? According to the well-known Rayleigh criterion [1], a diffraction-limited imaging system can resolve two point-like emitters as long as their transverse separation is larger than or comparable to the Rayleigh length $\mathrm {x_R} = \lambda D/R$, where $\lambda$ is the wavelength, $R$ is the size of the pupil, and $D$ the distance to the emitters. The criterion follows from the fact that, due to diffraction through the pupil, the image of a point-source is a spot of size $\mathrm {x_R}$, determined by the point-spread function (PSF) characterizing the optical system [2]. Ways to circumvent this fundamental limit have been developed over time; one of the most important being super-resolved fluorescence microscopy [3], which makes use of nonlinear optics to ensure that nearby emitters do not emit at the same time. Within quantum imaging [4], methods have been proposed that exploit the fact that a quantum state of light with exactly $N$ photons has an effective wavelength that is $N$ times shorter [5], which in turn implies a PSF $N$ times sharper [6].

Note that the above strategies rely on source engineering, which is not always an option, as for example in astronomic observations. Furthermore, they are based on direct imaging (DI), i.e., pixel-by-pixel measurement of the intensity of the field focused on the image plane. However, DI ignores the information carried by the phase of the field, which instead can be uncovered by interferometric measurements. The latter may allow us to demultiplex the field in the image plane and measure it in a non-local basis of optical modes obtained as a superposition of the field in different pixels. With this approach, one may hope to achieve the ultimate quantum limit of detection [7].

In 2015, Tsang, Nair, and Lu [8] framed the problem of finding the ultimate far-field optical resolution as a problem of quantum estimation [9,10]. They considered the estimation of the midpoint and transverse separation of two point-like emitters of equal intensity. For a finite number of photons detected, the statistical error is bounded from below by the Cramér–Rao bound, which in turn is expressed in terms of the Fisher information. One can show that for DI, the Fisher information approaches zero when the separation becomes smaller than the Rayleigh length, which implies large statistical errors: a manifestation of the Rayleigh resolution criterion. In contrast, interferometric measurements show a different behavior for the Fisher information. Spatial-mode demultiplexing (SPADE) of the field in the image plane can be optimized, depending on the particular form of the PSF, to yield a constant Fisher information, independently of the value of the transverse separation. As a matter of fact, in the regime of single-photon imaging, SPADE (as well as other interferometric techniques [11]) has been shown to be the globally optimal measurement, for this and for similar estimation problems [8,12,13], among the broader set of all measurements compatible with the principles of quantum mechanics. The optimality of SPADE is not confined to estimation theory though, and applications have been explored to the related problem of hypothesis testing [14–16] and beyond the case of two point-sources [17].

The work of Tsang et al. sparkled a renewed interest in quantum imaging, which in turn led to a number of theoretical and experimental results, see Ref. [18] and the website of Tsang’s group [19] for a list of contributions. In practice, interferometric measurements can be implemented experimentally in a number of ways, e.g., with spatial light modulators [20–23], by image inversion [24], with a photonic lantern [25], and through multi-plane light conversion [26–28]. While SPADE has been proven optimal (in an information-theoretic sense) in the regime of highly attenuated signals with at most one photon per detection [8], the same approach may yield sub-Rayleigh resolution in the regime of bright sources, suggesting quantum-inspired methods for imaging and precision measurement [26,27,29,30]. The advantage of SPADE over DI may persist, with some subtle modifications and if the transverse separation is not too small, even in the presence of noise and cross talk, as a long as they are not too strong [31–34]. Recently, a modern version of the Rayleigh criterion has been formulated based on this approach [35], which takes into account the information contained in the phase [36–38].

In this paper, we report an experimental demonstration of sub-Rayleigh precision measurement, obtained through SPADE in the single-photon regime. Our experiment simulates a pair of ultra-weak incoherent point-sources that are observed through a diffraction-limited optical system characterized by its Rayleigh length. Unlike previous experimental work [20–24,28], here we consider sources of different intensities. This is motivated by the fact that pairs of natural emitters typically have unequal brightness. In particular, this may have potential application for the observation of exoplanets [15,39,40]. We implement SPADE using a telecom-wavelength demultiplexer from Cailabs [41], which sorts the optical field in the image plan on the basis of Hermite–Gaussian (HG) modes. In an independent work by Rouvière et al. [28], this demultiplexer has been used to estimate the spatial separation between two equally bright sources [42]. In a previous experiment, Ansari et al. [43] applied SPADE for the multi-parameter estimation of the source separation, midpoint, and relative intensity, in a setup where sources where separated in time instead of space. As we discuss here, the error analysis is qualitatively different when the sources have unequal intensity, see also Refs. [34,39,44], and crucially depends on how the optical system is aligned. Unlike other works, in our table-top experiment, we align the optical system toward the brighter source. Our experimental approach can be as well applied to align the optical system toward the “center of mass” (obtained using the relative intensities as weights), which is more interesting for field applications. We expect these two approaches to be quantitatively similar in the limit where one source is much brighter than the other. Alignment to the brighter source yields results that are qualitatively different when compared with alignment toward the midpoint between the two sources, see in particular Ref. [44]. Our scheme exploits only two photon detectors, which are used to count photons in the lower-order modes $\text {HG}_{01}$ and $\text {HG}_{10}$. This is sufficient to estimate the distance $d$ between the sources or, alternatively, their relative intensity $\epsilon$. The experiment is supported by theoretical modeling. We compute theoretical bounds on the statistical errors and compare with DI. In particular, we show that, in principle, noiseless SPADE can improve the signal-to-noise ratio (SNR) by a factor $\epsilon ^{-1/2}$, which may be significant when the primary source is much brighter than the secondary one, as for example for imaging of exoplanets.

2. The Model

Our experimental setup consists of two sources separated by a distance $d$, with relative intensities $\epsilon$ and $1-\epsilon$. We consider the limit of faint signals, where, most of the time, no photon is detected and the field is in the vacuum state $|0\rangle$. The probability of detecting a photon is $\eta \ll 1$ and the probability of multiple photon events is negligibly small. In this regime, the state of a single spatio-temporal mode of the field is represented by a density matrix of the following form:

In our experimental setup, we align the optical system to the brighter source. This corresponds to centering the state $\psi _0$ in the origin of the reference system in the image plane, whereas the $\psi _d$ is centered at position $(d_x,d_y)$, with $d= \sqrt {d_x^2 + d_y^2}$. For a Gaussian PSF of width $w_0$, we have

2.1 Alignment

In our experimental setup, we align the optical system to the brighter source. This is in contrast to what is done in other works, where the alignment is on the midpoint or centroid [8,11,17,20,21,28]. To align the optical system, we maximize the signal in mode $\text {HG}_{00}$. In fact, the probability of detecting a photon in the lower HG mode is

2.2 Theoretical Error Bounds

Our main goal is either the estimation of the separation $d$ between the two sources (if the relative intensity is known) or their relative intensity $\epsilon$ (if the separation is known). To achieve this goal, we measure the field in the lower HG modes, $\text {HG}_{00}$, $\text {HG}_{01}$, $\text {HG}_{10}$. In addition to $p_{0}$ in Eq. (8), we have

2.3 Cross Talk

Ideally, SPADE would exactly project the field into the HG modes. In practice, demultiplexing is affected by cross talk. As our focus is on the lower order modes $\text {HG}_{00}$, $\text {HG}_{01}$, $\text {HG}_{10}$, the cross talk is characterized by the probability $\chi$ that a photon in the $\text {HG}_{00}$ mode (which is brighter) is detected in either the $\text {HG}_{01}$ or $\text {HG}_{10}$ modes. Therefore, the probabilities in Eqs. (19), (20) need to be replaced by

3. Comparison with Direct Imaging

Formally, DI is represented as a measurement of the operators $n(x,y) = a^\dagger (x,y) a(x,y)$ on the state in Eq. (1). For simplicity, we assume to know that the secondary source lays on the $x$ axis. In the limit of infinite number of pixels (i.e., infinite resolution), the probability density of detecting a photon at position $(x,y)$ in the image plane is

 

Fig. 1. Fisher information $F_d$ for the estimation of the source separation with DI, re-scaled by the factor $w_0^2/\epsilon /\eta$, plotted versus the dimensionless separation $d_a := d/w_0$. From bottom to top, the curves are obtained for $\epsilon = 0.01,0.1,0.5$. For small separation, the Fisher information is $F_d \simeq \epsilon ^2/w_0^2$. For large separation, when the sources are well separated, we have $F_d \simeq \epsilon /w_0^2$.

Download Full Size | PDF

Similarly, we can compute the Fisher information for the estimation of $\epsilon$,

From the Fisher information, we can compute the ultimate precision in the estimation of these parameters from DI. We have

In conclusion, for small $\epsilon$, noiseless SPADE allows for a reduced error compared with DI, yielding an improvement by a factor $1/\sqrt {\epsilon }$. However, this improvement is limited by noise, in particular by cross talk. As an example, Fig. 2 shows the effect of cross talk. The plot compares noiseless DI with noiseless and noisy SPADE. While noiseless SPADE outperforms noiseless DI by a factor up to $\sqrt {\epsilon }$, this advantage is quickly washed out by cross talk. As the figure shows, a cross talk factor $\chi$ of the order of $1{\% }$ already drastically reduces the performance of SPADE. As also discussed by Linowski et al. [34], in the presence of cross talk, SPADE outperforms DI only if the separation is not too small (note that our analysis of DI is different as we align the optical system to the brighter source and not to the midpoint). However, our comparison is not entirely fair as we are considering ideal DI with an infinite number of arbitrary small pixels. In practice, DI would have a finite number of pixels and would also be affected by cross talk between neighbor pixels.

 

Fig. 2. Statistical error (re-scaled by a factor $\sqrt {\eta n}/w_0$) in the estimation of the separation $d$, plotted versus the dimensionless parameter $d_a = d/w_0$. For $\epsilon =0.01$, the solid line shows the error from noiseless DI, computed using Eq. (33). The dashed lines show the error for SPADE, in the presence of cross talk, computed from Eq. (26), for negligible alignment error ($\sigma =0$); from bottom to top, the cross talk error is $\chi =0, 0.001,0.01$.

Download Full Size | PDF

Similar conclusions can be drawn for the problem of estimating the parameter $\epsilon$; a comparison of the statistical error for noiseless DI and noiseless and noisy SPADE is shown in Fig. 3.

 

Fig. 3. Statistical error (re-scaled by a factor $\sqrt {\eta n}/w_0$) in the estimation of the relative intensity $\epsilon$, plotted versus $\epsilon$. For $d=0.5 w_0$, the solid line shows the error from noiseless DI, computed using Eq. (34). The dashed lines show the error for SPADE, in the presence of cross talk, computed from Eq. (27), for negligible alignment error ($\sigma =0$); from bottom to top, the cross talk error is $\chi =0, 0.001,0.005$.

Download Full Size | PDF

4. Comparison with the Ultimate Quantum Limit

Our analysis can be compared with that of Řehaček et al. [39], who derived the quantum Fisher information (QFI) matrix for the estimation of the parameters $s_0,d,\epsilon$, where $s_0$ is the coordinate of the midpoint between the two sources. The QFI matrix, in turn, allows one to assess the ultimate quantum limit to parameter estimation. To use their result in our setting, we need to apply a change of variables from $(s_0,d,\epsilon )$ to $(x,d,\epsilon )$, where $x = s_0-d/2$ denotes the position of the brighter source. This is due to the fact that, in our experiment, we align the optical system to the brighter source and not to the midpoint. We also note that the alignment procedure is analogous to that of parameter estimation. After the change of variables, we obtain the QFI matrix of interest for our work:

From this expression, we see that the elements of the QFI matrix show different scaling laws with respect to $\epsilon$. The element $Q_{11}$, which refers to the estimation of $x$, is of order $1$, whereas the $Q_{22}$, which refers to the estimation of $d$, is of order $\epsilon$. The correlation terms $Q_{12}$ and $Q_{21}$ are also of order $\epsilon$. This suggests that despite correlations between $x$ and $d$, an error in the estimation of $x$ will have little impact on the estimation of $d$ if $\epsilon$ is small enough. Similar consideration applies to the estimation of $x$ and $\epsilon$.

To discuss the phenomenon of separation of scale in more detail, note that, in our work, we do not consider the joint estimation of the three parameters. When we estimate $d$, we assume that $\epsilon$ is known, and vice versa. In this scenario, it is sufficient to consider the $2 \times 2$ sub-matrices of the full QFI matrix. Consider first the estimation of $x$ and $d$, assuming $\epsilon$ is known. The sub-matrix of interest is

The same argument holds if $d$ is known and one considers the estimation of $x$ and $\epsilon$. In this case, the sub-matrix of interest is

In conclusion, the statistical correlations are small in the limit of small $\epsilon$. This implies that the alignment procedure and the estimation of the parameter (either the source separation of the relative intensity) can be approximated as independent. Moreover, as anticipated above, in our experimental setup, the impact of misalignment error is negligible compared with cross talk. This observation supports the theoretical analysis of our experimental approach, where we treated the estimation error as being effectively independent of the alignment error.

5. Experimental Setup

The experimental setup is shown in Fig. 4. We combine two collimated sources (at telecom wavelength) on a non-polarizing beam splitter (NPBS) to mimic two point-like sources. The brighter source is denoted as source A, the weaker one as source B, with intensity $I_A$ and $I_B$, respectively. The relative intensity is $I_B/I_A \simeq \epsilon$.

 

Fig. 4. Experimental setup. The weaker source (B source) is a heralded photon source that crosses a two-lens system (L) to match the beam waist with demultiplexer waist ($w_0=300\,\mathrm{\mu}$m) and then a polarizing beam splitter (1). The transmitted component, used for the experiment, crosses a polarizer (P) mounted on a motorized rotation stage (2), a second polarizer (P), two mirrors (M), reflects on a beam splitter (4) and impinges on the free-space input of the demultiplexer (5). The brighter source (A source) crosses a two-lens (L) system to match the beam waist with demultiplexer waist ($300\,\mathrm{\mu}$m), a fixed polarizer (P), two mirrors (M), the beam splitter (4) where it combines with B source and impinges on the input port of the demultiplexer (5). The intensity and position of source B can be tuned using the motorized rotation stage (2) and translation stage (3). The demultiplexer (5), PROTEUS-C from Cailabs, allows to perform intensity measurements on six HG modes, but just the $\text {HG}_{01}$ and $\text {HG}_{10}$ modes are detected through two free-running InGaAs/InP single-photon avalanche diodes cooled at $-90^{\circ }$ (6), whose electrical output signal is detected and counted by a time tagger (7).

Download Full Size | PDF

The source B is a heralded photon source. It consists of a $3$-mW continuous-wave laser at $775$ nm pumping a type-II waveguide of a periodically poled lithium niobate (PPLN) crystal stabilized at $33.90^\circ$ (within $0.01^\circ$) by a proportional integrative derivative (PID) controller. The pump, by means of spontaneous parametric downconversion, generates approximately $500 \times 10^{3}$ photon pairs (orthogonally polarized) per second at $1550$ nm wavelength. At the output, the beam crosses a two-lens system to match the beam waist with demultiplexer waist ($w_0=300 \, \mathrm{\mu}$m) and a polarizer beam splitter, where only the transmitted photons are used for the experiment. The other photons of the pairs, vertically polarized, are not measured. The intensity of the vertically polarized beam can be accurately tuned by means of two free-space polarizers. Note that, as the horizontal photon is not measured, the vertical one behaves as a thermal state. In detail, the beam from the source B crosses a film polarizer mounted on a motorized rotation stage and then another film polarizer (at fixed angle). In this way, by changing the angle of the rotation stage, it is possible to tune the B beam intensity while maintaining the polarization fixed (as the second polarizer is fixed). We align the angle of the fixed polarizer with direction of maximum detector responsivity to the electric field: in this way, we are confident that the polarization of the B beam is fixed for any angle of the first polarizer bypassing the issue of polarization-dependent detector efficiency. Finally, by means of two steering mirrors (the second being mounted on a micrometric translation stage to shift the beam position), the beam is coupled with the free-space input port of the demultiplexer after reflection on a beam splitter used to overlap B and A beams.

The source A is a few-mW light-emitting diode (LED) at $1550$ nm, attenuated by $90$ dB using fiber attenuators, and collimated using a two-lens system to match the beam waist with the demultiplexer waist. As for the source B, a polarizer is used to align with the direction of maximum detector efficiency and a pair of steering mirrors are used to couple the beam with demultiplexer input after crossing a beam splitter where A and B beams are overlapped. The position of the beam A remains fixed during the experiment, whereas the source B can be attenuated by changing the angle of rotation stage and shifted by moving the translator. With this scheme, it is possible to tune the separation $d$ and the intensity ratio $\epsilon$ independently.

The two beams are finally fed into a demultiplexer, PROTEUS-C model from Cailabs, which allows for intensity measurements on up to six HG modes. It accepts radiation from the free-space input port, and decomposes it in the lowest-order modes ($\text {HG}_{00}$, $\text {HG}_{01}$, $\text {HG}_{10}$, $\text {HG}_{11}$, $\text {HG}_{20}$, $\text {HG}_{02}$). The modes are coupled with six single-mode fibers following conversion into the $\text {HG}_{00}$ mode. Finally, the modes $\text {HG}_{01}$ and $\text {HG}_{10}$ are coupled, through a single-mode fiber, with free-running InGaAs/InP single-photon avalanche diodes cooled at $-90^{\circ } C$ and operating with dead time of $20 \, \mathrm{\mu}$s allowing a negligible dark count rate. The detectors generate an electric pulse that is recorded by a time-to-digital converter. When the beams separation is zero (i.e., the beams overlap completely) and the overall source (A+B beam) has a circular symmetry, the power leaked into the $\text {HG}_{01}$ and $\text {HG}_{10}$ modes reaches its minimum value.

The main experimental error in the setup is cross talk. The cross talk between $\text {HG}_{00}$ and $\text {HG}_{nm}$ may be quantified by the ratio $P_{nm} / P_{00}$, where $P_{nm}$ is the power on the $\text {HG}_{nm}$ output channel if only $\text {HG}_{00}$ is injected with a power $P_{00}$ from the input. The cross talk is a limitation that is due to the presence of light in high-order modes even if the incoming radiation is fully matched with the demultiplexer (in the ideal case of negligible cross talk, only $\text {HG}_{00}$ should be excited). The cross talk factor $\chi$ introduced in Section 2.3 is $\chi = P_{01} / P_{00} + P_{10} / P_{00}$.

In the experiment, we operate the demultiplexer in the single-photon regime for the estimation of either the beam intensity ratio $\epsilon = I_B/I_A$ or the beam separation $d$. The beams have to simulate two point-like sources with Gaussian PSF characterized by the demultiplexer waist $w_0$. To simulate a situation where we do not know whether there is a single source or there are two sources, we align the system by maximizing the $\text {HG}_{00}$ output when the two sources are completely superimposed. Data are collected by translating the B beam and keeping the A beam centered. In a practical scenario, one can only align the demultiplexer with the “center of mass” of light intensity, since the positions of the sources are unknown. Consequently, the proper procedure for device calibration would be to maximize the power in $\text {HG}_{00}$ mode as we did, but then, to move both the beams in opposite directions. However, as discussed in Section 2.1, the two procedures are equivalent up to an error of order $\epsilon d$.

5.1 Calibration and Measurement

We acquired the photon counts $C_{n}^{H_{01}}$ ($C_{n}^{H_{10}}$) in the $\text {HG}_{01}$ ($\text {HG}_{10}$) mode. The total photon counts impinging on detectors coupled to $\text {HG}_{01}$ and $\text {HG}_{10}$ modes is then:

For each value of $\epsilon$ and $d$, we detected photons in $\text {HG}_{10} + \text {HG}_{01}$ modes during $1$ s of acquisition time. This is obtained by combining $N_m=100$ independent measures $n_1^i$, for $i=1,\dots, N_m$, lasting $10$ ms each. Figure 5 shows the total counts $n_1(d_a)=\sum _{i=1}^{N_m} n_1^i(d_a)$ as a function of the dimensionless separation $d_a := d/w_0$, for four different values of $\epsilon$ (blue points), and the corresponding quadratic regressions $n_1^F(d_a)=a_1+b_1 d_a^2$ (red line), where $a_1$ and $b_1$ (related to cross talk and dark count) are fitting parameters. As discussed in Section 2.2, a quadratic law is expected for small separations, $d_a \ll 1$,

 

Fig. 5. Top of panels: counts $n_1$ as a function of $d_a$, for four values of the intensity ratio $\epsilon$; the blue points are from experimental data, with error bars estimated upon repeated measurements; the red line shows the result of a quadratic regression on the data points. Bottom of panels: residuals between experimental data and fitted curve. They appear randomly distributed indicating the quadratic function is a good approximation of HG in the range of measurements.

Download Full Size | PDF

We use the fitted curve $n_1^F(d_a)$ as a calibration curve to measure the distance set on the translator stage. We obtain the associated uncertainty $\delta d_a$ by error propagation from the uncertainty in the total photon count $n_1$,

Figure 6 shows the experimental uncertainty $\delta d_a$ (orange points) compared with the theoretical bounds for ideal DI computed in Section 3. The latter are computed using the estimated photon number $\eta n$ during $1$ s of acquisition time, which is in the range between $40 \times 10^{3}$ and $80 \times 10^{3}$. This shows that our setup, though affected by cross talk and dark counts (errors due to residual misalignment are negligible), allows us to beat ideal and noiseless DI if $d_a$ is not too small.

 

Fig. 6. Experimental uncertainty $\delta d_{a}$ (orange points) as function of $d_{a}$, for four different values of the intensity ratio $\epsilon$, and the theoretical error estimates (blue lines) for DI imaging calculated for $\eta n = 40 \times 10^{3}$ and $\eta n = 80 \times 10^{3}$ photons.

Download Full Size | PDF

We also note that our experimental errors are approximately 2–4 times larger than the SPADE theoretical bounds of Section 2.3 (computed using the experimental value $\chi = 0.0035$ [27]), and the gap increases at small separation $d_a \ll 1$. The discrepancy is likely due to finite-sampling and noise in the mechanical mounts and imperfections of the translation stage. In fact, as the measurement is performed by changing the position of translation stage, the imperfections of the translation generate increased uncertainty $\delta d_a$ that does not depend on the imaging system but just on the sources. In other words, in a real acquisition, these imperfections should not be present.

We carry out analogous calibration and measurement for the estimation of the parameter $\epsilon$. In Fig. 7, we represent the counts $n_1$ (blue points) as a function of $\epsilon$ ratio at fixed distance, and the result of a linear regression (red line) that we use as a calibration curve. The experimental uncertainty $\delta \epsilon$ is obtained as

 

Fig. 7. Top of panels: counts $n_1$ as a function of $\epsilon$ at four values of separation $d_a$ (blue points) with error bars obtained from $N_m$ repeated measurements; the red line show the result of a linear regression on the data points. Bottom of panels: residuals between experimental data and fitted curve. They appear randomly distributed confirming the proportionality between $\epsilon$ and $d_a$.

Download Full Size | PDF

In Fig. 8, the experimental uncertainty $\delta \epsilon$ is compared with theoretical bounds for ideal DI. As one would expect, SPADE performs better than DI for larger values of $d_a$ and small values of $\epsilon$. Also, the experimental errors are larger than the SPADE theoretical error bounds (less than a factor two). In this case, the discrepancy is smaller compared to the discrepancy on the measure of $d_a$, most likely because the measurements are performed by changing $\epsilon$ without touching the translation stage.

 

Fig. 8. Experimental uncertainty $\delta _\epsilon$ (orange points) as a function of $\epsilon$ at four values of separation $d_a$ and corresponding theoretical error estimates (blue lines) for DI imaging calculated for $\eta n = 40 \times 10^{3}$ and $\eta n = 80 \times 10^{3}$ photons.

Download Full Size | PDF

6. Conclusions

In a diffraction-limited optical system, the resolution of direct imaging (DI) is limited by the width of the point-spread function, according to the Rayleigh criterion. However, when interferometric measurements are employed, the information carried by the phase of the field may allow us to beat the Rayleigh resolution limit. This is particularly important in the regime of single-photon imaging, where interferometric measurements achieve the ultimate precision limit as set by the quantum Cramér–Rao bound.

In this work, we have demonstrated single-photon spatial-mode demultiplexing (SPADE) using multi-plane light conversion technology to sort the transverse field into its Hermite–Gaussian (HG) components. Working with two single-photon avalanche diodes, we have detected the lowest-order modes $\text {HG}_{01}$ and $\text {HG}_{10}$, which are sufficient to estimate either the separation between two point-like sources (if the relative intensity is known), or their relative intensity $\epsilon$ (if the separation is known). Unlike previous experimental works, here we have focused on sources of unequal intensity. This is motivated by the fact that pairs of natural sources typically have unequal intensities. Of particular interest is the situation when the primary source is much brighter than the secondary one, $\epsilon \ll 1$, which suggests application to exoplanet imaging and detection. In fact, it is in this regime that SPADE may outperform DI by a factor $\epsilon ^{-1/2}$.

Our experimental errors are dominated by cross talk and finite-sampling fluctuations are approximately three times larger than the asymptotic theoretical estimates. Nevertheless, we have been able to demonstrate an error in parameter estimation below what could be achieved by ideal, noiseless DI. By using larger samples and suitable opto-mechanical mounts, we may substantially reduce the discrepancy with the theoretical bounds, potentially allowing SPADE to beat noiseless DI by a greater extent and for smaller values of the source separation. In particular, the advantage of SPADE of DI may be more easily observed for smaller value of $\epsilon$; whereas in our experiment, the factor $\epsilon ^{-1/2}$ was always smaller than $10$. Moreover, there is room for experimental improvement by using faster and more efficient detectors (e.g., nanowires). Finally, the most crucial parameter is the cross talk value; this is not a fundamental limitation and may be improved by future advances in manufacturing demultiplexer devices.