Development of Computational 3 D MoM Algorithm for Nanoplasmonics

In this paper, we present an algorithm for full-wave electromagnetic analysis of nanoplasmonic structures. We use the three-dimensional Method of Moments to solve the electric field integral equation. The computational algorithm is developed in the language C. As examples of application of the code, the problems of scattering from a nanosphere and a rectangular nanorod are analyzed. The calculated characteristics are the near field distribution and the spectral response of these nanoparticles. The convergence of the method for different discretization sizes is also discussed.


I. INTRODUCTION
Nanoplasmonics studies interaction of optical fields with metallic nanostructures beyond the diffraction limit of light.At optical frequencies, metals exhibit electrical charge oscillations known as plasmons or surface plasmon resonances [1]- [3].The resonances of these electrical oscillations depend on the electrical properties of the metal element, on its dimensions and geometry and on polarization of the incident electromagnetic wave.In literature, one can find analysis of plasmonic nanoparticles with different geometries, such as spheres, rods [11], stars [12], nanoburguer [13], tetrahedral [14], circular [15] and triangular nanodisks [16], etc.The growing number of publications in this area in recent years can be explained by numerous possible applications such as superresolution microscopy [4], nanoantennas for nanophotonics [5]- [9], ultra-high-density optical data storage devices [10], among others.
In this paper, we present a computational algorithm developed for nanoplasmonics which is based on three-dimensional Method of Moments (3D MoM) [17].Using this algorithm, we analyze electromagnetic scattering of optical fields by metallic nanoparticles with different geometries.To characterize the complex dielectric constant of gold nanoparticles in optical frequencies, we use in the analysis the Lorentz-Drude model.The computational implementation of the method is performed using the language C. As examples of application of the code, two problems of scattering of a nanosphere and a rectangular nanorod are resolved.We calculate distribution of the near field, the spectral response and analyze the resonances of these particles.The convergence of the method for different discretization sizes is also investigated.To validate the developed code, we compare our results with simulations carried out by the commercial package Comsol Multiphysics and by the analytical Mie model for spherical particle.

A. Description of the Problem
In Fig. 1, the geometries considered in this paper are shown.The problems consist of the electromagnetic scattering from a single nanosphere (Fig. 1a) and from a rectangular nanorod (Fig. 1b) made of gold centered at the origin of the rectangular coordinate system in free space., where 0  is the permittivity of free space and the r  the relative permittivity which can be approximated by the Lorentz-Drude model for the gold with one interband term as follows: , c is the speed of light, and j is the imaginary unit.This model is a good approximation for the complex permittivity at wavelengths greater than 500nm.It characterizes the dispersion of metals at optical frequencies [2].1.075 10 s    

B. Integral Equation
The general electromagnetic scattering problem is presented in Fig. 2a.The total electric field E outside the volume V obj of the nanoparticle is given by the sum of the incident plane wave and the scattered from the nanoparticle wave: where superscripts i and s indicate the incident and scattered fields.The latter can be viewed as the field radiated by an equivalent polarization current density () ' eq J r , with ' r  obj V , as shown in Fig. 2b.The scattered fields obey the Maxwell's equations [19]: The equivalent polarization current density () ' eq J r exists only inside the material, and it is given by 0 - where r' is the vector that indicates the source point, and 3) and ( 4) we obtain the wave equation  k ω μ ε and  the wavelength.The solution of ( 6) is given by s ( ) , ) ( () 0 G is the dyadic Green's function for free space defined by 0 1 ,, ( ) ( ) , e 4 () Equation ( 7) is used to calculate the scattered field outside the volume V obj of the nanoparticle.
However, to calculate the scattered field inside the nanoparticle, where there is a singularity, one should modify (7) according to [20] s 'j ωμ '  GG r r r r , PV means the principal value of the integral and the second term is a correction factor.Substituting (10) in (2), we obtain the following integral equation

C. Solution by 3D MoM
This section presents a numerical solution of the integral equation ( 11) by 3D MoM.Firstly, we write the integral equation ( 11) in scalar form given by In this equation, we set 1 x = x , 2 x = y and 3 x = z .To solve the integral equation ( 12) by 3D MoM, we divide the volume V ob j into N subvolumes V m (m=1, …, N), where ()  r are constant in each subvolume.With m r as the point in the center of this subvolume, applying (12) to each subvolume m V , we obtain , , (14) An equivalent representation for ( 13) is This equation can be written as where [G] is a matrix of order 3N×3N, while [E] and [E i ] are vectors of dimension 3N.The elements of [G] are given by 0 where The total electric field at each point m r is determined by inverting the matrix [G].The elements of this matrix in (16) are calculated approximately by [19], [21]: p q mn mn mn n n mn x x mn mn pq x x mn mn mn j k V exp j G j cos cos j r (17) for ,  mn and Obviously, the higher the number of subvolumes N the better is approximation.The geometry of the subvolumes V n should be very close to a sphere of radius n a .However, good results can be obtained using cubic cells as we show in this work.
After finding the electric field inside the volume V obj , one can calculate electric field anywhere outside the nanoparticle by inversion of the system To do this, we use the solution [E] with ( 7) and ( 2) in points outside V obj .

III. NUMERICAL RESULTS
Based on the theoretical model presented in the previous section, we develop an algorithm in C language.To do this, we firstly define the constants and the variables of the problem.Then a 3D cubic domain of height and length 4a is created, where a=60nm is the radius of the sphere.This spherical volume is divided into cubic cells of dimensions dx=dy=dz.The spherical volume with radius a is created in the point (x 0 , y 0 , z 0 ) The central sphere is positioned at the origin of coordinate axis 0 0 0 ( 0) x y z    .For the case of the nanorod, we create a rectangular volume with the length L =60nm and square cross section with the size w =10nm.The axis of the nanorod is placed along the axis x, and it is centered at the origin of the coordinate system.The method does not require an Absorbing Boundary Conditions (ABCs) to simulate a free space radiation because the Green's function already takes into account the radiation condition.The ABCs are used in other numerical methods, such as the FDTD or professional programs such as software Comsol Multiphysics.The latter uses the Finite Element Method (FEM) and a special type of ABC, known as perfectly matched layer (PML), or radiation condition, for an artificial absorption of electromagnetic waves.

All the cells created inside the object are excited by a plane wave [E i ]. Subsequently, the elements of the matrix [ ] 3N X 3N G
are calculated for the volume of the object using ( 17) and ( 18) for mn  and mn  , respectively.The field inside the volume V obj is obtained by the inversion of the linear system This linear system is solved by the Gaussian elimination method for complex numbers.Then we calculate the total field in any point outside the volume V obj with (2) and with the fields inside V obj calculated in the previous step.

A. Nanosphere
Initially, for validation of the algorithm we compare our results with the classical solution for the sphere by analytical Mie model [18].We also compare our results with those obtained by the software Comsol Multiphysics.Subsequently, also we analyze a rectangular nanorod to confirm the possibility of implementing the algorithm for different geometries.All the simulations were realized in a core i7 computer with 16G of RAM.
The analyzed nanosphere has the radius a  60nm.In this case, we use the total number of cubic elements N=1791 with the dimensions dx=dy=dz=8×10 −9 .The results of calculations for the near fields in the axis x, y and z are shown in Figs.3-5.The amplitudes of these fields are normalized to the magnitude of the incident plane wave E i0 .The electromagnetic wavelength used in all simulations is =550nm, which is near the resonance.In Comsol simulation, we created a mesh of 131638 elements inside a spherical domain with diameter D=460nm.The spherical domain is limited by a PML absorbing boundary condition to simulate the open space.In Figs.3-5, we observe a good agreement between the results for the used discretizations.However, for the points near the surface of the sphere one can note a certain difference in the results.This can be explained by rapid variation of the field in this region.This means that we need a finer discretization in these regions to obtain a higher accuracy.We observe from this figure a good agreement between the calculated results.For the points near the sphere, electromagnetic fields vary rapidly with distance, therefore, it is necessary a finer discretization of the nanosphere, as it was observed also in Figs.3-5.

B. Rectangular Nanorod
The analyzed rectangular nanorod has length L= 60nm, and width w =10nm.For this case, we performed two simulations with different number of total cubic cells elements of the discretization N to estimate convergence of the method.In one simulation, we used N=750, cubic cells of dimensions dx=dy=dz=2×10 −9 , and in the other one we used N= 6000, with cubic cells of dimensions dx=dy=dz=1×10 −9 .We also simulated this nanorod using Comsol software.In the latter case we used a mesh with 371747 elements, and a spherical domain limited by a PML.All the simulations were realized in a core i7 computer with 16G of RAM.
We show in Fig. 7 the normalized electric field E x at the end of nanorod at the point (x = L/ 2 +d, y = 0, z = 0), where d =5nm and 20nm.Figs.8-10 show the spectral response of the normalized electric field in the axis z near the middle of the nanorod at the points (0,0, w / 2 +d), where w = 10nm is the side of the square cross section of the nanorod (Fig. 1b), and d =5nm, 10nm and 20nm, respectively.
The wavelength interval extends from 500nm to 1300nm, where the principal resonances of the nanorod occur.We observe from Figs. 7-10 a good agreement between the MoM and Comsol Multiphysics results.
These figures show that with increase of N, our results converge to the Comsol simulation.

IV. CONCLUSIONS
We presented in this paper a 3D MoM computational algorithm for efficient full-wave electromagnetic scattering analysis of plasmonic nanostructures.The method was codified in C language, and two gold nanoparticles were analyzed: nanosphere and rectangular nanorod.To validate the computational code, we compared our results with simulations carried out by the commercial package Comsol.In case of nanosphere, we also compared the obtained results with the classical analytical Mie model.In all cases, we observed a good agreement between the results obtained by our code, the analytical model and the commercial software.We also analyzed convergence of the method.For this purpose, we used different number of cubic cells N. For large values of N, our results approach to the Comsol simulation.The full wave method is quite general and can be used to analyze plasmonic nanostructures with different geometries and excitation sources.However, the method requires a good computational capacity in terms of memory and processing speed.

Fig. 3 .
Fig.3.Normalized near field distribution for =550nm along axis x for gold sphere with radius a=60 nm.

Fig. 4 .
Fig.4.Normalized near field distribution for =550 nm along axis y for gold sphere with radius a=60 nm.

Fig. 5 .
Fig.5.Normalized near field distribution for =550nm along axis z for gold sphere with radius a=60 nm.

Fig. 6
Fig.6 shows the variation of the x component of the normalized electric field versus wavelength for different points along the axis x.The range of the analyzed wavelengths is between 500nm and 1000nm.The resonance of the sphere occurs near =550nm.The points where the electric fields were calculated are positioned in (x=a+d, y=0, z=0), a is the radius of the sphere, and d runs through the values 20nm, 40nm, 80nm and 160nm.