1. Introduction

Near-field optics circumvents the diffraction limit by introducing a sharp probe tip to the proximity of the investigated sample specimen [1]. The high-momentum photons in the near-field zone is the necessary ingredient to gain deep subwavelength spatial resolution. In the modern scattering-type scanning near-field optical microscopy (s-SNOM) [2–4], the probe tip is polarized by the incident light, inducing a strongly enhanced field intensity at the apex of the tip, which interacts with the sample underneath. This interaction further modifies the amplitude and phase of the field scattered by the tip. The detected scattered field in the far-field zone can then be used to infer complex material properties with deep subwavelength spatial resolution.

In any imaging or spectroscopic experiment, establishing a concrete understanding of the signal contrast mechanism is the cornerstone for quantitative information extraction from the raw data. This is deemed a challenging task in s-SNOM due to the nontrivial electrodynamic interaction between the light, the nano-sized tip, and the sample material in the near-field zone [5]. Because of the difficulty in modeling the irregular shape of the tip, idealization is often inevitable. Based on this idea, numerous successful models with different complexities have emerged over the past two decades as solutions to explain and analyze the s-SNOM contrast. Prominent examples include the point-dipole model [6,7], the finite-dipole model [8], the lightning-rod model [9], and a few others [10,11]. Furthermore, extensions of the models to account for more general experimental conditions such as layered thin films [12–16] and anisotropic materials [17–20] have been proposed.

Recently we proposed a numerical framework that extends the original point-dipole model to surfaces with arbitrary topographic features and dielectric compositions [21]. Although the method demonstrates great merits and is time-efficient, it suffers from several drawbacks. For example, the point-dipole model is known to severely underestimate the penetration depth of the electric field due to the rapid decay of the dipole field [22]. This results in discrepancy between the predicted spectra and the experimental measurements on thin films [22]. Secondly, the point-dipole model does not represent the realistic tip shape. When strong resonances occur in the sample material, the predicted spectral response often significantly deviates from the experimental measurement [8]. Therefore, the previous method welcomes improvements that can address such issues. In this article, we demonstrate an alternative to our previously reported method that is more advanced and efficient in explaining and modeling s-SNOM data. While exhibiting similar time efficiency with the initial method, the new version considers more realistic tip geometries which enables a wide palette of additional use cases, augmenting as well the applications that we previously discussed in detail in our past work [21]. As a proof-of-concept demonstration, we employ an elongated spheroid tip [10,23] to investigate various near-field phenomena such as the near-field interaction with biaxial and uniaxial materials, the dependence of the spatial resolution on various experimental parameters, and several other important issues that preoccupy the s-SNOM community.

2. Quasi-electrostatic simulations including the extended tip geometry

Electrodynamic simulations have led to fruitful results in the context of s-SNOM signal modeling [24–31]. However, the computationally intensive nature of full-wave simulation prevents its wide use, with a limited number of applications having been reported. Simulating a near-field spectrum or image with a high spectral or spatial resolution is often impractical for commercial computers. On the other hand, it should be noticed that most analytical models are constructed under the quasi-electrostatic limit, where the incident wavelength is considered to be much larger than the tip and the probed region of interest on the sample. Although this approximation does not seem appropriate, especially in the visible and near-infrared (IR) regimes as the tip length is typically in the order of ten microns, it has been nonetheless shown that models under the quasi-electrostatic approximation capture the near-field physics very well especially in mid- and far-IR [6,8–10,32]. The prominent advantage of the quasi-electrostatic approximation is that it makes numerical simulations stable and time-efficient by only solving $- \nabla \cdot \nabla V = {\rho _f}$ with the appropriate boundary conditions [21]. Compared to its full-wave simulation counterpart, we find that our quasi-electrostatic simulation is at least one to two orders of magnitude faster. In addition, due to the absorption boundaries and small tip-sample distance used, full-wave simulations can often encounter convergence issue, which is not present in quasi-electrostatic simulations. We exploit these valuable features to set up a new s-SNOM simulation method that is described next in detail.

Various commercial and open-sourced Maxwell equations solvers such as COMSOL Multiphysics provide convenient platforms to conduct low-frequency electromagnetic simulations. To simulate the tip-sample near-field interaction in a quasi-static limit, the first step is to construct a uniform background excitation field ${E_i}$, which polarizes the tip and induces the field enhancement at the tip apex. As p-polarized light is predominantly used in the experiment and the elongated tip geometry yields a much larger polarizability in the tip shank direction, we only consider the vertical field and its dipole moment. The horizontal components can be considered in a similar fashion but are not explicitly included in our current simulations. As schematically shown in Fig. 1(a), the uniform excitation field ${E_i}$ can be generated by setting the boundary conditions for the top and bottom surfaces as constant voltage $V = {V_0}$ and $V = 0$, similar to a parallel plate capacitor connected to a fixed voltage source. The tip and sample can be subsequently added into the simulation space. Commonly used tip geometries in the analytical models include elongated spheroid [8,10], hyperboloid [9], cone [11], and so on. We adopt a metallic elongated spheroid tip as an example throughout this work, but other tip configurations can also be used depending on the application of interest (See Supplement 1). The length of the spheroid is set as $L = 600\; \textrm{nm}$ in this study as we find this value to yield a good fit to the experimental data. In practice, L can be used as an empirical fitting parameter to capture the characteristics of the tip used in the experiments. The apex radius of curvature a relates to $L\; $ and the width of the spheroid w by $a = \frac{{{w^2}}}{{2L}}$. By default, we use $a = 27\; \textrm{nm}$, a dimension that aligns well to the size of common commercially available tips used for s-SNOM imaging. The appropriate mesh setting can be determined by checking the convergence while iteratively increasing the mesh density. We find that the convergence is ensured by setting the mesh size near the spheroid apex comparable to the tip-sample distance $dz$.

 

Fig. 1. (a) Schematics of the simulation setup. (b) The dipole moment of the spheroid on a Si sample as a function of tip-sample distance. (c) The corresponding approach curves demodulated at second, third, and fourth harmonics. Higher demodulation order exhibits a faster decay rate. (d) Surface charge density along the surface of the spheroid when it is 10 nm and 500 nm away from the Si surface. The inset shows the surface charge distributions. (e) Simulated near-field spectra of SiO₂ thin films of different thicknesses normalized to Si. (f) The amplitude of the dipole moment as a function of frequency and tip-sample distance for a 300 nm thick SiO₂ film on Si substrate.

Download Full Size | PDF

Once the simulation is complete, we obtain the potential $V({x,y,z} )$ and subsequently the electric field ${\boldsymbol E}({x,y,z} )$. The dipole moment of the spheroid in the vertical direction can then be calculated by integrating the surface charge density $\sigma $ over the spheroid surface

Before performing elaborate simulations, we verify the validity of the method by means of simple examples. First, we simulate ${p_z}$ as a function of the tip-sample distance $dz$ on Si as shown in Fig. 1(b). The corresponding demodulated signals ${S_n}$ (approach curves) are shown in Fig. 1(c). Both ${p_z}$ and ${S_n}$ exhibit rapid decay as the tip pulls away from the sample surface as expected. A higher demodulation order yields a faster decay rate, which is a well-known phenomenon [36–39]. Here we take a look at the surface charge distribution in the spheroid in two cases – when it is interacting, or not, with the Si surface. The charge density $\sigma $ on the spheroid surface along the vertical direction when $dz = 10\; \textrm{nm}$ and $dz = 500\; \textrm{nm}$ are shown in Fig. 1(d). When the spheroid is far away from the Si surface, the charge distribution is purely induced by the background field ${E_i}$. Therefore, the positive and negative charges exhibit a symmetric distribution toward both ends of the spheroid. When the spheroid is brought in the vicinity of the Si surface, the near-field interaction breaks the symmetry and causes a drastically higher charge density on the apex, which is responsible for the strong field enhancement as well as for the modified dipole moment.

Next, we demonstrate the advantage of the proposed method compared to the previous point-dipole model framework [21]. In Fig. 1(e) we plot the broadband near-field amplitude spectra of SiO₂ thin films deposited on a Si substrate. The spectra are normalized to that on Si. The thickness t of the film varies from 300 nm to 2 nm. The simulated spectra show great consistency with the previously published experimental data (see also Supplement 1) while the point-dipole model prediction is noticeably worse [22]. This is because the point-dipole does not capture the effect of the extended tip shank very well while the elongated spheroid is a much superior approximation to the real conic or pyramidal tip. Since the measurement of thin films is ubiquitous in s-SNOM experiments [40–45], an accurate model for data analysis in this regard is indispensable. For reference, the simulation of one spectrum in Fig. 1(e) only took a few minutes using a common commercial desktop.

Besides s-SNOM, the proposed simulation method can be readily adapted to simulate other optical scanning-probe measurements. For example, in Fig. 1(f) we show the amplitude of ${p_z}$ as a function of $dz$ and frequency $\omega $ for a 300 nm SiO₂ film on a Si substrate. This result is representative of the recently developed peak-force SNOM technique, where the vertical interaction can be reconstructed without the Lock-in demodulation of the scattering signal [38,46,47]. Other optical scanning probe techniques such as photo-induced force microscopy (PiFM) [48,49], photo-thermal-based techniques like AFM-IR [50–52], tip-enhanced photo-luminescence [53], and tip-enhanced techniques based on non-linear optical effects [54–56] could potentially benefit from similar simulations as well. For example, to simulate PiFM signal, one can evaluate the electric force on the tip due to the charge accumulation on the sample. To simulate AFM-IR contrast, the electromagnetic simulation needs to be coupled with heat transport simulation to determine the thermal expansion of the sample under near-field heating.

It should be noted that in the simulation, the background electric field ${E_i}$ is not a constant when ${V_0}$ is set as a constant. This is due to the change of the sample's dielectric function or geometry. To illustrate this point let us consider an ideal parallel plate capacitor with a dielectric slab inside. The dielectric slab represents the sample in the simulation. Assume the spacing between the plates is ${d_1} + {d_2}$. The thickness of the slab is ${d_2}$. The relative dielectric function of the slab is $\varepsilon $. Continuity of ${D_ \bot }$ across the dielectric boundary leads to ${\varepsilon _0}{E_i}$ $= \varepsilon {E_d}$, where ${E_d}$ is the field inside the dielectrics and ${\varepsilon _0} = 1$ is the relative vacuum dielectric constant. Also, we have $V = {E_i}\; {d_1} + {E_d}\; {d_2}$. Combining the two equations gives

3. Near-field response of uniaxial and biaxial anisotropic materials

Most crystals naturally exhibit anisotropy due to the lattice symmetry. However, the incorporation of sample anisotropy in some models is difficult and mathematically cumbersome. This is the reason why the impact of sample anisotropy on the near-field scattering signal is still not thoroughly studied. One great advantage of our proposed method is that the full dielectric tensor can be straightforwardly included in the simulation. Next, we investigate the near-field response of uniaxial and biaxial anisotropic materials. An earlier experimental study investigated the dielectric response of mica flakes under 40 nm thin using broadband near-field spectroscopy (nano-FTIR). Using a multilayer finite-dipole model, the authors found that only the out-of-plane dielectric function contributes to the near-field spectra while the in-plane response in mica is not detectable. We attempt to test this finding using the spheroid tip in our simulation. The in-plane and out-of-plane dielectric functions of mica are directly obtained from the previous study [57] and input into the COMSOL simulations in the full tensor form. The simulated ${S_2}$ amplitude spectra for 40 nm, 7 nm, and 4 nm mica on Si substrate are consistent with the experimental data in [57] and shown in Fig. 2(a). Next, we closely examine the in-plane and out-of-plane responses individually. In Fig. 2(b) we show the calculated spectra of 40 nm mica on Si substrate. The three spectra are obtained using the full dielectric tensor, only the in-plane dielectric function, and only the out-of-plane dielectric function. Our result indicates that although only considering the out-of-plane component of the dielectric tensor yields a similar result compared to considering the full tensor, the effect of the in-plane response is non-negligible at this sample thickness. For thinner films such as 4 nm mica on Si substrate, we find that indeed the out-of-plane response is dominant as shown in Fig. 2(c).

 

Fig. 2. (a) Simulated ${S_2}\; $ amplitude spectra of 40 nm, 7 nm, and 4 nm mica on slab Si substrate. (b) Simulated ${S_2}\; $ amplitude spectra of 40 nm mica on Si substrate considering the full dielectric tensor, only the out-of-plane dielectric function, and only the in-plane dielectric function. (c) Similar to (b) but for 4 nm mica. (d) In-plane and out-of-plane dielectric functions of hBN, showing distinctive phonon resonances. (e) ${S_3}$ amplitude spectra simulated for 300 nm thick hBN on a Si substrate considering the full dielectric tensor (solid curve) and the isotropic effective dielectric function (dashed curve). (f) ${S_3}$ amplitude spectra simulated for 5 nm thick hBN on a Si substrate (black curve) and 5 nm suspended hBN (orange curve).

Download Full Size | PDF

This finding for mica is seemingly in conflict with the previous reports on graphene, where the in-plane response was found to be the dominant factor to the near-field signal contrast [15]. The noticeable difference between a thin layer of mica and graphene is that the optical conductivity of graphene, especially for doped variants, is orders of magnitude higher. Next we investigate how in-plane and out-of-plane responses contribute to the near-field signal using a prototypical material hexagonal boron nitride (hBN), where the in-plane and out-of-plane phonon resonances occur at distinctive frequencies, unlike mica or other anisotropic crystals such as sapphire. Figure 2(d) shows the real part of the dielectric functions for the in-plane and out-of-plane directions of hBN [58], where the in-plane resonance occurs at a higher frequency (∼1350 cm⁻¹) while the out-of-plane resonance at a lower frequency (∼750cm⁻¹) with comparable oscillator strength. It has been shown that an effective dielectric function ${\varepsilon _{eff}} = \sqrt {{\varepsilon _{in}}{\varepsilon _{out}}} $, which is the geometric mean of the in-plane and out-of-plane dielectric functions, provides a good approximation for the bulk sample [8,10]. This is the direct result of the method of image charge [59]. First we study the thick film case and simulate a 300 nm thick hBN on a Si substrate using the anisotropic dielectric tensor and compare that to the result using the effective dielectric function. The calculated ${S_3}$ amplitude spectra are shown in Fig. 2(e). Both the in-plane and the out-of-plane resonances show signatures in the near-field spectra. It is evident that considering the effective dielectric function indeed leads to a practically identical prediction compared to considering the full dielectric tensor. That is, for bulk samples or thick films, s-SNOM probing does not have a directional preference and the contributions from in-plane and out-of-plane dielectric responses are equal.

In the case of thin films, the effective dielectric function ${\varepsilon _{eff}} = \sqrt {{\varepsilon _{in}}{\varepsilon _{out}}} $ breaks down and the situation is more complicated. Not only the strengths of in-plane and out-of-plane responses of the film are important, but also the dielectric property of the substrate plays an important role as it can anisotropically screen the film response. This can be seen in Fig. 2(f), where ${S_3}$ amplitude spectra are simulated for 5 nm hBN on a Si substrate as well as for 5 nm suspended hBN. In the hBN/Si case, the in-plane resonance is suppressed while in the hBN/air case the out-of-plane resonance is diminished. This behavior has a qualitative interpretation in terms of substrate screening. The tip-scattered field is related to the total dipole moment on the thin film and the substrate. On a high dielectric substrate with isotropic dielectric function ${\varepsilon _s}$, the dipoles are screened in such a way that the horizontal dipoles are weakened by a factor of ${f_h} = \frac{2}{{1 + {\varepsilon _s}}}$ while the vertical dipoles are enhanced by a factor of ${f_v} = \frac{{2{\varepsilon _s}}}{{1 + {\varepsilon _s}}}$, as can be shown by the method of image charge. In the limit of $|{{\varepsilon_s}} |\gg 1$, one has ${f_h} = 0$ and ${f_v} = 2$. This is further verified by the form of the frequency- and momentum-dependent reflection coefficient of the structure [10,60,61]:

Furthermore, in models such as the point-dipole model [15,61] and the lightning-rod model [9], uniaxial materials with in-plane rotational symmetry can be incorporated by considering ${r_p}({\omega ,q} )$. However, with biaxial materials, the in-plane rotation symmetry is broken, which poses difficulty in the modeling efforts [19]. This is where our proposed method demonstrates superiority. As a concrete example, we consider the near-field response from the biaxial van der Waals material MoO₃. MoO₃ is known to host unique anisotropic surface phonon polaritons in the mid-IR and THz regimes, which have attracted tremendous research interest [23,62,63]. The dielectric function of MoO₃ in the relevant frequency range can be obtained from [64]. We first consider a 100 nm MoO₃ slab on a Si substrate and simulate the electric field distribution at various characteristic frequencies. As shown in Fig. 3(a), unidirectional polariton propagation in the crystal [100] and [001] directions are observed at 740 cm⁻¹ and 915 cm⁻¹, while elliptical polariton waveform emerges at 985 cm⁻¹. The biaxial anisotropy is responsible for this unique signature of MoO₃. Next we demonstrate the simulation of the near-field spectrum of MoO₃. 100 nm, 30 nm, and 10 nm MoO₃ slabs on a Si substrate are simulated. The amplitude and phase spectra are shown in Fig. 3(b) and (c). The ability to simulate subtle spectral responses of biaxial materials is the key to the extraction of quantitative information from near-field measurements. Besides spectral domain, the near-field response in real space is another interesting aspect especially when polariton occurs. Here we demonstrate that our method is suitable to study real-space phenomena such as the edge-launched polaritons. We simulate a line scan in the [100] direction across the edge of a 250 nm MoO₃ on a Si substrate at 930 cm⁻¹. ${p_z}$ as a function of distance to the edge X and $dz$ is shown in Fig. 3(d). The demodulated near-field response $|{{S_2}} |$ is shown as the black curve, where clear periodic oscillation is observed. Compared to the calculated phonon polariton dispersion (inset of Fig. 3(d)), we find that the periodicity exactly matches the phonon polariton wavelength, which indicates that the observed oscillation occurs due to the edge-launched polariton instead of the tip-launched one. This is due to the fact that the tip radius (27 nm) is about one order-of-magnitude smaller than the polariton wavelength thus the momentum matching is inefficient.

 

Fig. 3. (a) Electric field in the z-direction ${E_z}$ simulated at 740 cm⁻¹, 915 cm⁻¹, and 985 cm⁻¹, where different polariton propagation directions due to the biaxial anisotropy are observed. (b) and (c) Simulated ${S_3}$ amplitude and phase spectra for 100 nm, 30 nm, and 10 nm MoO₃ slabs on Si substrate. (d) Simulated p_z as a function of dz during a line scan across the edge of a 250 nm thick MoO₃ slab on Si substrate at 930 cm⁻¹. The black curve shows the corresponding amplitude of ${S_2}$, where oscillation with ${\lambda _p}$ period corresponding to the edge-launched phonon polariton is observed. The inset depicts the dispersion of the phonon polariton calculated by the imaginary part of the reflection coefficient $Im({{r_p}} )$. The star marks the observed frequency and momentum in the simulation. (e) Simulated ${S_3}$ amplitude spectra on two 10 nm MoO₃ slabs with a relative twist angle of $0^\circ ,\; 30^\circ ,\; 45^\circ ,\; 60^\circ ,\; $ and 90$^\circ $.

Download Full Size | PDF

In addition, topological transitions of polariton dispersion were recently observed in twisted MoO₃ slabs [65–67], making these a unique platform for investigating interesting physics and applications. This complex system again poses tremendous difficulties with respect to conventional modeling endeavors while our method can characterize the twisted layer in a facile manner. Next we demonstrate the spectra as a function of the twist angle $\theta $ in Fig. 3(e). The twisted layer can be characterized by a non-diagonal dielectric tensor ${\varepsilon _{twist}} = {R^T}{\varepsilon _{Mo{O_3}}}R$, where ${\varepsilon _{Mo{O_3}}}$ is the diagonal dielectric tensor of MoO₃, R is the rotation matrix, and ${R^T}$ is the transpose of R. Small yet systematic changes are observed for $\theta $ from 0$^\circ $ to 90$^\circ $. Note that this effect is intrinsically near-field, not to confuse with the far-field factor. The possibility to accurately simulate the spectral details is of paramount importance in delicate systems like this and other 2D material heterostructures. Our method provides a quantitative and versatile modeling approach for characterizing the spectrum of arbitrary materials, which is an important step for the quantitative analysis of near-field data.

It is worth noting again that in our simulation, the multiplicative far-field factor is not explicitly included. The spectral feature is purely due to the tip-sample near-field interaction. When an anisotropic sample is present, the anisotropy can also manifest itself via the far-field reflection off the sample surface, where the angle of incidence could be of significant importance and should be carefully addressed.

4. Experimental factors for the lateral resolution

Besides the flexibility it provides for incorporating different material properties, another major advantage of our numerical approach is the ability to consider the topographic and compositional inhomogeneities and simulate near-field images and line scans, a task usually very difficult with conventional analytical models. We demonstrate the utility of our method in this regard through the investigation of the lateral resolution in s-SNOM. A previous study employing both experimental works and full-wave simulations shows that, contrary to common beliefs, the resolution of THz s-SNOM does not simply scale with the tip apex radius linearly [27]. Instead, blunter tips yield a much larger scattering signal while still preserving a reasonable resolution. This makes blunter tips a better option than sharper tips considering the fact that THz light sources typically emit at low power. Here we shed additional light on this problem by employing our proposed methodology.

First we simulate the near-field signal line profile as the tip scans over a metal (Au)-dielectric (Si) interface. To avoid topography-related complications, the interface is designed to be flat as schematically shown in the inset of Fig. 4(a). Simulated ${S_n}$ amplitude profiles for $n = 1,2, \ldots 8$ are shown as the dots in Fig. 4(a). The profiles are fitted by the arctan function proposed in [27]. The best fits are shown as solid curves in Fig. 4(a) and the fitting parameters provide a quantitative estimate for the lateral resolution [27]. Figure 4(b) demonstrates the lateral resolution determined from fitting for different demodulation orders using different tip radii ranging from 10 nm to 200 nm. The tapping amplitude is kept constant at 50 nm in all simulations. The general trend points to the fact that higher demodulation leads to a sharper boundary profile. This is because the lateral extension of the field modulated at higher harmonics is increasingly more confined [24]. More importantly, comparing 200 nm tip radius to 10 nm tip radius, we can see that although the tip radius is 20 times larger, the resolution is only worsened by a factor of less than 3. Importantly, the relative signal contrast as well as the raw scattering intensity are observed to increase along with the tip radius (Fig. 4(c)). These observations are highly consistent with the findings of the previous reports involving both experimental work and sophisticated full-wave simulations [24,27,28], which further validates the efficiency and accuracy of our simulation method based on quasi-electrostatics.

 

Fig. 4. (a) Simulated line profiles of ${S_n}$ amplitude as the spheroid scans over the Au and Si interface. The tip radius is 27 nm. (b) The lateral resolution obtained from fitting for different tip radii and demodulation order. (c) Relative contrast and absolute signal intensity on Au for ${S_3}$ amplitude as a function of tip radius. (d) The lateral resolution obtained from fitting for different material boundaries. The results include Au-Si ($\varepsilon = 12$) boundary at 1000 cm⁻¹, Au-PMMA ($\varepsilon = 2$) boundary at 1500 cm⁻¹, and Au-SiO₂ ($\varepsilon ={-} 3.25 + 1.5i$) boundary at 1130 cm⁻¹. (e) Lateral resolution obtained from fitting for different tip oscillation amplitudes. (f) Relative contrast and absolute signal intensity on Au for ${S_3}$ amplitude as a function of the tip oscillation amplitude.

Download Full Size | PDF

The time-efficient nature of the simulations enables us to systematically study other contributing factors that account for the attainable resolution. Next we simulate the resolution dependence on material properties. More specifically, we study two additional interfaces: an Au-PMMA interface at 1500 cm⁻¹, where PMMA has a small dielectric constant (${\varepsilon _{PMMA}} = 2$) and an Au-SiO₂ interface at 1130 cm⁻¹, where the phonon polariton resonance occurs (${\varepsilon _{Si{O_2}}} ={-} 3.25 + 1.5i$). The determined resolutions are plotted in Fig. 4(d) along with that from the previous Au-Si interface (${\varepsilon _{Si}} = 12)$. Both tip radius (27 nm) and tapping amplitude (50 nm) are kept constant. Comparing the resolutions obtained for Au-Si and Au-PMMA, we can see that the resolution is slightly higher in the PMMA case due to the larger dielectric contrast. However, the resolution is drastically worse in the Au-SiO₂ case. We attribute the worsened resolution to the long-wavelength lossy phonon polariton in SiO₂, similar to the blurry boundary recently observed in WTe₂ due to the plasmon polaritons [68]. The oscillation amplitude of the tip is another important factor that determines the achievable resolution. In Fig. 4(e) we study the resolutions obtained using different oscillation amplitudes. Smaller oscillation amplitude is found to improve resolution. However, when the oscillation amplitude is small, demodulation of the raw scattering signal to high harmonics of the tip oscillation frequency results in noisy data, which reduces the reliability of the fitting. Therefore, we only determine the resolution values up to the 7th harmonic for 25 nm tapping amplitude and the 3rd harmonic for 10 nm tapping amplitude. As shown in Fig. 4(f), although the relative contrast is higher when a smaller tapping amplitude is used, the raw signal intensity is lower.

To summarize our findings, higher demodulation order, smaller tip radius, larger contrast in dielectric function, and small tapping amplitude generally result in higher lateral resolution. However, in practice, higher resolution is usually not the only consideration. Raw signal intensity and signal contrast are often equally important because of the finite signal-to-noise ratio limited by the light source power, stability, and detector noise. Therefore, there is a fine balance in choosing the ideal experimental parameters. Identifying this optimal balance experimentally is cumbersome due to the typical signal integration times used in s-SNOM, and hence the image acquisition duration. Thus, our simulation method can be a valuable tool to help identify optimal experimental parameters for a particular s-SNOM application of interest.

Further on, we notice that the resolution and contrast determined from the flat material interface might not reflect other sample situations well. We prove this point by investigating a case where a small Au disk is embedded in a PMMA host, as schematically shown in the inset of Fig. 5(a). A simulated near-field ${S_3}$ amplitude image is shown in Fig. 4(e) using a tip radius of 27 nm and tapping amplitude of 50 nm. The line profile across the red dashed line is plotted in Fig. 5(b) together with several other profiles obtained using different tip radii. The trend is obviously very different from the flat interface case discussed previously, which is evident by the resolution and the relative contrast shown in Fig. 5(c). In the embedded small disk case, the resolution is observed to scale mostly linearly to the tip radius while the relative contrast drastically decreases when a large tip radius is used. This observation makes the larger tip radius highly unfavorable when imaging samples with small feature sizes. Based on the simulation results, we conclude that the ideal tip radius for a specific experiment depends on the characteristic size of the feature in the sample. A tip radius comparable in size to the feature of interest is the most optimal if both the resolution and the relative contrast are considered.

 

Fig. 5. (a) Simulated ${S_3}$ amplitude image for a 20 nm radius Au disk embedded in PMMA. (b) Simulated line profiles of ${S_3}$ amplitude as the spheroid scans over the red dashed line in (a) for different tip radii. The vertical dashed line indicates the edge of the embedded disk. (c) Lateral resolution and relative contrast as a function of tip radius.

Download Full Size | PDF

5. Conclusion and discussions

In this article, we propose a numerical method that is highly capable of accurately simulating the tip-sample interaction and tip scattering in s-SNOM experiments with tip geometry on demand. This method significantly augments the previously proposed simulation framework based on the point-dipole model [21]. The ability to consider the elongate tip shank provides a realistic prediction of the penetration depth of the electric field, which is important for analyzing spectra or images taken on layered thin films.

We use the proposed method to study anisotropic near-field interaction and find that for thick samples ($t$>300 nm), s-SNOM is practically probing the geometric mean of the in-plane and out-of-plane dielectric functions. For thin samples ($t$ ∼10 nm), however, the near-field response is a complex mixture of in-plane and out-of-plane responses. Furthermore, our flexible method is able to account for biaxial anisotropy and even the presence of non-diagonal dielectric tensors, which is considered difficult in conventional analytical models.

The relations between the lateral resolution, signal contrast, raw signal intensity, and various external experimental parameters such as tip radius, demodulation orders, tapping amplitude are studied in greater depth. When the sample feature is much larger than the tip radius, a larger tip radius leads to stronger scattering intensity and larger material signal contrast while preserving a good lateral resolution. However, when the sample feature is small, e.g. nanoparticles, a larger tip radius becomes disadvantageous. Our simulations show that in this case the optimal tip radius is the one that is comparable to the feature size.

Due to the quasi-electrostatic nature of the simulation, retardation effects are ignored in our simulations. In certain situations these effects could play an important role [9,30], where full-wave simulations should be employed instead. A prominent example is when the incident light excites antenna resonances of the tip [69]. In this case, surface current forms standing wave, and a stronger field enhancement under the tip apex is realized. However, in most cases, our proposed method is demonstrated to be highly accurate, versatile, and easy to implement. It features several unique improvements over the previous methods and models. Therefore, we envision that it will become an important tool for future s-SNOM data analysis such as dielectric function extraction from spectroscopic or imaging measurements [14,70–72]. Given that the proposed method is highly time-efficient, it can also be regarded as a valuable tool to support future machine learning related research [73–75] focused on s-SNOM and other optical techniques based on tip-enhancement, where a fast and reliable way to generate training data matching the experimental conditions of interest is crucial.