Computational polarized Raman microscopy on sub-surface nanostructures with sub-diffraction-limit resolution

: Raman microscopy with resolution below the diffraction limit is demonstrated on sub-surface nanostructures. Unlike most other modalities for nanoscale measurements, our approach is able to image nanostructures buried several microns below the sample surface while still extracting details about the chemistry, strain, and temperature of the nanostructures. In this work, we demonstrate that combining polarized Raman microscopy adjusted to optimize edge enhancement effects and nanostructure contrast with fast computational deconvolution methods can improve the spatial resolution while preserving the flexibility of Raman microscopy. The cosine transform method demonstrated here enables significant computational speed-up from O ( N 3 ) to O ( N log N ) - resulting in computation times that are significantly below the image acquisition time. CMOS poly-Si nanostructures buried below 0.3 − 6 µ m of complex dielectrics are used to quantify the performance of the instrument and the algorithm. The relative errors of the feature sizes, the relative chemical concentrations and the fill factors of the deconvoluted images are all approximately 10% compared with the ground truth. For the smallest poly-Si feature of 230 nm, the absolute error is approximately 25 nm.


Introduction
Raman microscopy, or micro-Raman spectroscopy, is a powerful tool to investigate the fundamental properties of materials including the chemical composition, strain, and temperature distribution. [1,2] In the past two decades, multiple research fields have been significantly impacted by this technique. For example, for low-dimensional materials such as graphene and MoS 2 , Raman spectroscopy is used to study the phonon dispersion and the inter-layer interaction because Raman scattering directly probes the vibration modes. [3,4] In biological research, Raman microscopy is used to image living cells because it is non-invasive and label-free. [5,6] However, compared with other widely used characterization tools such as scanning electron microscopy (SEM), [7] atomic force microscopy (AFM) [8] and stochastic optical reconstruction microscopy (STORM), [9] Raman microscopy has poor spatial resolution due to the optical diffraction limit. This critical disadvantage limits the application of Raman microscopy when higher spatial resolution is required. In this work, we demonstrate that combining polarized Raman microscopy adjusted to optimize edge enhancement effects and nanostructure contrast with fast computational deconvolution methods can improve the spatial resolution while preserving the flexibility of Raman microscopy. The cosine transform method demonstrated here allows for significant computational speed-up from O(N 3 ) to O(N log N) and addresses the required constraints including non-differentialable regularization and non-negative pixel values.
In the literature, various techniques have been proposed to enhance the spatial resolution of Raman microscopy. For example, researchers have modified the instruments to couple Raman scattering with near-field enhancement effects. One successful example is the tip-enhanced Raman spectroscopy (TERS), which is based on the near-field enhancement from the metallic tip-sample interaction. [10] A typical TERS setup is not limited by the diffraction limit but the size of the tip apex. There are multiple reports on TERS with sub-10 nm spatial resolutions. [11][12][13][14][15] However, TERS can only resolve features on the sample surface and the near-field enhancement is sensitive to the tip morphology. [10,12] These drawbacks hinder the the accessibility and universality of TERS and other similar near-field instruments.
In the past decade, researchers also demonstrated resolution enhancement of Raman microscopy using far-field methods which were originally used for fluorescence microscopy. For example, Watanabe et al. reported structured line illumination (SLI) Raman microscopy of which the spatial resolution is 1.4× higher than the theoretical limit in wide-field microscopy. [16] Based on image scanning microscopy (ISM), Roider et al. reported a fiber-coupled scanning Raman microscope which has a hexagonal packed collection fiber bundle and a multiline detector. [17] By pixel reassignment, they demonstrated a 30% resolution enhancement. The idea of localization fluorescence microscopy [9] is also adapted with surface enhanced Raman scattering (SERS). [18][19][20] Moreover, similar to stimulated emission depletion (STED) microscopy, [21] nonlinear saturation effects have been exploited with stimulated Raman scattering (SRS) and coherent anti-Stokes Raman scattering (CARS). [22][23][24] . Among these techniques, as with TERS, the SERS based approaches are not suitable for the imaging of nanostructures located several microns below the surface. The main limit of the STED-like and other nonlinear methods is that the damage threshold of the sample has to be high because intense laser sources are often required to induce the nonlinearity. [25] SLI only has spatial resolution enhancement parallel to the illumination line while the resolution in the perpendicular direction is still limited by the confocal slit width. The ISM technique is promising and universal with the spatial resolution limited by the system point-spread-function (PSF). As can be seen in this report, in principle, our approach can be combined with the ISM technique to further enhance the spatial resolution.
Without intensive modification of the instruments, conventional microscopes are used alongside computational reconstruction of the ground truth images. In this scheme, similar to most optical image processing, the observed image is treated as the convolution of the PSF with the ground truth image. Multiple deconvolution algorithms have been reported for confocal Raman microscopy. In 2003 and 2006, Tomba et al. formalized the deconvolution problem as a quadratic programming problem with Tikonov regularization to resolve liquid-polymer interfaces of 4 µm. [26,27] In 2016, Cui et al. used a Markov-Poisson maximum a posterior (MPMAP) algorithm with an iterative kernal to resolve 200 nm poly methyl methacrylate (PMMA) patterns on a Si wafer while the diffraction limit was approximately 300 nm. [28] In 2019, Winterauer et al., used a primal-dual method for interior point least square (IPLS) with the non-negativity constraint to reconstruct the image of poly-(3,4 ethylenedioxythiophene) (PEDOT) nanowires on Si substrate. [29] They demonstrate a spatial resolution of 125 nm while the diffraction limit was approximately 320 nm. Recently, researchers have considered the hyperspectral nature of Raman images in the context of resolution enhancement. In 2015, Offroy et al. coupled deconvolution algorithm and chemometrics to characterize aerosols. [30] In their algorithm, the spectra and the concentration distribution of the pure chemicals are extracted iteratively by multivariate curve restorationalternating least squares (MCR-ALS). The image of the concentration is then divided into multiple lower resolution sub-images with pixel shift and optimization is performed on this multi-frame image set. They achieved 65% enhancement of the spatial resolution based on the 'step-edge' criterion. In 2020, Winterauer et al. reported the results of resolution enhancement on PEDOT nanowires using a hyperspectral-based pipeline. [31] They extracted the convolved images using principle components analysis (PCA) and the deconvolution was achieved using IPLS method similar to their previous report. The smallest resolvable separation of two nanowires was 14 nm. In general, the computational methods can approximately enhance the spatial resolution with a factor of 30%-100%.
Computational approaches for resolution enhancement of Raman microscopy are desirable because minimal instrument modification is required and they can also be used together with other advanced techniques as long as optical diffraction influences the imaging. However, these approaches are often not well specified and several critical issues are not clearly discussed in the literature. First, in the deconvolution algorithms, the symmetry of the PSF and the boundary conditions are usually assumed but not utilized to accelerate the computation. This is due to the difficulties of incorporating the regularization terms and the non-negativity constraint into the inverse problem. As we will present in Section 4, the complexity of a straightforward deconvolution algorithm is O(N 3 ) where N is the total pixel number. Meanwhile, the total acquisition time of a Raman image scales with N. Therefore, a fast deconvolution algorithm is specifically meaningful for Raman images to lower the computational overhead. Second, the dependence of the deconvolution results on the input parameters is not fully investigated. As an ill-posed problem, the deconvolution is usually solved with regularization terms and the proper choice of the regularization parameters is non-trivial. In the literature, these parameters are usually empirically selected. Another input parameter is the PSF which is the convolution kernel and is usually experimentally measured but the effect of the measurement errors on the results is unknown. Third, most authors do not quantify the error of the deconvolved images with the ground truth or only compare the feature sizes. Other critical properties, such as the chemical concentrations, are rarely discussed while one of the key advantage of Raman microscopy is the capability of chemical analysis.
The experimental aspect of sub-diffraction-limit Raman microscopy also requires further investigation. In terms of the samples under test, no result on a buried sample has been reported but these samples are of great interest in many fields. For example, the devices buried under a few hundred nm to several µm of dielectric materials are common in modern integrated circuit (IC) industry, [32,33] where Raman microscopy is used as a characterization tool. For these samples, the PSF measurements should be performed with care. As Everall pointed out in a series of reports, [34][35][36] for high NA objectives, which are typically used for Raman setup, the existence of hundreds nm of oxide can lead to significant spherical aberration and thus the PSF may be distorted. In terms of the instrumentation, in the literature, simple Raman microscopes are usually used and only the total intensity of Raman scattering is considered while other properties of Raman scattering are not utilized. For example, it is well known that the Raman scattering is highly polarization-sensitive. [37][38][39][40] As we will present in this work, the polarization selection can be used to enhance the contrast of a Raman image.
In the work, we report our computational polarized 785 nm Raman microscopy on poly-Si nanostructures buried 0.3 − 6 µm below the surface on two complementary metal-oxidesemiconductor (CMOS) chips. The deconvolution is solved iteratively using an algorithm based on the fast cosine transform, which has complexity of O(N log N). The dependence of the deconvolved images on the regularization parameter and the size of the PSF is discussed and we suggest that one train the regularization parameter on a separate pattern. The errors of the feature sizes, the relative chemical concentrations, and the fill factors of the deconvolved images are investigated. The report is organized as follows. In Section 2, we present our Raman microscope and the test samples. In Section 3, we discuss the measurement of the PSF. In Section 4, we review the algorithm used for the deconvolution. In Section 5, we present the deconvolution results and the related discussions.

Instrumentation and samples
In this section, we present our polarized Raman microscope and the test samples. In Fig. 1(a), the schematic plot of our setup is shown. We emphasize that the polarization configuration is built for two reasons. First, as will be elaborated in the next section, the polarization selection rule can be used to eliminate the signal from crystalline substrates. Second, the polarization can be aligned with the features to enhance the Raman signal. This is known as the edge enhancement effect due to the electric field localization, which is similar to the tip enhancement effect but does not require external confinement. Researchers have utilized this effect to improve the spatial resolution on Si and Ge stripes and measure the local strain fields. [37,38,41] It has been shown that the strength of the electric field localization depends on the polarization. For example, consider a Si/SiO 2 interface, the tangential electric field is continuous while the perpendicular electric field decreases on the Si side because the dielectric constant of Si is multiple times larger than SiO 2 . The effect is more significant when the interface is further reduced to a corner. We will also discuss this effect in the next section. The excitation source is a single frequency 785 nm distributed Bragg reflector (DBR) laser diode. Here the near infrared laser is used to minimize the parasitic heating. The diode is coupled with a polarization-maintaining fiber (PMF) after which the laser light is collimated, cleaned by a narrow bandpass filter, and expanded. The beam is then reflected by a dichroic mirror and passes through a polarization beam splitter (PBS) and an achromatic half-wave plate (AHWP). The PBS keeps the excitation and the back-scattered Raman light co-polarized and the AHWP rotates this polarization direction relative to the sample. A NA= 0.95 near-infrared objective is used to focus the excitation light onto the sample. A XYZ piezo stage is used to mount the sample. The minimum step size and the bidirectional repeatability are both approximately 10 nm. The back-scattered Raman light is collected by the same objective and passes through the AHWP, PBS, DM, a 785 nm longpass filter and a beam expander. The beam is then coupled into another PMF using an off-axis parabolic (OAP) mirror. The spectrum is analyzed by a dispersive spectrometer with a grating of 600 grooves per mm and a liquid nitrogen-cooled CCD camera. By inserting a beam splitter (BS) between the AWHP and the objective, an epi-illuminated bright-field microscope can be enabled. A CMOS camera records the bright-field images which are used to roughly locate the region of interest and the focal plane.
In Figs. 1 we present two test CMOS chip samples, named A and B. Sample A and B are realized in Micron 220 nm DRAM [42] and GLOBALFOUNDARIES 55BCDLite, respectivaly. As shown in Figs. 1(b) and 1(c), the layer stacks of A and B both consist of poly-Si, crystalline Si, oxide and dielectric materials, and the test patterns are all in the poly-Si layer. The main difference is that B is an intact CMOS chip with the poly-Si layer buried under the full dielectric stack (≈ 6 µm) while A was back etched and bonded onto a handle wafer and thus there is only 200-300 nm oxide on top of the poly-Si. Also, in A, the poly-Si is relatively far away (>2 µm) from the Si handle while the poly-Si is only 300 nm away from the Si substrate in B. The detailed fabrication process of A can be found in this paper. [42]

PSF measurements
We assume that the PSF has a Gaussian shape and the full-width half-maximum (FWHM) along x and y are measured by knife edge measurements. Here the knife edge measurement refers to scanning the focal spot across poly-Si edges along x and y with the assumption that the poly-Si edges are ideally sharp. It is critical that the poly-Si layer is thin (≈ 100 nm) since the knife edge measurements on thick materials can lead to severe errors. [36] In Figs. 2(a) and 2(b) , we plot the Si Raman peak intensity near 520 cm −1 versus positions on sample A. Here the polarization is 45 • relative to x. Note that the Raman signal is enhanced near the edge. Similar to the literature, we also observe that the strength of the enhancement depends on the polarization direction. (See Supplement 1) The PSF curve fitting with edge enhancement is more complicated than the situation in which the ground truth signal is considered as a step function. Since the localization is usually within tens of nm, we assume that the ground truth of the signal is step-like plus a delta function at the edge. The fit FWHM of the PSF are 422.1 ± 24.7 nm and 393.9 ± 24.3 nm in x and y, respectively. (Figs. 2(a) and 2(b)) The mean and the standard deviation are from 20 repeated measurements at different locations of the edges. If one fits the experimental data as a Gauss error function using the falling edge without considering edge enhancement, the fit FWHM will be approximately 50 − 100 nm smaller. According to the discussion in Section 5, this can introduce 10 − 20% errors in the final results.
The same measurement is performed on sample B. (Figs. 2(c) and 2(d)) Unlike sample A, the (100) Si substrate of sample B is only 300 nm away from the poly-Si, which is the typical shallow trench isolation (STI) thickness used in CMOS. Given the diffraction limit of the focal spot in z is approximately 2 µm, the substrate has significant contribution to the Raman signal and the contrast of the poly-Si layer is low even with edge enhancement. The contrast can be improved by aligning the co-polarization direction along the [100] of the substrate Si. In principle, this should eliminate the crystalline Si Raman associated with the optical modes at ≈ 520 cm −1 while only slightly lower the poly-Si Raman because the orientation of the poly-Si is random. [39,40] Here we align our laboratory coordinate, namely xy, with the [100] of the substrate and align the polarization along x. The fit FWHM of the PSF are 409.3 ± 15.6 nm and 401.7 ± 7.2 nm in x and y, respectively. Note that though the thickness of the dielectric stack is much thicker on sample B, the PSF does not vary significantly compared with sample A.

Deconvolution algorithm
In this section, we present the deconvolution algorithm. Let D be the 3D (2D spatial plus 1D spectral) data set of the hyperspectral Raman image. The Raman spectrum at each pixel is first compressed into the relative chemical concentration using PCA. The deconvolution is then performed on the image of the chemical concentration. As mentioned before, the observed image is the convolution of the PSF and the ground truth. In this work, the deconvolved image f is solved as an optimization problem with L1 regularization Here f and g are the column-major vectorization of the deconvolved and the raw images, respectively. H is the blurring matrix determined by the PSF and the boundary conditions. λ is a hyperparameter tuning the strength of the regularization. D x and D y are the derivative matrices along x and y, respectively. ∥·∥ i is the Li-norm. In this work, we use the Neumann boundary condition to minimize the ringing effects at the boundaries.
The L1 regularization is solved iteratively with the split Bregman method. [43,44] ( , which is the barrier parameter controlling the convergence of d to Df , is updated after several iterations. The update rule of µ l is specified later. b k is chosen in each iteration according to the Bergman distance. Equation (2) can be split and solved iteratively as and b k is updated as Equation (4) has a closed form solution Equation (3) is standard quadratic programming with non-negativity constraint and it is the most time-consuming sub-problem. A simple and efficient algorithm to solve Eq. (3) can be found in this report. [45] However, in order to incorporate the fast cosine transform, the non-negativity constraint is solved iteratively using the quadratic penalty method [46,47] as where η l >0 is the barrier parameter controlling the convergence of f to p. Equation (7) is now a quadratic programming problem without constraint and has a closed form solution Since the coefficient matrix is Hermitian and positive definite, Eq. (9) can be solved first using the Cholesky decomposition with O(N 3 ) complexity and then the forward and backward substitutions with O(N 2 ) complexity. [48] Here N is the total pixel number of the observed image. Without updating µ l and η l , the Cholesky decomposition only needs to be performed once during the iterations.
As pointed out by Ng et al., given that both the PSF and the boundary conditions have mirror symmetry with respect to x and y, H is a Toeplitz-plus-Hankel block matrix and can be diagonalized by the tensor product of two cosine transform matrices. [49] Here ⊗ represents the tensor product and C is the 1D cosine transform matrix with δ being the Kronecker delta. Λ is a diagonal matrix which can be computed using the standard basis e 1 as Note that the multiplication of (C ⊗ C) and a vector can be computed using the 2D fast cosine transform with O (N log N). Also note that though D x and D y do not have the mirror symmetry, D T D = D T x D x + D T y D y does have the symmetry and thus where Γ is diagonal and can be computed in the same way as Eq. (12).
Equation (9) can then be written as Here Λ, Γ and H T g only need to be computed once with O(N log N). D is sparse and thus In total, Eq. (15) has computational complexity of O (N log N). This is a significant speed-up compared with O(N 3 ).
After the convergence of f , the barrier parameter is updated as with α>0, β>0. The iteration continues until f ≥ 0 is satisfied and d converges to Df . The deconvolution algorithm can be summarized in Alg. 1. The parameters used in this report are µ 0 = η 0 = 10 −2 , α = β = 0.5, ϵ 0 = 10 −3 , ϵ d = ϵ f = 10 −4 . Turning these parameters changes the rate of convergence, but does not have significant effects on the deconvolution results. The selection of λ is discussed in the next section.
Finally, we point out that the fast cosine transform is used here because of the Neumann boundary condition. If periodic boundary conditions or anti-reflective boundary conditions are used, one can use the fast Fourier transform and the fast sine transform, respectively. [50] All the discussions above are still valid.

Deconvulution results and discussion
In this section, we present our deconvolution results on CMOS samples. The pipeline of our work can be summarized as follows. First, we approximate the PSF using the fitting approach as mentioned in Section 3. Then we perform deconvolution on a training pattern using the method in Section 4 with varying regularization parameter λ. In this step, the optimal λ is determined by matching the deconvolution results with the ground truth. Finally, with this trained λ, we perform deconvolution on test patterns different from the training pattern and evaluate the results. Note that both the PSF and λ are determined by experiments independent from the test patterns. Over-fitting is avoided in this pipeline.
Our training pattern is a 2D poly-Si grid on sample A. The bright-field image of this pattern is similar with the one in Fig. 1(d). The grid has a periodicity of 930 nm and the bar width is approximately 325 nm (35% duty cycle). The Raman image was acquired by scanning the sample at a step size of 50 nm and the total pixel number is 61 × 61 (3 × 3 µm 2 field of view). Here the polarization is 45 • relative to x. In Fig. 3(a), we compare the first loadings and a representative raw spectrum normalized by its first score. Given the two curves are approximately identical, we conclude that the chemical concentrations and the Raman spectra of the sample can be approximated by the first score and the associated loadings, respectively. Some authors suggest the use of MCR-ALS in this step, [30] but in our work we notice that the results have no significant difference between using PCA and MCR-ALS. In Fig. 3(b), the raw image of the relative chemical concentration is present. The image is vectorized as g in Eq. (1) and is normalized by the mean valueḡ of all the pixels. This normalization is critical since the raw Raman counts from different experiments can vary significantly. In Fig. 3(c), we present the ground truth of the 2D grid. Here we choose the shift in x and y to approximately match Fig. 3(b). It can be observed that the raw image is blurred and the bar width is not resolved.
In Fig. 3(d)-(f), we present the deconvolution results of Fig. 3(b) using Eq. (1) with various λs. As expected, L1 regularization preserves edges and generates images with piece-wise flat regions. We observe that the deconvoluted images capture the shape and the periodicity of the grid. Each individual bar can also be resolved since the chemical concentration vanishes between adjacent bars. As λ increases, the bar shape becomes wider and shorter. This is expected since the penalty on the derivative becomes larger. This is clear in Figs. 3(g) and 3(h) where we plot the data slices on the white dashed lines along x and y respectively. In the Supplement 1, we also present the situations that the deconvolution fails when λ is too small or too large.
To quantify the error of the deconvolution, we compare the fill factor (FF) of the deconvolved image with that of the ground truth image. Here the FF is defined as the proportion of the poly-Si area over the total area. For a perfect 2D grid, the ground truth of FF is where p, w are the periodicity and the width, respectively. d = w p is the duty cycle of the grid bar. In the deconvoluted image, FF is approximated by counting the pixel number with the chemical concentration higher than the averaged value FF ≈ # of pixels with f >f total # of pixels (19) where f is the relative chemical concentration in Eq. (1) andf is the mean value of f . In our results, we observe thatf ≈ḡ = 1. This is because that g is normalized byḡ and convolution does not change the mean value. In Fig. 3(i), we present FF of the decovoluted image with varying λ. It can be observed that in general FF increases with λ. This is consistent with Figs. 3(g) and 3(h) in which the bar width increases with λ. Given FF GT ≈ 0.58, the optimal λ * ≈ 0.0048. The deconvoluted image corresponding to λ * is Fig. 3(e). Though the FF is correct, it is also worthwhile evaluating the error of the deconvolved feature size (bar width) and the poly-Si chemical concentration of the bar in Fig. 3(e). Here we approximate the bar width by counting pixels with f >f in the data slices in Figs. 3(g) and 3(h). The averaged bar width is 340 ± 42 nm, which is approximately 5% larger than the ground truth (325 mm). The poly-Si chemical concentration is approximated by the average of the f with f >f (poly-Si region), which is 1.63 ± 0.33. Here the ground truth is 1/FF GT ≈ 1.73. The means of both the bar width and the chemical concentration are relatively closed to the ground truth but the standard deviations are approximately 15% − 20%. The relatively large standard deviation is likely from the imperfection of the sample and is also a side effect of the edge enhancement effect. For example, In Fig. 3(b), g has different values near the centers of different bars. The difference can be amplified by the deconvolution and reflected in f . Also, as we pointed out in the PSF measurements, the polarized Raman intensity will be enhanced by the local geometry. Both factors will contribute to the variance from the ground truth in which we assume perfect geometry and neglect the local effects. Nevertheless, by choosing the optimal λ based on FF, the deconvolution still generates reasonably precise feature sizes and relative chemical concentrations. However, it is unlikely that one can determine the optimal λ a priori because the ground truth is usually unknown. In these situations, the optimal selection of the regularization parameters has been discussed by many authors. [51][52][53] In our work, we train λ on a known pattern and present that this optimal parameter can be used for the deconvolution on unknown patterns. Specifically, we use the pattern in Fig. 3 as the training pattern and set λ = 0.0048. The first test pattern is also a 2D poly-Si grid on sample A with smaller features: p = 760 nm, w = 230 nm (d = 0.3). The raw image, the ground truth, and the deconvolution with λ = 0.0048 are present in Fig. 4(a)-(c). We also present the data slices on the dashed white lines in Fig. 4(d) and (e). It can be observed that the bars are resolved in Fig. 4(c)-(e). The average bar width in Figs. 4(d) and 4(e) is 207 ± 35 nm (ground truth 230 nm). The averaged relative chemical concentration of the poly-Si region in Fig. 4(c) is 1.96 ± 0.43 (ground truth 1/FF = 1.96 ). All the values are closed to the ground truth (≈ 10% mean error). The relative large standard derivations, as discussed, are mainly due to the inherent imperfection of the measurements. We emphasize that here the deconvolution does not require information on the unknown test pattern.
The second test pattern is on sample B, which is a full-stack CMOS chip. The detailed difference of the two samples are mentioned in the PSF measurement section. Here the polarization is along the Si <100> of the substrate (x of the image) to minimize the background Raman signal. The laser power is approximately 6× higher than the previous experiments. The test pattern is a 1D grating with p ≈ 650 nm, w ≈ 275 nm (d ≈ 0.42). The optical image is present in Fig. 1(e). Similar to Fig. 4, the raw image, the ground truth image, the deconvolution with the trained λ = 0.0048, and the data slices are present in Fig. 5. Note that for a 1D grating, FF GT = d. We find that the averaged bar width is 291 ± 38 nm (ground truth 275 nm). Here the average value is from all the 61 y slices. The averaged poly-Si chemical concentration in Fig. 5(c) is 2.10 ± 0.31 (ground truth 2.37). FF = 0.044 (ground truth 0.042). The errors are still approximately 10%. This experiment shows that the deconvolution on Raman images can be applied to intact CMOS samples and that the training process can be performed on a separate sample with different experiment conditions.
We now consider the dependence of the deconvolution results on the PSF. Here we intentionally change the size of the PSF to introduce systematic errors. The Raman image of sample B is used as the example. We set the FWHM of the PSF as FWHM ′ x/y = k FWHM B x/y (20) where FWHM B x/y is the mean measured FWHM on sample B. In Figs. 6(a)-6(d), we present the deconvoluted images and the data slice with k = 0.8, 1.0, and 1.2 and λ = 0.0048. We observe that if the PSF is small compared to k = 1 the contrast is lower but the feature size (bar width) does not vary significantly. Meanwhile, when the PSF is larger than k = 1, the feature size shrinks and the contrast increases simultaneously. In Figs. 6(e) and 6(f), we can verify these observation by sweeping k and plotting the relative chemical concentrations and the bar width. This trend can be understood by considering that the mean of the image does not change after deconvolution. First, the PSF size determines the strength of the blurring effect. Specifically, when the PSF is small the blurring effect is weak, and thus the deconvoled image is close to the raw image. When the PSF is large, the deconvoled image has to be sparse and compact to match the strong blurring. Therefore, with the same mean value, the contrast increases with the PSF size. Second, the asymmetrical behavior is a result of the non-negativity constraint. In Fig. 6(e), when k<1, as k increases, the relative chemical concentration increases in the poly-Si area and decreases in the no-poly-Si area. The magnitudes are close and thus the bar width does not vary significantly. When k>1, the relative chemical concentration in the no-poly-Si area drops to zero and cannot decrease anymore. The bar width then starts to decrease to compensate for the increasing chemical concentration in the poly-Si area. Here we do not discuss FF because the features of the 1D grating are relatively simple and the FF versus k curve is almost identical to the bar width curve with a linear scaling factor.
In our experiment, the PSF measurement is relatievly precise with the standard deviation ≈ 6% on sample A and ≈ 4% on sample B. According to Figs. 6(e) and 6(f), the measurement error that propagates into the deconvolued image should be <10%, which is close to the error of the deconvolution itself. Meanwhile, the systematic error of the PSF (non-Gaussian focus, non-perfect poly-Si edge, etc.) is relatively difficult to calibrate. The systematic error may contribute to the deviation of the deconvolution with k = 1 from the ground truth in Figs. 6(e) and 6(f). We minimize this error by measuring the PSF using thin (≈ 100 nm) poly-Si edge-like patterns on the same chips of the test patterns.
Finally, we test the computational speed-up of Alg. 1 compared with other methods. We use the same data set and the parameters in Fig. 4. The raw image is truncated into multiple sizes to show the scaling of the running time with the problem size. Two other algorithms are used for comparison. In the first algorithm, Eq. (9) is solved using the Cholesky decomposition. The algorithm is similar to Alg. 1, with the Cholesky decomposition computed in each outer loop and the forward and backward substitutions computed in each inner loop. In the second algorithm, Eq. (3) is directly solved using the quadprog function in the Optimization Toolbox of MATLAB. The results are present in Fig. 7, where we also plot the acquisition time of the truncated images. It can be observed that Alg. 1 is 1 − 2 orders of magnitude faster than the Cholesky method and 2 − 3 orders of magnitude faster than quadprog. The gaps keep increasing as the problem size increases. Also Alg. 1 is more that two orders of magnitude faster than the acquisition time and thus the computational overhead is negligible. Meanwhile, the running time of the other two methods will take 5% − 100% of the acquisition time. These observations clearly show the superior performance of Alg. 1. It can be expected that the advantage of Alg. 1 would be even more significant with problems with larger sizes.   Fig. 4(a) with the total pixel number 30 × 30, 40 × 40, 50 × 50, and 60 × 60. The parameters are consistent with Fig. 4. The algorithms are implemented in MATLAB R2021a using a PC with Intel Xeon CPU E5-1620 v2 @ 3.70GHz and 80GB RAM.

Conclusion
In this work, we present our work of computational polarized Raman microscopy with subdiffraction-limit resolution on two CMOS chips. The polarization configuration is used to eliminate the background Raman signal from the substrates as well as to induce the edge enhancement effect to improve the image contrast. We first determine the PSF with the existence of edge enhancement. The hyperspectral image is then compressed into the image of chemical concentrations using PCA. With the compressed image, the deconvolution is formulated as a regularized optimization problem and is solved based on the fast cosine transform with O (N log N) complexity. The dependence of the decovolution results on the regularization parameter is investigated on a pattern with 325 nm features. We show that if the raw image is normalized by the mean value, one can train the regularization parameter on a known pattern and apply it to unknown patterns even if the physical properties and the experiment conditions are different. In this way, we can resolve 230 nm features on a back-etched CMOS sample and 275 nm features on an intact CMOS sample. These spatial resolutions are approximately 50% − 70% of the diffraction limit. The errors of the feature size and the relative chemical concentration are approximately 10%. Some portion of these errors might come from the imperfection of the sample and the side effect of edge enhancement. We also investigate how the size of the PSF affects the deconvolution results and confirm that the error propagation from the PSF measurement is less than 10%. At the end, we benchmark our algorithm and show significant speed-up compared to two other methods.
It is worthwhile to point out that we use the prior information of the feature orientations when choosing the polarization direction. If this information is unknown or the local feature orientation is disordered, one can perform two scans with orthogonal polarization and the two images should have complementary contrast enhancement. (For example, see Fig. S1 in the Supplement 1.) If the acquisition time is not a constraint, one can generate images with multiple polarization angles from 0 to 90 • . This set of images can be averaged to achieve enhancement of features in different directions. If the scale of the disordered sample is much smaller than the PSF, the contrast enhancement may be difficult to observe since the signal is averaged over the PSF. One may still be able to retrieve the map of the local feature orientation by deconvolution. However, in this situation, the local feature orientation is a property to resolve rather than a property that can be utilized.
In conclusion, we set up a complete pipeline for the deconvolution with polarized Raman microscopy. This pipeline includes the PSF determination with the existence of edge enhancement, independent parameter training, and result testing. The performance of this approach is evaluated by discussing the parameter selection and quantifying the errors of various properties. In the future, we expect researchers to use our method on more complicated samples, for example, compounds with multiple chemicals. Since our approach has minimal limitation on the instruments, it is promising to combine it with other techniques such as ISM, SRS and CARS. We also hope that the discussion on the size of the PSF can contribute to the development of algorithms which require even less prior information, such as blind deconvolution. [53] Funding. National Research Foundation Singapore (NRF-CRP16-2015-04 to N.I.N); Singapore-MIT Alliance for Research and Technology Centre.