Analysis of superresolution via 3D structured illumination intensity correlation microscopy

Intensity correlation microscopy (ICM), which is prominently known through antibunching microscopy or super-resolution optical fluctuation imaging (SOFI), provides superresolution through a correlation analysis of antibunching of independent quantum emitters or temporal fluctuations of blinking fluorophores. For correlation order $m$ the PSF in the signal is effectively taken to the $m$th power, and is thus directly shrunk by the factor $\sqrt{m}$. Combined with deconvolution a close to linear resolution improvement of factor $m$ can be obtained. Yet, analysis of high correlation orders is challenging, what limits the achievable resolutions. Here we propose to use three dimensional structured illumination along with $m$th-order correlation analysis to obtain an enhanced scaling of up to $m+m=2m$. Including the stokes shift or plasmonic sub-wavelength illumination enhancements beyond $2m$ can be achieved. Hence, resolutions far below the diffraction limit in full 3D imaging can potentially be achieved already with low correlation orders. Since ICM operates in the linear regime our approach may be particularly promising for enhancing the resolution in biological imaging at low illumination levels.


Introduction
Since it was first shown that the resolution limit, posed by diffraction, can be overcome [1], a variety of superresolution microscopy methods have been developed.Yet, each technique comes with certain requirements and limitations, thus justifying an ongoing pursuit of novel methods.One group of methods relies on stimulated ground or excited state depletion and a non-linear fluorophore response to deterministically engineer the effective excitation point spread function (PSF) [1][2][3].Other methods stochastically localize single photoswitchable molecules via centroid fitting of the PSF [4][5][6][7].
Another branch of methods makes use of intensity correlations that are evaluated from an image series of a widefield microscope [8,9].For these correlation microscopy techniques, statistically blinking fluorophores [8] or quantum emitters that exhibit anti-bunching [9] can be used to enhance the resolution by shrinking the effective PSF by the factor √ m (with correlation order m).Especially the first approach, known as superresolution optical fluctuation imaging (SOFI), is widely applied due to its combination of resolution improvement with low complexity of use [10][11][12].Yet, in practice high correlation orders are not evaluated due to strong brightness skewing in the final SOFI images and long measurement times that are required for a reliable evaluation of high correlation orders [10].This problem together with the moderate √ m scaling currently prevents SOFI to obtain resolutions far below the diffraction limit.
In parallel, structured illumination microscopy (SIM) was developed, where by the use of spatial frequency mixing the resolution is doubled within the linear wave optics regime [13].The non-linear derivative saturated SIM (SSIM) leads to an in principle unlimited resolution, though at the cost of necessitating high intensities [14].Other derivatives combine SIM with the thirdorder process of CARS [15,16], or with surface plasmons [17][18][19] to access higher spatial frequency information.3D-SIM doubles both, the lateral and axial resolution [20].Standing wave fluorescence microscopy (SWFM) techniques [21,22] highly enhance the axial resolution via an axially structured illumination, but not the lateral one.Double-objective illumination and detection techniques [23,24] with 3D-SIM attain the axial resolution of SWFM and the lateral one of 2D-SIM [25].The SIM toolbox is now considered to be one of the most powerful and versatile superresolution techniques, due to its combination of resolution improvement with good acquisition speed and flexibility of use [26].
Recently it was shown that SIM and correlation microscopy based on antibunching can be combined to enhance the lateral resolution [27].For correlation order m the enhancement scales favorably as m + √ m, which is a large improvement over the moderate √ m scaling of antibunching microscopy.However, measurements of quantum correlations that exhibit antibunching require experimental efforts and long measurement times.Current demonstrations are thus limited to proofof-principle experiments.
Here, we propose to use 3D structured illumination [20] to enhance the resolution of SOFI, a new technique which we call structured illumination SOFI (SI-SOFI).The advantage is an easy implementation with the same m+ √ m scaling, valid for lateral and axial resolution enhancements.Already for second-order SI-SOFI the factor of 3.41 approximately equals the resolution enhancement √ 12 = 3.46 of 12th-order SOFI, and for third-order SI-SOFI the factor 4.73 even matches the resolution of 22th-order SOFI.We point out that both SOFI and SIM operate within the linear regime and are established techniques in the field of superresolution microscopy.SI-SOFI thus bears the potential for deep-subwavelength resolution at low illumination levels, especially relevant for imaging of biological specimen.We present the theory and show simulations that demonstrate the theoretically predicted resolution enhancement, which scales favorably with m + √ m.

Theory
In this section we show that 3D-SI-SOFI leads to a resolution enhancement of m + √ m in full 3D imaging.We assume R ≡ r for the coordinates in the object and image plane, respectively, i.e. a magnification of one.Let h(r) be the 3D PSF of a given widefield microscope, with r = (x, y, z).H(k) ≡ F T {h(r)} denotes the corresponding 3D optical transfer function (OTF) obtained by Fourier transform (F T ) of h(r), where k = (k x , k y , k z ) denotes the spatial frequency in reciprocal space.The lateral and axial width of the PSF determine the resolution power a microscope provides to discern individual close-by emitters.A shrunk PSF equals an increased resolution.In reciprocal space the resolution is equally determined via the size of the OTF, where a larger OTF support equals a higher resolution.
The lateral width of the PSF is usually smaller than the axial one.Moreover the axial resolution can not properly be defined in widefield microscopy, which is due to the missing z-cone in Fourier space (see Fig. 1a) and out-of-focus light being captured by a planar 2D detector placed in the image plane [28].Optical sectioning capability in z-direction can however be retrieved by the measurement of z-stacks and deconvolution, using a pinhole as in confocal microscopy, or a variety of other means.In SOFI, for example, the missing z-cone is intrinsically removed [29], as will also be shown below (after Eq. ( 5)).
The 3D PSF and 3D OTF of real imaging systems are generally computed numerically.Analytical expressions can only be derived for a few idealized systems.Here, as an approximation, we consider a 3D Gaussian PSF of the form [8] where w ρ and w z denote the lateral and axial width of the PSF, and ρ = (x 2 + y 2 ) −1/2 .The OTF is then also a 3D Gaussian in reciprocal space (see Fig. 1b).We choose the ratio w z /w ρ = 3.0 to mimic a typical widefield microscopy PSF, where the resolvable distances along ρ and z between two close-by emitters are given by [30] ∆ρ min = 0.61 In Eq. ( 2), A is the numerical aperture of the microscope objective, n the refractive index and λ the wavelength of the emitted fluorescence light.For simplification purposes we assume equal wavelengths for excitation and emission, i.e. λ ex = λ em ≡ λ, throughout this paper.∆ρ min is the well-known Rayleigh limit [31,32], whereas the axial resolution limit ∆z min is rather a formally assigned quantity, due to reasons discussed above.The fluorophores are driven by monochromatic light of intensity I str (r) and we assume a linear response, i.e. a driving intensity (far) below saturation.In widefield microscopy and SOFI plane-wave illumination leads to the flat excitation intensity I str (r) = I 0 .
We describe the model system to be imaged by the fluorophore density distribution n(r, t) ∝ N i=1 δ(r − r i )s i (t), which is comprised of individual point-like emitters at positions r i that emit incoherent fluorescence radiation with time-dependent intensity s i (t).In classical microscopy the emitters can be assumed to have a constant intensity s i (t) = si .In SOFI the emitters are statistically blinking since they randomly and independently enter a bright or a dark state.The timedependent intensity of each emitter can thus be described by s i (t) = si + ∆s i (t), where si is the average intensity and ∆s i (t) describes the fluctuations with zero-mean expectation value.
Considering the convolution by the PSF (expressed through the operation " * "), each time-resolved image taken in the image plane reads The image average takes the form I(r, t) = Ī(r) = N i=1 h(r − r i )s i , which is the tradionally measured quantity.The corresponding second-order equal-time spatial intensity auto correlation function (r 1 , t 1 = r 2 , t 2 = r, t), calculated by classical means, is given by The first summand is the squared mean intensity.In the second summand, however, the individual emitter positions are convolved with the squared PSF and the contributions are scaled by the variance of the emitter intensity fluctuations.In classical microscopy one would assume ∆s i (t) 2 = 0, what highlights the importance of blinking.The cross terms ∆s i (t)∆s j (t) in the second line drop out for i = j due to statistical independence of the emitters.To separate the superresolving signal of the second summand from the first summand, it is necessary to evaluate the so-called cumulant, which for second order reads Due to the squared PSF the resolution in the final image is enhanced by the factor √ 2. In Fourier space, vice versa, the new OTF F T {h(r) 2 } = H(k) * H(k) is enlarged and the missing z-cone is filled, what results in true optical sectioning capability [29].Cumulants of higher order m are evaluated in the same manner by a combination of intensity correlation functions up to the same order m [8].The final SOFI images will have the form N i=1 h(r − r i ) m ∆s i (t) m , where the PSF is taken to the mth power and scaled by ∆s i (t) m .Correspondingly the resolution is enhanced by the factor √ m.Even though high correlation orders can in principle be evaluated, the moderate √ m scaling prevents resolutions far below the diffraction limit.In addition different molecular brightnesses and/or blinking ratios ∆s i (t) become more and more pronounced with rising order m, what skews the final image [10].While this problem is mitigated by use of balanced cumulants [33], in practice still only low correlation orders are utilized.
Since the absolute scaling is not relevant in the mathematical treatment we will consider equal average intensities si ≡ 1 and cumulant scalings ∆s i (t) m ≡ 1 for all i = 1, . . .N emitters and m = 2, 3, . . .correlation orders.[34] The final SOFI images then read with the PSF h m (r) ≡ h(r) m .The corresponding OTF is H m (k) ≡ F T {h m (r)} and will be required below.
Another means to enhance the resolution is to use structured illumination I str (r), which for the 2D case reads I str (r . By this the individual emitter intensities are scaled by I str (r i ), what for SIM results in measured images of the form [13] where α = tan −1 (k y /k x ) denotes the orientation of the structured illumination and ϕ is an adjustable phase.
The multiplication n(r) × I str (r, α, ϕ) in real space corresponds to a mixing of the object's spatial frequencies k with the spatial frequency k 0 of the known standing wave illumination pattern in Fourier space.Hence, initially unobservable spatial frequencies k > k max are encoded in the microscope's OTF support.Taking a set of linearly independent images and applying computational post-processing allows for the retrieval of this information.The resolution enhancement reads (k 0 + k max )/k max and reaches two for diffraction limited illumination with k 0 = k max .In 2D-SIM two plane waves are superposed which stem from the ±1 diffraction orders of a grating, what leads to the electric field distribution E = e ikxx+ikzz + e −ikxx+ikzz .The corresponding intensity pattern reads

see the previous definition). In 3D-SIM the central beam is added, what results in the electric field
E = e i(kxx+kyy)+ikzz + e ikz + e −i(kxx+kyy)+ikzz .(8) Note that k = 2π/λ is the wavenumber of the illumination and the spatial frequencies k = (k x , k y , k z ) have to fulfill (k 2 x + k 2 y + k 2 z ) −1/2 = k to obtain a propagating wave.The resulting 3D intensity pattern reads [20] where we added independent variable phases along the lateral and axial direction.This pattern contains seven spatial frequency components.Taking the Fourier transform yields F T {I(r)} = I 0 7 j=1 e iϕj δ(k − k j ), with seven delta peaks in Fourier space and phases ϕ j (which contain combinations of ϕ r and ϕ z ).For ±1 diffraction orders propagating at the angles ±60 • the positions read x + k 2 y ) −1/2 .These are shown in Fig. 1c, where the axes are normalized to the wavenumber k = 2π/λ.The central (blue) dot corresponds to plane wave illumination while the outer green dots arise due to the spatial frequencies of the standing wave pattern.
Using this 3D illumination, taking the Fourier transform of the expression in the first line of Eq. ( 7) and Image f) depicts the Fourier transform of the squared structured illumination, where the outer (red) dots represent the contributions from the first higher harmonics.Again, combining images e) and f) yields the OTF supports displayed in g), i.e., of second-order SOFI, second-order SOFI with 2D-SIM for a single α, second-order SOFI with 3D SIM for a single α and second-order SOFI with 3D-SIM for four orientations α = 0, π 4 , 2π 4 , 3π 4 .The total support for this case is already enhanced by the factor 3.41 along all axes.
utilizing convolution theorems yields [20] where k j , ϕ j and c j , respectively, represent the center positions, the corresponding phase factors and weights [see prefactors of cosines in Eq. ( 9)].Note that ñ(k) = F T {n(r)} contains the sought-after spatial frequency information of the unknown object under investigation.In a single image all seven components are superposed and can not be disentangled.Similarly to 2D-SIM, at least seven independent 3D images (e.g. by measurement of z-stacks and deconvolution) are required.This is established by varying the phases ϕ r and ϕ z and creates the linear system An = G, where the matrix elements of A are given by the phase terms e iϕj of Eq. ( 10) and n denotes a vector with entries ñ(k − k j ).The vector G possesses Eq. ( 10) as entries, which represent the Fourier transforms of (experimentally) measured data.The system is solved by applying the inverse matrix n = A −1 G.
So far a single orientation α was considered.To sufficiently cover the enlarged OTF support the procedure needs to be repeated for three orientations α = tan −1 (k y /k x ) = 0, 1π 3 , 2π 3 (see the final image in Fig. 1d).In the next step all disentangled components are shifted back to their true positions in Fourier space and superposed, leading to an effectively doubled OTF support size.In addition, one can optimize the Fourier components by rescaling them via the application of a Wiener filter, leading to ñnew (k) = ñ(k)/(H 2 (k) + γ).This establishes a mostly constant scaling over the OTF sup-port.The constant γ > 0 prevents division by zero and should be chosen noise dependently.For high noise levels a larger γ should be chosen.Finally a triangular apodization can be applied to reduce ringing in the final image [13].The final image of 3D-SIM is obtained by the inverse Fourier transform of the post-processed and merged Fourier components.Alternatively the assembly can be conducted in real space by first applying the inverse Fourier transform to each individual component and then multiplying by the complex wave e ikj r .This approach avoids the need for near-integer numbers (k x , k y , k z ) j in Fourier space [13].Note that in 3D-SIM the missing z-cone shown in Fig. 1a can be filled, leading to true 3D imaging capability [20], yet still limited by diffraction to a two-fold improvement.
SI-SOFI overcomes the limitations of linear SIM and SOFI by fruitfully combining both.That is, the structured illumination encodes information from outside the original OTF support via spatial frequency mixing and and the correlation analysis effectively raises all signals to the mth power.In mathematical terms this is established by the combination of Eqs. ( 6) and ( 7) which results in [27] SI-SOFI m (r) = Now, h(r) and I str (r i , α, ϕ) are taken to the mth power, such that higher harmonics up to cos(mk 0 r) arise and the individual OTF H m (k) in Fourier space is enlarged by √ m.For 2D-SIM combined with correlation analysis the lateral resolution was shown to be enhanced by the factor m + √ m [27].To illustrate the outcome of Eq. (11) we use the 3D structured illumination of Eq. ( 9) and consider secondorder SI-SOFI.Taking the Fourier transform of Eq. ( 11), we thus obtain (for m = 2) c j e iϕj ñ(k−k j ) , (12) where 19 spatial frequency components arise.Higher harmonics now arise along the lateral and axial direction, as illustrated in Fourier space by the added most outer (red) dots in Fig. 1f.Moreover, the OTF is enlarged by √ 2 (see Fig. 1e).The total improvement thus reaches 2 + √ 2 ≈ 3.41 along all axes (see the 3D volume in Fig. 1g).The larger OTF H 2 (k) also reduces the need for many orientations α such that four α = 0, 1π 4 , 2π 4 , 3π 4 are sufficient [27].The image reconstruction process follows the one introduced for 3D-SIM, where now 19 independent 3D images need to be acquired.

Simulations
To illustrate the mathematical description of section 2, we performed a couple of simulations, where we chose the 3D PSF of Eq. ( 1).Since we are interested in the resolution power of SIM, SOFI and SI-SOFI relative to wide-field microscopy we use dimensionless units and merely normalize r by the Rayleigh limit ∆ρ min of Eq. ( 2).We chose a setup with three close-by emitters distributed in the 3D object space within a sub-diffraction limited volume (see Fig. 2a).The coordinates r 1 = (−0.16,0.16, 0.05), r 2 = (0.26, −0.26, 0.57) and r 3 = (0.26, −0.26, 0.68) with pair-wise lateral or axial separations where chosen to demonstrate the full 3D resolution capabilites of SI-SOFI.We calculated the 3D data as it would be obtained by use of a pixelated CCD detector in combination with the measurement of z-stacks and deconvolution.For illustration purposes we used theoretical data without noise.For a discussion and details on the noise analysis we refer the reader to refs.[8,9,27].For the post-processing we chose γ = 0.05, both for SIM and SI-SOFI.For SIM the phases ϕ r = 0, 2π 5 , 4π 5 , 6π 5 , 8π 5 and ϕ z = 0, 2π 3 are required per orientation α.In SI-SOFI the required phases are ϕ r = 2π 9 j (j = 0, . . ., 8) and ϕ z = 0, 2π 3 , 4π 3 .Fig. 2 shows the simulations results for six different signals.These are a) widefield microscopy, b) secondorder SOFI, c) 3D-SIM, d) fourth-order SOFI, e) secondorder SI-SOFI and f) twelfth-order SOFI.Images d) and f) were chosen for comparison purposes since they provide improvements of √ 4 = 2.0 and √ 12 = 3.46, respectively, which are close to the ones expected from 3D-SIM (see Fig. 2c) and second-order SI-SOFI (see Fig. 2e).The simulation results are in good agreement with the theoretical predictions.It can be seen that SI-SOFI equally enhances the lateral and axial resolution and thus provides full 3D superresolution.

Conclusion
In this paper we proposed to enhance the resolution power of SOFI by adding 3D structured illumination to the setup.We presented a mathematical treatment that predicts for SI-SOFI a resolution enhancement of m + √ m for correlation order m, both along the lateral and axial direction.We confirmed the results via simulations that matched the theoretical predictions.Moreover we pointed out that SIM in combination with SOFI is easy to implement and clearly outperforms each method on its own.We note that since SI-SOFI fully operates within the linear regime it bears the potential to increase the resolution in particular in imaging biological specimen at low illumination levels, i.e., especially in cases where other methods can not be utilized.
Further enhancements along the axial direction can be achieved by combining SOFI with the double-objective 3D-SIM technique known as I 5 S [25].The three added coherent beams from the second objective lead to very fast modulations along the axial direction.Additionally, the OTF is enlarged along the axial direction as in 4Pi-and I 5 -microscopy [23,24].Again, SI-SOFI would square the PSF and the excitation pattern.
Another promising future route may be to combine SI-SOFI with plasmonic SIM (PSIM) techniques [17][18][19].Even though these techniques are limited to 2D, they al- low the exploitaton of spatial frequencies k 0 > k max of the standing wave pattern.In PSIM k 0 = 2k max should not be exceeded to prevent gaps in the OTF support coverage.The resolution enhancement is thus limited to the factor 3 [18,19].SI-SOFI, however, would highly benefit from spatial frequencies k 0 > 2k max since the enlarged OTF H m (k) prevents an early formation of gaps and the higher harmonics cos(mk 0 r) reach out to very high spatial frequencies.

Fig. 1 .
Fig. 1.Illustrations of total OTF supports of 3D-SIM (left column) and second-order 3D SI-SOFI (right column).Image a) depicts the OTF of a regular widefield microscope and image b) the approximation H(k) in form of a 3D Gaussian.Image c) shows the Fourier transform of the structured illumination of Eq. (9) with the center positions kj = (kρ, kz)j (j = 1, . . ., 7) (see the central blue and outer green dots) given in the main text.Combining a) and c) yields the images in d), where the OTF support of widefield microscopy, 2D SIM for a single orientation α, 3D-SIM for a single orientation α and 3D-SIM for three orientations α = 0, π 3 , 2π 3 are shown.The final image of d) provides a two-fold enlarged support along all axes.Image e) shows the OTF H2(k) of second-order SOFI, which is enlarged by the factor √ 2 along all axes.Image f) depicts the Fourier transform of the squared structured illumination, where the outer (red) dots represent the contributions from the first higher harmonics.Again, combining images e) and f) yields the OTF supports displayed in g), i.e., of second-order SOFI, second-order SOFI with 2D-SIM for a single α, second-order SOFI with 3D SIM for a single α and second-order SOFI with 3D-SIM for four orientations α = 0, π 4 , 2π 4 , 3π 4 .The total support for this case is already enhanced by the factor 3.41 along all axes.
A.C. and J.v.Z.gratefully acknowledge funding by the Erlangen Graduate School in Advanced Optical Technologies (SAOT) by the German Research Foundation (DFG) in the framework of the German excellence initiative.G.S.A. thanks the Air Force Office of Scientific Research (Award number: FA 9550-18-1-0141) for supporting this work.