Rigorous modeling and physical interpretation of terahertz near-field imaging

Apertureless scanning near-field optical microscopy (SNOM) operating with terahertz (THz) laser pulses is a subject of great research interest. The Mie scattering theory is commonly used to explain ...


INTRODUCTION
Terahertz (THz) imaging attracts much attention for modern applications, for instance, in security systems [1], packaging inspection [2], and biochemical engineering [3,4]. THz waves are electromagnetic radiation in the frequency band about from 0.1 to 10 THz, corresponding to wavelengths in the submillimeter range. According to the Rayleigh criterion, the conventional spatial resolution of THz imaging microscopy cannot be much smaller than 30 µm (corresponding to about 10 THz), which is a big problem when the size of the sample shrinks to (sub)micron scales. An idea to overcome the diffraction limit and improve the resolution in THz imaging is to use near-field scanning microscopy (SNOM). The SNOM combines the approach of collecting the optical field in the near zone (the near-field) of the sample surface with a scanning probe microscopy technique. In the apertureless scattering SNOM (or s-SNOM), a metallic probe that interacts with the incident radiation and enhances the field near the probe tip is brought to a sub-wavelength proximity of the sample, giving the s-SNOMs a high resolution performance [5,6].
The experiments recently reported by [7] demonstrated that the resolution for an s-SNOM imaging system at 2 THz can go down to 150 nm, which is much smaller than the corresponding wavelength of the incident signal (about 150 µm). These authors also point out [7]- [9] that the experimental observations are not consistent with the predictions of the Mie scattering theory, which traditionally has been used to interpret the optical data obtained with s-SNOMs. A simple circuit model based on antenna theory was found to agree with the THz scattering measurements [7,8]. In our paper, we describe a rigorous numerical model for this s-SNOM THz imaging system using the finite element method (FEM). Our approach, though only two-dimensional, deals with the electromagnetic interactions on the basis of Maxwell's equations and therefore overcomes the problems associated with the Mie scattering method [9]. Our results are consistent with the three-dimensional experimental observations reported previously [8]. We also show that an intuitive, slightly modified antenna radiation model is fully capable of explaining the explicit THz scattering results.

PHYSICAL MODEL
In the experiment reported by Chen and coworkers [8], a tungsten probe is brought into the proximity of a gold surface evaporated onto a silicon substrate. A p-polarized pulse with a central frequency of 2 THz and an incidence angle of 60 • is used to illuminate the probe tip-surface system. A p-polarized field couples more efficiently into the s-SNOM system than an s-polarized one, since the electric field vector of the s-polarized signal is perpendicular to the direction of the probe and therefore has a limited capability of exciting plasmons on the tip surface. The specular reflection of the incident signals is detected in the far zone of the tip-sample structure Figure 1(a). The altitude of the probe tip is changed during this experiment [8]. The measurements at different tip heights showed that for separations between the tip and the gold surface smaller than about 1.5 µm the specular reflection of the incident pulse grows as the tip altitude decreases; this reflection also increases as the tip altitude becomes higher beyond 1.5 µm [8]. However, such a dependence of the detected THz signal cannot be explained using the Mie scattering model [7,9].
To elucidate such a behavior, a detailed physical model involving the electric field-material interactions is required. The main material property to be specified in this study is the electric permittivity for the 2 THz radiation band. The relative permittivity of gold (Au) is obtained from the measurements of the normalized surface resistance of gold samples at the ambient temperature with a non-resonant cavity [10]. The relative permittivity of tungsten (W) is calculated by using the Lorentz-Drude (LD) model [11]. The LD model, which gives a better approximation of the permittivity than the more commonly used Drude model, considers not only the interactions of free electrons (intra-band transitions) but also takes into account the effects of electrons that are promoted from lowerlying bands into the conduction band (inter-band transitions). The relative electric permittivity constants at 2 THz obtained in this way are: tungsten ε r = −8.6 × 10 3 + i 6.7 × 10 4 , and silicon ε r = 11.6 + i 2.5 × 10 −3 . These values are employed in the numerical simulations below.

MATHEMATICAL MODEL
A real s-SNOM scattering experiment takes place in three dimensions (3D), and it is clear that a two-dimensional (2D) model cannot capture all the phenomena, such as the effects of the radiation passing by the side of the probe tip. However, as demonstrated in the following sections, the 2D model incorporates most of the salient features of the s-SNOM system and shows good physical agreement with the experimental data [8]. It is also considerably simpler and demands much less computer resources compared to a full-scale 3D model [12]. For the numerical modeling in a 2D geometry, we apply the finite element method (FEM) implemented in the RF module of the COMSOL Multiphysics simulation package [13].
In the 2D model, the calculation area is separated into two regions: the near-field zone and far-field zone Figure 1(a). The near-field zone is the core of the s-SNOM imaging system that includes both the tungsten tip and the sample surface. The size of this zone is about 1.4 mm × 1.4 mm, which holds about 10 × 10 full wavelengths of the 2 THz radiation band (λ ≈ 150 µm). The detector is placed in the far-field zone, about 2 m away from the probe tip. The specular reflection registered by the detector is calculated according to the Huygens-Fresnel principle and the data for the field along the integration line in the near-field zone Figure 1(b). The amplitude and phase of the scattered electric wave on this inte- gration line are obtained as the exact solution of the Maxwell equations in the near-field zone Figure 1(b). This area is surrounded by perfectly matched layers (PMLs) [14], which are employed to simulate open boundaries, i.e., no boundaries at all, to avoid backscattering.
The tungsten probe in the simulation has a height of 1.2 mm; this corresponds to the part of the probe that is illuminated by the incident plane wave signal. The probe tip has a conical shape with a vertex angle of 21 • and the end of the tip consists of a sphere of radius 100 nm Figure 2(a). A gold layer with a thickness of 200 nm is placed under the tip. In the simulation, the altitude of the tip (i.e., the distance from the bottom of the tip to the surface of the gold layer) is varied from 200 nm to 10 µm, and the intensity of the scattered field in the specular direction is calculated for these tip altitudes.
Although the near-field zone is relatively small in comparison to the entire s-SNOM system, its accurate simulation would still require an enormous amount of mesh elements (as well as computer memory and time) since the skin depths of tungsten and gold at 2 THz are very small (140 nm and 70 nm, respectively). In the conventional approach of meshing, the mesh size in a metal should be much smaller than the skin depth. In our design, the working area including the tungsten probe is divided into two sections Figure 2(b): the lower part (approximately as shown in Figure 2(b) and the upper part which is the rest of the working area (top section outside Figure 2(a). In the upper part of the probe and in the gold layer above the substrate, the impedance boundary condition can be applied [13]. Such a boundary condition simulates a real metal only on its boundaries, and allows avoiding the meshing inside the metal. This greatly decreases the number of mesh points in the calculation. In the upper part the tungsten probe is relatively thick and therefore it transmits hardly any radiation. The field in the narrow gap between the top of the tungsten probe and the PML is very weak and it has no effect on the detected specular radiation. In the lower part of the probe we still make use of the normal (physical) boundary conditions, since the size of that part is comparable to the skin depth. The integration line in Figure 1(b) is far from this lower interaction part and the field on the integration line is effectively an outgoing wave. The specular reflection inten-sity is readily obtained from this field in the form of far-field diffraction.

RESULTS AND DISCUSSION
Applying the model described above, we have calculated the specular reflection as a function of the tip altitude. The experimental results of this dependence contradict the Mie scattering model when the probe tip is brought close to the gold surface [8]. The simulations (triangle curve in Figure 3) show that when the altitude is higher than 1.5 µm, approaching the tip towards the gold surface leads to a decrease of the THz signal. This can be explained by the shadow area of the probe which increases (or equivalently, the directly illuminated area decreases) as the gap between the probe and the gold layer becomes smaller. However, for the distances smaller than 1.5-2 µm, the detected signal starts to increase while the shadow is still increasing. This result is hard to explain involving Mie scattering theory. Chen et al. believe that the capacitance of the tip-surface system plays an important role in this phenomenon [8]. The authors have proposed an antenna model which takes into account the geometric shape of the entire tungsten tip rather than just the terminating sphere at the end of the tip [7]. In this model, the probe tip-surface configuration is represented by an antenna-like circuit. The incident signal is dissipated in ohmic heating and re-radiation, corresponding to resistance R. Additionally the circuit is characterized by inductance L and capacitance C. Such a simplified, yet physical, approach demonstrates good agreement with the experimental data [8]. According to this model, the detected specular reflection signal P sp can be calculated as where Z 0 ≈ 377 ohm is the vacuum impedance, R = 6000 ohm, L = 0.3 nH, and the radiation frequency ω = 4π × 10 12 Hz [7]. Eq. (1) takes into account the linear dependence of the scattered signal on a geometrical shadow of the tip. We have to slightly modify the original antenna 09007-3 Journal of the European Optical Society -Rapid Publications 4, 09007 (2009) Y. Li, et. al. model [7,8], to specify the relationship between the capacitance C and the tip altitude in our geometry. The underlying idea is that the system's capacitance changes when the tip altitude is varied. The capacitance can be approximated by an expression resembling that of the usual parallel-plate capacitor [15], viz., where A e f f is an effective area of the capacitor, ε 0 is the vacuum permittivity, and h is the tip altitude. However, Eq. (2) must be appropriately fitted since our 2D capacitor differs from the conventional 3D two-plate capacitor. In particular, the effective area A e f f accounts for the specific shapes of the 2D capacitor "plates" (i.e., a triangular-shaped tip and a straight line), and an additional term C bg , identified as a background capacitance, is introduced. Such a capacitance can be treated as the influence of the shaft on the antenna properties [16]. It is physically reasonable assume that contribution from relatively large conical tip to capacitance of the tipsample system can not change with the rate proportional to small change of the tip altitude. Most efficient coupling is expected at the tip altitudes where antenna dissipates maximum of the incident field resulting in the resonance drop of the scattered field [17]. To match this physical antenna model with the explicit simulation data, we find that A e f f = 4 × 10 −12 m 2 and C bg = 0.015 fF. The specular reflection intensity, calculated with these fitting parameters in Eqs. (1) and (2), is shown by the solid-dot curve in Figure 3. Figure 3 clearly manifests good matching of the numerical simulations and theoretical predictions obtained with the antenna model both in the "shadow" (the tip altitude higher than about 2 µm) and the "resonance" regions (the tip altitude lower than 2 µm). The results are also qualitatively in good agreement with the experimental data as regards the shapes of the curves and the position of the "resonance" [8]. The Mie scattering model failed to explain this resonant region. One cannot expect quantitatively exact agreement between 2D calculations and 3D measurements, since in a real experiment the majority of specularly reflected light comes from an area on the side of the probe tip. Therefore the influence of the tipsample interaction on the specular reflection is rather weak, especially at higher tip altitudes. In 2D studies, on the other hand, the specular reflection can originate only from directly underneath the tip and so the specular reflection intensity depends more strongly on the tip-sample separation. However, this model contains most important physical features to describe properly the tip-sample system. For example, substituting semi-isolating Si for the gold layer removes the resonance behavior of the scattered signal. Low concentration of free electrons in the new layer (which leads to relevant change of the material permittivity) is not enough to generate plasmons responsible for dipole coupling with the tip. Thus, only "traditional" specular scattering behind the tip shadow is obtained with the modeling Figure 4 that is expectable result for real physical systems. An additional factor limiting the rigorous comparison with the experimental data is a single frequency field used in the simulation. One can assume that the broadband spectrum of a real pulse results in larger range of the tip altitudes where the tip-sample system is affected by the resonance response. This leads to flattening the resonance be-haviour and, consequently, to slower increase of the reflected signal than it occurs with monochromatic field. The effect can be clearly noticed from comparison with the experimental results [8]. The important point is that in the 2D geometries the s-SNOM interactions at THz frequencies can be rigorously computed (without recourse to particular theories such as Mie scattering) and the results can be accurately represented by a simple antenna model with appropriately fitted parameters.

CONCLUSIONS
We have demonstrated a simplified 2D model of THz SNOM imaging system. The near-field scattering calculations involving the metallic probe tip and the sample surface are performed using the COMSOL Multiphysics package. The modeling and numerical computations make use of standard FEM algorithms, as well as the technique of perfectly matched layers and the impedance boundary conditions for metals. The important results are that in 2D geometries the tip-enhanced s-SNOM interactions at THz frequencies can be rigorously computed, without recourse to particular theories such as Mie scattering, and that the SNOM system can be accurately represented by a simple antenna model with appropriately fitted parameters.