Prediction of radiation pressure force exerted on moving particles by the two-level skeletonization

Xiao-Min Pan; Ming-Jiang Gou; Xin-Qing Sheng

doi:10.1364/OE.22.010032

1. Introduction

The micro-manipulation [1] based on optical tweezers [2–4] offers many benefits over its competing technologies, most notably the ability to work with delicate or potentially harmful samples in sealed environments without direct mechanical access. Some of the more sophisticated optically actuated devices make use of interactions among multiple rotating objects. The rotation control and rotational interactions have not only been demonstrated in two dimensions but also in three dimensions [5]. Precise control of the motion of a particle in tweezers relies on the prediction of RPF exerted on the particle. Obviously, the more accurately the RPF can be predicted, the more reliable the complex manipulation can be reached. Therefore, an accurate and efficient computational tool to predict RPF exerted on moving particles is badly required to implement the optical tweezers devices with increasing complexity. From the simulation point of view, we can either simulate the moving particle in the fixed optical tweezers directly or change the state of the optical tweezers by varying parameters of the light beam while assuming the particle is static. To ease the implementation, the latter is employed in our study.

Different types of methods can be employed to predict RPF because the behavior of optical tweezers depends upon the size of the trapped particle relative to the wavelength of light used to trap it. In cases where the dimensions of the particle are much greater than the wavelength, a simple ray optics treatment is sufficient. If the wavelength of light far exceeds the particle dimensions, the particles can be treated as electric dipoles in an electric field. For optical trapping of dielectric objects of dimensions within an order of magnitude of the trapping beam wavelength, the only accurate models involve the treatment of either time dependent or time harmonic Maxwell equations using appropriate boundary conditions. These models can be generally grouped into analytical approaches or full-wave algorithms. Analytical methods are always limited to particles with canonical shapes while full-wave ones are always expensive. In [6–9], analytical algorithms focused on the spherical particles have been reported. It is a much more difficult task to predict the RPF exerted on nonspherical particles. Many efforts have been reported on this topic [10–16]. The sizes of the arbitrarily shaped particles in these methods are so severely limited that the beam profile effect on RPF is difficultly exploited. Most recently, multilevel fast multipole algorithm (MLFMA), one of the most widely used full-wave method in computational electromagnetic community, has been employed to predict RPF on arbitrarily shaped particles [17]. Suppose N is the number of unknowns, MLFMA can reduce both the computational complexity of O(N² ∼ N³) and the storage complexity of O(N²) for surface integral equation (SIE) to O(NlogN) [18–21]. Additionally, the use of SIE with triangular patches reduces considerably the unknowns since only the discretization on the particle surface is necessary. However, to model a moving particle in optical tweezers, a set of shaped light beams should be used. Numerically, discretizing different excitation beams gives rise to multiple right hand sides (RHS’s) for the matrix system. As each RHS must be solved separately, the time cost becomes O(N_bNlogN) with N_b the number of excitation beams. When N_b becomes large, the simulation can be prohibitively expensive. For multiple incidents electromagnetic scattering problem, some methods based on model-based parameter estimation (MBPE) [22], asymptotic waveform evaluation (AWE) [23], singular value decomposition (SVD)/QR factorization [24] or other techniques [25, 26] are developed. Nonetheless, their application is either limited by the lack of rigorous error control schemes or confined to the cases when N_b is moderate due to the efficiency issue.

Encouraged by the skeletonization framework in [27], a full-wave algorithm is proposed to predict RPF exerted on moving particles in this work. The algorithm is based on the interpolative decomposition (ID) [28–30] and skeleton concept [31–34]. The algorithm introduces skeleton beams to accelerate the MLFMA computation, which are selected by conducting ID skeletonization on the excitation matrix consisting of all RHS vectors. A bottle-neck in this process is that the peak memory usage of skeletonization may exceed what hardware can offer. By dividing the whole parameter space into serval intervals, the peak memory usage can be decreased. However, the efficiency of the computation would degrade because the number of skeletons would increase as the number of intervals grows. To solve this difficulty, a two-level skeletonization scheme is proposed.

The remainder of the paper is organized as follows. Section 2 outlines combined tangential formulation (CTF) of SIE and the method of moments (MoM) system for a moving particle in optical tweezers. Section 3 details the proposed method, including the skeletonization framework, the cost analysis and the two-level scheme of the algorithm. Section 4 presents numerical experiments to show performance of the proposed algorithm. Section 5 concludes the paper.

2. CTF and MoM

Figure 1 shows the typical time-harmonic scattering (with a time factor e^−jωt of the angular frequency ω) by the arbitrarily shaped homogeneous particle with permittivity ε_i, and permeability μ_i (i denotes internal) embedded in an homogeneous lossless medium with constitutive parameters ε_e; μ_e (e denotes external). Ω_i and Ω_e, respectively, denote the inner and outer regions, which are assumed to be open sets and are separated by a boundary surface S, with a unit normal n̂ (directed from S toward Ω_e). E^inc and H^inc are the impressed electric and magnetic fields in Ω_e. By invoking Love’s equivalent principle for the exterior medium, the region Ω_i is filled up with the same material of the region Ω_e, and the sources in Ω_i are removed. The equivalent currents $J_{e} = \hat{n} \times H_{e}^{tot}$ and $M_{e} = E_{e}^{tot} \times \hat{n}$ , positioned on the external face S_e of the surface S, generate the original scattered fields $E_{e}^{sca}$ and $H_{e}^{sca}$ in the region Ω_e, and null fields in the region Ω_i, i.e.,

E_{e}^{sca} (r, J_{e}, M_{e}) + E^{inc} = {\begin{array}{l} 0 & if r \in Ω_{i} \\ E_{e}^{tot}, & if r \in Ω_{e} \end{array}

H_{e}^{sca} (r, J_{e}, M_{e}) + H^{inc} = {\begin{array}{l} 0 & if r \in Ω_{i} \\ H_{e}^{tot}, & if r \in Ω_{e} \end{array}

An equivalent problem can be established for Ω_i as well, where Ω_e is filled up with the same material of Ω_i, and the sources in Ω_e are removed. The equivalent currents (

J_{i} = - \hat{n} \times H_{i}^{tot}

,

M_{i} = - E_{i}^{tot} \times \hat{n}

), defined on the external face S_i of the surface S, produce the original scattered fields

E_{i}^{sca}

and

H_{i}^{sca}

in Ω_i, and null fields in Ω_e:

E_{i}^{sca} (r, J_{i}, M_{i}) = {\begin{array}{l} E_{i}^{tot}, & if r \in Ω_{i} \\ 0, & if r \in Ω_{e} \end{array}

H_{i}^{sca} (r, J_{i}, M_{i}) = {\begin{array}{l} H_{i}^{tot}, & if r \in Ω_{i} \\ 0 & if r \in Ω_{e} \end{array}

The scattered electric and magnetic fields in Eqs. (1)–(4) can be evaluated as

E_{q}^{sca} (r, J_{q}, M_{q}) = η_{q} ℒ_{q} (J_{q}) - 𝒦_{q} (M_{q}),

H_{q}^{sca} (r, J_{q}, M_{q}) = η_{q}^{- 1} ℒ_{q} (M_{q}) + 𝒦_{q} (J_{q}),

where q = i, e; and

η_{q} = \sqrt{μ_{q} / ε_{q}}

is the the wave impedance for the region q. The operators ℒ and 𝒦 corresponding to the region q can be written as

ℒ_{q} {X} (r) = - j k_{q} \int_{S} d S^{'} [I + \frac{1}{k_{q}^{2}} \nabla \nabla^{'} \cdot] X (r^{'}) g_{q} (r, r^{'}),

𝒦_{q} {X} (r) = - \frac{Ω (r)}{4 π} X (r^{'}) + \nabla \times {\int ―}_{S} d S^{'} g_{q} (r, r^{'}) X (r^{'}),

where ∫ stands for the principal value integration; 0 ≤ Ω(r) ≤ 4π is the internal solid angle;

g_{q} (r, r^{'}) = \frac{e^{- j k_{q} | r - r^{'} |}}{4 π | r - r^{'} |}

, denotes the homogeneous-space Green’s function with

k_{q} = ω \sqrt{μ_{q} ε_{q}}

the wavenumber for the region q. It holds true from the boundary condition and the definition of equivalent currents that J_i = −J_e and M_i = −M_e. The system consisting of the four equations (1)–(4) is overdetermined for the two unknowns (J_e, M_e). To solve (J_e, M_e), different SIE formulations are generally constructed by combining tangential (T) and/or normal (N) forms of the four equations: tangential electric field integral equation (T-EFIE), tangential magnetic field integral equation (T-MFIE), normal electric field integral equation (N-EFIE) and normal magnetic field integral equation (N-MFIE). The CTF of SIE is the combination of [35, 36],

η_{e}^{- 1} {(T-EFIE)}_{e} + η_{i}^{- 1} {(T-EFIE)}_{i},

η_{e} {(T-MFIE)}_{e} + η_{i} {(T-MFIE)}_{i} .

where η_e/i is the wave impedance for the region e or i. For the numerical solution of SIEs, the surfaces are discretized by small planar triangles and (J_e, M_e) are expanded by Rao-Wilton-Glisson (RWG) [37] basis functions. Using a Galerkin scheme with the RWG testing functions, the CTF can be recast in the matrix form as

Z \cdot x = b, or (\begin{matrix} Z^{J, J} & Z^{J, M} \\ Z^{M, J} & Z^{M, M} \end{matrix}) \cdot {\begin{matrix} x^{J} \\ x^{M} \end{matrix}} = {\begin{matrix} b^{J} \\ b^{M} \end{matrix}},

where Z is the impedance matrix of size 2N × 2N with N the number of coefficients of J_e or M_e. b is the RHS vector obtained by discretizing the incident field. The detailed discretization procedure can be found in literature, such as [38, 39]. After the coefficient vector of effective currents x is obtained, the scattered fields E^sca and H^sca at any position and thus the associated RPF can be obtained by the manner described in [40].

Fig. 1: The red blood cell (RBC) model.

Download Full Size | PDF

To investigate a moving particle in optical tweezers, we equivalently fix the particle at the origin of the coordinate system while shifting and tilting the excitation light beams. Therefore, the RPF exerted on a moving particle can be studied by changing the parameters characterizing the beams. Without loss of generality, this work uses Gaussian light beams as the excitation. The parameters to describe a Gaussian beam include: θ, φ, r_o, w_o, where (θ, φ) denotes the beam direction, r_o is the beam center, w_o is the beam waist. Using MoM to study the optical tweezers generated by the beams within the parameter space of κ ∈ [κ_min, κ_max] (κ = θ, φ, r_o, w_o, or their combinations) amounts to solving a set of systems of linear algebraic equations that are of the matrix form as:

Z \cdot X (κ) = B (κ),

where B(κ) is an 2N × N_b known excitation matrix consisting of N_b RHS vectors, X(κ) is the unknown solution matrix. For simplicity, κ and (θ, φ, r_o, w_o) are mutually used or both dropped when no confusion would arise. If Z⁻¹ is available, (11) can be solved efficiently. However, when N is large, direct evaluation of Z⁻¹ is always unaffordable, both in memory usage and CPU time. The simulation is therefore conventionally limited to electrically small dielectric objects. With the aid of MLFMA, both computational and storage complexities are reduced to O(NlogN) for the iterative solution of each RHS. The total cost of obtaining the solution matrix X is thus O(N_bNlogN). For a prescribed resolution on κ, the exact value of N_b is hardly estimated in advance. People always use a large N_b empirically to satisfy the resolution requirement. Consequently, the computation may become expensive even for the small particle. The excitation matrix B is always rank deficient and admits low-rank decomposition, e.g., SVD, QR, to reduce the cost. However, SVD or QR is only efficient for small matrix due to the high computational complexity of the decomposition itself.

3. The proposed fast algorithm

3.1. The skeletonization framework

Most recently, an efficient rank revealing algorithm, interpolative decomposition (ID), has been developed in [28]. ID has been proved much more efficient than SVD or QR in conducting low-rank decomposition [29, 34]. Suppose C is a complex m × n matrix of rank r with r ≤ m and r ≤ n, ID states that there exists a complex m × r matrix S whose columns consist of a subset of the columns of C and a complex r × n matrix R such that

C^{m \times n} \approx S^{m \times r} \cdot R^{r \times n},

when the exact rank of C^m×n is greater than r, but the (r + 1)-st greatest singular value of C^m×n is small. Before conducting the decomposition as shown in (12), a threshold ε_ID is often prescribed to control the error of the approximation. In the form of L₂-norm, ε_ID measures the difference between C and S · R [28]. The column vectors in S are generally called as skeleton vectors [29]. By applying ID to B(κ), we can reach

B (κ) = B^{S} (κ) \cdot R (κ),

where B^S(κ) and R(κ) are, respectively, the N × N_skel skeletonized excitation matrix and the N_skel × N_b projection matrix, with N_skel the number of skeleton light beams. Substituting (13) into (11) yields

Z \cdot X (κ) = B^{S} (κ) \cdot R (κ) .

Suppose

X^{S} (κ) = Z^{- 1} \cdot B^{S} (κ),

where X^S(κ) is the solution matrix corresponding to B^S(κ). According to (14) and (15), X(κ) can be obtained as

X (κ) = X^{S} (κ) \cdot R (κ) .

Equation (16) states that the complete solution matrix X can be obtained after X^S is obtained.

Except for the error arising from the ID skeletonization, no approximation is introduced during the deduction from (11) to (16). Therefore, the accuracy of proposed algorithm can be strictly manipulated by ε_ID since ID is error controllable [28].

The direct implementation of ID in (12) requires [28]

C_{ID} = l \cdot C_{H} + O (r \cdot m + r \cdot l \cdot n)

floating-point operations, where C_H is the cost of applying C^H to a vector, the superscript “H” denotes transpose conjugate operation. As shown in [28], l = r + 5 or l = r + 10 is sufficient. From the theoretical point of view, r is a constant for a given physical system. If C is dense but rank deficient, the cost can be decreased to (mnlog(r) + r²(m + n)) by making use of randomness [28, 29]. Since the matrix to be decomposed in this work is always sufficiently rank deficient, r ≪ n. The time for the ID skeletonization is thus dominated by the term O(mn), which is negligible in comparison with the total solution time as shown in Section 4. Additionally, it increases only linearly with n. When applying ID to B, m is equal to 2N, n is identical to N_b and r is N_skel. The relationship is used by default throughout the paper if it is not pointed out specifically. The peak memory usage of ID (in byte) is about

M_{ID} = 16 (m \cdot n + 2 n \cdot m_{2} + 17 m) \approx 16 (2 m \cdot n + 17 m),

where m₂ is the greatest integer less than or equal to m, such that m₂ is a positive integer power of 2. To remain the high accuracy of ID, double precision is generally employed. Although the memory required by ID can be released after skeleton beams are figured out, the peak memory consumption can be a bottle-neck as it may exceed what the computer can offer when both N and N_b are large. To fit the memory limit, a two-level scheme is developed.

3.2. The two-level scheme

It is not suggested to decrease M_ID by reducing N_b at the risk of violating the resolution. Here, a two-level scheme is proposed as the remedy. It is based on the multi-interval variation of the skeletonization algorithm, where n = N_b beams are usually uniformly distributed into multiple intervals. Suppose N_intv is the number of intervals, the number of beams in the i-th interval is n_i = n/N_intv. According to (18), the peak memory usage of ID skeletonization decreases linearly with the growth of N_intv. Additionally, skeletonization is conducted independently for each interval, still in the error controllable fashion. Thus, the multi-interval strategy will give rise to no loss of accuracy. The total time used by the skeletonization (to say, O(mn)) does not increase with N_intv because $n = \sum_{i}^{N_{intv}} n_{i}$ . However, its side-effect is that the total number of skeleton beams, N_skel, would increase since skeletonization for each single interval can not access the global information in B. Obviously, many intervals would be generated if mn is large. The growth of N_skel will degrade the efficiency of the proposed method. The side-effect can be eliminated by the proposed two-level scheme.

The two-level scheme consists of two stages as shown in Algorithm 1 and 2. The former describes how to carry out the skeletonization, while the latter details the procedure to recover the complete solution matrix X. In both Algorithm 1 and 2, κ denotes one of the parameters (κ = θ, φ, r_o, w_o) or their combinations. Algorithm 1 begins with dividing the whole parameter space into several intervals. The criterion for the division is usually the memory limit. And then, a loop is performed for the i-th interval to find out the N_i,skel skeleton beams { $κ_{i}^{S}$ } within the associated parameter sub-space, and to construct the skeletonized RHS matrix $B_{i}^{S}$ and projection matrix R_i. Next, a new RHS matrix H is produced by collecting all $B_{i}^{S}$ . Through ID skeletonization on H, a set of skeleton beams {κ^S} are generated and stored in H^S. The skeletons {κ^S} are named as two-level skeletons because they are selected from { $κ_{i}^{S}$ }. Algorithm 2 starts from obtaining the skeletonized solution matrix X^S by the standard MLFMA computation (other fast methods are also applicable). Next, the solution matrix X^H corresponding to { $κ_{i}^{S}$ }(i = 1, 2,··· ,N_intv) is computed through the ID procedure. After the skeletonized solution matrix $X_{i}^{S}$ is extracted from X^H, it is used to recover X_i for the i-th interval. At last, the complete solution matrix is generated by assembling all X_i’s together.

Algorithm 1:. the two-level multi-interval strategy–skeletonization stage

View Table | View all tables in this article

Algorithm 2:. The two-level multi-interval strategy–recovering stage

View Table | View all tables in this article

4. Numerical experiments

The numerical experiments are performed on an IBM server configured by two 6-core X5650 CPUs and 64 gigabytes (GB) memory. In all computations, 12 OpenMP threads are forked to accelerate the MLFMA as one CPU has 6 cores. The generalized minimum residual (GMRES) iteration process is terminated when the L₂-norm of the residual vector is reduced to 10⁻⁴. ε_ID = 10⁻⁴ and κ is sampled uniformly. The Davis-Barton fifth-order approximation electromagnetic fields [40] are used. The wavelength in the background medium is noted by λ in the following. To investigate the error arising from skeletonization, we define the difference of RPF as $δ_{RPF}^{p} (κ) = \frac{| F_{bf}^{p} (κ) - F_{skel}^{p} (κ) |}{| F_{bf}^{p} (κ) |}$ , where $F_{bf}^{p}$ and $F_{skel}^{p}$ are, respectively, the RPF obtained from the brute-force approach and the proposed fast algorithm, and p denotes one of the Cartesian components. The speed-up is computed by $C_{spdp} = \frac{T_{bf}}{T_{skel}}$ , where T_skel and T_bf are, respectively, the time with and without the skeleton-based algorithm; if brute-force computation is not conducted, T_bf is estimated by assuming the time for the solution of each RHS identical. T_skel and T_bf do not include the time for computing RPF from (J_e, M_e). It should be noted that C_spdp heavily depends on N_b. In the computations, RPF and incident angle, respectively, is defined in terms of (x, y, z) and (ξ, γ, ζ) coordinate, as shown in Fig. 1 and 2. The two coordinates are such defined that x = −ξ, y = γ and z = −ζ.

Fig. 2: Comparison of z components of RPF by Gaussian beams (N_b = 401, λ = 514.5nm and w_o = 2λ) on a prolate slightly volatile silicone oil particle (ε = 2.25) computed with and without the skeletonization. The transversal radius of the particle is b = 1μm while the radius a along z axis varies. The step of ζ_o is 0.5μm. ”BF” denotes brute-force.

Download Full Size | PDF

4.1. Accuracy and efficiency

Accuracy and efficiency of the proposed two-level scheme is validated by numerical experiments on a set of oil micro-particle models. Specifically, spheroidal particles of slightly volatile silicone oil (relative permittivity ε = 2.25) illuminated by Gaussian beams are considered, as shown in Fig. 2. Being same as those in [17], the ratio of a/b for the three oil particles are, respectively, 1.00, 1.05 and 1.10, where b is fixed at 1μm. All the three particles are discretized by 30924 triangle patches, leading to 92772 unknowns. That is, m in (18) is 92772. The center of the particle is fixed at the coordinate origin and λ = 514.5nm. By setting the incident angle q = 0° and varying z_o from −100μm to 100μm with the step of 0.5μm (ξ_o = 0 and γ_o = 0), RPF’s are computed by the brute-force manner and by the proposed algorithm. In the brute-force computation, 401 RHS’s have to be solved separately. As a result, T_bf’s for the three cases are, respectively, 280.7, 401.0 and 441.1 minutes. The corresponding T_skel’s are 3.5, 5.0 and 5.5 minutes because only 5 skeleton beams are selected for all the three cases, with the C_spdp of about 80. To investigate the error caused by skeletonization, the relative errors are presented in Fig. 3. As it is shown, the maximum $δ_{RPF}^{z}$ is less than 2.0 × 10⁻⁴.

Fig. 3: The relative error of z components of RPF by Gaussian beams (N_b = 401, λ = 514.5nm and w_o = 2λ) on a prolate slightly volatile silicone oil particle (ε = 2.25) computed with and without the skeletonization. The transversal radius of the particle is b = 1μm while the radius a along z axis varies. The step of ζ_o is 0.5μm.

Download Full Size | PDF

M_ID reaches 13986 MB for the computations above because only one interval is used to find the skeletons. In the following, the performance of the two-level scheme is revealed by calculations on the oil particle with a/b = 1.10 through manually setting N_intv to be 2, 3, 4 and 5. The associated statistics is listed in table 1. It can be seen that M_ID decreases linearly with N_intv at the expense of the monotonic growth of $\sum_{i}^{N_{intv}} N_{i, skel}$ . Definitely, the efficiency would degrade in the original multi-interval computations. However, the degradation is avoided when the two-level scheme is applied because N_skel remains a constant. Theoretically, the number of skeleton beams is a constant if the excitation matrix B delivers sufficient resolution. When multiple intervals are employed to figure out skeletons, $N_{skel} = \sum_{i}^{N_{intv}} N_{i, skel}$ increases with the number of intervals because the global information is not available in each interval. Since these RHS’s are required to be solved separately, the efficiency of the multi-interval variation decreases with the increasing of N_intv. In the two-level scheme, N_skel becomes independent of N_intv and remains a constant because the global information is successfully recovered during the second level skeletonization. The additional cost of the two-level skeletonization scheme is that the time used by the second level skeletonization. It is always negligible in comparison with those for the N_intv intervals because $\sum_{i}^{N_{intv}} N_{i, skel}$ is much less than N_b, as shown by the statistics. Consequently, the cost of the two-level strategy is insensitive to N_intv and the decrease of efficiency is avoided. The two-level scheme is still error controllable since the ID skeletonization is well manipulated. As shown in Fig. 4, the relative errors for most of the RHS’s are less than 10⁻⁴. The maximum $δ_{RPF}^{z}$ is only about 8.0×10⁻⁴. The two-level scheme makes the choice of N_b easy. The exact N_b is hardly estimated precisely for a prescribed resolution. In practice, it is chosen heuristically. If it is too small, the resolution requirement can hardly be satisfied. In contrast, the efficiency would become a problematic issue if it is unnecessarily large. This problem is not completely overcome by the multi-interval strategy because $N_{skel} = \sum_{i}^{N_{intv}} N_{i, skel}$ grows with N_intv. However, it is substantially solved because N_skel is independent of N_b as well as N_intv. To validate the analysis, we increase N_b to 801 for the computation on the particle with a/b = 1.10. The brute-force computation requires 884.3 minutes to solve all the RHS’s while the computation with the two-level scheme remains 5.5 minutes, the same as the case with N_b = 401, since N_skel is still 5. The accuracy is either under control as shown in Fig. 5.

Table 1:. Statistics of computations on the spheroid oil particle with a/b = 1.10 when different number of intervals are employed in the two-level scheme.

View Table | View all tables in this article

Fig. 4: The relative error of z components of RPF when different number of intervals are employed to figure out the skeleton beams (b = 1.00μm and a = 1.10μm; the other parameters are the same as those for the computations in Fig. 2).

Download Full Size | PDF

Fig. 5: z components of RPF (N_b = 801, the parameters are the same as those for Fig. 2).

Download Full Size | PDF

4.2. Capability

As discussed in Section 1, the problem of evaluating RPF exerted on a moving particle by a single excitation beam can be converted into that of computing RPF’s exerted on a static particle by multiple beams. The latter is much easier to be implemented than the former. With the conversion, numerical experiments on a moving red blood cell (RBC) model in water is conducted to demonstrate the capability of the proposed algorithm, as shown in Fig. 1. An ordinary RBC is always modeled by a rotationally-symmetric biconcave surface defined similar to [20, 41],

r (Θ, Φ) = a {sin}^{s} Θ + b

where s = 5, a = 3.30μm, and b = 0.55μm. Different from the model depicted by (19), the dimensions along x and y of our model are, respectively, scaled by the factors of 1.8 and 0.5. Additionally, the angle Θ_c is equal to 10°, as shown in Fig. 1. The volume of the model remains 80.9μm³. The relative permittivity ε is 1.3 for the water and 1.4 for the RBC model. For all Gaussian beams, λ = 532nm, w_o = 2λ, ξ_o = 0 and γ_o = 0. Multiple Gaussian beams are produced by varying θ, φ and ζ_o. To be specific, the ranges of θ and φ are [45°, 55°] and [85°, 95°]. The angle step is 0.2° along both θ and φ. ζ_o varies within the range of [−2λ, 2λ] with the step of 0.1λ. As a result, 106 641 = 51 × 51 × 41 RHS’s are involved in the computation. For such a large target, iterative solution of each RHS costs about half an hour for the MLFMA. Consequently, the brute-force computation requires more than 53 320 hours to complete the iteration. However, the time can be substantially decreased by the skeletonization algorithm because N_skel is only 63. The time used to figure out the skeletons is only about 12 hours for such a challengeable problem. Figure 6 and 7 presents the RPF varies with φ and ζ_o when θ is fixed at 45°, 50° and 55°. In the figure, the unit of force is 10⁻⁹N/W.

Fig. 6: x and y components of RPF (denoted by F_x and F_y) exerted on the RBC model when θ, φ and ζ_o of the beam vary within [45°, 55°], [85°, 95°] and [−2λ, 2λ]. The step of angle and distance is, respectively, 0.2° and 0.1λ. The unit of force is 10⁻⁹N/W.

Download Full Size | PDF

Fig. 7: z components of RPF (F_z) exerted on the RBC model (The other parameters are the same as those in 6).

Download Full Size | PDF

Due to the complex shape of the particle, x, y and z components of RPF behave quite differently. It is interesting to see that strong inverse RPF appears within some particular attitudes. With both positive and negative RPF’s, the complicated motion/rotation of the particle can be realized and manipulated. Theoretically, positive force is caused by reflected fields while negative force results from refracted fields. Inverse RPF is obtained if negative force is stronger than positive one [9,42]. However, the property of the fields and thus the RPF is dependent on many factors including the constitutive parameters of the particle and the background, the shape of the particle, the profile of the incident beam, the relative position between the particle and the beam, etc. Detailed analysis and explanation on RPF for moving particles with complex shapes require a large volume of numerical experiments. Indeed, it is one of our undergoing research work.

5. Conclusion

The proposed algorithm conducts low-rank decomposition on the excitation matrix consisting of all RHS’s to figure out the so-called skeleton light beams by interpolative decomposition (ID). A two-level skeletonization scheme is proposed to overcome the bottleneck associated with the peak memory usage during the ID skeletonization. Numerical experiments show that the proposed algorithm is efficient and error controllable. Although not demonstrated here, the performance of the algorithm is not sensitive to shape of moving particles or the excitation light beams because the ID skeletonization is a pure algebraic tool.

Acknowledgments

This work was partly supported by Program for New Century Excellent Talents in University under Grant NCET-12-0045, by the 973 Program under Grant 2012CB720702, by the NSFC under Grants 60901005 and 61371002, and by the Excellent Scholars Support Fund of Beijing under Grant 2012D009011000002.

References and links

1. K. C. Neuman and S. M. Block, “Optical trapping,” Rev. Sci. Instrum. 75, 2787–2809 (2004). [CrossRef]

2. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, “Observation of a single-beam gradient force optical trap for dielectric particles,” Opt. Lett. 11, 288–290 (1986). [CrossRef] [PubMed]

3. A. Ashkin, J. M. Dziedzic, and T. Yamane, “Optical trapping and manipulation of single cells using infrared-laser beams,” Nature 330, 769–771 (1987). [CrossRef] [PubMed]

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

5. H. Shpaisman, D. B. Ruffner, and D. G. Grier, “Light-driven three-dimensional rotational motion of dandelion-shaped microparticles,” Appl. Phys. Lett. 102, 071103 (2013). [CrossRef]

6. G. Roosen and C. Imbert, “Optical levitation by means of two horizontal laser beams: A theoretical and experimental study,” Phys. Lett. A 59, 6–8 (1976). [CrossRef]

7. J. S. Kim and S. S. Lee, “Radiation pressure on a dielectric sphere in a gaussian laser beam,” Opt. Acta 29, 801–806 (1982). [CrossRef]

8. K. F. Ren, G. Greha, and G. Gouesbet, “Radiation pressure forces exerted on a particle arbitrarily located in a gaussian beam by using the generalized lorenz-mie theory, and associated resonance effects,” Opt. Commun. 108, 343–354 (1994). [CrossRef]

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

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

11. B. T. Draine and P. J. Flatau, “Discrete-dipole approximation for scattering calculations,” J. Opt. Soc. Am. A 11, 1491–1499 (1994). [CrossRef]

12. S. H. Simpson and S. Hanna, “Computational study of the optical trapping of ellipsoidal particles,” Phys. Rev. A 84, 053808 (2011). [CrossRef]

13. C.-F. Kuo and S.-C. Chu, “Numerical study of the properties of optical vortex array laser tweezers,” Opt. Express 21, 26418–26431 (2013). [CrossRef] [PubMed]

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

15. F. Borghese, P. Denti, R. Saija, and A. l. Maria, “Optical trapping of nonspherical particles in the t-matrix formalism,” Opt. Express 15, 11984–11998 (2007). [CrossRef] [PubMed]

16. L. Bi and P. Yang, “Modeling of light scattering by biconcave and deformed red blood cells with the invariant imbedding T-matrix method,” J. Bio. Opt. 18, 055001 (2013). [CrossRef]

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

18. R. Coifman, V. Rokhlin, and S. Wandzura, “The fast multipole method for the wave equation: a pedestrian prescription,” IEEE Antennas Propag. Mag. 35, 7–12 (1993). [CrossRef]

19. X. M. Pan, W. 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]

20. O. Ergul, A. Arslan-Ergul, and L. Gurel, “Computational study of scattering from healthy and diseased red blood cells,” J. Bio. Opt. 15, 045004(2010). [CrossRef]

21. X. Q. Sheng, J. M. Jin, J. Song, W. C. Chew, and C. C. Lu, “Solution of combined-field integral equation using multilevel fast multipole algorithm for scattering by homogeneous bodies,” IEEE Trans. Antennas Propag. 46, 1718–1726 (1998). [CrossRef]

22. X. Wang and D. H. Werner, “Improved model-based parameter estimation approach for accelerated periodic method of moments solutions with application to the analysis of convoluted frequency selected surfaces and metamaterials,” IEEE Trans. Antennas Propag. 58, 122–131 (2010). [CrossRef]

23. B.-Y. Wu and X. Q. Sheng, “Application of asymptotic waveform evaluation to hybrid FE-BI-MLFMA for fast RCS computation over a frequency band,” IEEE Trans. Antennas Propag. 61, 2597–2604 (2013). [CrossRef]

24. A. Schroder, H. D. Bruxns, and C. Schuster, “A hybrid approach for rapid computation of two-dimensional monostatic radar cross section problems with the multilevel fast multipole algorithm,” IEEE Trans. Antennas Propag. 60, 6058–6061 (2012). [CrossRef]

25. P. D. Ledger and K. Morgan, “An adjoint enhanced reduced-order model for monostatic RCS computation,” Electromagnetics 28, 54–76 (2008). [CrossRef]

26. P. Zhen, M. B. Stephanson, and J. F. Lee, “Fast computation of angular responses of large-scale three-dimensional electromagnetic wave scattering,” IEEE Trans. Antennas Propag. 58, 3004–3012 (2010). [CrossRef]

27. X. M. Pan and X. Q. Sheng, “Fast computation of two-dimensional spatial electromagnetic scattering from large-scale targets,” Computational Electromagnetics Workshop (CEM), 2013 pp. 1–3 (2013). [CrossRef]

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

29. N. Halko, P. G. Martinsson, and J. A. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM Rev. 53, 72 (2011). [CrossRef]

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

31. K. L. Ho and L. Greengard, “A fast direct solver for structured linear systems by recursive skeletonization,” SIAM J. Sci. Comput. 34, A2507–A2532 (2012). [CrossRef]

32. X. M. Pan and X. Q. Sheng, “Hierarchical interpolative decomposition multilevel fast multipole algorithm for dynamic electromagnetic simulations,” Progr. Electromagn. Res. 134, 79–94 (2013). [CrossRef]

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

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

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

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

37. S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” IEEE Trans. Antennas Propag. 30, 409–418 (1982). [CrossRef]

38. P. Yla-Oijala, M. Taskinen, and S. Jarvenpaa, “Surface integral equation formulations for solving electromagnetic scattering problems with iterative methods,” Radio Sci. 40, RS6002 (2005). [CrossRef]

39. O. Ergul and L. Gurel, “Comparison of integral-equation formulations for the fast and accurate solution of scattering problems involving dielectric objects with the multilevel fast multipole algorithm,” IEEE Trans. Antennas Propag. 57, 176–187 (2009). [CrossRef]

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

41. J. Q. Lu, P. Yang, and X.-H. Hu, “Simulations of light scattering from a biconcave red blood cell using the finite-difference time-domain method,” J. Bio. Opt. 10, 024022 (2005). [CrossRef]

42. T. C. B. Schut, G. Hesselink, B. G. De Grooth, and J. Greve, “Experimental and theoretical investigations on the validity of the geometrical optics model for calculating the stability of optical traps,” Cytometry 12, 479–485 (1991). [CrossRef] [PubMed]

Prediction of radiation pressure force exerted on moving particles by the two-level skeletonization

Abstract

1. Introduction

2. CTF and MoM

3. The proposed fast algorithm

3.1. The skeletonization framework

3.2. The two-level scheme

4. Numerical experiments

4.1. Accuracy and efficiency

4.2. Capability

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Tables (3)

Equations (20)

Optics Express

N _intv	1	2	3	4	5
M_ID(MB)	13986	6993	4462	3497	2797
$\sum_{i}^{N_{intv}} N_{i, skel}$	5	8	10	12	15
C_ID(sec)	44	44	44	45	45
N _skel	5	5	5	5	5
T_skel(min)	5.5	5.5	5.5	5.5	5.5