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.
© 2010 Optical Society of America
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  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)  and finite element method (FEM) . 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 with polarizabilities located at positions . Each dipole has polarization3], and is the time harmonic E-field amplitude at each dipole location due to the incident field plus contributions from dipoles:8]. Defining the diagonal tensors as , and substituting into Eqs. (1, 2),2 ):2); thus A is a square matrix. Once 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 .
The original implementation of the DDA  used the Clausius–Mossotti polarizability, given by6), α 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 :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 
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 . The single and double precision interaction matrices consume and , 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],
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 , including the phase shifts, in the z direction,13, 14) are applied to the planar component in Eq. (11), the reflected dipole fields become11) on the first term, with the integration variable here being λ, the above equation is simplified to3).
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,14], which can be calculated using13, 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),24)] and [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 , faster convergence can be achieved by evaluating the Hankel function formulation (Fig. 6 ),
The essential integrals and 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,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,
dipole:5), and here we explicitly define37, 38, 39, 40, 41, 42, 43, 44, 45), is46) is the same as it would be for the structure in free space. Note that we have omitted Coulomb’s constant in every term in Eq. (46) and henceforth for computational reasons and to be consistent with ; 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 to a medium with a lower refractive index , the incident and transmitted angles (from normal to the media interface) are governed by Snell’s law,57)— 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 by60), in Eq. (58) is complex, and that results in the component of the field in the z direction, , being exponentially decaying.
The Fresnel formulas for the TM and TE modes,61, 62)] to give the corresponding transmitted fields in the second medium,
5. SCATTERED FAR FIELD
The far-field zone can be defined at any given point that is much farther than the intra-dipole distances . The scattered far field  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  with those from the image (reflected) dipoles,8 ). The radial component of the scattered field approaches zero in the far field, and thus 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,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  we can employ the scheme for calculating the T-matrix for an imaging system discussed in .
We tested our MATLAB functions that calculate Eqs. (32, 33, 34, 35), which contain the Sommerfeld integrals, against the FORTRAN subroutines from . Input coordinates were taken from the upper half-space , 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 , 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 plane, i.e., . The incident light was an s-polarized plane wave with a wavelength of , at which the refractive index of the Si substrate is , 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 . The scattered far field at a given point of detection can be calculated using Eq. (67). The phase function of the normalized intensity is depicted in Fig. 10 . Next, we changed the diameter of the sphere to 300 nm and , and plotted the phase function in Fig. 11 from 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 . 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  and silicon  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 .
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 , 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 under 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 , 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 [Fig. 16a], the effect of the AFM probe appears to produce resonance modes toward the longer wavelength end of the spectrum, peaking at , while the resonance of the sphere itself [Fig. 16c] is still evident at . As the AFM tip is drawn away, e.g., , its effect is lessened, and the resonance peak approaches that for the sphere by itself.
The AFM probe also produces a sharp resonance peak, for both cases and (Fig. 16), in the short wavelength end of the spectrum at , 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 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  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  where the conversion formulas for the four permutations of the source (dipole) and observer locations (above/below surface) are listed.
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 ) 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.
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.
2. E. Purcell and C. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J. 186, 705–714 (1973). [CrossRef]
3. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). [CrossRef]
4. M. A. Taubenblatt and T. K. Tran, “Calculation of light scattering from particles and structures on a surface by the coupled-dipole method,” J. Opt. Soc. Am. A 10, 912–919 (1993). [CrossRef]
5. R. Schmehl, B. M. Nebeker, and E. D. Hirleman, “Discrete-dipole approximation for scattering by features on surfaces by means of a two dimensional fast Fourier transform technique,” J. Opt. Soc. Am. A 14, 3026–3036 (1997). [CrossRef]
6. G. H. Yuan, P. Wang, Y. H. Lu, Y. Cao, D. G. Zhang, H. Ming, and W. D. Xu, “A large-area photolithography technique based on surface plasmons leakage modes,” Opt. Commun. 281, 2680–2684 (2008). [CrossRef]
7. R. Fikri, D. Barchiesi, F. H’Dhili, R. Bachelot, A. Vial, and P. Royer, “Modeling recent experiments of apertureless near-field optical microscopy using 2d finite element method,” Opt. Commun. 221, 13–22 (2003). [CrossRef]
8. J. D. Jackson, Classical Electrodynamics, 3rd ed. (Wiley, New York, 1998).
9. B. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. 333, 848–872 (1988). [CrossRef]
10. H. Kimura, “Light-scattering properties of fractal aggregates: numerical calculations by a superposition technique and the discrete-dipole approximation,” J. Quant. Spectrosc. Radiat. Transf. 70, 581–594 (2001). [CrossRef]
11. A. Sommerfeld, “Über die Ausbreitung der Wellen in der drahtlosen Telegraphie,” Ann. Phys. (Leipzig) 28, 665–737 (1909).
12. W. C. Chew, Waves and Fields in Inhomogeneous Media (Van Nostrand Reinhold, 1990).
13. A. Baños, Dipole Radiation in the Presence of a Conducting Half-Space (Pergamon, 1966).
14. E. K. Burke and G. J. Miller, “Modeling antennas near to and penetrating a lossy interface,” IEEE Trans. Antennas Propag. AP-32, 1040–1049 (1984). [CrossRef]
15. G. J. Burke and A. J. Poggio, Numerical Electromagnetics Code (NEC-4)—Method of Moments, Part I: Program Description—Theory, Tech. Rep. UCID-18834 (Lawrence Livermore Laboratory, 1981).
16. R. J. Lytle and D. L. Lager, Numerical Evaluation of Sommerfeld Integrals, Tech. Rep. UCRL-51688 (Lawrence Livermore Laboratory, 1974).
17. D. L. Lager and R. J. Lytle, Fotran Subroutines for the Numerical Evaluation of Sommerfeld Integrals Unter Anderem, Tech. Rep. UCRL-51821 (Lawrence Livermore Laboratory, 1975).
18. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge U. Press, 1999).
19. S. Tojo and M. Hasuo, “Oscillator-strength enhancement of electric-dipole-forbidden transitions in evanescent light at total reflection,” Phys. Rev. A 71, 012508 (2005). [CrossRef]
20. P. C. Waterman, “Matrix formulation of electromagnetic scattering,” Proc. IEEE 53, 805–812 (1965). [CrossRef]
21. M. I. Mischenko, J. W. Hovenier, and L. D. Travis, Light Scattering by Nonspherical Particles: Theory, Measurements and Applications (Academic, 2000).
22. D. W. Mackowski, “Discrete dipole moment method for calculation of the T-matrix for nonspherical particles,” J. Opt. Soc. Am. A 19, 881–893 (2002). [CrossRef]
23. V. L. Y. Loke, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “T-matrix calculation via discrete-dipole approximation, point matching and exploiting symmetry,” J. Quant. Spectrosc. Radiat. Transf. 110, 1460–1471 (2009). [CrossRef]
25. V. L. Y. Loke, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Modelling of high numerical aperture imaging of complex scatterers using T-matrix method,” in Proceedings of ELS XII Helsinki (2010), pp. 138–141.
26. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370–4379 (1972). [CrossRef]
27. E. Palik, ed., Handbook of Optical Constants of Solids (Academic, 1985).
28. L. N. Aksyutov, “Temperature dependence of the optical constants of tungsten and gold,” J. Appl. Spectrosc. 26, 656–660 (1977). [CrossRef]
29. B. J. Frey, D. B. Leviton, and T. J. Madison, “Temperature-dependent refractive index of silicon and germanium,” Proc. SPIE 6273, 62732J (2006). [CrossRef]
31. F. N. Dönmezer, M. P. Mengüç, and T. Okutucu, “Dependent absorption and scattering by interacting nanoparticles,” in Proceedings of the Sixth International Symposium on Radiative Transfer (2010).
32. H. Zhang and E. D. Hirleman, “Prediction of light scattering from particles on a filmed surface using discrete-dipole approximation,” Proc. SPIE 4692, 38–45 (2002). [CrossRef]
33. Y. Harada and T. Asakura, “Radiation forces on a dielectric sphere in the Rayleigh scattering regime,” Opt. Commun. 124, 529–541 (1996). [CrossRef]
34. A. G. Hoekstra, M. Frijlink, L. B. F. M. Waters, and P. M. A. Sloot, “Radiation forces in the discrete-dipole approximation,” J. Opt. Soc. Am. A 18, 1944–1953 (2001). [CrossRef]
35. T. A. Nieminen, V. L. Y. Loke, A. B. Stilgoe, G. Knöner, A. M. Brańczyk, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Optical tweezers computational toolbox,” J. Opt. A, Pure Appl. Opt. 9, S196–S203 (2007). [CrossRef]