Open Access
Add to BibSonomyAbstract
The hybrid of finite element and boundary integral (FE-BI) method is employed to predict nano-optical trapping forces of arbitrarily shaped metallic nanostructures. A preconditioning strategy is proposed to improve the convergence of the iterative solution. Skeletonization is employed to speed up the design and optimization where iteration has to be repeated for each beam configuration. The radiation pressure force (RPF) is computed by vector flux of the Maxwell’s stress tensor. Numerical simulations are performed to validate the developed method in analyzing the plasmonic effects as well as the optical trapping forces. It is shown that the proposed method is capable of predicting the trapping forces of complex metallic nanostructures accurately and efficiently.
© 2015 Optical Society of America
1. Introduction
Optical trapping of metal nanoparticles has been widely investigated experimentally since Svoboda and Block [1] stably trapped Au nanospheres. At optical frequencies, the metal’s free-electron gas can sustain surface and volume charge-density oscillations, called plasmons, with distinct resonance frequencies. 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. Making using of these effects, new approaches have been developed to stably trap and manipulate a variety of objects [2]. Various parameters have to be optimized in order to achieve nice performance of the trapping. These parameters include not only geometry-dependent parameters and the material composition of the object, but also incident beam configurations: orientation, positions, beam waist radius, etc. Modeling study is thus essential for designing and optimizing advanced optical trapping because it prevents waste of sources and time required for building and testing prototypes of novel designs by delivering the detailed analysis without their actual realizations. For a simple structure, such as a circular cylinder or a sphere, the trapping force can be calculated analytically [3, 4]. However, a numerical solution of the governing Maxwell’s equations is indispensable for structures with complicated configurations, for example, the target consisting of multiple arbitrarily shaped nanostructures.
A variety of numerical approaches have been developed which can simulate the plasmonic effects. They can be generally classified into time domain methods and frequency domain methods. Time domain method, e.g., finite-difference time-domain (FDTD), demands a dispersive model to handle the frequency-dependent material properties [5]. The accuracy of the method depends highly on the approximation of the dispersive model in comparison to the actual material properties over the bandwidth of interest. In addition, the structured mesh used in FDTD usually leads to staircase error at curved surfaces. In contrast, frequency domain methods permit direct use of experimental data describing constitutive relations in the computation. They provide the flexibility to handle complex materials with no appropriate dispersive model. Besides, some frequency domain methods allow the use of unstructured mesh that is capable of approximating curved surfaces more accurately. Among the frequency domain methods, semi-analytical methods, e.g., T-matrix method [6], discrete dipole approximation (DDA) [7] and multiple multipole method (MMP) [8], allow the computation of plasmon resonances for complex geometrical configurations. Nevertheless, most of the techniques are specific for one type of geometrical configuration, or they are not flexible enough for all the above-mentioned geometrical arrangements. Besides, the methods based on the fundamental multipole expansion potentially lead to an ill-conditioned system when modeling arbitrarily shaped nanostructures.
Recently, full-wave electromagnetic solvers in frequency domain have gained a significant interest in the optical community due to their capability and versatility. These solvers can be designed according to the integral equations (IEs) or to the differential equations (DEs). The approaches based on IEs include volume integral equation (VIE) methods [9, 10] and the surface integral equation (SIE) ones [11–21]. SIE is also named as boundary integral (BI) or boundary element (BE). Integral formulations are very attractive because the radiation condition at infinity are naturally satisfied. SIE requires only the surface mesh of the target. Methods based on SIE are limited to piecewise-constant materials. VIE formulations can handle inhomogeneous materials as well as homogeneous ones at the expense of discretizing the whole spatial domain occupied by the target. IE formulations generate a fully populated impedance matrix, which can be solved directly or iteratively. The complexity of computing the inverse of the impedance matrix reaches O(N3) in the traditional manner, where N is the number of unknowns. For applications with many unknowns, direct solution becomes infeasible. The iterative solver with the aid of a fast method to accelerate the matrix-vector multiplication (MVM) is always a nice alternative for many engineering or scientific applications. However, the convergence of the iterative solution can be a troublesome issue for the applications of interest in this work because the magnitude of the permittivity of metal can be very large at optical frequencies. The large permittivity leads to a very ill-conditioned impedance matrix due to the imbalance among the matrix elements [22]. Consequently, the iterative solution often converges slow or fails to converge. Additionally, many beam configurations should be investigated to reach an optimized trapping. The iterative solution has to be repeated many times. The finite element (FE) method [23, 24] is one of the formulations based on DEs, which is proven very accurate, flexible to handle realistic structures. It requires volume mesh for 3D targets as VIE does. Because it produces a sparse matrix, the cost of the storage and CPU time is at O(N) per iteration in iterative solvers. Although the FE matrix is always ill-conditioned, efficient solution is often available thanks to the sparsity of the matrix. However, other than the target, the surrounding medium must be discretized for the absorption boundary condition (ABC) setups [23]. In addition, the spurious reflections from the ABC would degrade the accuracy of the simulation. Actually, all of the aforementioned algorithms have their inherent pros and cons, especially for simulations of metallic plasmonic nanostructures.
In this work, we describe a hybrid of FE and BI (FE-BI) solution to take advantage of both the FE and SIE methods. FE-BI uses BI rather than ABC for the boundary condition [26–29]. The disadvantages of ABC are totally avoided while the nice properties of FE remain. The traditional FE-BI is tailored with respect to the applications of interest. Specifically, a preconditioning strategy is proposed to accelerate the convergence. What’s more, a skeletonization strategy is developed to accelerate the repetitive computation for different excitations resulting from the optimization process. The performance of FE-BI on predicting radiation pressure force (RPF) exerted on nanostructures is investigated.
The remainder of the paper is organized as follows. Section 2 discusses the FE-BI formulation to compute the RPF. After a brief introduction of the iterative approach to solve the FE-BI matrix system, Section 3 details the proposed preconditioning technique and the skeletonization strategy to accelerate the repetitive computation for different excitations. Section 4 presents numerical simulations to show the validity and performance of the proposed algorithm. Section 5 concludes the paper.
2. Evaluation of RPF by FE-BI
2.1. Maxwell’s stress tensor and RPF
The optical force exerted on a given target is evaluated by integrating the Maxwell’s stress tensor over a spherical surface S enclosing that particular target, which can be written as
where is the outward normal of S and is the Maxwell’s stress tensor. Suppose the background medium has the permittivity ε0 and permeability μ0, E(r) and H(r) are the total field on S, can be written as where * stands for the conjugate operation and is the identity matrix. According to Eq. (1) and Eq. (2), the evaluation of RPF is converted into the problem of solving E(r) and H(r) from the Maxwell’s equations. In this work, we employ the full wave numerical method, FE-BI, to solve the Maxwell’s equations.2.2. Formulation of FE-BI
Figure 1 shows the typical configuration of many arbitrarily shaped nanostructures with the relative permittivity εi and permeability μi (i = 1,2,3,⋯) that are embedded in the host region V0 characterized by the permittivity ε0 and permeability μ0. Vi denotes the inner region of structure i. Vi and V0 are assumed to be open sets and are separated by the boundary surface Si for nanostructure i, with the unit normal (directed from Vi toward V0). The structures are illuminated by the optical beam Einc and Hinc at the angular frequency ω with a time factor e−jωt. Our aim is to find the field in the host medium (E0,H0) due to the impressed field (Einc,Hinc). According to the equivalence theorem, the original problem can be decomposed into two parts: the exterior and interior equivalent problems. The BI method and the FE method, respectively, is employed to analyze them. The two equivalent problems are related by forcing the continuity of the tangential fields at the interfaces.
By invoking Love’s equivalent principle for the exterior medium, regions ∪Vi are filled up with the same material of the region V0, and the sources in Vi are removed. The equivalent currents and , positioned on the external face of the surface Si for the i-th nanostructure, generate the original scattered fields Esca and Hsca in the region V0, and null fields in the region Vi, i.e.,
Equation (3) and Eq. (4) is the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE). The scattered electric and magnetic fields can be evaluated as
with the the wave impedance for the region V0. Operator ℒ and is defined as where ⨍ stands for the principal value integration; 0 ≤ Ω(r) ≤ 4π is the internal solid angle; , denotes the homogeneous-space Green’s function with the wavenumber for the region V0.For the interior equivalent problem, the fields inside each Vi can be formulated into the variational problem with the function given by
where . The tangential fields satisfy the continuity condition at the interfaces of the exterior and interior regions. It can be written asIn accordance with the boundary conditions in Eq. (9) and the definition of (Ji,Mi), the FE variational equation Eq. (8) can be combined with the boundary integral equation of the EFIE formulation Eq. (3) or that of the MFIE formulation Eq. (4) to solve the fields (E0,H0). Since EFIE and MFIE both suffer from the interior resonance, the combined field integral equation (CFIE) is often employed as an alternative. With the combination coefficient α, it has the form
Meshing the interior regions by tetrahedrons, electric field E can be written as
where and denote the number of edges inside region Vi and on the boundary Si. and are the commonly used vector basis functions, namely, edge element function [23]. The superscripts ”I” and ”S” denote the interior and exterior regions. (xe)k is the unknown coefficient associated with . The normalized magnetic field on the boundary of interior and exterior regions has the similar form where (xh)k is the unknown coefficient associated with . Substituting Eq. (11) and Eq. (12) into Eq. (8) for all the interior regions, we arrive at the FE matrix, where T denotes the transpose, and . The matrices KII, KIS, KSI, KSS and B are the standard FE matrices. Their elements are given by with i = 1,2,3,⋯. Obviously, for edge k, spans over several neighboring elements that share it. The FE matrices in Eq. (14) are all sparse. Furthermore, KII and KSS are symmetric, and KIS = (KSI)T. Equation (13) cannot be solved unless a relation between and is established. Such a relation is provided by Eq. (10). By adopting the same FE expansions that were used in Eq. (11) and Eq. (12), the expansions for (Ji,Mi) has the form of, with r′ ∈ Si, i = 1,2,⋯. The function gk in Eq. (15) is also known as the RWG function [25]. For all regions, substituting Eq. (15) into Eq. (10) and using Galerkin’ method, we have where [P] and [Q] are two full matrices whose elements are given byThe elements in the right hand side (RHS) vector y are given by
Combining Eq. (13) with Eq. (16), we yield
3. Fast solution of FE-BI
The iterative solution strategy developed in [29] is outlined in subsection 3.1. The proposed acceleration techniques are discussed in subsection 3.2 and subsection 3.3.
3.1. The iterative solution strategy
The total number of unknowns of the FE-BI matrix system is N = NI + NS. N can be considerably large as the number of nanostructures increases, although and are always moderate. It is generally infeasible to compute the inverse of the matrix in Eq. (19). If solved iteratively, Equation (19) may fail to reach the convergence because the FE matrices are always very ill-conditioned. As revealed in [29], the combination of direct methods with iterative methods always performs nice. In that method, a matrix M is conceptually defined as,
Combining Eq. (20) with Eq. (19) yields,
Equation (21) is always well-conditioned. Iterative solver, e.g., the generalized minimum residual (GMRES), works well for it. Additionally, fast methods such as multilevel fast multiple algorithm (MLFMA) can be employed to accelerate the MVM [11, 16, 17]. In comparison with the SIE, Eq. (21) always converges much faster for nanostructures, as would be indicated in Section 4.
3.2. The preconditioning strategy
It is worthy pointing out that Nk spans over several neighboring elements that share the common edge k. The vector basis function Nk in the region Vi finds no tetrahedrons inside Vl that share edge k if l ≠ i. The FE matrices for the region i, i.e., , , , and Bi, are completely decoupled from those for the region l, , , , and Bl. For the application with 3 nanostructures, the FE-BI matrix system Eq. (19) can be reformulated as
where Pil and Qil (i,l = 1,2,3) represents the coupling between the boundary Si and Sl arising from BI. yi (i = 1,2,3) is the part of excitation vector in Eq. (18) when r locates on Si. Equation (22) can be written as withThe first term in the left hand side of Eq. (23) accounts for interactions within each region. The resultant sub-system, i.e., Zii · xi, is a standard FE-BI system that can be treated by the method depicted in Eq. (20) and Eq. (21). The second term in the left of Eq. (23) denotes the coupling among different regions. The decomposition in Eq. (23) gives rise to an efficient preconditioning strategy. Specifically, with the sub-systems’ solution, an preconditioned system can be reached as,
Obviously, the extension of the preconditioning technique in Eq. (24) to many nanostructures is quite straightforward. It should be noted that we always use the iterative solution of Eq. (21) to realize the preconditioning operation of . Additionally, the preconditioning strategy of Eq. (24) is not applicable to the application with one single nanoparticle/nanostructure.
3.3. The acceleration of optimization process
As pointed out in Section 1, various parameters have to be optimized in order to achieve nice performance of the trapping. For the Gaussian beam [30], these parameters include θ, ϕ,ro or wo, where (θ,ϕ) gives the orientation of the incident beam in the spherical coordinate; ro and wo is, respectively, the beam center and waist. Suppose κ is a set of controlling parameters uniquely describing a given light beam, Eq. (19) becomes a system with multiple excitations when κ varying. For simplicity, the notation (θ,ϕ,ro,w0) and κ are mutually used or both dropped in the following when no confusion would arise. The associated system can be written as
where Y(κ) consists of Ninc excitation column vector defined in Eq. (19). X(κ) is the corresponding unknown solution matrix. During the optimization process, κ is always over-sampled because it is hard to precisely find the suitable sampling rate. The obtained RPF thus varies smoothly with κ. Consequently, the associated Y(κ) and X(κ) are both rank deficient in the numerical point of view. Unfortunately, the rank is hard to be analytically estimated. To exploit the rank deficient property, we conduct the interpolative decomposition (ID) on Y(κ), where Ys(κ) and R(κ) are, respectively, the skeletonized excitation matrix and the projection matrix, with the number of skeleton incident beams [31]. Substituting Eq. (26) into Eq. (25) yields where Xs(κ) is the skeletonized solution matrix corresponding to Ys(κ). Equation (27) states that the complete solution matrix X can be recovered after Xs is computed. Except for the error arising from the ID skeletonization, no approximation is introduced during the above process. Therefore, the accuracy of the algorithm can be well manipulated with the aid of the error control scheme in the ID [32]. Because ID is numerically proved error controllable [32], the skeletonization process is error controllable [18, 19, 35, 36]. After applying the skeletonization scheme to Eq. (25), the number of light beams to be handled reduces from Ninc to . The acceleration rate almost reaches as the operation of Xs(κ) · R(κ) is generally quite cheap [31,33]. The acceleration rate can be very high because is always much smaller than Ninc. Since Z−1 is always expensive to compute, Xs(κ) is often solved iteratively.For a general problem, the skeletonization consists of three steps: assembling the excitation matrix Y and finding skeleton beams, obtaining the skeletonized solution matrix Xs and recovering the complete solution matrix X. The skeletonization strategy had been employed to accelerate the SIE solution in [33] and the radar cross section evaluation in [34]. Here, we extend it to handle the optimization of the nano-optical trapping forces. It is worthy noting that the actual dimension of Y is NS × Ninc since y has only NS non-zero entries.
4. Numerical simulations
The numerical simulations are performed on an IBM server configured by X5650 CPU. GM-RES iteration process is terminated when the L2-norm of the residual vector is reduced to 10−4. κ is sampled uniformly. The Davis-Barton fifth-order approximation electromagnetic fields [37] are used. All light beams are vertical polarized. In the following simulations, the system with multiple excitations is always generated through illuminating the object by the same beam with different configurations: orientation or/and position. All the constitutive parameters of metal are taken from [38]. The wavelength in the background medium is denoted by λ.
4.1. Convergence and Accuracy
Convergence and accuracy analysis is carried out by numerical simulations on the Au sphere with radius 74.712 nm where the Mie analytical solution exists. The analysis includes both cases where the sphere is on-resonance and out of resonance. Since the sphere is resonant around 548.6 nm, we conduct the computations at 6 frequency points, where the wavelength is, respectively, 500.0 nm, 547.6 nm, 548.1 nm, 548.6 nm, 549.1 nm and 549.6 nm. The computation of the 500.0nm case is out of resonance. The comparison study in [16] revealed that different SIE formulations behaved distinctly on iterative convergence as well as accuracy. Particularly, SIE of the Poggio-Miller-Chang-Harrington-Wu (PMCHW) formulation suffers from the slow convergence but would always outperform its competitors on accuracy if convergence is reached. In contrast, SIE of the electric and magnetic current combined field integral equation (JMCFIE) formulation always converges faster than other SIE formulations. Interested readers can find the detailed discussions on JMCFIE and PMCHW in [16]. Here, JMCFIE is employed for SIE simulation because the associated PMHCW computation dose not converge.
The scattering results obtained by both FE-BI and JMCFIE are compared with the Mie solution. All FE-BI computations use the mesh of 91386 tetrahedrons with 5098 boundary triangle patches between Au and the background air. The boundary patches are employed by BI for the FE-BI computation. The same set of triangle patches is either utilized by JMCFIE for the SIE computation. Table 1 lists the statistics of computations. As shown by the table, FE-BI converges faster than JMCFIE in all cases. The accuracy of FE-BI is good for both resonance and non-resonance cases. It should be noted that no preconditioner is employed for JMCFIE and neither the preconditioning strategy depicted in subsection 3.2 is turned on.
Table 1. Statistics of computations on the Au sphere with radius 74.712 nm. The error is computed by , where EMie and Ec is the far fields corresponding to Mie’s series and each computation. L, the number of samplings along θ, is equal to 181 for the angle range of 0°–180°.
An Au brick, as shown in Fig. 2, is employed to study the accuracy convergence of the FE-BI in terms of mesh density. The brick is illuminated by the plane wave of the wavelength λ = 1064 nm at the orientation (θ = 135°,ϕ = 270°). The size of the brick is (0.4, 0.2, 0.6) λ. The center of the block is at the origin of the coordinate. Four sets of meshes are employed, where the average mesh size is, respectively, 0.16, 0.12, 0.05 and 0.03 dielectric wavelength. As it is known, the sharp edges always pose a challenge to the simulation, especially when a high near field accuracy is required. To investigate the accuracy on the near field, we present the distribution of the equivalent surface electric current on the brick in Fig. 2. Due to the symmetry of the geometry, the current should be symmetric with respect to the Y-Z plane. The data in Fig. 2 indicates that with the accuracy improves when the mesh becomes dense. The FE-BI seems to require a little dense mesh than SIE, in comparison with the study in [16] and our numerical simulations not presented here.
Fig. 2 Distribution of equivalent surface electric current on the brick when illuminated by the plane wave with the wavelength λ = 1024 nm at the orientation (θ = 135°,ϕ = 270°). The size of the brick is (0.4, 0.2, 0.6) λ. The center of the block is at the origin of the coordinate. The magnitude of the current is scaled for a nice view. (a): the full view; (b)-(e): viewed from the −y direction for the cases with different meshes.
4.2. Simulation of Trapping Forces
4.2.1. A single Au particle
An Au particle as shown in Fig. 3a is employed to show the accuracy and efficiency of the proposed FE-BI scheme in predicting optical trapping forces. The geometry center of the particle is at the origin of the coordinate. The wavelength and the waist of the incident beam is 1064 nm and 1.029 λ. Without the loss of generality, the mesh is generated in such a way that the geometry fidelity is remained and the common discretization requirement of BI is satisfied. Specifically, the target is meshed by 27899 tetrahedrons with 3382 boundary patches between Au and air. The beam propagates along −y direction. In the computation, the beam center moves along y-axis from −21280 nm to 21280 nm at the step of 106.4 nm. A system with 401 excitation beams is produced. Figure 3a presents the iteration history of the FE-BI, JMCFIE and PMCHW when the center of the incident beam is at (0, 0, 0). Figure 3b gives the variation of the y-component of the RPF when the beam center moving.
Fig. 3 The iteration history and RPF data for the Au particle simulation. The wavelength (λ) and the waist (wo) of the incident beam is 1064 nm and 1.029 λ. The beam propagates along −y direction. (a): iteration history when the bean center is at (0, 0, 0). The residual, an indicator of the error, is evaluated by computing the L2-norm of the residual vector in GMRES. (b): The y-component of the RPF. ”Fy-FE-BI” and ”Fy-no-skel”, respectively, denotes the data obtained by FE-BI with and without skeletonization. ”Fy-JMCFIE” represents the data obtained by JMCFIE. The beam center moves along y-axis from -21280 nm to 21280 nm at the step of 106.4 nm. The unit of force is 10−9N/W.
To reach a fair comparison, the SIE computations employ the same surface mesh as that employed by the BI in FE-BI computations. No preconditioner is employed by JMCFIE and PMCHW. Neither the preconditioning strategy depicted in subsection 3.2 is employed. It is indicated that FE-BI converges in about 100 steps while JMCFIE requires over 700 steps to reach convergence. PMCHW fails to converge in 2000 steps.
Due to the symmetry of the particle, it is expected that the y-component dominates the net RPF. Figure 3b compares the y-component of the RPF obtained by FE-BI and JMCFIE. The RPF data with respect to PMCHW is not available due to the failure of the convergence. As shown in the figure, two results agree well. Table 2 details the computational resources. Although FEM requires volume meshing, the total memory requirement by the FE-BI is much less than that consumed by SIE because the FEM matrix is sparse. Since FE-BI solution converges much faster than JMCFIE, FE-BI is efficient than SIE in memory usage as well as CPU time. The skeletonization scheme works well in reducing the total solution time. In this computation, only 3 skeleton Gaussian beams are selected. The total iteration time is reduced by a factor of over 130. Figure 3b indicates that the difference between the data obtained by FE-BI with and without skeletonization is indistinguishable.
Table 2. Statistics of computations on the Au particle as shown Fig. 3a. Iteration counts and time are averaged for 3 skeleton incident beams.
4.2.2. Three Au spheres
The metal plasmonic effects can greatly enhance the optical trapping forces. To show the capability of the proposed method on analyzing the plasmonic effects, the application with multiple nanostructures is simulated. As shown in Fig. 4, the system consists of three identical Au spheres of the diameter 106.4nm. The spheres are, respectively, placed at (d,0,0), (0, 0, 0) and (−d,0,0). By varying d, we compute the RPF exerted on the sphere located at (0, 0, 0) by the Gaussian beam with the wavelength of 532nm. The waist of the beam is fixed at 1.029λ. The light illuminates the spheres along −z direction. In the computation, the beam center varies along z-axis from −10640nm to 10640nm at the step of 106.4nm. A system with 201 incident beams is produced. Figure 4 gives the dominant component, i.e., z-component, of the RPF when the beam center varying. To indicate the plasmonic effects clearly, we calculate the RPF exerted on the system where only one single sphere is at (0, 0, 0). It is denoted by the one-sphere case in Fig. 4.
Fig. 4 The variation of the z-component of the RPF exerted on the Au sphere located at (0, 0, 0) when the center of the vertical polarized incident Gaussian beam varying. The wavelength and waist of the Gaussian is, respectively, 532 nm and 1.029 λ. The beam propagates along −z direction. The unit of force is 10−9N/W.
Fig. 5 The field distribution on the X-Z plane for different d’s when the center of the vertical polarized Gaussian beam is at (0, 0, 0). (a): the one-sphere case; (b-d): the the three-sphere cases with different d’s. The wavelength and waist of the Gaussian is, respectively, 532 nm and 1.029 λ. The beam propagates along −z direction. The unit of the electric field is 10−9V/m.
It is found that RPF fluctuates markedly with d. For all cases, z-component of RPF reaches its peak when the beam center arrives at z = 0. The RPF is almost identical to that of the one-sphere case when d > 532nm or d =218nm. Since the waist of the Gaussian beam is 1.029 λ, only the sphere at (0, 0, 0) is effectively illuminated by the incident light. This explains why RPF becomes identical to that of the one-sphere case when d > 532nm. The RPF exerted on the sphere at (0, 0, 0) is weaker than that in the one-sphere case when 218 nm < d < 532nm. Two facts, namely, the energy and momentum conservation and the interactions among spheres, contribute to the phenomenon. Mutual interactions among spheres would change the behavior of the forces. Due to the loss of Au at the working frequency, part of the energy is absorbed by the particles. Consequently, less energy acts on the sphere at (0, 0, 0). If we replace Au with the lossless dielectric, the force would be monotonically enhanced by the mutual interactions with the decrease of d. When d < 218nm, the RPF is significantly enhanced by the strong near fields resulting from the plasmonic effects. For the case of d = 125nm, the peak RPF becomes about 1.5 times larger than that of the one-sphere case. Figure 5 gives the field distribution on the X-Z plane when the center of the Gaussian beam is at (0, 0, 0). As shown in the figure, near fields become stronger when d decreases. Field enhancement and hot spots can be obviously observed in the d = 125nm case.
4.2.3. Many nanostructures
In system design, beam configurations should be optimized. To demonstrate the capability of the proposed method in handling the optimization process, a system consists of many nanoparticles are considered. As shown in Fig. 6, an Au sphere is placed above an Au rod array. The center of the sphere is at the origin of the coordinate. The wavelength and waist of the Gaussian is, respectively, 879.98 nm and 1.029 λ. The center of the Gaussian beam moves from −5279.88 nm to 5279.88 nm along x-axis at the step of 105.6 nm. Simultaneously, the orientation of the incident light varies from 0° to 90° along θ with the step of 1° while ϕ = 180°. The total number of incident beams reaches over 9000. Figure 7 presents the variation of the x-component of the RPF exerted on the sphere when the incident direction and the center of Gaussian beam varying simultaneously. To show the fluctuation clearly, the figure only presents a small part of the data. The result indicates that both the orientation and the center of the beam can influence the trapping forces greatly.
Fig. 7 Variation of the x-component of the RPF exerted on the Au sphere (see Fig. 6) when the orientation and center of the vertical polarized Gaussian beam varying simultaneously. The wavelength and waist of the Gaussian is, respectively, 879.88 nm and 1.029 λ. The unit of force is 10−9N/W.
15 skeleton beams are selected from over 9000 beams. 190 steps are averagely required by the iterative solution without the preconditioning strategy in subsection 3.2, whereas less than 80 steps are needed for the convergence with the preconditioner. As a result, each skeleton beam can be handled in about 2 minutes. The solution of the equivalent current can be finished in about 30 minutes.
5. Conclusion
The optical trapping forces exerted on metal nano-system with complex configuration are simulated by the hybrid of finite element and boundary integral (FE-BI) method. The convergence difficulty in surface integral equation (SIE) is greatly alleviated. The skeletonization scheme performs well in improving the simulation efficiency for the optimizing process with many beam configurations: orientation, positions, beam waist radius, etc. Numerical results indicate that the proposed method can accurately and efficiently predict the plasmonic effects as well as the optical trapping forces. The proposed method offers a viable alternative to cope with complex metallic nanostructures and has clear advantages such as the fast convergence and reduced computational demand.
Acknowledgments
This work was partly supported by Program for New Century Excellent Talents in University under Grant NCET-12-0045, by the 111 Project of China under grant B14010, by the NSFC under grant 61421001 and 61371002.
References and links
1. K. Svoboda and S. M. Block, “Optical trapping of metallic rayleigh particles,” Opt. Lett. 19, 930–932 (1994). [CrossRef] [PubMed]
2. O. M. Marago, P. H. Jones, P. G. Gucciardi, G. Volpe, and A. C. Ferrari, “Optical trapping and manipulation of nanostructures,” Nature Nanotech. 8, 807–819 (2013). [CrossRef]
3. K. F. Ren, G. Grehan, and G. Gouesbet, “Prediction of reverse radiation pressure by generalized lorenz-mie theory,” Appl. Opt. 35, 2702–2710 (1996). [CrossRef] [PubMed]
4. F. Xu, K.-F. Ren, G. Gouesbet, X.-S. Cai, and G. Grehan, “Theoretical prediction of radiation pressure force exerted on a spheroid by an arbitrarily shaped beam,” Phys. Rev. E 75, 026613 (2007). [CrossRef]
5. Stephen K. Gray and T. Kupka, “Propagation of light in metallic nanowire arrays: Finite-difference time-domain studies of silver cylinders,” Phys. Rev. B 68, 045,415 (2003). [CrossRef]
6. M. I. Mishchenko, “Radiation force caused by scattering, absorption, and emission of light by nonspherical particles,” J. Quant. Spectrosc. Radiat. Transfer 70, 811–816 (2001). [CrossRef]
7. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). [CrossRef]
8. E. Moreno, D. Erni, C. Hafner, and R. Vahldieck, “Multiple multipole method with automatic multipole setting applied to the simulation of surface plasmons in metallic nanostructures,” J. Opt. Soc. Am. A 19, 101–111 (2002). [CrossRef]
9. W.-B. Ewe, H.-S. Chu, and E.-P. Li, “Volume integral equation analysis of surface plasmon resonance of nanoparticles,” Opt. Express 15, 18200–18208 (2007). [CrossRef] [PubMed]
10. G. Miano, G. Rubinacci, and A. Tamburrino, “Numerical modeling for the analysis of plasmon oscillations in metallic nanoparticles,” IEEE Trans. Antennas Propag. 58, 2920–2933 (2010). [CrossRef]
11. X. M. Pan, W. C. Pi, M. L. Yang, Z. Peng, and X. Q. Sheng, “Solving problems with over one billion unknowns by the MLFMA,” IEEE Trans. Antennas Propag. 60, 2571–2574 (2012). [CrossRef]
12. K. Sendur, “An integral equation based numerical solution for nanoparticles illuminated with collimated and focused light,” Opt. Express 17, 7419–7430 (2009). [CrossRef] [PubMed]
13. A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3d simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A 26, 732–740 (2009). [CrossRef]
14. C. Forestiere, G. Iadarola, G. Rubinacci, A. Tamburrino, L. D. Negro, and G. Miano, “Surface integral formulations for the design of plasmonic nanostructures,” J. Opt. Soc. Am. A 29, 2314–2327 (2012). [CrossRef]
15. G. Benjamin and O. J. Martin, “Scattering on plasmonic nanostructures arrays modeled with a surface integral formulation,” Photonics and Nanostructures - Fundamentals and Applications 8, 278–284 (2010). [CrossRef]
16. M. G. Araujo, J. M. Taboada, D. M. Solis, J. Rivero, L. Landesa, and F. Obelleiro, “Comparison of surface integral equation formulations for electromagnetic analysis of plasmonic nanoscatterers,” Opt. Express 20, 9161–9171 (2012). [CrossRef] [PubMed]
17. O. Ergul and L. Gurel, “Fast and accurate analysis of large-scale composite structures with the parallel multilevel fast multipole algorithm,” J. Opt. Soc. Am. A 30, 509–517 (2013). [CrossRef]
18. X. M. Pan, L. Cai, and X. Q. Sheng, “An efficient high order multilevel fast multipole algorithm for electromagnetic scattering analysis,” Progr. Electromagn. Res. 126, 85–100 (2012). [CrossRef]
19. X. M. Pan and X. Q. Sheng, “Preconditioning technique in the interpolative decomposition multilevel fast multipole algorithm,” IEEE Trans. Antennas Propag. 61, 3373–3377 (2013). [CrossRef]
20. D. M. Solis, J. M. Taboada, F. Obelleiro, L. M. Liz-Marzan, and F. J. Garcia de Abajo, “Toward ultimate nanoplasmonics modeling,” ACS Nano 8, 7559–7570 (2014). [CrossRef] [PubMed]
21. A. Ji, T. V. Raziman, J. Butet, R. P. Sharma, and O. J. F. Martin, “Optical forces and torques on realistic plasmonic nanostructures: a surface integral approach,” Opt. Lett. 39, 4699–4702 (2014). [CrossRef] [PubMed]
22. L. Landesa, M. G. Araujo, J. M. Taboada, L. Bote, and F. Obelleiro, “Improving condition number and convergence of the surface integral-equation method of moments for penetrable bodies,” Opt. Express 20, 17237–17249 (2012). [CrossRef]
23. J. M. Jin, The Finite Element Method in Electromagnetics (Wiley-IEEE Press, New York, 2002), 2nd ed.
24. G. K. Christopher, J. N. Stephen, and V.-D. Tuan, “Investigating the plasmonics of a dipole-excited silver nanoshell: Mie theory versus finite element method,” Nanotechnology 21, 315203 (2010). [CrossRef]
25. W. C. Chew, J. M. Jin, E. Michielssen, and J. Song, Fast Efficient Algorithms in Computational Electromagnetics (Artech House, Boston, MA, 2001).
26. X. Q. Sheng, J. M. Jin, J. Song, C. C. Lu, and W. C. Chew, “On the formulation of hybrid finite-element and boundary-integral method for 3d scattering,” IEEE Trans. Antennas Propag. 46, 303–311 (1998). [CrossRef]
27. Y. Ji and T. H. Hubing, “Emap5: A 3d hybrid FEM/MoM code,” Appl. Computat. Electromagn. Soc. J. 15, 1–12 (2000).
28. J. M. Jin, Theory and Computation of Electromagnetic Fields (Wiley-IEEE Press, 2010). [CrossRef]
29. X. Q. Sheng and E. K.-N. Yung, “Implementation and experiments of a hybrid algorithm of the mlfma-enhanced FE-BI method for open-region inhomogeneous electromagnetic problems,” IEEE Trans. Antennas Propag. 50, 163–167 (2002). [CrossRef]
30. M. L. Yang, K. F. Ren, M. J. Gou, and X. Q. Sheng, “Computation of radiation pressure force on arbitrary shaped homogenous particles by multilevel fast multipole algorithm,” Opt. Lett. 38, 1784–1786 (2013). [CrossRef] [PubMed]
31. X. M. Pan and X. Q. Sheng, “Accurate and efficient evaluation of spatial electromagnetic responses of large scale targets,” IEEE Trans. Antennas Propag. 62, 4746–4753 (2014). [CrossRef]
32. E. Liberty, F. Woolfe, P. G. Martinsson, V. Rokhlin, and M. Tygert, “Randomized algorithms for the low-rank approximation of matrices,” Proc. Natl. Acad. Sci. USA 104, 20167–20172 (2007). [CrossRef] [PubMed]
33. X. M. Pan, M. J. Gou, and X. Q. Sheng, “Prediction of radiation pressure force exerted on moving particles by the two-level skeletonization,” Opt. Express 22, 10032–10045 (2014). [CrossRef] [PubMed]
34. H. W. Gao, J. W. Hao, X. M. Pan, and X. Q. Sheng, “Application of interpolative decomposition to FE-BI-MLFMA for fast computation of monostatic scattering from 3-d complex composite objects,” Antennas and Wireless Propagation Letters”, IEEE 13, 1490–1493 (2014).
35. X. M. Pan, J. G. Wei, Z. Peng, and X. Q. Sheng, “A fast algorithm for multiscale electromagnetic problems using interpolative decomposition and multilevel fast multipole algorithm,” Radio Sci. 47, RS1011 (2012). [CrossRef]
36. X. M. Pan and X. Q. Sheng, “Improved algebraic preconditioning for MoM solutions of large-scale electromagnetic problems,” IEEE Antennas Wireless Propag. Lett. 13, 106–109 (2014). [CrossRef]
37. J. P. Barton and D. R. Alexander, “Fifth-order corrected electromagnetic field components for a fundamental Gaussian beam,” J. Appl. Phys. 66, 2800–2802 (1989). [CrossRef]
38. P. B. Johnson and R. W. Christy, “Optical constants of the noble metals,” Phys. Rev. B 6, 4370–4379 (1972). [CrossRef]
