## Abstract

To address the large number of parameters involved in nano-optical problems, a more efficient computational method is necessary. An integral equation based numerical solution is developed when the particles are illuminated with collimated and focused incident beams. The solution procedure uses the method of weighted residuals, in which the integral equation is reduced to a matrix equation and then solved for the unknown electric field distribution. In the solution procedure, the effects of the surrounding medium and boundaries are taken into account using a Green’s function formulation. Therefore, there is no additional error due to artificial boundary conditions unlike differential equation based techniques, such as finite difference time domain and finite element method. In this formulation, only the scattering nano-particle is discretized. Such an approach results in a lesser number of unknowns in the resulting matrix equation. The results are compared to the analytical Mie series solution for spherical particles, as well as to the finite element method for rectangular metallic particles. The Richards-Wolf vector field equations are combined with the integral equation based formulation to model the interaction of nanoparticles with linearly and radially polarized incident focused beams.

© 2009 Optical Society of America

## 1. Introduction

Nano-optics is a rapidly growing field with a diverse set of existing and emerging practical applications. Near-field optical techniques that enhance localized surface plasmons may obtain intense optical spots beyond the diffraction limit for optical data storage [1]. The magnetic storage industry is also interested in sub-wavelength optical spots for heat assisted magnetic recording to overcome the superparamagnetic limit [2-4]. The interaction of light with nanostructures reveals unique information about the structural and dynamic properties of matter, and is of great importance for biological and solid-state applications. Nano-optical transducers have been widely used in near-field scanning optical microscopy [5,6]. Hartschuh et al. [7] obtained 20 nm resolution images of carbon nanotubes using an apertureless configuration. The resolution and scanning time of the scanning near-field optical microscopes, however, are limited by the spot size and transmission efficiency of the nano-optical systems. Therefore, advances in nano-optical transducers benefit scanning near-field optical microscopes and applications that utilize these microscopes. In addition, intense sub-wavelength optical spots have potential applications in nanolithography [8] and bio-chemical sensing [9]. All of these applications benefit from small optical spots. The transmission efficiency of nano-optical systems should also be maximized for practical applications since transmission efficiency of the nano-optical system will determine the data transfer rate of storage devices and scan times of near-field scanning microscopes.

Various parameters have to be optimized in order to achieve large transmission efficiency while keeping the optical spot size well below the diffraction limit. These parameters include not only geometry-dependent parameters and the material composition of the nano-optical transducer, but also source-dependent parameters, such as operational wavelength and the numerical aperture of the incident beam. Selecting an optimum set of parameters for a nano-optical transducer is important in achieving small spots and large transmission efficiencies. Optimizing the performance of nano-optical parameters requires modeling and simulation of these structures through 3-D full-wave solutions of Maxwell’s equations. An extensive parametric study of the aforementioned transducers requires efficient and accurate solutions of Maxwell’s equations.

Due to the large number of geometry, material composition, and source-related parameters, the development of efficient and accurate modeling and simulation tools for near-field optical systems is necessary. In this study, an integral equation based numerical solution is developed for nano-optical particles when they are illuminated with collimated and focused incident beams. The numerical technique developed in this study requires only the discretization of the nano-optical transducer, rather than the entire structure. Therefore, it results in a fewer number of unknowns than the numerical algorithms currently being utilized for solutions of nano-optical systems, such as finite difference time domain and finite element method. This study is organized as follows: In Sec. 2, a detailed formulation of the integral equation based numerical solution is provided. The integral equation is discretized into a matrix equation using the method of weighted residuals. In Sec. 3, the results of the numerical technique are compared to the results of the analytical Mie series solution for spherical particles and the finite element method for rectangular metallic particles. In Sec. 4, the formulation is extended to incident focused beam excitation. Richards-Wolf vector field equations are combined with the integral equation based formulation to model linearly and radially polarized focused beams. Concluding remarks appear in Sec. 5.

## 2. Method of weighted residuals

In this section, we provide a formulation for an integral equation based modeling and design tool for nano-optical systems. Similar tools have been successfully used for the analysis and design of other nano-optical systems in the literature. Nano-optical system modeling studies in the literature utilize differential equation based approaches, such as finite difference time domain (FDTD) [10-15] and finite element method (FEM) [4,15], as well as integral equation based techniques [16-20]. Previous integral equation based techniques have not presented three-dimensional results when the incidence excitation is composed of linearly and radially polarized tightly focused beams. A tightly focused beam of incident light provides a large incident electric field onto nanoparticles, improving the near-field radiation in the vicinity of the particle. Therefore, it is desirable to obtain integral equation based solutions when the incidence excitation is composed of linearly and radially polarized tightly focused beams. In this study, a three-dimensional integral equation based solution is obtained when the incidence excitation is composed of linearly and radially polarized tightly focused beams.

A full-wave implementation of the method of weighted residuals (MWR) [21-25], which is also known as the method of moments (MoM), has a number of advantages over FDTD and FEM for nano-optical system analysis. In MWR, the effects of the surrounding medium and boundaries are taken into account using a Green’s function formulation. Therefore, MWR requires only the discretization of the nano-optical transducer, whereas FDTD and FEM require the discretization of the entire computational space. Therefore, the resulting matrix equations of the MWR are smaller in size. An additional advantage of an integral equation based approach is the reduction of the additional error due to the discretization of the boundaries. In an integral equation based approach, the boundary conditions are handled in Green's function formulation; therefore, there is no additional error due to the discretization of the boundaries. In a differential equation based approach, such as FDTD and FEM, however, there is additional error introduced into the solution due to artificial boundary conditions. In addition, the integration of complicated excitation functions, such as focused beams in a dense medium, is easier in an integral equation based MWR compared to FDTD.

In this study, an integral equation based full-wave solution of Maxwell’s equation is developed. To discretize the integral equation into a matrix equation, we employ MWR in this work. In the rest of this section, the details of the solution are given when the incident beam is a collimated beam, which is approximated by a plane wave. In Sec. 4, the formulation will be extended to linearly and radially polarized focused incident beams, which are represented by the Richards-Wolf vector field equations.

The total electric field is a result of the interaction of an incident optical beam with a nanoparticle. The total electric field *E*⃗* _{tot}*(

*r*⃗) is composed of two components

where *E*⃗* _{inc}* (

*r*⃗) and

*E*⃗

*(*

_{scat}*r*⃗) are the incident and scattered electric field components, respectively. The incident electric field can be defined as the electric field propagating in space in the absence of a scattering object. The scattered electric field

*E*⃗

*(*

_{scat}*r*⃗) in Eq. (1) represents the fields resulting from the interaction of the incident field

*E*⃗

*(*

_{inc}*r*⃗) with the particles. In three-dimensional space, the scattered field

*E*⃗

*(*

_{scat}*r*⃗) can be written as

where *J*⃗(*r*⃗) is the induced current over the particle, *ω* is the angular frequency, *μ* is the permeability, and

is the dyadic Green’s function in free space at point *r*⃗ due to a point source at point *r*⃗′ . By applying the boundary conditions at the surface of a conducting metal, the electric field integral equation is obtained as

In Eq. (4), *E*⃗* _{inc}^{tng}* (

*r*⃗) is the tangential component of the known incident electric field on the particle and the

*J*⃗(

*r*⃗) is the unknown induced current on the nanoparticle. This is a Fredholm’s-type integral equation of the first kind since the unknown appears inside the integral. To solve Eq. (4) for

*J*⃗(

*r*⃗), we will expand it into a summation

where *b*⃗* _{j}* (

*r*⃗) represents known basis functions with unknown coefficients

*I*.

_{j}In this work, triangular rooftop basis functions are used to discretize the induced current over the nanoparticle. These basis functions are originally proposed by Glisson and Wilton [26] on rectangular domains and used on triangular domains by Rao et al. [27]. Triangular rooftop basis functions have been very popular due to their ability to model realistic geometries. Particle geometry is discretized in order to expand the induced current with triangular basis functions. Discretization examples for particles are shown in Figs. 1 (a) and (b) for spherical and rectangular particles, respectively. Over these triangular domains the induced current can be discretized using the triangular rooftop basis functions, which is illustrated in Fig. 2.

Mathematical representation of the triangular basis function associated with the n^{th} edge is given as

where *l _{n}* is the length of the edge,

*A*

_{n}^{+}is the area of the triangle

*T*

_{n}^{+}, and

*A*

_{n}^{-}is the area of the triangle

*T*

_{n}^{-}. In Fig. 2, two triangles

*T*

_{n}^{+}and

*T*

_{n}^{-}are the two triangles associated with the n

^{th}edge of the discretized particle.

Substituting the expansion given in Eq. (5) back into the Eq. (4), and changing the order of the integration and summation we obtain

Due to the approximation of the induced current with the summation in Eq. (5), there is a residual error in Eq. (7). The residual error in space can be written as

In the method of weighted residuals the error is distributed so that it is minimized in the minimum mean square sense. For this purpose, a new set of functions, known as weighting functions *W⃗ _{i}*(

*r*⃗) are used. The residual error

*ℜ*(

*r*⃗) is distributed in space by equating the inner product of the residual error

*ℜ*(

*r*⃗) with the weighting function

*W*⃗

*(*

_{i}*r*⃗) to zero

The weighting function *w*⃗* _{i}*(

*r*⃗) in Eq. (9) can be selected in a number of different ways [23]. In this study, Galerkin’s weighting method is chosen, in which weighting functions are selected as identical to basis functions. Such a selection yields the best result in the minimum mean square sense.

By placing the weighting functions into Eq. (9) we can obtain the resulting equations for the unknown coefficients of the basis functions. After mathematical manipulations, the result can be expressed as a system of linear equations as

where *Z _{i,j}* is the impedance matrix element on the i

^{th}row and j

^{th}column which is given as

and *V _{i}* is the excitation source element on the i

^{th}row given as

By solving the matrix equation in Eq. (10), we obtain the unknown coefficients of the basis functions in the induced current expansion in Eq. (5).

The matrix and vector elements in Eq. (11) are obtained using numerical integration techniques over triangular domains. An important issue in evaluating Eqs. (11) and (12) is the singularity in the kernel of the integrals. For the diagonal elements *Z _{i,i}*, the observation point and source point can be very close to each other or even coincide. In such instances, the numerical integration diverges, even though the integrals in Eq. (11) are integrable. To avoid this numerical problem, the singularity extraction technique is applied in Eq. (11). The integrals in Eq. (11) are divided into two parts: (1) the part that can be treated using the numerical integration, and (2) the part that is evaluated analytically. For example, the first term on the right hand side of Eq. (11), which has a first order singularity, can be separated into numerically- and analytically-treatable parts as

In this study, the singularity extraction technique for triangular domains [28] is employed to avoid numerical inaccuracy.

The singularity in the second term on the right hand side of Eq. (11) is third order, which is more difficult to handle analytically compared to a first order singularity. However, for the triangular rooftop basis functions used in this study, the second term on the right hand side of Eq. (11) can be simplified as

which has a first order singularity and is handled with the formulations given in the literature [28].

## 3. The interaction of metallic nanoparticles with a collimated beam

Using the integral equation based formulation given in the previous section, the interactions of a collimated beam with both a conducting metallic sphere and cube are studied. The collimated beam is modeled as a linearly polarized plane wave propagating in the *z* direction, which is mathematically represented as

The results of the integral equation based solution are compared to the results of the analytical Mie series solution for spherical particles and the finite element method for rectangular metallic particles [29].

To compare the results, the radar cross sections of particles are calculated for different cross sections of the far-field. The radar cross section is defined as

where *E*⃗* _{inc}*(

*r*⃗) and

*E*⃗

*(*

_{scat}*r*⃗) are the incident and scattered fields, respectively. The incident field for this problem is defined by Eq. (15) and the scattered field is obtained using a far-field approximation of Eq. (2). In the far-field region, the distance between the source and the observation point can be approximated by

Substituting Eq. (17) back into Eq. (2), the scattered field in the far-zone can be written as

Substituting the coefficients obtained from the solution of the matrix equation into the induced current expansion in Eqs. (5) and (18), and changing the order of the integration and summation, the expression for the far-zone scattered field can be simplified as

Using Eqs. (15), (16), and (19), the scattering cross section of various particles now can be obtained.

In Fig. 3, the radar cross section of a sphere with a radius of 140 nm is presented to compare MWR results with the analytical Mie series solution. The operating wavelength of the laser source is 700 nm. A comparison of the MWR results with the analytical Mie series solution shows a good agreement between the results. The percent relative error between the MWR results and Mie series solution is presented for the 140 nm particle in Fig. 4. The results suggest that the error is smaller than 1.7 % on various cuts. The main source of this error is the discretization of the nanoparticle. This error can be further reduced by increasing the number of unknowns. A similar comparison is provided in Fig. 5 for a larger sphere with a radius of 350 nm. The results in Fig. 5 also show a good agreement between the results. The relative error between the results is presented in Fig. 6. The results suggest that the maximum error for this case is about 2.5 %.

In Fig. 7, the scattering cross section of a conducting metallic cube with a side length of 200 nm is obtained on various cuts in the far-field. There is no analytical solution for a cube, therefore, we utilized an FEM solution as a reference. Similar to the previous calculations, a linearly polarized plane wave is utilized. The operating wavelength is 700 nm. In Fig. 7 (a) and (b) the θ component of the radar cross section is plotted on ϕ=0° and ϕ=90° cuts. The MWR and FEM results show a good agreement.

The method in this study is capable of addressing the near-field computations. Once the unknown coefficients in Eq. (5) are calculated, these equations can be substituted back into the electric field integral to calculate the near-field distributions.

## 4. Linearly and radially polarized focused beam

It is also very desirable to obtain solutions when the incidence excitation is composed of linearly and radially polarized focused beams. In the previous sections, the integral equation based solutions are provided when the incident beam is a plane wave. In this section, the formulation is extended to the case where the incident beam is a focused beam. In the previous formulation, the incident electric field in Eq. (15) was substituted in Eq. (12) to find the excitation vector elements. The main difference in the formulation in this section is that instead of Eq. (15), we utilize the electric field distribution for a focused beam obtained from a lens system.

Richards and Wolf developed a method for calculating the electric field semi-analytically near the focus of an aplanatic lens [30, 31]. Using the Richards-Wolf method, we can obtain both transverse and longitudinal components near the focus for both linear and radial polarizations. The total electric field in the vicinity of the focus is given as

where *α* is the half angle of the beam. In Eq. (20), *a*⃗(*θ*,*ϕ*) is the weighting vector for a plane wave incident from the *θ*, *ϕ* direction. Here it should be noted that *a*⃗(*θ*, *ϕ*) is a polarization dependent quantity, which is given as

for linear and radial polarizations, respectively. The interaction of a focused beam with linear or radial polarization can be obtained by substituting Eq. (20) back into Eq. (12).

In Fig. 8, various components of the near-field radiation from a sphere are plotted when the incident beam is a linearly polarized focused beam obtained from an optical lens system with a numerical aperture of 0.85. The operating frequency is 700 nm. The results are plotted for spherical particles with radii 70 and 140 nm. The *E _{x}* and E

_{z}components are plotted on the

*ϕ*=

*π*/2 cut as a function of

*θ*. For small spheres, the

*E*component has a maximum at

_{x}*θ*=

*π*/2. As the spherical particle gets larger, we observe a shift of the location at which the

*E*component has a maximum field. This is due to the increased interaction between a larger sphere and a wider range of angular components of a focused beam. As the size of the spherical particle gets larger, the particle interacts more with components that are incident to larger angles. A similar shift is also observed in the

_{x}*E*component, as shown in Fig. 8.

_{z}In Fig. 9, various components of the near-field electric field are plotted for a radially polarized incidence beam. The results are plotted for spherical particles with radii 70 and 140 nm. *E _{x}* and

*E*components are plotted on the

_{z}*ϕ*=

*π*/2 cut as a function of

*θ*. The incident beam parameters are identical to the previous set of results with the exception that a radial polarization is used instead of a linear polarization. Contrary to the results in Fig. 8,

*E*shows a minimum at

_{x}*θ*=

*π*/2 in Fig. 9. This is due to the difference in the strength of various components of the linearly and radially polarized incident focused beams. For the linearly polarized focused wave, the

*x*-component of the electric field is much stronger than the other two components. The radially polarized wave, on the other hand, has a strong

*z*-component in the focal region.

## 5. Conclusion

In this work, an integral equation based numerical solution was developed. The formulations for both plane waves and focused beams were given. For focused beams, the Richards-Wolf vector field equations were combined with the integral equation based formulation to model both linearly and radially polarized focused beams. The results of the integral equation based solution were compared to the results of the analytical Mie series solution for spherical particles and the finite element method for rectangular metallic particles. The methods showed a good agreement.

## Acknowledgments

This work was performed with the support of the European Community Marie Curie International Reintegration Grant (IRG) Agreement Number MIRG-CT-2007-203690.

## References and links

**1. **T. D. Milster, “Horizons for optical data storage,” Optics and Photonics News **16**, 28–32 (2005). [CrossRef]

**2. **R. E. Rottmayer, S. Batra, D. Buechel, W. A. Challener, J. Hohlfeld, Y. Kubota, L. Li, B. Lu, C. Mihalcea, K. Mountfield, K. Pelhos, C. Peng, T. Rausch, M. A. Seigler, D. Weller, and X. Yang, “Heat-assisted magnetic recording,” IEEE Trans. Magn. **42**, 2417–2421 (2006). [CrossRef]

**3. **T. W. McDaniel, W. A. Challener, and K. Sendur, “Issues in heat-assisted perpendicular recording,” IEEE Trans. Magn. **39**, 1972–1979 (2003). [CrossRef]

**04. **K. Sendur, W. Challener, and C. Peng, “Ridge waveguide as a near field aperture for high density data storage,” J. Appl. Phys. **96**, 2743–2752 (2004). [CrossRef]

**5. **D. W. Pohl, W. Denk, and M. Lanz, “Optical stethoscopy: image recording with resolution *λ*/20,” Appl. Phys. Lett. **44**, 651–653 (1984). [CrossRef]

**6. **A. Lewis, M. Isaacson, A. Harootunian, and A. Muray, “Development of a 500 Å spatial resolution light microscope,” Ultramicroscopy **13**, 227–231 (1984). [CrossRef]

**7. **A. Hartschuh, E. J. Sanchez, X. S. Xie, and L. Novonty, “High-resolution near-field Raman microscopy of single-walled carbon nanotubes,” Phys. Rev. Lett. **90**, 095503 (2003). [CrossRef] [PubMed]

**8. **L. Wang and X. Xu, “Numerical study of optical nanolithography using nanoscale bow-tie-shaped nano-apertures,” J. Microsc. **229**, 483–489 (2008). [CrossRef] [PubMed]

**9. **B. Liedberg, C. Nylander, and I. Lundstroem, “Surface plasmon resonance for gas detection and biosensing,” Sens. Actuators **4**, 299–304 (1983). [CrossRef]

**10. **W. A. Challener, I. K. Sendur, and C. Peng, “Scattered field formulation of finite difference time domain for a focused light beam in a dense media with lossy materials,” Opt. Express **11**, 3160–3170 (2003). [CrossRef] [PubMed]

**11. **J. T. II Krug, E. J. Sánchez, and X. S. Xie, “Design of near-field probes with optimal field enhancement by finite difference time domain electromagnetic simulation,” J. Chem. Phys. **116**, 10895 (2002). [CrossRef]

**12. **T. Yamaguchi, “Finite-difference time-domain analysis of hemi-teardrop-shaped near-field optical probe,” Electron. Lett. **44**, 4455427 (2008). [CrossRef]

**13. **T. Yamaguchi and T. Hinata, “Optical near-field analysis of spherical metals: Application of the FDTD method combined with the ADE method,” Opt. Express **15**, 11481–11491 (2007). [CrossRef] [PubMed]

**14. **L. Liu and S. He, “Design of metal-cladded near-field fiber probes with a dispersive body-of-revolution finite-difference time-domain method,” Appl. Opt. **44**, 3429–3437 (2005). [CrossRef] [PubMed]

**15. **T. Grosges, A. Vial, and D. Barchiesi, “Models of near-field spectroscopic studies: comparison between finite-element and finite-difference methods,” Opt. Express **13**, 8483–8497 (2005). [CrossRef] [PubMed]

**16. **J. P. Kottmann and O. J. F. Martin, “Plasmon resonant coupling in metallic nanowires,” Opt. Express **8**, 655–663 (2001). [CrossRef] [PubMed]

**17. **J. P. Kottmann and O. J. F. Martin, “Accurate Solution of the Volume Integral Equation for High-Permittivity Scatterers,” IEEE Trans. Antennas Propag. **48**, 1719–1726 (2000). [CrossRef]

**18. **J. P. Kottmann, O. J. F. Martin, D. R. Smith, and S. Schultz, “Dramatic localized electromagnetic enhancement in plasmon resonant nanowires,” Chem. Phys. Lett. **341**, 1–6 (2001). [CrossRef]

**19. **J. P. Kottmann, O. J. F. Martin, D. R. Smith, and S. Schultz, “Field polarization and polarization charge distributions in plasmon resonant nanoparticles,” New J. Phys. **2**, 27 (2000). [CrossRef]

**20. **J. Jung and T. Sondergaard, “Green’s function surface integral equation method for theoretical analysis of scatterers close to a metal interface,” Phys. Rev. B **77**, 245310 (2008). [CrossRef]

**21. **J. H. Richmond, “Digital computer solutions of the rigorous equations for scattering problems,” Proc. IEEE **53**, 796–804 (1965). [CrossRef]

**22. **R. F. Harrington, “Matrix methods for field problems,” Proc. IEEE **55**, 136–149 (1967). [CrossRef]

**23. **R. F. Harrington, *Field Computation by Moment Methods*, (IEEE Press, New York, NY, 1993). [CrossRef]

**24. **E. K. Miller, L. Medgyesi-Mitschang, and E. H. Newman, Eds., *Computational Electromagnetics* (IEEE Press, New York, NY, 1992).

**25. **R. C. Hansen, Ed., *Moment Methods in Antennas and Scattering*, (Artech, Boston, MA, 1990).

**26. **A. W. Glisson and D. R. Wilton, “Simple and efficient numerical methods for problems of electromagnetic radiation and scattering from surfaces,” IEEE Trans. Antennas Propag. **28**, 593–603 (1982). [CrossRef]

**27. **S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag. **30**, 409–418 (1982). [CrossRef]

**28. **D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler “Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains,” IEEE Trans. Antennas Propag. **33**, 276–281 (1984). [CrossRef]

**29. **
All the FEM calculations in this report are performed with High Frequency Structure Simulator (HFSSTM) from Ansoft Inc.

**30. **E. Wolf, “Electromagnetic diffraction in optical systems I. An integral representation of the image field,” Philos. Trans. R. Soc. London Ser. A **253**, 349–357 (1959).

**31. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. tructure of the image field in an aplanatic system,” Philos. Trans. R. Soc. London Ser. A **253**, 358–379 (1959).