1. Introduction

Surface plasmons can be excited when nanoparticles are illuminated at an appropriate optical frequency. They are a collective oscillation of free electrons along the surface of metallic structure. A strong enhancement of near-field amplitude and a large scattering cross section for a narrow wavelength band can be observed at the surface plasmon resonance. For a simple structure, such as a circular cylinder or a sphere, the resonance can be calculated analytically [1, 2]. However, for a finite structure with a complicated geometry, a solution of the governing Maxwell’s equations is needed to obtain its optical properties at resonance.

In addition to the geometry of the structure, the computation of surface plasmon resonance of metals such as silver and gold is difficult due to their large refractive index at optical wavelengths. Such material properties produce a very short effective wavelength inside the metals and strong field discontinuities at the boundary. Various numerical methods have been proposed to characterize the surface plasmons of metallic structures. These numerical methods can be generally categorized into time domain method and frequency domain method. Time domain method such as finite-difference time-domain requires a dispersive model to incorporate the frequency-dependent material properties [3, 4]. The accuracy of the method is highly dependent on the approximation of the dispersive model in comparison to the actual material properties over the simulation bandwidth. Besides, the finite-difference time-domain method usually uses strucutred mesh and this produces staircase error at curved surfaces.

Various frequency domain methods, such as multiple multipole method [5], discrete dipole approximation [6], boundary element method [7, 8, 9] and volume integral equation (VIE) method [10, 11, 12], have been developed to study the surface plasmon problem of piecewise homogeneous nanoparticles. In contrast to the time domain method, these frequency domain methods do not require any dispersive model for the frequency-dependent material properties. They permit direct use of experimental permittivity data in the computation and thus provide the flexibility to solve scattering problems of material with no appropriate dispersive model. In addition, some frequency domain methods allow the use of unstructured mesh that is able to approximate curved surfaces more accurately.

The VIE approach presented in this paper offers the flexibility to model arbitrarily shaped geometry as well as an easy approximation of inhomogeneous material properties. Although the boundary element method generates less unknown, the VIE is more attractive for electromagnetic problem involving complicated material properties since a single VIE is sufficient to handle the complex problem. In this paper, we use the VIE to characterize the surface plasmon resonance behaviors of arbitrarily shaped nanoparticles. The geometry of the nanoparticles is discretized by an unstructured mesh and the material properties of the nanoparticles are obtained directly from experimental data. Different from previous work, we have adopted divergence-conformingvector basis functions in the discretization of the resultant integral equation. By using divergence-conforming vector basis functions, the continuity of fields is ensured and no artificial charges are generated during the simulation. These are important features for solving VIE accurately, especially for nanoparticles with complicated geometry and material properties. Several numerical examples will be presented to demonstrate the applicability of the proposed method for the analysis of surface plasmon resonance of complicated nanoparticles.

Our paper is organized as follows. In the next section, we give a brief description on the VIE formulation and the divergence-conforming vector basis functions. In Section 3, we demonstrate the accuracy and efficiency of our computational scheme, follow by the results of different shaped nanoparticles such as hollow and star-shaped cylinders. The effects of illumination directions and surrounding media to surface plasmon resonance will also be discussed. Finally, we draw some conclusions in Section 4.

2. Formulation

Consider an arbitrarily shaped two dimensional (2D) scatterer which is a structure with a translation symmetry along the z-axis and consists of inhomogeneous dielectric material properties. The object is embedded in an isotropic homogeneous background medium with permeability µ_b=µ₀ and permittivity ε_b. The scatterer is illuminated by an incident transverse-electric TEz wave (electric field in xy-plane), which is excited by impressed sources in the background medium. The scatterer is assumed non-magnetic (µ_r=1) and has permittivity ε(r) at location r.

By invoking the volume equivalence principle, the dielectric material can be removed and replaced by equivalent volume electric current densities J _V [13, 14]. The equivalent current densities are radiating in the unbounded background medium. By using the equivalent current densities, the total electric field E in the general dielectric material region can be written as the sum of the incident electric field and the scattered electric field due to J _V. Thus the VIE can be expressed as

where D=εE is the electric flux density, ${\eta }_b=\sqrt{{\mu }_b⁄{\epsilon }_b}$ is the intrinsic impedance of the background medium, $k_b=\omega \sqrt{{\mu }_b{\epsilon }_b}$ is the wavenumber of the background medium and $G=\frac{1}{4j}{H}_0^{\left(2\right)}\left(k_b\mid r-r\prime \mid \right)$ is the two dimensional scalar Green’s function in the background medium. It is also noted that the equivalent current densities are related to the electric flux density through

where κ_e is the contrast ratios of the permittivity and is defined as

The resultant VIE can be solved numerically by using method of moments [15]. The 2D cross section of the structure is discretized by using triangular elements that permit flexible modeling of arbitrarily shaped cross section. For a highly curved surface, finer triangular elements are needed to approximate the curvature. It is also possible to use parametric triangular elements that conform to the curved surface [14]. The dielectric property in each individual triangular element is assumed constant for inhomogeneous materials. When solving integral equations using the method of moments, the equivalent current density is normally used as the unknown quantity to be determined. However, for solving VIE, the equivalent current density is not directly used as the unknown quantity. Instead, the electric flux density D is used as the continuity of the normal components of D can be ensured by using proper basis functions. For triangular elements, the suitable basis functions are the triangular rooftop basis functions [16, 17]. A triangular rooftop basis function is defined on the common edge of a pair of triangular patches as

and its divergence is given as

where l_n is the length of the n^th edge, A ^± represents the area of the triangles T ^±, and ρ ⁺ and ρ ^- are vectors pointing away and toward the free vertex respectively. An example of triangular rooftop basis function is shown in Fig. 1.

 

Fig. 1. A triangular rooftop basis function.

Download Full Size | PDF

On every element, only a maximum of three basis functions, corresponding to the three edges, will exist. The features of triangular rooftop basis functions, such as divergence conforming and continuity of the normal component of flux density to the interior face, make them suitable to be implemented in the volume integral equation [17]. The electric flux density is then expanded using the triangular rooftop basis functions f_n

where D_n is the coefficients to be determined. For the n^th edge located at the exterior boundary of the object, an auxiliary triangle is introduced in the exterior region where the free vertex of the auxiliary triangle coincides with the center of the n^th edge. From (2), J_V can be expressed as

 

Fig. 2. SCS as a function of wavelength for a silver circular cylinder with radius =50 nm. The plasmon resonance at λ=347 nm is well reproduced by our method.

Download Full Size | PDF

Subsequently, we substitute (7) into (1), and test (1) with f_m. This results in a N×N linear system which can be written in a matrix form as

where matrix Z̄ represents the interaction between equivalent sources and column vectors V and I are respectively the excitation vector and unknown coefficients vector to be determined. The elements of matrix Z̄ can be computed using

and the elements of vector V can be computed using

Once all the elements are obtained, the matrix equation (8) can be solved either by a direct method like Gaussian Elimination or by an iterative matrix solver like Conjugate Gradient. For a large matrix problem, the solution can be accelerated using fast algorithms such as fast multipole method [18] and adaptive integral method [19].

3. Numerical results

In this section, four numerical examples are presented to demonstrate the capability of the proposed method for analyzing the surface plasmon resonance properties of nanoparticles. All the frequency-dependent material properties used in the computation are directly obtained and interpolated from experimental data [20]. In all examples, the light incident angle is the angle between the incident light and the horizontal axis. The first example is a silver circular cylinder with a radius of 50 nm. The scattering cross section of the circular cylinder is computed and shown in Fig. 2. The results obtained by using Mie series are also plotted for comparison and a very good agreement with the results obtained by using our method is observed. The plasmon resonance at 347 nm has been well reproduced by using our method.

In the second example we consider the surface plasmon resonance of elliptical, rectangular and square cylinders which have the same cross section area as the circular cylinder described in example 1. The eccentricity of the ellipse is 0.8 and the height-to-width ratio of the rectangle is 0.6. The scattering cross sections of these cylinders with light incident from angle 180° and 90° are computed and shown in Fig. 3. The scattering cross sections of the cylinders obtained using boundary element method are also computed and shown in square symbols for comparison. A very good agreement has been observed between the results computed by VIE and boundary element method. When the light incident angle varies from 180° to 90°, the resonance for the elliptical cylinder is redshifted from 343 nm to 375 nm. We also found that the full width at half maximum (FWHM) of the resonance for the elliptical cylinder has been expanded from 36 nm to 156 nm when the incident angle varies from 180° to 90°. The rectangular cylinder also exhibits same optical behaviors as elliptical cylinder, i.e. redshift of the resonance from 373 nm to 378 nm and increment of the FWHM from 28 nm to 153 nm, when light incident angle varies from 180° to 90°. The near-field distributions of these cylinders at their respective resonant frequencies are plotted in Figs. 4 and 5 for light incident angle 180° and 90°, respectively. For both incident angles, strong fields can be observed at the surface of the nanoparticles, in particular at the corners of the rectangular and square cylinders.

 

Fig. 3. SCS as a function of wavelength for different shaped cylinders computed by VIE (solid lines) and boundary element method (square symbols). For both elliptical and rectangular cylinders, their surface plasmon resonances are redshifted and their FWHM become wider when the light incident angle varies from 180° to 90°.

Download Full Size | PDF

 

Fig. 4. The normalized near-field distribution of nanoparticles at their corresponding resonant wavelengths with light incident from incident angle 180° as shown in inset of Fig. 3(a).

Download Full Size | PDF

 

Fig. 5. The normalized near-field distribution of nanoparticles at their corresponding resonant wavelengths with light incident from incident angle 90° as shown in the inset of Fig. 3(b).

Download Full Size | PDF

The third example we consider is a hollow cylinder as shown in Fig. 6(a). The inner and outer radii of the cylinder are 75 nm and 50 nm, respectively. The circular void is shifted 10 nm to the right of the center. The hollow cylinder is illuminated by light propagating along the x-axis as indicated in Fig. 6(a). The scattering cross sections of the hollow cylinder with background and void filled with different materials, i.e. air, water and glycerol, are computed and shown in Fig. 6(b). We also computed and plotted the scattering cross section of a solid cylinder in Fig. 6(b) by filling the void with silver. Generally, the solid cylinder produces a higher magnitude of scattering cross section than hollow cylinders in all three different background materials. For the solid cylinder immersed in different materials, we only recorded a single resonance in the computation and observed a redshift of the resonance when the refractive index of the background medium increases. However for the hollow cylinder, we obtained multiple resonances in the computation and observed a redshift of the resonances when the cylinder is placed in a higher refractive index background. The multiple resonances are probably due to the resonances of the void and the coupling between the void and the core of the cylinder. The near-field distributions of the solid and hollow cylinders immersed in different background materials are computed at their respective principal resonances and are shown in Figs. 7 and 8, respectively. In the near-field plots of the hollow cylinders, we found that a strong field enhancement can be obtained inside the void for all three different background materials but the field strength is weakened when higher refractive index background materials are used.

The last example we considered is five- and six-pointed star-shaped cylinders as shown in Fig. 9(a) [21]. The length d of these two star-shaped cylinders are set equal to 100 nm and the corners of the stars are filleted with radius r. Two values of r, i.e. 1 nm and 5 nm have been considered in this example. The stars are illuminated by light propagating along -y-axis and their scattering cross sections are computed and shown in Fig. 9(b). The surface plasmon resonances of the five- and six-pointed stars with fillet radius r=1 nm are 435 nm and 400 nm, respectively. As for the stars with fillet radius r=5 nm, the resonances of the five- and six-pointed stars are blueshifted to 393 nm and 383 nm, respectively. The near-field distributions of the five- and six-pointed stars are computed at their respective resonant frequencies and are shown in Fig. 10. In the near-field plots for the stars with r=1 nm, we observed a very strong field near the corners transverse to the illumination direction. However, a more evenly distributed field surrounding the corners are observed in the near-field plots for the stars with r=5 nm.

 

Fig. 6. Geometry and scattering cross section as a function of wavelength for a hollow (solid) silver cylinder. The cylinder is considered solid when n ₁=n _silver. The solid and hollow cylinders are immersed in different background materials.

Download Full Size | PDF

 

Fig. 7. The normalized near-field distribution of a solid silver cylinder (n ₁=n _silver) immersed in different background materials.

Download Full Size | PDF

 

Fig. 8. The normalized near-field distribution of a hollow silver cylinder (n ₁≠n _silver) immersed in different background materials.

Download Full Size | PDF

 

Fig. 9. Geometry and scattering cross section as a function of wavelength for silver five-pointed and six-pointed stars. The corners of the stars are filleted with different radii r.

Download Full Size | PDF

 

Fig. 10. The normalized near-field distribution of five- and six-pointed stars with different fillet radii r.

Download Full Size | PDF

4. Conclusion

In this paper, an efficient method based on the volume integral method has been presented for the analysis of surface plasmon resonance of nanoparticles. This method uses directly the experimental frequency-dependent material data and is applicable to modeling of arbitrarily shaped object. The resultant integral equation is subsequently solved by using method of moments with the use of divergence-conforming vector basis functions. Four numerical examples have been presented in this paper which show that the surface plasmon resonance of arbitrarily shaped cylinders can be captured by using the proposed method. The results demonstrated that the VIE approach is very good for fast and accurate analysis of the surface plasmon behaviors of homogeneous and inhomogeneous arbitrarily shaped metallic nanoparticles.