1. Introduction

Controlling the nano-sized metal and/or dielectric particles by light [1] is of great interest to realize complicated nano-particle composite photonic crystals [2] and to enable biological single cell or molecule manipulation [3]. Optical binding effects have been investigated experimentally for various dielectric nano- and micro-particles separate or in close proximity [4, 5]. To evaluate the optically induced forces on particles, methods such as generalized Mie theory [6], coupled dipole approximation [7, 8], Maxwell’s stress tensor [2, 9, 10] and that combined with bispherical coordinates [11–13] have been applied to bi-spheres and tri-spheres. However, it has been shown that the resonance states of metal nano-spheres shift drastically from visible to infrared wavelengths when the spheres approach and then touch, and that enhanced field is observed at the contact points of spheres [14,15]. At the fundamental resonance states of clusters, the optically induced forces may act differently from those observed at non-resonance or high-order resonance states. To the best of author’s knowledge, systematic numerical analysis of the infrared-induced forces on various topologies of metal nano-clusters has not been performed, which may be due to the difficulty in the treatment of plasmonic states and induced forces on metal sphere configurations, as well as the numerical expenditure of solving such problems in three-dimensions.

To analyze particles of irregular shapes for which no analytical method is available, some versatile numerical methods such as the finite-difference methods are necessary. For the analysis of optical forces on dielectric objects, not metal objects, the finite-difference time-domain (FDTD) method [16] with the Lorentz force formulation has been reported in [17]. Recently the FDTD method with Maxwell’s stress tensor has been applied to metal spheres and bi-pyramidal particles in [18]; they address the advantage of Maxwell’s stress tensor formulation over the Lorentz force formulation concerning that the surface terms of the Lorentz and the Coulomb forces on metal objects may affect the accuracy of the volume integral of the force density. In this paper, however, we investigate using the FDTD method with the Lorentz force formulation the optically induced forces on various metal nano-sphere clusters. We show first the accuracy of the numerical method by comparing to rigorous analysis of the Mie theory [19, 20] for the radiation pressure exerted on single metal spheres of various radii, as well as by comparing to published results for two separate metal spheres at close distances. For the present method the force exerted on a metal particle is calculated by a volume integral of the Lorentz and the Coulomb forces obtained from the dynamic electromagnetic fields of the FDTD method; we show later that the volume integral retains the accuracy also for the metal objects. An advantage of the volume integral of the Lorentz and the Coulomb forces would be that the integral is defined solely by the volume of the object no matter how complicated the shape is. In contrast, for Maxwell’s stress tensor formulation the surface and its normal vector must be defined for the integral, for which some uncertainty may arise when the object has any irregular or discontinuous surfaces.

We plan to investigate further a nonlinear problem of second harmonic generation in a system including various metal particles and clusters, where the principle of superposition for the typical field expansion analysis methods no longer holds because of the nonlinearity, and the time-domain analysis method would be a powerful tool. We have already examined the accuracy of the FDTD method for the analysis of surface plasmon waves on thin metal structures [21, 22]. As a next step toward more comprehensive analysis, we investigate in this paper the optically induced forces in similar, but linear, systems. We show for topologically simple straight and triangular clusters that the fundamental mode of the plasmonic resonance induced by an infrared light generally exhibits strong binding forces on the spheres. The force originates from the polarized electromagnetic field associated with the local plasmonic oscillation, and it behaves like the van der Waals attractive force between the induced multipoles of metal spheres [23]; however, the difference is that the van der Waals force acts from a distance as long as 10 nm, whereas the force investigated in this paper is induced by an infrared light and acts only when particles are touching. The mode and the frequency of the plasmon resonance on clusters are identified and the force magnitude and the direction exerted on the particles are calculated.

2. Analysis method of electromagnetic field and induced force

We use Yee’s FDTD leapfrog scheme [24] with auxiliary differential equations (ADE) [25,26], which allows effective implementation of the frequency dependent metal permittivity in the FDTD analysis. In brief, a complex permittivity of classical oscillators of the Lorentz and the Drude models is expressed by equations of motion, and they are solved together with Maxwell’s equations. The complex permittivity of various metals has been compiled as fitting equations by Rakić et.al. [27]. We consider silver spheres based on their results of the Lorentz-Drude (LD) model for the implementation into the FDTD method. No quantum size effect is taken into consideration because it has been shown appropriate if the radius of the metal sphere is larger than approximately 20 nm [28]. The LD model of the complex metal permittivity is given as a function of angular frequency ω by

To analyze electromagnetic forces, the formulation based on Maxwell’s stress tensor and that on the Lorentz and the Coulomb forces are equivalent [31]. The essential difference is the surface integral of the stress tensor and the volume integral of the force density. We choose, as described in Introduction, the Lorentz force formulation for aiming at applications to more general and complicated structures. The induced Lorentz and the Coulomb forces are computed instantaneously from the electric (E) and the magnetic (H) field distributions of the FDTD analysis by

Careful consideration is needed for the integral of Eq. (5); due to the finite spatial discretization cell size Δ, the electric field on the surface of metal objects varies rapidly over a single cell inside and outside the surface, so does the magnetic field. As a consequence, the force induced on the surface also distributes over the same region. Therefore, the integral of the force density is taken over the region that is slightly outside the objects, e.g. for a metal sphere of radius r, the integral is taken over the spherical volume of radius r + 1.5Δ, which gives convergent force values as Δ approaches infinitesimal limit as shown later. In fact if the surrounding medium is uniform, the bound charge density Eq. (3) and the bound current density Eq. (4) outside the metal object vanish, hence the integral outside the metal object does not affect the actual force integral. The numerical results of the force density distribution, which will be discussed in the later section (Fig. 7), support this fact.

 

Fig. 7 Normalized field E_x (a) and force density F_x (b) distributions for two silver spheres of r = 30 nm and g = 3 nm. Discretization is Δ = 0.375 nm. The thin (green) dashed line is the outline of the sphere. In (c) E_x is shown for idealized lossless silver-like metal sphere. The same distributions of (a) and (c) along the x-axis are compared in (d). Only the limited region near the gap for x ≥ 0 and z ≥ 0 is shown.

Download Full Size | PDF

In our analyses, the fundamental resonance frequency is first solved by the FDTD analysis using a broad-band (i.e. short-pulse) excitation, the resonance mode and the induced optical force are then analyzed by the subsequent run of the FDTD analysis with a band-limited (i.e. long-duration-pulse) excitation or with a smoothly-increasing continuous wave (CW) excitation at the resonance frequency.

The optical force is the response to the incident excitation field, and is composed of a double-frequency Fourier content and a zero-frequency content in the positive (for a repulsive force) or negative (for an attractive force) direction. Thus the time averaged value of the curve indicates the magnitude and the direction of the induced force. If a CW source excitation is used, after the response of the system has settled into a steady state, the time averaged value of the force is obtained straightforwardly. As in Eqs. (2) to (4) the optically induced force is proportional to the intensity of the illumination, the force value is evaluated by the force per unit intensity pN/(W/μm²).

To reduce the computational region and ensure the excitation of proper mode fields, symmetry boundary conditions are employed; a perfect electric conductor (PEC, i.e. odd-symmetry tangential E-fields and even-symmetry normal E-fields) and a perfect magnetic conductor (PMC, i.e. even-symmetry tangential E-fields and odd-symmetry normal E-fields) are imposed according to the structural symmetry and the incident field polarization. The perfectly matched layer (PML) absorbing boundary condition [32] is used when necessary to realize an infinitely large space beyond the boundary.

3. Numerical verification

3.1. Radiation pressure on single metal sphere

Our preliminary calculations of the TM mode resonance frequency of single silver spheres for r = 10 nm to 200 nm have shown that a discretization of Δ = r/40 gives an error within 6% of theoretical values. This Δ has been chosen as small as possible for the necessary accuracy and the available computational resource. For verification of the present force analysis method, we have already confirmed the two-dimensional (2D) analysis of forces on dielectric cylinders [17]. We analyze in this paper the radiation pressure on single silver and gold spheres of various radii for r = 20 nm to 100 nm situated in air and illuminated by a plane wave at wavelength λ = 632.8 nm as in [9]. The three-dimensional (3D) analysis of the radiation pressure on a sphere is difficult because electromagnetic fields exist in the directions both perpendicular and tangential to the sphere surface, and the convergence is not as fast as the 2D cases. In addition, to derive an accurate force value a large analysis space is required to ensure the plane wave incident field as accurately as possible, for which the computation becomes extensive. A distortion from the plane wave incidence results in a large error in the induced force as discussed later.

The analysis configuration of the radiation pressure is shown in Fig. 1. A quarter sphere is analyzed by using PEC at y = 0 and PMC at z = 0. The light field is excited by a plane wave with a y-polarized electric field. A plane wave propagates in the x-direction, which is absorbed by PML at x = 0 and x = X_max. The incident time variation was chosen to be a smoothly increasing CW in order to observe the steady state clearly. For the analysis of smaller spheres of radius r = 20 nm and 30 nm, the analysis space larger than (X_max, Y_max, Z_max) = (10r, 5r, 5r) and the outer boundary conditions of PEC at y = Y_max and PMC at z = Z_max is needed to facilitate more suitable plane wave field distribution, while reducing the influence from the mirror images of the sphere across the symmetry boundaries. The use of absorbing boundary in this case results in an incident field that is more like a beam than a plane wave. The force values are significantly small for spheres much smaller than the wavelength; if the incident field deviates from the plane wave the resultant force values are affected sensitively. However, the analysis space should be smaller than approximately a quarter wavelength in order to reduce the multiple reflection of the scattered field by the outer boundaries.

 

Fig. 1 Top and side views of the configuration for the radiation pressure analysis. Only a quarter structure is analyzed using symmetry conditions.

Download Full Size | PDF

For larger spheres of radius r = 80 nm and 100 nm, the analysis space is chosen as large as (X_max, Y_max, Z_max) = (10r, 5r, 5r), which is larger than the wavelength; to reduce the multiple reflection by the outer boundaries, PML is used at z = Z_max. The multiple reflection by the outer boundaries will cause temporal fluctuation in the time signal, which gives a substantial error in estimating the time averaged force values. The force values in these cases are large enough not to be affected seriously by the unideal plane wave field, but easily affected by the multiple reflections between the analysis boundaries. A space discretization Δ = r/20 is chosen for our available computational capacity, and only for some cases Δ = r/40 has been tested to confirm the numerical convergence.

The results are compared in Table 1 to the rigorous analytical values of the Mie theory, which have been obtained by the analysis code “BHMIE” in [20] modified to calculate the radiation pressure efficiency and verified carefully also with data in [19]. Note here that the force values per unit intensity calculated from [9] in MKSA unit are erroneously an order of magnitude smaller than the actual values, presumably because the illumination intensity in MKSA unit has been described mistakenly an order of magnitude larger. These results indicate reasonable accuracy of 3 to 10 % errors from the Mie theory. The results for Δ = r/20 and r/40 are close and their numerical convergence is confirmed. From these results the numerical volume integral of the Lorentz force density is considered to be accurate. In Figs. 2 and 3 the time variation of the incident field (upper curves) and the force received from the plane wave illumination (lower curves) are shown; they show clearly that the system has reached a steady state. At first glance the time response of the force in Fig. 2 seems to be a sinusoidal wave. However, it is slightly shifted to the positive direction, and the time average of this signal gives a positive value, which corresponds to the radiation pressure in the direction of the wave propagation. Therefore, as described previously, the distortion of the incident wave from an ideal plane wave significantly affects the time signal and eventually the force value, especially when the sphere is much smaller than the wavelength. For larger spheres comparable to the wavelength, as shown in Fig. 3, the large positive shift is clearly observed in the time signal, the relative error in the force value becomes small.

 

Fig. 2 Time evolution of the electric field and the radiation pressure exerted on a silver sphere of r = 20 nm for the analysis space 16r×12r×12r.

Download Full Size | PDF

 

Fig. 3 Time evolution of the electric field and the radiation pressure exerted on a silver sphere of r = 100 nm for the analysis space 8r×4r×4r.

Download Full Size | PDF

Table 1. Radiation pressure per unit intensity of a plane wave illumination at wavelength λ = 632.8 nm on single silver and gold spheres of various radii with parameters of the analysis space X_max×Y_max×Z_max. Unit is in [pN/(W/μm²)]. The last line in each section shows the exact value obtained by the Mie theory. For finer discretization (Δ = r/40), only a part of the table has been filled due to too extensive computation

View Table | View all tables in this article

To examine the effects of the surface components of the Lorentz and the Coulomb forces to the volume integral of the force density, which is suggested in [18] as a source of error, we have compared the results of an identical code executed both with single precision and double precision arithmetic. The results were identical to five to six significant figures, which imply the numerical round-off error in the volume integral can be omitted.

3.2. Force on two separate spheres

For further verification of the method regarding the applicability to irregular configurations, we analyze systems of separate two silver spheres of Fig. 4 for r = 30 nm with g = 1.5 nm and 3 nm in air, which have been analyzed in [11] with a static approximation. The resonance frequency was first detected by a broad-band pulse excitation, and then at the particular resonance frequency the force was calculated by a continuous wave excitation. The light field is excited by a plane wave which propagates in the y-direction (perpendicular to the surface of this paper) with an x-polarized electric field. Symmetry boundary conditions are used, i.e. PEC at x = 0 and PMC at y = 0 and at z = 0 to reduce the analysis space into an eighth. The size of the analysis space is then (X_max, Y_max, Z_max) = (4r, 2r, 3r), the outer boundaries are imposed by PEC at x = X_max, PML absorbing boundary at y = Y_max and PMC at z = Z_max. For these cases the force exerted on the sphere is mainly induced by the strong field in the gap region, which is even enhanced by the resonance and is much stronger than the radiation pressure by the incident wave. The strong gap field reduces the influence of the analysis boundaries to the calculated force values. Various discretization sizes were tested to see the numerical convergence as far as our computational capacity allows.

 

Fig. 4 Configuration of the separate two spheres.

Download Full Size | PDF

The typical time responses of the field at the gap (normalized value) and the force induced on the upper one of the spheres are shown in Fig. 5 for r =30 nm and g=1.5 nm. In these results the field and the force settle by 20 fs. The force is downward (attractive force) and the time average of the curve f_x (lower curve) gives a negative value. The analysis parameters and the results are summarized in Table 2. The resonance frequencies, the field enhancement factors and the force values are plotted in Fig.6. The convergence of the calculation is excellent as the discretization becomes small. In particular, the force values asymptotically approach to a certain value, which indicates the appropriateness of the numerical modeling. It is found from these results that Δ needs to be smaller than sixth to eighth of the gap (g/6 to g/8) to obtain an accurate force value, as well as smaller than r/20 to r/40.

 

Fig. 5 Time evolution of the electric field and the radiation pressure exerted on the upper one of the two separate silver spheres of r = 30 nm, g = 1.5 nm and Δ = 0.375 nm with a continuous wave excitation.

Download Full Size | PDF

 

Fig. 6 Numerical convergence of variables for the r = 30 nm spheres in Table 2 for (a) g = 1.5 nm and (b) g = 3.0 nm.

Download Full Size | PDF

Table 2. Resonance frequency ν_r [THz], field enhancement η_E and induced force per unit incident intensity f_x [pN/(W/μm²)] for separate two silver spheres of radius r = 30 nm with the separation g = 1.5 nm and 3 nm in air. Discretization levels are refined from Δ = r/20 = 1.5 nm to r/120 = 0.25 nm. The directions of the forces are all binding

View Table | View all tables in this article

The force values obtained by the static approximation [11] vary with the incident wavelength. However, we can find results that correspond to our calculation by reading the peak values of the resonance in their plots, and we note some appreciable differences of 10 to 20 % in the resonance frequency and one to two orders of magnitude difference in the force values. Note that the force and the field in the gap region of the metal spheres change drastically with the gap distance, thus the discrepancy may appear even larger. The reason for the discrepancy is considered as follows. The optical constants used in [11] are adopted from [33], of which complex relative permittivity has an imaginary part smaller in magnitude than that of (1)(see also [34]), e.g. (1) gives ɛ_r = −7.2 −i0.71 at 610 THz while [33] gives ɛ_r = −6.6 −i0.21 at 680 THz comparing at the resonance frequencies in Table 2 for the case of g = 1.5 nm. The static approximation in [11] assumes the field distribution obtained by solving Laplace’s equation of real numbers, for which the field is based on a static and scalar field approximation. The static approximation is valid for cases where field variation is confined in a region much smaller than the wavelength of illumination. In the present cases, however, the wavelength is 450 nm in air while the sphere diameter is 2r = 60 nm, i.e. sphere size is as large as one eighth of the wavelength, for which the approximation may not be sufficiently accurate. The field enhancement η_E in [11] depends largely on the resonance position from a local minimum (off-resonance) to a local maximum (on-resonance) from 50 to 500; if the resonance position differs the resulting force would have as much variation as ${\eta }_E^2$. Similar large discrepancy between the static approximation [11] and an electrodynamic calculation is also reported in [35].

Another reason would be that the present system of metal spheres is highly lossy with a large imaginary part of ɛ_r, which may well lead to the field penetration into metal that is different from those of the static approximation. The field distributions in an actual silver sphere and a tentatively idealized lossless silver-like metal sphere are calculated and shown in Fig. 7 (a) and (c), and their values are compared in (d). Although the qualitative distribution is similar in (a) and (c), the quantitative difference is clearly seen in (d). Consequently, it is anticipated that this penetration effect leads to the strong damping of the field, hence to the smaller field enhancement. For such a highly lossy system with a strong field variation, the full-vector field analysis is known to be more accurate. The FDTD method has some errors in the analysis of curved surfaces when using the approximation of stair casing with rectangular unit cells. However, the convergence analysis with different discretizations indicates that the numerical artifact due to the approximation of the curved surface is small. In fact, once a convergence at the infinitesimal limit of discretization is obtained, the finite-difference scheme is of slightly better than first-order in accuracy [36]. From the force distribution plot in Fig. 7(b), it is clear that the major contribution of the force locates along the surface of the sphere almost limited in the region inside and outside by Δ. The possible artifacts of the limited analysis time have been investigated using the CW excitation and found that the system has reached the steady state.

4. Analysis of metal nano-clusters

The clusters shown in Fig. 8 are analyzed in this paper: straight 2-, 3- and 4-sphere clusters and a triangular cluster of r = 20nm all surrounded by air. The detailed model of the contact point of spheres is also shown in Fig. 8. The resonance frequency of touching spheres is known to depend largely on the area of the contact of spheres [14,15]. In our numerical model, due to the fixed discretization Δ = r/40, a point contact between two spheres can not be modeled. Instead the spheres are assumed to be in contact in a finite circular area of radius d ≈ 6Δ = 0.15 r. We consider this contact area as a slight deformation of the spheres when they touch; it helps to avoid the field singularity at the contact point and to allow stable calculation.

 

Fig. 8 Configurations of the analyzed clusters and the detail of the contact point between two spheres. The y-axis is in the direction from the rear to the front face of the paper.

Download Full Size | PDF

4.1. Two-sphere clusters

The numerical convergence is tested first with the analysis of two touching spheres. The discretization cell size Δ is refined while the size of the contact area is retained, and the analysis results are compared. The calculated fundamental resonance frequency and the force exerted on one of the spheres are 252 THz, 2730 pN/(W/μm²) for Δ = r/40 = 0.5 nm, and 249 THz, 2317 pN/(W/μm²) for Δ =r/80 = 0.25 nm, respectively. Therefore, by referring to the previous convergence analysis of Fig. 6, the convergence is considered to be obtained and the discretization is sufficiently small to model the touching spheres. The computational time on an eight-core workstation was several hours for the case of Δ = 0.5 nm, and several ten hours for Δ = 0.25 nm. The settling time of the field after the smoothly increasing CW excitation is 30 fs to 50 fs as shown in Fig. 9, which depends on the cluster topology and is slightly slower than the cases of the single sphere and the two separate spheres because the resonance frequency is lower. In Fig. 9 a slight overshoot is observed because the raise time of the CW excitation was short. However after the settling time has passed, the fields and the forces exhibit stable single mode oscillation.

 

Fig. 9 Time evolution of the electric field and the radiation pressure exerted on the upper one of the two touching silver sphere cluster of r = 20 nm, and Δ = 0.5 nm with a continuous wave excitation.

Download Full Size | PDF

The optical field and the force density distributions are plotted in Fig. 10 for the fundamental resonance mode of the 2-sphere cluster. Due to the finite conductivity of silver, the field penetrates into the sphere strongly at the contact region and a strong field is observed in the narrow gap region right beside the contact area. The force distribution F_x is observed mainly along the surface of the sphere in the attractive direction in the vicinity of the contact region, whereas the forces on the top and the bottom ends of the cluster are exerted in the repulsive directions. Our calculation results indicate the maximum field enhancement η_E is approximately 70 to 90 at the contact point, thus the intensity enhancement is ${\eta }_I={\eta }_E^2\approx 5\times {10}^3\text{to}\hspace{0.17em}8\times {10}^3$. Since numerical convergence is expected and considering available computational capacity, the spatial discretization is chosen to be Δ = r/40 for the rest of the analyses.

 

Fig. 10 Normalized E_x field and force F_x distribution of the fundamental resonance mode of the 2-sphere cluster of r = 20 nm. The total force integrated over a sphere is indicated by white arrows; f_x = −2317 pN/(W/μm²).The resonance frequency is ν_r = 249 THz.

Download Full Size | PDF

4.2. Clusters of various shapes

The field and the force distributions of the 3- and 4-sphere straight clusters and the triangular cluster are shown in Figs. 11, 12 and 13, respectively. These field and force plots indicate an analogous distribution of strong fields at the contact region and relatively weaker fields at the top and the bottom of the clusters. As anticipated, clusters composed of more spheres have lower resonance frequencies. Interestingly, the total force exerted on each sphere is also attractive for the fundamental modes of these clusters. As in Fig. 13, under E_x-polarized illumination, the 3-sphere triangular cluster exhibits attractive forces that bind all the 3 spheres; namely the direction of the force exerted on the lower-left sphere is θ = tan⁻¹(f_2x/f_2z) = 51.7° (angle from the z-axis), which directs slightly downward than to the center of the upper sphere and slightly upward than to the center of the cluster. Moreover the force in the x-direction satisfies f_1x ≈ −2f_2x as the spheres are kinetically balanced. It is interesting to note that, although these forces are obtained solely from the electromagnetic field analysis, they satisfy surprisingly well the geometrically expected relations of the forces.

 

Fig. 11 E_x and F_x distribution for the 3-sphere straight cluster. r = 20 nm, f_x = −2750 pN/(W/μm²). ν_r = 206 THz.

Download Full Size | PDF

 

Fig. 12 E_x and F_x distribution for the 4-sphere straight cluster. r = 20 nm, f_1x = −887 pN/(W/μm²), f_2x = −1490 pN/(W/μm²), ν_r = 176 THz. The sizes of the arrows are not to scale.

Download Full Size | PDF

 

Fig. 13 E_x field and F_x, F_z distribution for the triangular cluster. r = 20 nm, f_1x = −91.2 pN/(W/μm²), f_2x = 45.4 pN/(W /μm²), f_2z = 35.9 pN(W/μm²), ν_r = 224 THz. The sizes of the arrows are not to scale.

Download Full Size | PDF

The resonance spectra of the various silver clusters investigated in this paper are shown in Fig. 14. The resonance peaks are in the infrared range at around 1500 nm (200 THz), with their widths 100 nm to 200 nm (approximately 10 THz to 20 THz, which agree with the results in [10]) whereas the peak positions are apart more than or comparable to the peak widths. However, the peak widths can increase as an ensemble average if the size and the shape of the spheres vary. For better-controlled particles the identification of each resonance mode will be more feasible. The advantage of considering the fundamental resonance of the whole cluster is that the peak exists in the infrared range, being far away from the higher-order resonances that are composed of the superposition of many resonances of individual spheres. Additional calculations for spheres immersed in water (modeled by a constant ɛ_r = 1.77) show that the silver clusters exhibit similar effects, and that the force values are approximately 40 % weaker than those in air. With available laser intensity of 100 mW per 100 μm², the typical force value in our results e.g. 1000 pN/(W/μm²) gives an induced force of 1 pN. Other possible forces such as thermal fluctuation and gravitational forces are found much smaller than the optical forces. The van der Waals force becomes significantly strong for touching spheres, which may exceed the present optically induced forces [11, 23]. However, these binding forces for the fundamental resonance are strong and we will further investigate their influence to the equilibrium of a specific cluster formation under a careful selection of the illumination wavelength. The infrared induced attractive forces on metal particles will have interesting effects on the clustering phenomena.

 

Fig. 14 Frequency spectra for the silver sphere clusters of r = 20 nm in air obtained by the Fourier transform of the FDTD time signals with the Hanning time window.

Download Full Size | PDF

5. Conclusion

The finite-difference analysis of optically induced forces on metal nano-spheres has been performed by using the Lorentz force formulation, which has been verified by comparison to rigorous theoretical analysis results. The finite-difference method allows flexible and accurate analysis of fields and forces on complicated metal nano-particle systems. Through the numerical analysis it has been found that the optical forces on the fundamental plasmonic resonances of the metal nano-spherical clusters of r = 20 nm are strong in the attractive direction. Simple clusters composed of up to several metal nano-spheres have the fundamental resonance in the infrared range and their modes are analogous for all the cluster topologies. For the resonance states significantly enhanced fields are observed at the contact points of the spheres, as reported in the previous papers, for which the induced binding forces have been estimated quantitatively. This will be the starting point of the further investigation of more complicated plasmonic particle systems with both linear and nonlinear optical materials.