Superresolution via Structured Illumination Quantum Correlation Microscopy (SIQCM)

We propose to use intensity correlation microscopy in combination with structured illumination to image quantum emitters that exhibit antibunching with a spatial resolution reaching far beyond the Rayleigh limit. Combining intensity measurements and intensity auto correlations up to order $m$ creates an effective PSF with FWHM shrunk by the factor $\sqrt{m}$. Structured Illumination microscopy on the other hand introduces a resolution improvement of factor 2 by use of the principle of moir\'e fringes. Here, we show that for linear low-intensity excitation and linear optical detection the simultaneous use of both techniques leads to an in theory unlimited resolution power with the improvement scaling favorably as $m + \sqrt{m}$ in dependence of the correlation order $m$. Hence, yielding this technique to be of interest in microscopy for imaging a variety of samples including biological ones. We present the underlying theory and simulations demonstrating the highly increased spatial superresolution, and point out requirements for an experimental implementation.


INTRODUCTION
Superresolution optical far-field microscopy has undergone a tremendous evolution since roughly two decades ago it was shown that the classical resolution limit [1,2] posed by diffraction can be overcome [3,4], resulting in the development of a large variety of methods achieving superresolution. One group of methods relies on stimulated ground or excited state depletion and a non-linear response of fluorescence markers to given excitation intensities to deterministically engineer the effective excitation point spread function (PSF) [3][4][5][6]. Other methods stochastically localize single photoswitchable molecules with an accuracy of a few ten nanometers via centroid fitting of the PSF [7][8][9][10]. Another branch of methods makes use of higher-order intensity cross correlations in the Fourier plane [11][12][13] or auto correlations in the image plane of a microscope [14][15][16][17]. For the latter group of correlation microscopy techniques (CM), either super-poissonian bunched light emission due to statistical fluctuations [14] or sub-poissonian anti-bunched light emission of fluorescence markers can be used to enhance the resolution, both in widefield [16] and confocal microscopy [17]. Finally, structured illumination microscopy (SIM) leads to a doubled resolution by use of the principle of moiré fringes and linear wave optics [18], and its non-linear derivative saturated SIM (SSIM) leads to further improvements and in principle unlimited resolution, though at the cost of necessitating high intensities [19]. Other derivatives combine SIM with the third-order process of CARS or with graphene plasmons to access more higher spatial frequency information than ordinary SIM [20][21][22]. Note that sub-wavelength phenomena can also be found in other fields of physics, for example in sub-wavelength atom localization due to the non-linear behavior of coherent population trapping (CPT) [23][24][25] and sub-wavelength lithography via Rabi-oscillations [26,27]. CPT was also proposed to highly increase the resolution in a microscopy themed derivative [28].
Here we report on a novel superresolution method that relies on intensity correlation measurements in the image plane of a microscope in combination with structured illumination to image fluorophores that exhibit anti-bunching. We therefore term it Structured Illumination Quantum Correlation Microscopy (SIQCM). Linear low-intensity excitation and linear detection suffice such that the technique holds promise to highly enhance the resolution in biological imaging. Detrimental effects due to high intensities that are required by many superresolution techniques, leading to phototoxicity and photobleaching in fluorophores, do not arise. We demonstrate that already very low correlation orders m provide highly enhanced superresolution, that scales favorably as m + √ m. The present manuscript focuses on the highly enhanced lateral resolution using a simple widefield microscopic setup, however one can easily extend the scheme as CM as well as SIM each on their own already provide optical sectioning capability for 3D imaging [16,29,30].
The technique makes use of mth-order correlations and antibunched photon emission, inherently present in most common fluorophores, even at room temperature [31][32][33][34]. Hence, the required quantum emitters are already broadly in use in fluorescence microscopy. Furthermore, SIM is a well established technique in biological imaging with commercial microscopes, attaining the theoretically arXiv:1702.04319v1 [physics.optics] 14 Feb 2017 predicted resolution enhancement, being widely spread.

THEORY
Let h(r) be the PSF of a given microscope, where r denotes the position in the image plane and H(k) ≡ F T {h(r)} is the corresponding optical transfer function (OTF) obtained by Fourier transform (F T ), where k denotes the spatial frequency in reciprocal space. Later, in the mth-order correlation microscopy signals CM m the effective PSF reads h m (r) ≡ (h(r)) m and its corresponding OTF shall be defined as H m (k). In general h m (r) gets narrower for increasing correlation order m and its full width half maximum (FWHM) approximately scales as 1/ √ m. Vice versa, the observable region in reciprocal space is increased for H m (k) by √ m. Microscopes usually possess the circularly symmetric Airy disk (2J 1 (r)/r) 2 with r = |r| as PSF [35], what allows to resolve individual incoherent emitters as individual sources of radiation as long as their separation d is at least on the order of d ≥ λ/2 or more precisely d ≥ 0.61 λ/A [1,2]. We denote the Rayleigh limit as d R ≡ 0.61 λ/A, with A the numerical aperture of the microscope objective and λ the wavelength of the emitted fluorescence light. W.l.o.g we assume a magnification of one (or rather minus one) throughout our theoretical treatment such that the coordinates in the object and image plane R and r, respectively, can be regarded as equal, i.e. R ≡ r.
To measure fluorescence photons in the image plane the fluorophores in the object plane need to be driven by an excitation light field. In classical linear optics the fluorophores respond linearly to a given excitation intensity I 0 . Treating the fluorophores quantum mechanically as a two-level system with ground |g and excited state |e , however, this is only the case for intensities I 0 I sat , where the saturation intensity I sat ∝ 1 τ 2 l depends on the exited states lifetime τ l . The general expression for the intensity emitted by a two-level system driven by a given excitation intensity (in units of the I sat ) reads [36] I ∝ 1 2 In ordinary classical microscopy with fluorophores that possess lifetimes on the order of a few ns or below intensities usually remain in the linear regime. To induce nonlinear responses, e.g. required by STED microscopy [3] or SSIM [19], very high intensities are necessary that are accompanied by detrimental effects to biological imaging. In contrast, our approach contents with low intensity and linear response of fluorophores to achieve highly increased superresolution. Let us first assume a continuous and spatially uniform excitation illumination in the object plane with intensity I str (r, t) = I 0 and the fluorophore density distribution n(r) ∝ N i=1 δ(r − r i ) to be comprised of individual point-like sources at positions r i that emit statistically independent, i.e. incoherent radiation. Note that we can also assign (relative) weights to the independent emitters in case their photon emission rates differ. Differences in (relative) emission rates would be enhanced in the (higher-order) intensity auto correlations. However, usually fluorophores emit sufficiently uniform and our technique does not require very high correlation orders to achieve highly enhanced superresolution, in contrast to SOFI [14]. Further, in SOFI this problem is resolved by using balanced cumulants [37] and our higher order correlation signals can be adapted accordingly. Therefore, and to keep the analysis illustrative we consider uniform emission rates here.
Considering a linear response, i.e. I 0 I sat , the intensity in the image plane reads i is the positive frequency part of the electric field operator and σ − i is the lowering operator acting on the fluorophore at r i , which can be approximated by a two-level system with ground and excited states |g i and |e i . The phases φ i are varying randomly and independently on time scales larger than the excited states lifetime τ l and introduce the incoherence as the expectation value e iφi e −iφj = 0 for i = j. Note that the intensity I(r) ≡ G (1) (r) can be recognized as Glauber's first-order equal-time intensity correlation function G (1) (r 1 , t 1 ; , assuming an ergodic system [38].
Taking the square ( we obtain an incoherent sum of narrowed PSFs h 2 (r − r i ), however in addition also the detrimental cross terms. These cross terms can be removed by subtracting the second-order intensity auto correlation function G (2) Here, the squared terms h 2 (r) vanish as each twolevel system can emit only one photon simultaneously, what is the sought-after anti-bunching CM 2 signal [15][16][17]. Higher-order CM m signals are derived analogously taking into account higher-order correlation functions up to G (m) (r). The resolution enhancement of this signal moderately scales as √ m with the correlation order, what is also illustrated in reciprocal space by the central (blue) circles in Fig. 2 that define the observable regions for ordinary intensity measurements, CM 2 and CM 3 (from left to right). Now, considering a two-dimensional structured illumination I str (r, t) = I 0 [ 1 2 + 1 2 cos(k 0 r + ϕ)], a linear response of the fluorophores and ordinary intensity measurements one obtains a doubled resolution by the principle of moiré fringes. The illumination pattern and the investigated sample produce beat patterns in the object and the image plane such that initially unobservable spatial frequencies in reciprocal space are shifted by the amount k 0 = |k 0 | = k 2 x + k 2 y into the observable region and thus can be accessed (cf. left side in Fig. 2). In general it is useful to define I str = I str (r, α, ϕ), where α = tan(k y /k x ) is the orientation and ϕ is the adjustable phase of the pattern. Note that larger k 0 effectively enlarge the observable region in reciprocal space and thus the resolution by a higher amount, however k 0 is limited by diffraction and the given numerical aperture A of the microscope objective. Hence, by use of far field wave optics infinitesimally dense fringe spacings in the source plane can not be produced. Following the derivation for Eq. (2) with adjusted I str (r), the resulting signal reads (see also Ref. [18]) Rewriting this expression into h(r) * [n(r) · I str (r, α, ϕ)] and taking the Fourier transform yields where we used the convolution theorem and the identity F T {e ik0r g(r)} =g(k − k 0 ). The density in reciprocal space is denoted byñ(k) which arises together with its shifted versions, offset by ±k 0 . One image does not allow to separate the three individual components, such that three images with three different phases ϕ = 0, 2π 3 , 4π 3 are required, creating the linear system A n = G, where the matrix A describes the resulting system and n denotes a vector with entriesñ(k),ñ(k − k 0 ) andñ(k + k 0 ). The vector G on the right hand side of the system possesses the entriesĨ(k, α, 0),Ĩ(k, α, 2π 3 ) andĨ(k, α, 4π 3 ) which represent the Fourier transforms of the (experimentally) measured data. The system is solved by applying the inverse matrix n = A −1 G. To sufficiently cover the enlarged area in reciprocal space it is necessary to chose at least three orientations α = 0, 1π 3 , 2π 3 (cf. left side in Fig. 2) resulting in a total of 9 measurements.
Taking a non-linear fluorophore response into account higher harmonics of cos(k 0 r + ϕ) arise enabling access to higher spatial frequencies in reciprocal space via SSIM. The arising higher harmonics can be read out easily when plugging I 0 cos(k 0 r) into Eq. (1) as the excitation illumination and compiling the Fourier cosine series which reads I ∝ n b n cos(n · k 0 r) .
However, this comes at the cost of necessitating high intensities that lead to phototoxicity and photobleaching in most biological samples and fluorophores. Furthermore, the Fourier coefficients b n rapidly decrease with increasing n, such that only a limited number of higher harmonics surpasses noise inherently present in every (experimental) signal. Another drawback is the necessity for a very high number of images as each higher harmonic requires two additional phases ϕ and more orientations α are needed to cover the enlarged observable region in reciprocal space [19].
Our new approach combines the strength of both methods to enhance the already superresolving signals tremendously within the linear low-intensity regime. A schematic sketch is displayed in Fig. 1. Considering linear SI I str (r, α, ϕ) and the CM 2 signal we obtain the SIQCM signal and by Fourier transform where the Fourier components corresponding to the first higher harmonic arise and the individual disks in Fourier space, governed by H 2 (k), are enlarged by ≈ √ 2, leading to an overall resolution improvement of 2 + √ 2 ≈ 3.41. Eq. (8) is depicted in the middle of Fig. 2, where due to the additional higher harmonic images with five different phases ϕ = 0, 2π 5 , 4π 5 , 6π 5 , 8π 5 per orientation are required. resulting in a total of 20 images to sufficiently cover the highly enlarged observable area. When speed is not the major goal one can chose more orientations α to obtain a higher quality what is also considered in regular SIM. After obtaining the individual Fourier components n(k),ñ(k±k 0 ) andñ(k±2k 0 ) they need to be assembled properly in reciprocal space, that is applying the same procedure which is conducted in SIM. The extracted raw components are so far scaled by the circularly symmetric OTF H 2 (k) or by its shifted versions H 2 (k ± k 0 ) and H 2 (k ± 2k 0 ). To obtain an approximately homogeneous disk we divide the enlarged observable area (cf. Fig. 2) into subregions and rescale the components by use of a Wiener filterñ new (k) =ñ(k)/(H 2 (k) + γ), where the constant γ > 0 prevents division by zero. In general the modulus of γ depends on the signal-to-noise ratio (SNR) a given measurement provides. After assembly we apply a triangular apodization, resembling the Fourier transform of an Airy disk, to the homogeneous disk to reduce ringing in the final image [19]. The final image is obtained by taking the modulus of the inverse Fourier transform of the assembled (and post-processed) disk in reciprocal space. Using more advanced deconvolution methods proposed and applied in SIM would result in an even further enhanced resolution [39,40].

SIMULATIONS
For the simulations we chose masks with point-like emitters and calculated data as it would be detected by a CCD with discrete and finite pixels. Note that, here we are assuming perfect data, i.e. discrete intensity values matching theoretical calculations without noise. Experimental requirements to obtain the sought after SIQCM signals with preferably high SNR, i.e. sufficient statistics for the second and higher-order correlations, will be discussed later. For rescaling by use of the Wiener filter we used γ = 0.05. The pixel size and post-processed area where chosen in such a way that the offsets in reciprocal space approximately match integer numbers as we considered a real valued sinusoidal modulation. To remove the necessity for integer numbers in reciprocal space (challenging to realize in a real experiment) one can use a complex wave vector in real space [19].  Simulations illustrating the resolution power of ordinary intensity measurements G (1) (r), CM 2 (r), SIM and SIQCM 2 (r) are presented in Fig. 3, where a 3 x 3 array of independent emitters with separations d = 1.0 d R , d = 0.5 d R and d = 0.29 d R is imaged by use of the enlisted techniques. The first array is resolved by every method, as the chosen distance corresponds to the classical resolution limit, however with G (1) (r) barely resolving individual emitters. CM 2 provides a moderately increased resolution and SIM the second best resolution power. We want to point out that even though CM 2 and SIM already provide superresolved images, SIQCM 2 outperforms both methods by far and provides the highest resolution power. Reducing the source separation to d = 0.5 d R only SIM and SIQCM 2 can resolve the individual emitters and finally, for d = 0.29 d R , only our new method resolves the array. Resolving the last array corresponds to a resolution improvement of 3.45, exactly matching the theoretical prediction.
To show the resolution power of our technique that scales very favorably as m + √ m compared to ordinary CM that merely scales as √ m we also present simulations for third-order SIQCM (see illustration on the right hand side in Fig. 2). We chose six orientations α = 0, 1π 6 , 2π 6 , 3π 6 , 4π 6 , 5π 6 resulting in a total of 42 images as seven phases ϕ are required per orientation. In Fig. 4 the resulting final images for the same 3 x 3 array with d = 0.29 d R imaged by use of G (1) (r), CM 3 (r) and with three different reconstruction approaches for SIQCM 3 (r) are presented together with the corresponding observable region in reciprocal space. The source distribution that was previously just resolved by SIQCM 2 is not resolved by CM 3 but clearly resolved by the SIQCM 3 signal. For the first reconstruction approach we subdivided the observable regions in sections and rescaled the Fourier components by the Wiener filter. The resulting image (third column) is simply the modulus of the Inverse Fourier transform of the homogenous disk. To remove the ringing we further omitted small imaginary parts acquired throughout the numerical evaluation (which should be zero in theory) and cropped negative values (fourth column). Note though that this approach might not be used this easily on real data with noise. The reconstruction method presented in the last column shows the triangular apodization applied to the homogenous disk. This method has also been used to produce the images in the last column of Fig. 3. The final image is again simply the modulus of the Inverse Fourier transform. In the given simulation the approach that crops negative values to remove ringing performs best and provides the smallest FWHM of the effective PSF.
Resulting images for another mask are depicted in Figs. 5 and 6 to show that the method can be applied to arbitrary emitter distributions. In contrast to the first mask three emitters have been omitted to obtain an irregular array. The SIQCM signal again outperforms the CM and SIM signal by far, as was discussed in detail in the previous section. show the corresponding reciprocal space. Images show from left to right: G (1) (r), CM3(r) and SIQCM 3 (r) with three different reconstruction approaches using a homogenous disk, a homogenous disk with cropped negative values and with a triangular apodization.

CONCLUSION AND OUTLOOK
We introduced a new quantum imaging technique we call SIQCM which is based on the profitable merger of linear SIM with anti-bunching CM. For a linear lowintensity standing wave illumination pattern and linear detection of photon auto correlations in the image plane of a microscope our technique provides in theory unlimited superresolution with improvement scaling favorably as m + √ m with the correlation order m. Hence, it has the potential to increase the spatial resolution in imaging a variety of samples and in particular biological ones. Further, we anticipate the SIQCM concept to be applicable to super-poissonian bunched light emission (e.g. used in SOFI, due to on-off blinking of fluorophores), where auto correlations in the image plane can be combined into cumulants that equally lead to a signal with narrowed PSF. Adding structured illumination will not only introduce offsets by ±k 0 but also higher harmonics with offsets up to ±mk 0 . Optical sectioning capability provided by CM as well as SIM can also be implemented enabling three-dimensional imaging with most likely increased axial resolution, due to higher harmonics in zdirection compared to 3D-SIM [30].
Our new SIQCM approach would bring similar benefits to two-photon microscopy [41], and vice versa. Considering a standing wave excitation pattern with wavelength within the red or near infrared part of the spectrum, short wavelength photons from the UV or blue part of the spectrum are emitted by fluorophores due to two-photon absorption. Since the absorption cross section is inherently dependent on the squared excitation intensity the resulting effective illumination structure is of the form I str (r, t) = I 0 [ 1 2 + 1 2 cos( k0 2 r + ϕ)] 2 , where k0 2 corresponds to the near infrared illumination wavelength and shifts by ±k 0 (corresponding to the fluorescence wavelength) as used in regular SIM already appear in the fluorescence intensity signal. Evaluating correlations additionally would result in taking the 2m-th power of the structured illumination resulting in higher harmonics with offsets up to ±mk 0 . The well-known advantages of two-photon microscopy, high penetration depth, energy deposition (and thus photobleaching) only within the vicinity of the focal plane and inherent optical sectioning capability would be added to our highly improved superresolution.
Using bunched light emission our approach should be applicable with state of the art technology and reasonable speed as SOFI already provides acquisition times of a few seconds. To obtain the sought-after CM and SIQCM signals auto correlations in the image plane can also be determined by evaluating cross correlations of neighboring pixels, what reduces experimental requirements and introduces an effectively denser sampling in the image plane [15,16]. The latter fact is of practical importance as resolutions achievable with SIQCM often exceed the sampling density of CCD cameras in use and thus interpolation can be circumvented.