Optical quantum super-resolution imaging and hypothesis testing

Estimating the angular separation between two incoherent thermal sources is a challenging task for direct imaging, especially at lengths within the diffraction limit. Moreover, detecting the presence of multiple sources of different brightness is an even more severe challenge. We experimentally demonstrate two tasks for super-resolution imaging based on hypothesis testing and quantum metrology techniques. We can significantly reduce the error probability for detecting a weak secondary source, even for small separations. We reduce the experimental complexity to a simple interferometer: we show (1) our set-up is optimal for the state discrimination task, and (2) if the two sources are equally bright, then this measurement can super-resolve their angular separation. Using a collection baseline of 5.3 mm, we resolve the angular separation of two sources placed 15 μm apart at a distance of 1.0 m with a 1.7% accuracy - an almost 3-orders-of-magnitude improvement over shot-noise limited direct imaging.

SUPPLEMENTARY FIG. 2. Experimental g (2) (τ ) autocorrelation function of the pseudo thermal source. The central peak at zero delay time shows a g (2) (0) = 1.977 ± 0.003 which is in close agreement with a theoretical value of a true thermal state. The coincidence peaks are separated by 1 µs which is consistent with the clock repetition rate of 1 MHz. All results have been collected in real-time via the QuCoa (PicoQuant) software with a total integration time of 60 s in order to limit evaluation errors.

II. INTERFEROMETRIC CALIBRATION
This section shows the diagram of a one-shot phase modulation electrical signal used to generate a thermal state via an electro-optic modulator (EOM). See subsection B of the Method section in the main manuscript for a more detailed description.
... th SUPPLEMENTARY FIG. 3. Diagram of a one-shot phase modulation signal to the EOM. The multiplexed reference signal is repeated every two electrical pulses (orange rectangles) effectively halving the clock repetition rate to 500 KHz. The thermal signal instead (blue rectangles) is selected uniformly at random within the range [0, 2π) and similarly repeated every other electrical signal. The duration of all pulses is matched to the clock repetition rate of 1 MHz, i.e. 1 µs and the final pattern is then repeated every 200 ms.

III. OPTICAL MASK FABRICATION AND CHARACTERISATION
The optical masks were formed using patterned etching of thin chromium layers on a fused silica substrate. The 1.5 mm thick fused silica substrates were coated with 90 nm thickness of chromium using an electron-beam vacuum evaporation process forming a layer sufficiently thick to be fully opaque to the near-infrared radiation used in this experiment. The substrate was then coated in 5 nm positive photoresist (AZ 1505), and the pattern (pinholes, reference and alignment markers) inscribed using a Heidelberg DW66+ laser-writer. Once developed, the chromium was removed using a chemical etchant (TechniEtch Cr01). Fig. 4 displays pictures of the masks taken with a Leica DMRM microscope (15x magnification).  (15x magnification). Image (a) shows a mask with 30 µm wide pinholes separated by 90 µm while image (b) shows a mask with pinholes of the same size but a 35 µm separation. Pictures' contrast has been saturated to better resolve the two pinholes against the background.
A single-photon sensitive CCD camera (Rolera EM-C 2 Bio-Imaging Microscopy Camera) was used to extract profile intensity images of the transmitted light by the two circular pinholes at different distances. Fig. 5 shows the normalised intensities at imaging distances z = 2, 3 and 20 cm for a mask with 30 µm wide pinholes separated by 1 mm using thermal radiation. At short distances, the intensity profiles depict the classical Airy diffraction pattern expected from circular apertures where faint secondary rings are visible. In these configurations, the two sources can still be resolved, however, as the distance between the masks and the camera increases, diffraction prevails and the profiles merge together removing any knowledge of the initial sources. The masks were also tested using coherent radiation to ensure that the thermal generation process successfully removed any spatial correlation that could potentially disrupt the interferometer's mode sorting mechanism. Fig. 6 shows the resulting image for a mask with 30 µm wide pinholes separated by 150 µm placed at a distance z = 15 cm from the camera. Interferometric fringes are clearly recognisable showing a fringe separation of ≈ 0.821 mm which is in good agreement with the theoretical value expected for this imaging system. y (mm) x (mm) Without any phase and intensity modulation provided by the EOMs, the two point-like sources undergo "classical" interference as demonstrated by the visible interferometric fringes. The mask used had 30 µm wide pinholes separated by 150 µm placed at a distance of 15 cm from the camera. The separation of the fringes is ≈ 0.821 mm which is consistent with the expected theoretical value for the imaging system used.

IV. QUANTUM STATE DISCRIMINATION
In a previous work [1], we have shown that a two-mode interferometer has the same sensitivity in estimating the separation between two sources as SPADE, given comparable numerical apertures. Here we show that the same two-mode interferometer is also optimal in our discrimination problem.
x 0 x 0 +s photon counter 50:50 BS collector SUPPLEMENTARY FIG. 7. Schematic of two sources with a separation of s in the object plane, at a distance z0 from the collectors. Two collectors at d1 and d2 direct light into a two-mode interferometer consisting of a phase shift of α and a 50:50 beam splitter, followed by photon counters.
Consider the set-up in Fig. 7 where two collectors are placed at positions d 1 and d 2 orthogonal to the optic axis; the collectors are at a distance z 0 from the sources. One source (the star) is positioned at x 0 , and the planet, if present, is positioned at x 0 + s.
Assuming we are in the paraxial regime, the optical path difference of the planet between the two collectors is The optical path difference of the star to the two collectors is where k is the wavenumber of the radiation. The states to discriminate between are: We define the angular separation θ = s/z 0 . In the basis of ρ 0 , the two density matrices are The QRE between ρ 0 and ρ 1 is approximately Now, we apply the measurement in Fig. 7. We put the collected light at d 1 , d 2 through a phase shift e iα , followed by a 50:50 BS. We assume the operators transform as For H 0 , the measurement outcomes are For H 1 , they are Define we set α to This is analogous to the SPADE method, where we align the apparatus to the weighted center. That is, the optical system will be aligned towards the weighted center (centroid). If only the star is present, then the centroid corresponds to the position of the star. Otherwise, it is somewhere between the star and the planet, namely, x 0 (1 − ϵ) + (x 0 + s)ϵ.
As an example, we could observe a blurred image without being able to discriminate the two point-sources, yet we can still identify the point of maximum apparent intensity and point the telescope in its direction. In other words, we are assuming a preliminary phase in which the position of the centroid has been estimated. The classical relative entropy of this measurement is which is optimal in the limit that θ, ϵ ≪ 1.
As an example, we now show how the relative entropies are computed experimentally. For a particular value of ϵ, we measured the detection statistics for H 0 and H 1 at the value of α that maximises the output probability at one of the detectors. This value of α corresponds to the phase shift that maximises the relative entropy. For example for H 0 , α = −ϕ maximises P H0 (a), and the detection probabilities are P H0 (a) = 1/2 (1 + ν) , where ν is the interferometer's visibility. For the statistics for H 0 , we set ϵ = 0. During the experiment, detector a registered 59461 photons, and detector b registered 176 photons.
Therefore we take The goal of this section is to calculate the statistics of the coherent state at the output of detectors a and b, since this is used for calibrating the adjustable phase shift α (Fig. 8). In the experiment, the laser is input into a multimode fibre. On exiting the multimode fibre, the light is highly divergent, behaving as a point source as it propagates towards the two slits. Since there is no randomness introduced in the process, we model the output at the two slits as a plane wave of coherent states. The reference laser provides a calibration for the applied phase α. The electric fields at d 1 and d 2 are given by the combination of the two fields of the two sources, whose modes we label as b 1 and b 2 : where η is the transmission parameter.
Since the two slits are of the same size (equal intensity), we model the coherent states at b 1 and b 2 as where β is the amplitude of the coherent state, and κ describes a potential phase difference between the output at b 1 and b 2 . The optical path differences are Using the relationships in Eqs. (25)-(28) and after some simplifications, the photon numbers at d 1 and d 2 are: We need to calculate the expectation value of this observable, O = Re E † d1 E d2 e iα to obtain the correlations: The expectation value of the operator in Eq. (31) is After the collectors at d 1 and d 2 , we have a phase-shifter and 50:50 beam splitter. This transformation gives We can now calculate the statistics of the coherent state at the two detectors, A and B. We use capitalised letters to distinguish the operators from the detectors' labels a and b. Inverting the above gives If the interferometer is imperfect, where some of the signal is replaced by noise, we model this as By combining the results in Eqs. (29)-(41), we can calculate the normalised contrast between the two detectors R. This is defined as the difference in photon counts, divided by the total: ν cos(α) cos(κ) + cos(ϕ) cos(κ) cos(ϕ) + 1 = ν cos(α) when κ = 0 . (42) For κ = κ 2 − κ 1 = 0, which is our case here, since the two slits are equidistant from the output of the multimode fibre, the above expression reduces to cos α. That is, the reference laser behaves almost like a single point source.
The parameter κ being non-zero will only reduce the visibility of the measurement, and the effect is almost negligible. We measure the parameter ν cos(α) directly from experimental data, which is then used to update the probability distribution in the maximum likelihood method.
Here there is a normalisation factor that easily dealt with. We know the update distributions P(µ | ϕ), these are P(a | ϕ, α, ν) = 1 2 1 + ν cos(α) cos (ϕ) , Given a detection event µ = a, b where the adjusted phase was α, the probability for ϕ can be updated via Therefore, after the detection event, Eq. (48) is updated using Eqs. (46)-(47), depending on whether a or b occurred. After m detections, we have the vector of detection events, ⃗ µ m = (a, b, b, a, ...), for a given vector of adjustable phases ⃗ α = (α 1 , α 2 , ..., α m ). In the experiment, α is constant for each data point. We have Now, since all the functions we deal with here are sinusoids, we can conveniently express them as a Fourier series. After m clicks, the probabilities can be expressed as a Fourier series where m here corresponds to the higher order of the Fourier coefficient. The coefficient of the term e ikϕ is denoted a k , which depends on ⃗ µ m and α m . Normalising α 0 to 1/2π will keep the entire distribution normalised. We equivalently write the update events in this Fourier form. For example, if detector b fires, then We write Eq. (51) as therefore, the update coefficients are a 0 = π, a 1 = a −1 = π 2 ν cos(α). This example is particularly relevant, because the factor ν cos(α) is equal to R in Eq. (42), and is directly measured in the experiment using the calibration laser. Before the first detection, Eq. (50) only contains one term, a 0 = 1. After each detection event given by the probabilities in Eqs. (46)-(47), the number of Fourier coefficients grows by 2; a k are updated using Eq. (49), which once again uses Eqs. (46)-(47). After the coefficients are obtained, P (ϕ) can be evaluated for each value of ϕ, and we extract the maximum. There are two peaks for ϕ of equal intensity, since cos(ϕ) = cos(−ϕ). We show an example in Fig. 9 with ν cos α = 0.981 after 12740 detection events where 1478 were output at detector b.

VI. PRACTICAL EXAMPLE WITH AN IMAGING INSTRUMENT
As an example, the exoplanet LkCa 15 c [2] was observed with a Large Binocular Telescope (LBT) with variable baselines ranging from 1.4 to 7.0 m, at wavelengths 2.18 µm and 3.8 µm.
In imaging mode, the PSF of a circular lens is the Airy function, which we approximate with a Gaussian: for which the standard deviation is well-approximated by where k is the wavenumber and r is the diameter of the combined mirrors.
If the instrument is diffraction-limited, then the smallest PSF obtainable has σ ≈ 29 milliarcseconds (with a wavelength of 2.18 µm and a 7 m baseline). In Fig. 10 we show the precision limits of such an instrument operating in direct imaging (DI) mode, as well as using our method. The orange shaded region denotes the achievable precision if the visibility is between 96% and 98%. The red points on the plot correspond to the selected values of angular separation θ = 0.1σ, 0.5σ, 1σ and 1.5σ respectively, with a visibility parameter ν = 0.98 (for comparison, the largest optical interferometer, CHARA, has a visibility value ν > 99%). As we see, using our method, an instrument with a PSF of ≈ 30 milliarcseconds can resolve two binary stars with separations θ < 10 milliarcseconds with reasonable accuracy.