An Integral Equation based Numerical Solution for Nanoparticles Illuminated with Collimated and Focused Light

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. 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.


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][3][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 nanooptical systems. Therefore, advances in nano-optical transducers benefit scanning near-field optical microscopes and applications that utilize these microscopes. In addition, intense subwavelength 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 nanooptical 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 nearfield 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.

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][11][12][13][14][15] and finite element method (FEM) [4,15], as well as integral equation based techniques [16][17][18][19][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][22][23][24][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 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 represents the fields resulting from the interaction of the incident field ) (r E inc with the particles. In three-dimensional space, the scattered field ) (r E scat can be written as (2) where ) (r J is the induced current over the particle, ω is the angular frequency, µ is the permeability, and (3) 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 (4) In Eq. (4), ) (r E tng inc is the tangential component of the known incident electric field on the particle and the ) (r J 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 ) (r J , we will expand it into a summation (5) where ( ) r b j represents known basis functions with unknown coefficients j I .
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 Substituting the expansion given in Eq. (5) back into the Eq. (4), and changing the order of the integration and summation we obtain (7) 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 (8) 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 The weighting function ( ) r w i 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 (10) where j i Z , is the impedance matrix element on the i th row and j th column which is given as (11) and i V is the excitation source element on the i th row given as (12) 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 i i Z , , 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  (13) 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 (14) which has a first order singularity and is handled with the formulations given in the literature [28].

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 (15) 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 (16) where ) (r E inc and ) (r E scat 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 (17) 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 (19) 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  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.

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.   9. Electric field components when a radially polarized focused beam of light interacts with spheres of various sizes. The radially polarized focused beam is obtained from an optical lens system with a numerical aperture of 0.85. The operating frequency is 700 nm.

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.