Super-resolution imaging using the spatial-frequency filtered intensity fluctuation correlation

We report an experimental demonstration of a nonclassical imaging mechanism with super-resolving power beyond the Rayleigh limit. When the classical image is completely blurred out due to the use of a small imaging lens, by taking advantage of the intensity fluctuation correlation of thermal light, the demonstrated camera recovered the image of the resolution testing gauge. This method could be adapted to long distance imaging, such as satellite imaging, which requires large diameter camera lenses to achieve high image resolution.


Experimental setup
illustrates the laboratory-demonstrated camera setup. We used a standard narrow spectral bandwidth pseudo-thermal light source consisting of a 10 mm diameter 532 nm wavelength laser beam scattered by millions of tiny diffusers on the surface of a rotating ground glass. The object imaged was the 5-3 element of a 1951 USAF Resolution Testing Gauge. Like a traditional camera, the imaging lens, L I , had an aperture limited by an adjustable pinhole to approximately 1.36 mm diameter. However, this imaging device has two optical arms behind its imaging lens L I . The light transmitted by the object falls on a single-mode 0.005 mm diameter fiber tip that scans the image plane of arm one and is then interfaced with a photon counting detector D 1 . The combination of the scanning fiber tip, D 1 , and photon counting detector acts as a CCD array. In arm two, the second lens L B is placed behind the image plane and performs a Fourier transform of the field distribution of the image plane of arm two. D 2 , a scannable multimode 0.105 mm diameter fiber interfaced with a photon counting detector, is placed in the Fourier transform plane of L B and integrates only the higher spatial frequencies while filtering out the lower spatial frequency modes; we emphasize that the placement of this "spatial filter" does not depend on knowledge of the Fourier transform of the image.
The intensity fluctuations in this experiment were recorded by a Photon Number Fluctuation Correlation (PNFC) circuit 9,10 , which independently records the arrival time of each photo-detection event at D 1 or D 2 . The intensity, measured by the number of photons detected per second, is divided into a sequence of short time windows, each of which needs to be less than the second-order coherence time of the light; it is important that the width of the time window not be too long. The software first calculates the average intensity per short time window, n s , where s = 1, 2 indicates the detector, and then the difference or fluctuation term for each time window: . The corresponding statistical average of 〈 Δ n 1 Δ n 2 〉 is thus It should be emphasized that when the fluctuation correlation is calculated between the two detectors, Δ n 1 and Δ n 2 have not yet been time-averaged. The time averaging is performed after the correlation, as indicated by the multiplication appearing inside the sum over j.

Experimental results
Typical experimental results are presented in Fig. 2. In this measurement, the 5-3 element of a 1951 USAF Resolution Test Target was imaged in one dimension by scanning D 1 in the x-direction along the slits. Figure 2(a) shows a completely unresolved classical image of the three slits, I 1 (x 1 ), that was directly measured by the scanning detector D 1 . For reference the gray shading indicates the location of the slits. Figure 2(b) shows two results: the black dots plot the autocorrelation [I 1 (x 1 )] 2 , while the blue triangles show the fluctuation autocorrelation 〈 [Δ I 1 (x 1 )] 2 〉 at each point. 〈 [Δ I 1 (x 1 )] 2 〉 was calculated from the intensity fluctuation autocorrelation of D 1 . The measurements in Fig. 2(b) have a 2 resolution gain, similar to that of Oh et al. 6 . Using the Rayleigh limit 11,12 δx = 1.22λs I /D 11,12 , where s I is the distance from lens L I to the image plane and λ is the wavelength of illumination, the expected resolution of the autocorrelation in the image plane is approximately δ = .
x/ 2 0 13 mm which, as seen in 2(b), is not enough to resolve the three slits which have a slit-to-slit separation of about 0.13 mm. However, by spatially filtering arm two, the three slits of the 5-3 element of the gauge were clearly separated when correlated with arm 1, as seen in Fig. 2(c). The error bars in Fig. 2(a) and (b) are quite small, especially when compared to those in Fig. 2(c). This is a typical negative feature of second-order measurements; compared to first-order classical imaging, in order to achieve the same level of statistics the reported imaging mechanism needs a longer exposure time. How much longer the measurement takes depends on the power of the light source and other experimental parameters. Figure 2(c) is the sum of two measurements obtained by placing the bucket fiber tip at two points in the Fourier transform plane: x 2+ = 0.22 mm and x 2− = − 0.24 mm; this selects the higher spatial frequency modes which fall onto the two fiber tips and "blocks" all other spatial frequency modes. We represent this mathematically with the filter function F(x 2 ) = Π (x 2 − x 2+ , D F ) + Π (x 2 − x 2− , D F ) to simulate the physical "spatial filtering", where Π (x, w) is a rectangle function of width w, D F = 0.105 mm is the fiber diameter, and x 2 is measured from the central maximum. Then in one dimension for κ ∝ kx 2 /f B , where f B is the focal length of the bucket lens, Δ R c (x 1 ) = 〈 Δ I 1 (x 1 )∫ dx 2 F(x 2 )Δ I 2 (x 2 )〉 , and only integrates the higher spatial frequencies collected by the bucket Figure 1. Experimental setup: a 10 mm diameter 532 nm wavelength laser beam scatters from a rotating ground glass (i) and strikes a 1951 USAF Resolution Testing Gauge (ii), then the imaging lens L I (iii) and a pinhole of ≈1.36 mm diameter (iv). The light splits to arm one, where it is collected by a scanning point detector D 1 , and arm two, where it is collected by a spatially filtered bucket detector, consisting of a lens L B located on the image plane and a multimode fiber tip in the Fourier transform plane. The image is then calculated using the Photon Number Fluctuation Correlation (PNFC) circuit.
detector in the neighborhood of x 2 = x 2+ , x 2− . Again, this "spatial filtering" does not require any knowledge of the object or its Fourier transform function.
One way to improve these results is to replace D 2 with a CCD array; the CCD would still be in the Fourier transform plane, but with the central pixels blocked. This would allow D 2 to collect more light of higher spatial frequencies. Although the limits of our equipment, software data storage, and time constraints prevented the authors from making such improvements, in this reported measurement all three slits of the resolution gauge are certainly well-resolved, while both the classical imaging and the autocorrelation mechanisms could not resolve it.

Discussion and theory
In the experiment, we use the spatial correlation of the noise, 〈 Δ I 1 (r 1 , t 1 )Δ I 2 (r 2 , t 2 )〉 , to produce an image from the joint photo-detection of two independent and spatially separated photodetectors, D 1 and D 2 . In the following, we outline the theory behind our experiment. First we briefly consider how a first-order or classical camera produces an image in its image plane, 〈 I 1 (ρ 1 )〉 .
The experiment was performed using a pseudothermal light source, created by placing a rotating ground glass in the path of a laser beam. The ground glass contains a large number of tiny scattering diffusers which act as sub-sources. The wavepackets of scattered light play the role of subfields; each diffuser scatters a subfield to all possible directions, resulting in the subfields acquiring random phases. Each sub-field propagates from the source plane to the image plane by means of a propagator or Green's function, E m (ρ 1 ) = E m g m (ρ 1 ), where E m is the initial phase and amplitude of the field emitted by sub-source m and g m (ρ 1 ) is the Green's function which propagates the light from the m th sub-source located at ρ m to the point ρ 1 at some distance z from the source plane. To simplify the problem, we assume the fields are monochromatic and ignore the temporal part of the propagator. Then the light measured at coordinate (r, t) is the result of the superposition of a large number of subfields, where 〈 I(r, t)〉 , the mean intensity, is the result of the mth subfield interfering with itself; Δ I(r, t), the intensity fluctuation, is the result of the mth subfield interfering with the nth subfield, m ≠ n, and is usually considered noise because 〈 Δ I(r, t)〉 = 0 when taking into account all possible random phases of the subfields. A classical imaging system measures the mean intensity distribution on the image plane, 〈 I 1 (ρ 1 )〉 , where we have assumed a point detector D 1 is placed at coordinate ρ 1 , the transverse coordinate of the image plane. In an ideal imaging system, the self-interference of subfields produces a perfect point-to-point image-forming function. The ideal classical image assuming an infinite lens is the convolution between the aperture function of the object |A(ρ O )| 2 and the image-forming δ-function which characterizes the point-to-point relationship between the object plane and the image plane 11-13 .
is a Green's function propagating the m th subfield from the source plane to the object plane over a distance z. and g O (ρ 1 ) is a function propagating the subfield from the object plane to the detection plane over a distance s O + s I , and including the imaging lens. A(ρ O ) is an arbitrary function describing the object aperture. In reality, due to the finite size of the imaging system, we rarely have a perfect point-to-point correspondence. Incomplete constructive-destructive interference blurs the point-to-point correspondence to point-to-spot correspondence. The δ-function in the convolution of Eq. 3 is then replaced by a point-to-spot image-forming function, or a point-spread function which is determined by the shape and size of the lens. For a lens with a finite diameter, one common model describes the shape or pupil of the lens as a disk of diameter D: where the sombrero-like point-spread function is defined as somb(x) ≡ 2J 1 (x)/x; J 1 (x) is the first-order Bessel function. The image resolution is determined by the width of the somb-function: the narrower the higher. A larger diameter lens results in a narrower somb-function and thus produces images with higher spatial resolution.
To simplify the mathematics, it is common to approximate a finite lens as a Gaussian ρ − e D ( /( /2) ) L 2 with diameter D, but a smoother falloff than the disk approximation. This leads to a Gaussian imaging-forming function: This Gaussian version of the imaging equation will be used later in numerical calculations to simplify the mathematical evaluation.
It is clear from Eqs 4 and 5 that for a chosen value of distance s O , a larger imaging lens and shorter wavelength will result in a narrower point-spread function, and thus a higher spatial resolution of the image. Now we consider the noise produced image that is observed from the measurement of Fig. 1 by means of 〈 Δ I 1 (ρ 1 )Δ I 2 (ρ 2 )〉 . To make the explanation of the experimental results easier to follow, first we examine the case where two point scanning detectors D 1 and D 2 are placed in the image planes of arm one and arm two: The calculation of Next, we complete the summation over m in terms of the subfields, or the sub-sources, by means of an integral over the entire source plane. This integral results in the well-known Hanbury-Brown Twiss (HBT) correlation: where Δ θ is the angular diameter of the light source relative to the object plane. To simplify further calculations, we assume a large value of Δ θ and approximate the somb-function to a δ-function evaluated at ρ O = ρ O′ , κ = κ′ . 〈 Δ I 1 (ρ 1 )Δ I 2 (ρ 2 )〉 is therefore approximately equal to: It is clear that when ρ 1 = ρ 2 in Eq. 8, the measurement of 〈 Δ I 1 (ρ 1 )Δ I 2 (ρ 1 )〉 produces an image with a 2 resolution gain, with an imaging resolution due to the image-forming somb-functions, i.e., ρ ρ ρ . When the lens is large enough to resolve the object, the result is a point-to-point reproduction of the image only when ρ 1 = ρ 2 ; otherwise for small lens apertures Eq. 8 forms a point-to-spot image when |ρ O + ρ 1 /μ| < λs O /D and |ρ O + ρ 2 /μ| < λs O /D. Now we move D 2 to the Fourier transform plane of L B of arm two, i.e., to its focal plane, effectively performing a Fourier transform of the field distribution of the image plane. In addition, D 2 is placed off-center relative to the Scientific RepoRts | 6:38077 | DOI: 10.1038/srep38077 optic axis of the lens to select part of the spatial frequencies on the Fourier transform plane, acting as a spatial frequency filter. Mathematically, As a result of the spatial filter function F(κ 2 ), the imaging resolution of Eq. 9 is much narrower than that of a first-order image; however, it is difficult to simplify this equation further in this form.
To get a better understanding of the physics behind Eq. 9, instead of modeling the finite radius of the lens as a disk, which results in the somb-function, we approximate the finite radius of the lens as the Gaussian function with a half-width D/2, and evaluate in one dimension. This leads to a Gaussian imaging-forming function instead of the somb-function. Working in one dimension, we change ρ O to x O ; ρ 1 to x 1 ; ρ 2 to x 2 , etc. Then Δ R c (x 1 ) simplifies to: Corresponding to the experimental measurement, where D 2 was placed at two off-center points in the Fourier transform plane, we model the filter function in one dimension by two rectangle functions: and κ 2 (x F ) are integrated from ±∞ , the resulting equation is an analytic expression. Define Then Eq. 10 is, after evaluating the κ 2 and x 2 integrals, Then it is easy to see that restricting the allowed spatial frequencies of the Erfi functions constrains the values x O is allowed to take, which, together with f(x 1 , x O ), improves the ability to resolve different points on the object plane. However, without evaluation Eq. 13 may still not be clear enough to show exactly how the resolution is affected, so we have included the following figures which plot some informative values to support our experimental observation. Figure 3(a) compares the theoretical first-order unresolved image of three slits with the second-order fluctuation correlation image calculated using Eq. 13. The filter function is calculated for a fiber diameter D F of 0.105 mm and varying distances x F from the center of the Fourier transform plane. Note that the plot for x F = 0.18 mm demonstrates similar behavior to the observed experimental data, including the shift of the left and right peaks away from center. It is clear that, for a lens diameter of 1.36 mm, the gold transparent plot is completely unresolved. However, as the cutoff frequency κ 2 = kx F /f B increases, the second-order resolution also increases as seen in the increasing separation of the peaks in Fig. 3. This is more clearly illustrated in Fig. 3(b) and (c); in (b) the imaging function of the first-order image (black) is plotted with the second-order imaging function (teal) for x F = 0.18 mm. It is clear that the second-order imaging function is much narrower. In (c) the half-width of the imaging function in Eq. 13 is calculated at x 1 = 0 for increasing values of x F . Using the estimated experimental parameters, the second-order imaging resolution starts equal to the first order at x F = 0 and increases to the experimental setup's limit of about 0.002 mm.
It is evident from the experiment and theoretical calculations that the increase in spatial resolution is strongly dependent on the chosen spatial filter. It is, in effect, applying a high-pass spatial filter to one arm, producing an edge-sharpening effect 11,12 . The interesting part is that the correlation of the spatially filtered intensity fluctuations with arm one produces a resolved image, especially since neither arm "sees" a resolved image. This correlation filters out the lower spatial frequencies of the unresolved image of arm one, yielding a resolved image in the intensity fluctuation correlation rather than the intensity.

Conclusion
In summary, by using a high-pass spatial filter in the non-resolving side of a two-arm camera, the measurement of the intensity fluctuation correlation 〈 Δ I 1 (ρ 1 )Δ I 2 (ρ 2 )〉 was able to resolve an object that could not be resolved by a traditional camera. This imaging method would be particularly useful for long-distance imaging in situations where it is impractical to have large lenses but high resolution is still desired, as it could take advantage of the large angular size of the sun, 0.5° relative to the earth, and the correspondingly small coherence length, on the order of 0.2 mm. In addition, since the thermal light image in 〈 Δ I 1 (ρ 1 )Δ I 2 (ρ 2 )〉 is in general turbulence-free 14 , this method would be particularly attractive for satellite cameras taking high resolution images of objects on the ground. Technically more complicated optics or electronics for practical sunlight imaging will be discussed separately.