Surface waves and atomic force microscope probe-particle near-field coupling: discrete dipole approximation with surface interaction

Vincent L. Y. Loke; M. Pinar Mengüç

doi:10.1364/JOSAA.27.002293

Journal of the Optical Society of America A
Vol. 27,
Issue 10,
pp. 2293-2303
(2010)
•https://doi.org/10.1364/JOSAA.27.002293

Surface waves and atomic force microscope probe-particle near-field coupling: discrete dipole approximation with surface interaction

Vincent L. Y. Loke and M. Pinar Mengüç

Open Access

Get PDF
Email
Share
Get Citation
Copy Citation Text
Vincent L. Y. Loke and M. Pinar Mengüç, "Surface waves and atomic force microscope probe-particle near-field coupling: discrete dipole approximation with surface interaction," J. Opt. Soc. Am. A 27, 2293-2303 (2010)

Export Citation
- BibTex
- Endnote (RIS)
- HTML
- Plain Text
Citation alert
Save article
Spotlight Summary

More Like This

Radiation forces in the discrete-dipole approximation
A. G. Hoekstra, et al.
J. Opt. Soc. Am. A 18(8) 1944-1953 (2001)

Second harmonic scattering from small particles using Discrete Dipole Approximation
Naveen K. Balla, et al.
Opt. Express 18(21) 21603-21611 (2010)

Discrete-dipole approximation for periodic targets: theory and tests
Bruce T. Draine, et al.
J. Opt. Soc. Am. A 25(11) 2693-2703 (2008)

Related Topics
Table of Contents Category
- Optics at Surfaces
Optics & Photonics Topics
?

The topics in this list come from the Optics and Photonics Topics applied to this article.

About this Article
History
- Original Manuscript: June 25, 2010
- Revised Manuscript: August 9, 2010
- Manuscript Accepted: August 9, 2010
- Published: September 28, 2010
Virtual Issues
Virtual Journal for Biomedical Optics Vol. 5, Iss. 14

September 24, 2010 Spotlight on Optics

Abstract

Evanescent waves on a surface form due to the collective motion of charges within the medium. They do not carry any energy away from the surface and decay exponentially as a function of the distance. However, if there is any object within the evanescent field, electromagnetic energy within the medium is tunneled away and either absorbed or scattered. In this case, the absorption is localized, and potentially it can be used for selective diagnosis or nanopatterning applications. On the other hand, scattering of evanescent waves can be employed for characterization of nanoscale structures and particles on the surface. In this paper we present a numerical methodology to study the physics of such absorption and scattering mechanisms. We developed a MATLAB implementation of discrete dipole approximation with surface interaction (DDA-SI) in combination with evanescent wave illumination to investigate the near-field coupling between particles on the surface and a probe. This method can be used to explore the effects of a number of physical, geometrical, and material properties for problems involving nanostructures on or in the proximity of a substrate under arbitrary illumination.

1. INTRODUCTION

According to geometrical optics, when an incident wave undergoes total internal reflection off the boundary between two different media, no energy is transmitted from the first to the second medium. In reality, while this remains true, an evanescent wave exists along the boundary surface. There will be a real wave vector component parallel to the surface and an imaginary component normal to the surface; in other words, the evanescent wave propagates along the surface, but the field intensity decays exponentially away from the surface. However, energy can be drawn from the surface if an object, in the second medium, with a relatively higher refractive index is brought within the near-field distance, i.e., a fraction of a wavelength from the surface.

As a practical example, we investigate the effects of an atomic force microscope (AFM) probe in the vicinity of a nanoparticle on a planar substrate with surface waves (Fig. 1 ). Our experimental collaborators [1] have demonstrated that the AFM probe can selectively heat and melt metallic nanoparticles by exploiting the phenomenon of plasmon resonances. This nano-writing capability may be further developed for nanolithography, including possible fabrication of nanocircuits. To fully account for the heating of a nanoparticle, the field intensity, radiative transfer, spontaneous emission, and other relevant phenomena need to be considered. As a starting point, we provide a methodology to calculate the electromagnetic field in nanostructures on a substrate and investigate conditions that contribute to higher field intensity.

For this purpose, we employ a modified discrete dipole approximation (DDA) method where the surface interactions are accounted for in the formalism. In the DDA [2, 3], the scattering object is represented as an array of point dipoles that interact with the incident field, as well as with one another. The DDA is extensively used for modeling light interaction with particles with arbitrary shapes, and that can be made from inhomogeneous and anisotropic materials. This method can be applied to objects ranging from Rayleigh particles to relatively large structures (of sizes many times the incident wavelength); the object size is usually limited by computational time and available random-access memory (RAM). The other limitation is that the DDA does not account for reflections when the structure is in the proximity of a surface. To account for surface effects, the dipole interactions with their image counterparts need to be considered. Using information from prior work [4, 5] we implemented what we believe to be a new formalism of the DDA with surface interaction, called DDA-SI, in MATLAB.

There are other viable methods for modeling the AFM tip-surface system, e.g., finite difference time domain (FDTD) [6] and finite element method (FEM) [7]. What advantage one method holds over the others depend on many factors such as the symmetry of the system, the sparseness of the scattering objects, the need for multiple frequencies, the volume of the required computational domain, etc. A thorough comparison between the methods (DDA, FDTD, and FEM) for this type of problem is warranted; however, it is outside the scope of this paper.

2. DDA FORMALISM

The scattering object is represented by point dipoles (Fig. 1), which are essentially Rayleigh scatterers, numbered $j = 1, \dots, N$ with polarizabilities $α_{j}$ located at positions $r_{j}$ . Each dipole has polarization

P_{j} = α_{j} E_{j},

where

α_{j}

is the polarizability tensor [3], and

E_{j}

is the time harmonic E-field amplitude at each dipole location

r_{j}

due to the incident field

E_{inc, j} = E_{0} exp (i k r_{j} i ω t)

plus contributions from

N - 1

dipoles:

E_{j} = E_{inc, j} - \sum_{k \neq j} A_{j k} P_{j},

where the off-diagonal

3 \times 3

tensors of the interaction matrix,

A_{j k} = \frac{exp (i k r_{j k})}{r_{j k}} [k^{2} ({\hat{r}}_{j k} {\hat{r}}_{j k} - 1_{3}) + \frac{i k r_{j k} - 1}{r_{j k}^{2}} (3 {\hat{r}}_{j k} {\hat{r}}_{j k} - I_{3})],

where

r_{j k}

is the distance from points

r_{j}

r_{k}

, and

{\hat{r}}_{j k}

is the unit vector in the direction from points

r_{j}

r_{k}

. The dipole interaction tensor

A_{j k} = k^{2} G_{j k}

, where

G_{j k}

is the Green’s function of the electric field from the radiating dipole, derived from the vector Helmholtz equation [8]. Defining the diagonal tensors as

A_{j j} = α_{j}^{- 1}

, and substituting into Eqs. (1, 2),

E_{inc, j} = A_{j j} P_{j} + \sum_{k \neq j} A_{j k} P_{j},

the problem is reduced to solving

3 N

unknown polarizations

P_{j}

in the following system of

3 N

linear equations (Fig. 2 ):

\sum_{k = 1}^{N} A_{j k} P_{j} = E_{inc, j} .

The system of equations is exactly determined (Fig. 2); thus A is a square matrix. Once

P_{j}

is solved (usually using an iterative method), the field at each dipole (2), the scattered field, the extinction coefficient, the scattering coefficient, and other quantities of interest can be calculated [3].

The original implementation of the DDA [2] used the Clausius–Mossotti polarizability, given by

α_{j}^{CM} = \frac{3 d^{3}}{4 π} (\frac{m_{j}^{2} - 1}{m_{j}^{2} + 2}) = \frac{3 d^{3}}{4 π} (\frac{ϵ_{j} - 1}{ϵ_{j} + 2}),

which only becomes valid approaching the direct current (zero-frequency) limit

(k_{0} d \to 0)

. Consider that for a non-absorbing material, m (refractive index) will be real and, according to Eq. (6), α will also be real. We know that non-absorbing materials scatter light; therefore the individual dipoles would require an imaginary component in α. The dipole radiation would otherwise be in phase with the incident beam, meaning no scattering would occur.

By considering that an infinite lattice of polarizable points has the same dispersion relation as a continuum of refractive index m, the lattice dispersion relation (LDR) was derived by Draine [9]:

α_{j}^{LDR} = \frac{α_{j}^{CM}}{1 + \frac{α_{j}^{CM}}{d^{3}} [(b_{1} + m^{2} b_{2} + m^{2} b_{3} S) {(k d)}^{2} - \frac{2}{3} i {(k d)}^{3}]},

where

\hat{a}

and

\hat{e}

are unit vectors defining the direction and polarization of the incident light. There are recent developments in the polarizability calculation for specific applications [10], but the use of the LDR for the dipole polarizability has proven satisfactory in most cases.

The two main validity criteria for the DDA are the number of dipoles required to represent the scatterer and the minimal lattice spacing; their relationship is defined by [3]

N d^{3} = \frac{4}{3} π r^{3} .

The left-hand side is the cubic volume of the space occupied by the N-dipole lattice with spacing d. The right-hand side is the volume of the equivalent sphere with radius r. Thus, the effective radius for a sphere with N dipoles is

r = {(\frac{3 N}{4 π})}^{1 / 3} d .

The lattice spacing, relative to the wavelength of the incident light, has to be sufficiently small such that

d \leq \frac{1}{k | m |} .

For a glass in air,

d = 0.1

would suffice. However, lattice spacings on the order of

d = 0.01

are required for noble metals due to high imaginary components of their refractive indices especially at shorter wavelengths.

The interaction matrix, which has complex elements, is the component that takes up most of the RAM during the calculation; its memory requirements scale with ${(3 N)}^{2}$ . The single and double precision interaction matrices consume $8 \times {(3 N)}^{2}$ and $16 \times {(3 N)}^{2} bytes$ , respectively.

3. DDA-SI FORMALISM

The standard DDA is only valid for free space scattering; in the presence of a surface, we need to consider reflection. Firstly, any incident field originating from above the substrate will be reflected off the surface. Thus, the incident field will be a sum of the incident and the reflected waves. However, in the scenario where an incident plane wave originates from below the surface, the field above the surface is just the transmitted field, or in the case when the incident plane wave undergoes total internal reflection (Fig. 1), only an evanescent incident field exists above the surface.

Secondly, in addition to the direct dipole-to-dipole interactions, we have to account for the interactions from the reflected “image” dipoles (Fig. 3 ). The calculation of the radiated field from each reflected dipole does not involve a straight forward calculation of the image Green’s function. This is because the spherical wave from each dipole is only partially reflected. Consider that the spherical wave can be decomposed into cylindrical and planar components (Fig. 4 ), as per the following Sommerfeld identities [11, 12],

\frac{e^{i k r}}{4 π r} = \frac{i}{4 π} \int_{0}^{\infty} \frac{k_{ρ}}{k_{z}} J_{0} (k_{ρ} ρ) e^{i k_{z} | z |} d k_{ρ},

\frac{e^{i k r}}{4 π r} = \frac{i}{8 π} \int_{- \infty}^{\infty} \frac{k_{ρ}}{k_{z}} H_{0}^{(1)} (k_{ρ} ρ) e^{i k_{z} | z |} d k_{ρ},

where

J_{0}

and

H_{0}^{(1)}

, the Bessel and Hankel functions, respectively, are the cylindrical wave components of the expansions, whereas

e^{i k_{z} | z |}

is the planar component.

Now, consider the presence of the surface; the cylindrical waves expand in parallel to the surface and are thus unaltered, whereas the planar waves undergo partial reflection as the wave amplitude diminishes by factors given by the Fresnel coefficients [12], including the phase shifts, in the z direction,

R^{TE} = \frac{μ_{r} k_{2 z} - k_{1 z}}{μ_{r} k_{2 z} + k_{1 z}} e^{i k_{2 z} z_{k}},

R^{TM} = \frac{ϵ_{r} k_{2 z} - k_{1 z}}{ϵ_{r} k_{2 z} + k_{1 z}} e^{i k_{2 z} z_{k}},

for the TE and TM polarizations, respectively. Thus, when Eqs. (13, 14) are applied to the planar component in Eq. (11), the reflected dipole fields become

{(\frac{e^{i k r}}{4 π r})}_{TE}^{I} = \frac{i}{4 π} \int_{0}^{\infty} \frac{k_{2 ρ}}{k_{2 z}} J_{0} (k_{2 ρ} ρ) \frac{μ_{r} k_{2 z} - k_{1 z}}{μ_{r} k_{2 z} + k_{1 z}} e^{i k_{2 z} (z_{j} + z_{k})} d k_{ρ},

{(\frac{e^{i k r}}{4 π r})}_{TM}^{I} = \frac{i}{4 π} \int_{0}^{\infty} \frac{k_{2 ρ}}{k_{2 z}} J_{0} (k_{2 ρ} ρ) \frac{ϵ_{r} k_{2 z} - k_{1 z}}{ϵ_{r} k_{2 z} + k_{1 z}} e^{i k_{2 z} (z_{j} + z_{k})} d k_{ρ} .

The oscillation of an electric dipole above a surface can be treated as two separate components, i.e., the vertical electric dipole (VED) and the horizontal electric dipole (HED). For the VED, the ϕ component is zero, and we are only concerned with the z and ρ components, i.e., the TM field. Here, we will discuss the solution of the electric field for the VED in greater detail than that of the HED, which is solved in a similar manner.

In Eq. (16), we are left with a complicated integral; however, it can be rearranged into known forms by means of variable substitutions [13], $γ_{1} = {(λ^{2} - k_{1}^{2})}^{1 / 2} \underset{λ \to 0}{\to} - i k_{1}$ , $γ_{2} = {(λ^{2} - k_{2}^{2})}^{1 / 2} \underset{λ \to 0}{\to} - i k_{2}$ , and $ϵ_{r} = k_{1}^{2} / k_{2}^{2}$ , and also noting that $k_{ρ}^{2} = k - k_{z}^{2}$ ,

{(\frac{e^{i k r}}{4 π r})}_{TM}^{I} = \frac{1}{4 π} \int_{0}^{\infty} [\frac{- 1}{γ_{2}} J_{0} (λ ρ) e^{- γ_{2} (z_{j} + z_{k})} + \frac{2 k_{1}^{2}}{k_{1}^{2} γ_{2} + k_{2}^{2} γ_{1}} J_{0} (λ ρ) e^{- γ_{2} (z_{j} + z_{k})}] λ d λ .

Using the Sommerfeld identity (11) on the first term, with the integration variable here being λ, the above equation is simplified to

{(\frac{e^{i k r}}{4 π r})}_{TM}^{I} = - \frac{e^{i k r}}{4 π r_{I}} + \frac{k_{1}^{2}}{2 π} \int_{0}^{\infty} [\frac{1}{k_{1}^{2} γ_{2} + k_{2}^{2} γ_{1}} J_{0} (λ ρ) e^{- γ_{2} (z_{j} + z_{k})}] λ d λ,

where

r_{I} = \sqrt{ρ^{2} + {(z_{j} + z_{k})}^{2}}

is the distance between the image of the radiating dipole and the receiving point (Fig. 3).

To complete the DDA model for the scatterer in the presence of a surface, we need to include the electric field contribution from the image dipoles,

E_{j}^{I} = \frac{k_{2}^{2}}{ϵ_{2}} G_{j, k}^{I} P_{k},

where

G_{j, k}^{I}

is the Green’s function for the reflected field [14], which can be calculated using

G_{j, k}^{I} = (I + \frac{\nabla \nabla}{ϵ_{2} k_{2}^{2}}) (\frac{e^{i k_{2} r_{I, j k}}}{4 π r_{I, j k}}) = {I + (\frac{1}{ϵ_{2} k_{2}^{2}}) [\begin{matrix} \frac{\partial^{2}}{\partial ρ^{2}} & \frac{\partial^{2}}{ρ \partial ρ \partial ϕ} & \frac{\partial^{2}}{\partial ρ \partial z} \\ \frac{\partial^{2}}{ρ \partial ρ \partial ϕ} & \frac{\partial^{2}}{ρ^{2} \partial ϕ^{2}} & \frac{\partial^{2}}{ρ \partial ϕ \partial z} \\ \frac{\partial^{2}}{\partial ρ \partial z} & \frac{\partial^{2}}{ρ \partial ϕ \partial z} & \frac{\partial^{2}}{\partial z^{2}} \end{matrix}]} (\frac{e^{i k_{2} r_{I, j k}}}{4 π r_{I, j k}}) .

The electric field components [13, 15] of the VED, in cylindrical coordinates, for the source and evaluation points above the surface are obtained by applying Eq. (20) to Eq. (18),

E_{ρ}^{V} = \frac{- P_{ρ}}{4 π ϵ_{2}} (\frac{\partial^{2}}{\partial ρ \partial z} k_{1}^{2} V_{22}^{'} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{\partial^{2}}{\partial ρ \partial z} \frac{e^{i k_{2} r_{I}}}{r_{I}}),

E_{z}^{V} = \frac{- P_{z}}{4 π ϵ_{2}} [(\frac{\partial^{2}}{\partial z^{2}} + k_{2}^{2}) k_{1}^{2} V_{22}^{'} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} (\frac{\partial^{2}}{\partial z^{2}} + k_{2}^{2}) \frac{e^{i k_{2} r_{I}}}{r_{I}}],

where

V_{22}^{'} = V_{22} - \frac{2}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r_{I}}}{r_{I}},

V_{22} = 2 \int_{0}^{\infty} \frac{e^{- γ_{2} (z_{j} + z_{k})}}{k_{1}^{2} γ_{2} + k_{2}^{2} γ_{1}} J_{0} (λ ρ) λ d λ .

A similar process can be applied to obtain the electric field components of the HED,

E_{ρ}^{H} = \frac{- P_{ρ}}{4 π ϵ_{2}} cos ϕ [\frac{\partial^{2}}{\partial ρ^{2}} k_{1}^{2} V_{22}^{'} + k_{2}^{2} U_{22}^{'} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} (\frac{\partial^{2}}{\partial ρ^{2}} + k_{2}^{2}) \frac{e^{i k_{2} r_{I}}}{r_{I}}],

E_{ϕ}^{H} = \frac{- P_{ρ}}{4 π ϵ_{2}} sin ϕ [\frac{\partial^{2}}{\partial ρ^{2}} k_{1}^{2} V_{22}^{'} + k_{2}^{2} U_{22}^{'} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} (\frac{\partial^{2}}{\partial ρ^{2}} + k_{2}^{2}) \frac{e^{i k_{2} r_{I}}}{r_{I}}],

E_{z}^{H} = - cos ϕ E_{2, ρ}^{V},

where

U_{22}^{'} = U_{22} - \frac{2 k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r_{I}}}{r_{I}},

U_{22} = 2 \int_{0}^{\infty} \frac{e^{- γ_{2} (z_{j} + z_{k})}}{γ_{2} + γ_{1}} J_{0} (λ ρ) λ d λ .

Here

V_{22}

[Eq. (24)] and

U_{22}

[Eq. (29)] are called the essential Sommerfeld integrals in the zeroth order Bessel function form. They are calculated by using appropriate branch cuts and contour integration paths (Fig. 5 ) to ensure convergence [15, 16]. When

z_{j} + z_{k} < ρ

, faster convergence can be achieved by evaluating the Hankel function formulation (Fig. 6 ),

V_{22} = \int_{0}^{\infty} \frac{e^{- γ_{2} (z_{j} + z_{k})}}{k_{1}^{2} γ_{2} + k_{2}^{2} γ_{1}} H_{0}^{2} (λ ρ) λ d λ,

U_{22} = \int_{0}^{\infty} \frac{e^{- γ_{2} (z_{j} + z_{k})}}{γ_{2} + γ_{1}} H_{0}^{2} (λ ρ) λ d λ .

The essential integrals $V_{22}$ and $U_{22}$ reoccur in Eqs. (21, 22, 23, 24, 25); the motivation behind the extensive algebraic manipulation is to minimize unnecessary recalculation and to further optimize the code [15, 17]. They are further grouped in the following set of parameters,

I_{ρ}^{V} = \frac{\partial^{2}}{\partial ρ \partial z} k_{1}^{2} V_{22},

I_{z}^{V} = (\frac{\partial^{2}}{\partial z^{2}} + k_{2}^{2}) k_{1}^{2} V_{22},

I_{ρ}^{H} = (\frac{\partial^{2}}{\partial ρ^{2}} k_{2}^{2} V_{22} + k_{2}^{2} U_{22}),

I_{ϕ}^{H} = - (\frac{1}{ρ} \frac{\partial}{\partial ρ} k_{2}^{2} V_{22} + k_{2}^{2} U_{22}),

and need only be calculated once at any given evaluation point with unique

(z, ρ)

coordinates; in other words, they do not have to be calculated for every dipole. We translated the FORTRAN subroutines from [17] that calculate Eqs. (32, 33, 34, 35) for our DDA-SI MATLAB implementation.

The reflection contribution will be added to the interaction matrix of the standard DDA implementation that uses the Cartesian coordinate system; thus the following cylindrical-to-Cartesian coordinate transformation is performed,

[\begin{matrix} E_{x} \\ E_{y} \\ E_{z} \end{matrix}] = [\begin{matrix} \frac{x^{2}}{ρ^{2}} E_{ρ}^{H} - \frac{y^{2}}{ρ^{2}} E_{ϕ}^{H} & \frac{x y}{ρ^{2}} (E_{ρ}^{H} + E_{ϕ}^{H}) & \frac{x}{ρ} E_{ρ}^{V} \\ \frac{x y}{ρ^{2}} (E_{ρ}^{H} + E_{ϕ}^{H}) & \frac{y^{2}}{ρ^{2}} E_{ρ}^{H} - \frac{x^{2}}{ρ^{2}} E_{ϕ}^{H} & \frac{y}{ρ} E_{ρ}^{V} \\ \frac{x}{ρ} E_{z}^{H} & \frac{y}{ρ} E_{z}^{H} & E_{z}^{V} \end{matrix}] .

Note that, to speed up the calculation, the fractional forms of the trigonometric functions,

sin ϕ = y / ρ

and

cos ϕ = x / ρ

, are used instead; the electric field components in Cartesian coordinates are therefore as follows:

$P_{x}$ dipole:

E_{x}^{(x)} = \frac{1}{4 π ϵ_{2}} {[{(\frac{x}{ρ})}^{2} I_{ρ}^{H} - {(\frac{y}{ρ})}^{2} I_{ϕ}^{H}] - \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} (\frac{\partial^{2}}{\partial x^{2}} + k_{2}^{2})} P_{x},

E_{y}^{(x)} = \frac{1}{4 π ϵ_{2}} {\frac{x y}{ρ^{2}} [I_{ρ}^{H} + I_{ϕ}^{H}] - \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial x \partial y}} P_{x},

E_{z}^{(x)} = \frac{1}{4 π ϵ_{2}} {\frac{x}{ρ} I_{ρ}^{V} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial x \partial z}} P_{x},

$P_{y}$ dipole:

E_{x}^{(y)} = \frac{1}{4 π ϵ_{2}} {\frac{x y}{ρ^{2}} [I_{ρ}^{H} + I_{ϕ}^{H}] - \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial x \partial y}} P_{y},

E_{y}^{(y)} = \frac{1}{4 π ϵ_{2}} {[{(\frac{y}{ρ})}^{2} I_{ρ}^{H} - {(\frac{x}{ρ})}^{2} I_{ϕ}^{H}] - \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} (\frac{\partial^{2}}{\partial x^{2}} + k_{2}^{2})} P_{y},

E_{z}^{(y)} = \frac{1}{4 π ϵ_{2}} {\frac{y}{ρ} I_{ρ}^{V} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial y \partial z}} P_{y},

$P_{z}$ dipole:

E_{z}^{(z)} = \frac{1}{4 π ϵ_{2}} {\frac{x}{ρ} I_{ρ}^{V} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial x \partial z}} P_{z},

E_{y}^{(z)} = \frac{1}{4 π ϵ_{2}} {\frac{y}{ρ} I_{ρ}^{V} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} \frac{\partial^{2}}{\partial y \partial z}} P_{z},

E_{z}^{(z)} = \frac{1}{4 π ϵ_{2}} {I_{z}^{V} + \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{e^{i k_{2} r}}{r} (\frac{\partial^{2}}{\partial z^{2}} + k_{2}^{2})} P_{z} .

With the above field contributions from the image dipoles, we can complete the system of equations for the DDA-SI to solve for

P_{k}

\sum_{k = 1}^{N} A_{j k}^{SI} P_{k} = \sum_{k = 1}^{N} (B_{j} + k_{0}^{2} G_{j k} + R_{j k}) P_{k} = E_{inc, j},

where the standard terms

B_{j} + k_{0}^{2} G_{j k}

are identical to those in Eq. (5), and here we explicitly define

G_{j k} = - \frac{exp (i k_{0} r_{j k})}{r_{j k}} [\begin{matrix} β_{j k} + γ_{j k} {\hat{r}}_{j k, x}^{2} & γ_{j k} {\hat{r}}_{j k, x} {\hat{r}}_{j k, y} & γ_{j k} {\hat{r}}_{j k, x} {\hat{r}}_{j k, z} \\ γ_{j k} {\hat{r}}_{j k, y} {\hat{r}}_{j k, x} & β_{j k} + γ_{j k} {\hat{r}}_{j k, y}^{2} & γ_{j k} {\hat{r}}_{j k, y} {\hat{r}}_{j k, z} \\ γ_{j k} {\hat{r}}_{j k, z} {\hat{r}}_{j k, x} & γ_{j k} {\hat{r}}_{j k, z} {\hat{r}}_{j k, y} & β_{j k} + γ_{j k} {\hat{r}}_{j k, z}^{2} \end{matrix}],

where

r_{j k} = {[{(x_{j} - x_{k})}^{2} + {(y_{j} - y_{k})}^{2} + {(z_{j} - z_{k})}^{2}]}^{1 / 2},

{\hat{r}}_{j k, x} = \frac{r_{j k, x}}{r_{j k}}, {\hat{r}}_{j k, y} = \frac{r_{j k, y}}{r_{j k}}, {\hat{r}}_{j k, z} = \frac{r_{j k, z}}{r_{j k}},

β_{j k} = [1 - {(k_{0} r_{j k})}^{- 2} + i {(k_{0} r_{j k})}^{- 1}],

γ_{j k} = - [1 - 3 {(k_{0} r_{j k})}^{- 2} + 3 i {(k_{0} r_{j k})}^{- 1}] .

The surface reflection term

R_{j k}

, formatted from Eqs. (37, 38, 39, 40, 41, 42, 43, 44, 45), is

R_{j k} = - [\begin{matrix} {\hat{r}}_{j k x}^{I 2} I_{ρ}^{H} - {\hat{r}}_{j k y}^{I 2} I_{ϕ}^{H} & {\hat{r}}_{j k x}^{I 2} {\hat{r}}_{j k y}^{I 2} (I_{ρ}^{H} + I_{ϕ}^{H}) & {\hat{r}}_{j k x}^{I 2} I_{ρ}^{V} \\ {\hat{r}}_{j k x}^{I 2} {\hat{r}}_{j k y}^{I 2} (I_{ρ}^{H} + I_{ϕ}^{H}) & {\hat{r}}_{j k x}^{I 2} I_{ρ}^{H} - {\hat{r}}_{j k y}^{I 2} I_{ϕ}^{H} & {\hat{r}}_{j k y}^{I 2} I_{ρ}^{V} \\ - {\hat{r}}_{j k x}^{I 2} I_{ρ}^{V} & - {\hat{r}}_{j k y}^{I 2} I_{ρ}^{V} & I_{z}^{V} \end{matrix}] - \frac{k_{1}^{2} - k_{2}^{2}}{k_{1}^{2} + k_{2}^{2}} \frac{exp (i k_{0} r_{I, j k})}{r_{I, j k}} [\begin{matrix} - (β_{j k}^{I} + γ_{j k}^{I} {\hat{r}}_{j k x}^{I 2}) & - γ_{j k}^{I} {\hat{r}}_{j k x}^{I 2} {\hat{r}}_{j k y}^{I 2} & γ_{j k}^{I} {\hat{r}}_{j k x}^{I 2} {\hat{r}}_{j k z}^{I 2} \\ - γ_{j k}^{I} {\hat{r}}_{j k y}^{I 2} {\hat{r}}_{j k x}^{I 2} & - (β_{j k}^{I} + γ_{j k}^{I} {\hat{r}}_{j k y}^{I 2}) & γ_{j k}^{I} {\hat{r}}_{j k y}^{I 2} {\hat{r}}_{j k z}^{I 2} \\ - γ_{j k}^{I} {\hat{r}}_{j k z}^{I 2} {\hat{r}}_{j k x}^{I 2} & - γ_{j k}^{I} {\hat{r}}_{j k z}^{I 2} {\hat{r}}_{j k y}^{I 2} & β_{j k}^{I} + γ_{j k}^{I} {\hat{r}}_{j k z}^{I 2} \end{matrix}],

where

r_{j k}^{I} = {[{(x_{j} - x_{k})}^{2} + {(y_{j} - y_{k})}^{2} + {(z_{j} + z_{k})}^{2}]}^{1 / 2},

{\hat{r}}_{j k x}^{I} = \frac{r_{j k x}^{I}}{r_{j k}^{I}}, {\hat{r}}_{j k y}^{I} = \frac{r_{j k y}^{I}}{r_{j k}^{I}}, {\hat{r}}_{j k z}^{I} = \frac{r_{j k z}^{I}}{r_{j k}^{I}},

β_{j k}^{I} = [1 - {(k_{0} r_{j k}^{I})}^{- 2} + i {(k_{0} r_{j k}^{I})}^{- 1}],

γ_{j k}^{I} = - [1 - 3 {(k_{0} r_{j k}^{I})}^{- 2} + 3 i {(k_{0} r_{j k}^{I})}^{- 1}] .

The size of the interaction matrix in Eq. (46) is the same as it would be for the structure in free space. Note that we have omitted Coulomb’s constant

1 / 4 π ϵ_{0}

in every term in Eq. (46) and henceforth for computational reasons and to be consistent with [3]; this factor can be re-introduced, if required, when calculating the total field. The interaction matrix (46) is a square, which gives us a choice of many available iterative methods for solving the linear system; we use the generalized residual method (gmres function in MATLAB).

4. EVANESCENT FIELD

When a plane wave propagates from a medium with a higher refractive index $n_{1}$ to a medium with a lower refractive index $n_{2}$ , the incident $θ_{1}$ and transmitted $θ_{2}$ angles (from normal to the media interface) are governed by Snell’s law,

n_{1} sin θ_{1} = n_{2} sin θ_{2} .

Beyond the critical angle,

θ_{c} = {sin}^{- 1} (n_{2} / n_{1})

, total internal reflection occurs. Even then, Snell’s law holds true, although—by solving Eq. (57)—

θ_{2}

will be complex, and consequently the wave vector and the polarization vector will also be complex [18, 19]. While it remains consistent with geometrical optics that no energy flows across the boundary, there exists an evanescent field at the boundary (Fig. 7 ) that has the wave vector (in Cartesian components) defined by

k = (0, n_{2} k_{0} sin θ_{2}, n_{2} k_{0} cos θ_{2}),

where

k_{0}

is the wave number of the propagating light in vacuum, and

sin θ_{2} = (n_{2} / n_{1}) sin θ_{1},

cos θ_{2} = i \sqrt{{(n_{2} / n_{1})}^{2} {sin}^{2} θ_{1} - 1} .

The corresponding TM and TE polarization vectors are thus

e_{p} = (0, i \sqrt{{(n_{2} / n_{1})}^{2} {sin}^{2} θ_{1} - 1}, (n_{2} / n_{1}) sin θ_{1}),

e_{s} = (1, 0, 0),

respectively. Note that due to Eq. (60),

k_{z}

in Eq. (58) is complex, and that results in the component of the field in the z direction,

E_{z} = E_{0} exp (i k_{z} z)

, being exponentially decaying.

The Fresnel formulas for the TM and TE modes,

T_{2, p} = \frac{2 n_{1} cos θ_{1}}{n_{1} cos θ_{1} + (n_{1} / n_{2}) \sqrt{n_{2}^{2} - n_{1}^{2} {sin}^{2} θ_{1}}} | E_{1, p} |,

T_{2, s} = \frac{2 n_{1} cos θ_{1}}{n_{1} cos θ_{1} + \sqrt{n_{2}^{2} - n_{1}^{2} {sin}^{2} θ_{1}}} | E_{1, s} |,

can be applied to the polarization vectors [Eqs. (61, 62)] to give the corresponding transmitted fields in the second medium,

E_{2, p} = e_{p} | T_{2, p} |,

E_{2, s} = e_{s} | T_{2, s} | .

5. SCATTERED FAR FIELD

The far-field zone can be defined at any given point that is much farther than the intra-dipole distances $(r ⪢ r_{j})$ . The scattered far field [5] in the upper half-space is the sum of the field contributions as a result of the dipole moments of every dipole as per the conventional DDA [3] with those from the image (reflected) dipoles,

E_{sca} (r) = k_{0}^{2} \frac{e^{i k_{0} r}}{4 π r} \sum_{j = 1}^{N} {e^{- i k_{sca} r_{j}} [(p_{j} \cdot {\hat{e}}_{1}) {\hat{e}}_{1} + (p_{j} \cdot {\hat{e}}_{2}) {\hat{e}}_{2}] + e^{- i k_{sca} r_{I, j}} [R^{TM} (p_{j} \cdot {\hat{e}}_{1}) {\hat{e}}_{1} + R^{TE} (p_{j} \cdot {\hat{e}}_{2}) {\hat{e}}_{2}]},

where the unit vectors

e_{1}

and

e_{2}

are the Cartesian expressions of the

e_{θ}

and

e_{ϕ}

vectors in spherical coordinates, respectively, in the scattering frame (Fig. 8 ). The radial component of the scattered field approaches zero in the far field, and thus

e_{r}

is irrelevant. The reflected terms in Eq. (67) are subject to the Fresnel reflection coefficients (14, 13).

The parallel and perpendicular components of the scattered field can be expressed as linear functions of the incident field,

[\begin{matrix} E_{p, sca} \\ E_{s, sca} \end{matrix}] = \frac{e^{i k (r - z)}}{- i k r} [\begin{matrix} S_{2} & S_{3} \\ S_{4} & S_{1} \end{matrix}] [\begin{matrix} E_{p, inc} \\ E_{s, inc} \end{matrix}],

where the components of the amplitude scattering matrix are

S_{1} = \frac{- i k_{0}^{3}}{4 π ϵ_{0}} \sum_{j = 1}^{N} {e^{- i k_{sca} r_{j}} + R^{TE} e^{- i k_{sca} r_{I, j}}} p_{j}^{(2)} \cdot {\hat{e}}_{2},

S_{2} = \frac{- i k_{0}^{3}}{4 π ϵ_{0}} \sum_{j = 1}^{N} {e^{- i k_{sca} r_{j}} + R^{TM} e^{- i k_{sca} r_{I, j}}} p_{j}^{(1)} \cdot {\hat{e}}_{1},

S_{3} = \frac{- i k_{0}^{3}}{4 π ϵ_{0}} \sum_{j = 1}^{N} {e^{- i k_{sca} r_{j}} + R^{TM} e^{- i k_{sca} r_{I, j}}} p_{j}^{(2)} \cdot {\hat{e}}_{1},

S_{4} = \frac{- i k_{0}^{3}}{4 π ϵ_{0}} \sum_{j = 1}^{N} {e^{- i k_{sca} r_{j}} + R^{TE} e^{- i k_{sca} r_{I, j}}} p_{j}^{(1)} \cdot {\hat{e}}_{2} .

It may be useful to calculate the T-matrix [20, 21] when repeated calculations of the scattered field are required for different incident conditions. The T-matrix contains the scattering properties of a particle or scattering system for a particular wavelength, and it needs to be calculated only once. Subsequently, given any illumination, the scattered field can be quickly calculated. There are various methods to calculate the T-matrix via the DDA; they are discussed in [22, 23]. In the case of the particle(s) and substrate scenario, we would calculate the T-matrix for the upper half-space since that is where detectors in the laboratory are practically situated. In situations where the scattered field or image is captured by an objective lens [24] we can employ the scheme for calculating the T-matrix for an imaging system discussed in [25].

6. BENCHMARKING

We tested our MATLAB functions that calculate Eqs. (32, 33, 34, 35), which contain the Sommerfeld integrals, against the FORTRAN subroutines from [17]. Input coordinates $r (ρ, z)$ were taken from the upper half-space $(z > 0)$ , and the maximum difference between the two implementations was approximately 0.02%.

Although our intended application uses evanescent field illumination as opposed to a plane wave from the upper half-space, the full implementation was benchmarked against prior work in [4], where the scattered far field from objects (of a few 100 nm in size) was modeled. This is valid because the model allows for arbitrary illumination. The scattering geometry is shown in Fig. 9 ; in the examples to follow, the detector is always on the incident $(y z)$ plane, i.e., $ϕ = π / 2$ . The incident light was an s-polarized plane wave with a wavelength of $λ = 632.8 nm$ , at which the refractive index of the Si substrate is $3.85 + 0.018 i$ , and the variable incident angle is γ.

A 540 nm diameter polystyrene sphere, which has a refractive index of 1.59, is illuminated by a plane wave incident at $γ = 0 °$ . The scattered far field at a given point of detection $r_{d}$ can be calculated using Eq. (67). The phase function $(θ = 0 ° - 90 °)$ of the normalized intensity $[(I / I_{0}) k^{2} r_{d}^{2}]$ is depicted in Fig. 10 . Next, we changed the diameter of the sphere to 300 nm and $γ = - 65 °$ , and plotted the phase function in Fig. 11 from $θ = - 90 °$ to 90° due to the asymmetry. Considering that the scattering profile is very sensitive to the diameter of the refractive index of the object on the surface, the results are sufficiently accurate for our AFM model where experimental errors will be appreciably higher.

7. AFM PROBE MODEL RESULTS AND DISCUSSION

The model involving an AFM probe, particle, and substrate system is depicted in Fig. 1. The particle is a 20 nm gold (Au) sphere that rests on a silicon (Si) surface. The AFM probe is gold and pencil-shaped. Its conical tip has a 50 nm base and a 20 nm hemispherical end. It has three degrees of freedom (x, y, and z directions), and the shaft length is varied for the study, but the other dimensions are otherwise unchanged. Simulation results for lateral displacements are not included in this paper.

The system is illuminated by an incident plane wave that undergoes total internal reflection, of either parallel (TM) or perpendicular (senkrecht, TE) polarization, from below the media interface, and the resultant evanescent field above the surface (Fig. 7) is calculated using Eq. (65) or Eq. (66), respectively. The angle of incidence below the Si substrate surface is 45° with the field intensity set arbitrarily to 1. For the simulations involving dimensional variations, the incident wavelength is 632 nm in vacuum, at which the refractive index for Si is $3.88 + 0.02 i$ . For spectral simulations, the wavelength range is 200–900 nm, and the corresponding refractive indices for Au and Si are obtained from [26, 27], respectively.

The refractive index data quoted for the temperature of 300 K is used in the simulations. The optical properties of gold [28] and silicon [29] are temperature dependent; and, in practice, the temperatures that need to be considered would range up to the melting point of gold. Until we account for temperature dependence, we can only approximate that the initial conditions are indicative of what happens at higher temperatures. Another factor to consider is that the refractive index of gold is also size and shape dependent [30].

In this study, we are particularly interested in the internal field of the nanostructures rather than the scattered field. This involves calculating the electric field at a given dipole using Eq. (2), which requires that we first solve for the dipole moments using Eq. (46). The nano-sphere is fashioned using 32 dipoles in a cubic lattice array with a lattice spacing of 4 nm. We first compare the effects of the TM and TE evanescent illumination; the field intensity was expectedly higher for the former (Fig. 12 ). This is because the polarization vector for the TE field (62) lacks the vertical component when the vertical coupling between the AFM tip and the sphere would be dominant given the geometry of the configuration.

Motivated by the effects of the length of the AFM probe from the DDA calculations [31], for free space interactions, we investigated the equivalent effects of the shaft length in our pencil-shaped AFM probe. Although, existing AFM probe tips are mostly pyramidal or conical, we attached a variable-length cylindrical shaft to study the suspected effects of reflection on the field intensity and thus the pencil shape. Figure 13 shows that the maximum intensity was achieved for the shaft length of $L = 60 nm$ under $λ = 632 nm$ TM evanescent field illumination. The field intensity was significantly lower under TE evanescent field illumination as shown in Fig. 14 , showing very little effect, and virtually no effect past $L > 50 nm$ , from the varying shaft length. This is consistent with the lack of vertical coupling resulting from the TE illumination.

Figure 15 shows the effect of the tip separation on intensity. Naturally, the system would have a different response if the wavelength of the incident field is changed. We calculated the intensity of the scattered far field for wavelengths from 200 to 900 nm for the tip separations of 1 and 14 nm (Fig. 16 ). For $d = 1 nm$ [Fig. 16a], the effect of the AFM probe appears to produce resonance modes toward the longer wavelength end of the spectrum, peaking at $λ \approx 810 nm$ , while the resonance of the sphere itself [Fig. 16c] is still evident at $λ \approx 630 nm$ . As the AFM tip is drawn away, e.g., $d = 14 nm$ , its effect is lessened, and the resonance peak $(λ \approx 680 nm)$ approaches that for the sphere by itself.

The AFM probe also produces a sharp resonance peak, for both cases $d = 1 nm$ and $d = 14 nm$ (Fig. 16), in the short wavelength end of the spectrum at $λ = 270 nm$ , although its height is a few orders smaller than those for longer wavelengths. In all cases, including the absence of the AFM probe, the field intensity in the particle is diminished for $λ < 500 nm$ due to the imaginary component of the refractive index of Si (Fig. 17 ) that increases toward the shorter wavelength range, giving rise to absorption.

The near field can also be calculated [23] for any given point outside the structures in the upper half-space. This can provide information on the field distribution surrounding the AFM probe and the particle. It is also possible to calculate the field under the surface, i.e., within the substrate; the equations for field evaluation points below the surface are available in [15] where the conversion formulas for the four permutations of the source (dipole) and observer locations (above/below surface) are listed.

8. CONCLUSION

We implemented the DDA-SI, which comprises the DDA method with the inclusion of surface interaction. It is a flexible tool for studying the light interaction with structures, ranging from a few nanometers to a few micrometers in size, on a planar surface under arbitrary illumination. As a demonstration of the methodology, we presented a preliminary study of near-field coupling between the AFM probe, a particle, and a surface, which is by no means exhaustive, to investigate and demonstrate the applicability of the method to the problem.

The proximity of an AFM probe clearly enhances the near-field coupling, increasing the field intensity in the particle on the substrate; intensities higher than that of the incident field is possible. Further study would involve experimentally realistic shapes that mirror industrial requirements, different materials, arbitrary arrangements of nanoparticles, multiple probes, and other combinations. The model will also be extended to filmed surfaces (as per DDFILM [32]) and stratified media; extra terms, i.e., Sommerfeld-like integrals that represent reflections from multiple surfaces, will need to be evaluated. It is also possible to model a localized rough surface by fashioning the “terrain” using dipoles with the same optical properties as the substrate, available RAM permitting.

The inclusion of other phenomena such as radiative heat transfer and emission will be the focus of ongoing work on the model. Optical forces and torques can also be calculated from the DDA model by calculating the gradient and scattering forces experienced by individual dipoles [33, 34] or a whole object [23, 35]. Casimir or van der Waals forces are only significant factors if we are working with non-illuminated or vacuum environments, in which case, they may need to be considered.

ACKNOWLEDGMENTS

This project is funded by TUBITAK-1001 (Grant No. 109M170) and FP-7-PEOPLE-IRG-2008 (Grant No. 239382 NF-RAD) at Özyeğin University, Istanbul, Turkey.

Fig. 1 Dipole model of an AFM probe, a particle, and an infinite surface. The evanescent wave above the surface results from the internally reflected plane wave from below the surface.

Abstract

1. INTRODUCTION

2. DDA FORMALISM

3. DDA-SI FORMALISM

4. EVANESCENT FIELD

5. SCATTERED FAR FIELD

6. BENCHMARKING

7. AFM PROBE MODEL RESULTS AND DISCUSSION

8. CONCLUSION

ACKNOWLEDGMENTS

Cited By

Figures (17)

Equations (75)

Journal of the Optical Society of America A