August 15, 2014 / Vol. 39, No. 16 / OPTICS LETTERS

4827

Comparison of T-matrix calculation methods for scattering by cylinders in optical tweezers Xiaoqiong Qi,1,2,* Timo A. Nieminen,1 Alexander B. Stilgoe,1 Vincent L. Y. Loke,1 and Halina Rubinsztein-Dunlop1 1

2

Quantum Science Laboratory, School of Mathematics and Physics, The University of Queensland, Queensland 4072, Australia State Key Laboratory on Integrated Optoelectronics, Institute of Semiconductors, Chinese Academy of Sciences, Beijing 100083, China *Corresponding author: [email protected] Received May 2, 2014; revised July 1, 2014; accepted July 3, 2014; posted July 7, 2014 (Doc. ID 211318); published August 12, 2014 The T-matrix method, or the T-matrix formulation of scattering, is a framework for mathematically describing the scattering properties of an object as a linear relationship between expansion coefficients of the incident and scattering fields in a basis of vector spherical wave functions (VSWFs). A variety of methods can be used to calculate the Tmatrix. We explore the applicability of the extended boundary condition method (EBCM) and point matching (PM) method to calculate the T-matrix for scattering by cylinders in optical tweezers and hence the optical force acting on them. We compare both methods with the discrete-dipole approximation (DDA) to measure their accuracy for different sizes and aspect ratios (ARs) for Rayleigh and wavelength-size cylinders. We determine range of sizes and ARs giving errors below 1% and 10%. These results can help researchers choose the most efficient method to calculate the T-matrix for nonspherical particles with acceptable accuracy. © 2014 Optical Society of America OCIS codes: (350.4855) Optical tweezers or optical manipulation; (290.5825) Scattering theory; (290.5850) Scattering, particles. http://dx.doi.org/10.1364/OL.39.004827

In optical tweezers, a high numerical aperture lens is used to tightly focus a laser beam; transparent dielectric particles can then be trapped at the focus [1]. The nearinfrared wavelength region of 800–1100 nm is commonly used since the light is poorly absorbed by most living matter [2], and the particles are typically of sizes of the order of the trapping wavelength [3]. It is often very useful to be able to computationally model optical trapping [4]. Unfortunately, the particles at such sizes are on one hand too large for the Rayleigh approximation to be valid, and on the other, too small for the geometric optics approximation to be applied. Thus, a full electromagnetic wavescattering calculation is required. The T-matrix method is a useful method for calculating scattering by particles of sizes on the order of the illuminating wavelength [5–9]. For spherical particles, the T-matrix is given by the analytical Lorenz–Mie solution for scattering by a sphere [10–13]. However, there is much practical interest in the trapping of nonspherical particles, for example, for optically trapped scanning-probe microscopy [14]. In this case, it is necessary to resort to numerical methods to calculate the T-matrix [15–17]. As a general description of the electromagnetic scattering by a particle, many methods could be used to calculate the T-matrix [18–25], which include the extended boundary condition method (EBCM) [19–21], the point matching (PM) method [6,22], and the discrete-dipole approximation (DDA) [17,22–25]. The most commonly used method is the EBCM, where only surface integrals are required so that it reduces to a 2D method, and even 1D for axially symmetric particles. However, the computation suffers from ill conditioning for particles that become too nonspherical, such as highly elongated or flattened particles with aspect ratios (ARs) very different from 1. This is because of the values of the basis functions used; vector spherical wave functions (VSWFs) vary by many orders of magnitude over the surface of the particle. Similar problems apply to PM. The sizes and shapes for which the different methods become too inaccurate will 0146-9592/14/164827-04$15.00/0

vary with the different methods, since it is integrals of products of VSWFs that appear in EBCM, and the VSWFs themselves in PM. In addition, the linear system in PM is overdetermined. Thus, a study of the accuracy of the different methods can allow a suitable method to be chosen. The computational efficiency of the methods matters as well. In principle, DDA will work for particles of all ARs, but as it is a volume-based method, it can be significantly slower than the surface-based methods (EBCM and PM). In turn, the PM method can outperform the EBCM in terms of speed, by avoiding the time-consuming surface integrals required by the EBCM. Thus, knowledge of the relative accuracies of the methods allows the fastest and sufficiently accurate method to be chosen. In this Letter, we compare the accuracy and efficiency of these three methods to calculate the T-matrix (EBCM, PM, and DDA). We compare the calculated axial optical forces acting on different sizes and shapes of cylindrical particles with their main axis lying along the trappingbeam axis. First, we use the convergence of DDA results with decreasing dipole spacing (DS) to obtain a reliable calculation of the optical force. We then proceed to compare forces calculated using the EBCM and PM. We use the axial optical force to test the accuracy of T-matrix methods, since it depends on both the gradient and scattering forces. The scattering force is sensitive to errors affecting the power because errors in power can have a large effect on the change in axial momentum of the beam, since this depends on the difference between the final momentum and initial momentum. The axial optical force F z experienced by the particle in the optical trap provides a suitable parameter for comparison. We define the normalized deviation (ND) of F z as R jF DDA z − F T zj2 dz R : (1) ND  jF DDA zj2 dz DDA does not result in ill-conditioned linear systems for extreme ARs and so provides a suitable baseline © 2014 Optical Society of America

OPTICS LETTERS / Vol. 39, No. 16 / August 15, 2014 −4

1

x 10

(b)

0.8 0.7 0.6 0.5 0.4 0.3 0.2

0.12

0.3

x 10

For spheres with d=200nm

Normalized Deviation (ND)

Normalized Deviation (ND)

6

(a)

0.9

(a)

0.25

For spheres with d=2000nm

DDA method EBCM method PM method

0.2

5 4 3 2

0.15 0.1 0.05

10

15 20 25 30 35 40 Dipole spacing of DDA method (nm)

45

50

0

50

55

60 65 70 75 80 85 90 Dipole spacing of DDA method (nm)

0

95 100

0.1

−2

−1

0

1

2

3

(c)

−3

−2

−1

0

1

2

3

z (m)

4 x 10 −6

DDA method EBCM method PM method d = 1000nm L = 2000nm

0.06

F (N)/pN

−0.04 −4

4 −6

x 10

0.08

Fig. 1. Normalized deviation (ND) of the DDA method with various DS from the Lorenz–Mie theory: (a) for a small sphere with d  200 nm and (b) for a large sphere with d  2000 nm.

0.02

−0.02 −3

d = 1000nm L = 1000nm

0.04

−0.05

z (m)

0

DDA method EBCM method PM method

(b)

0.06

0

−0.1 −4

1

0.1

0.1 0.08

d = 1000nm L = 150nm

F (N)/pN

−4

F (N)/pN

4828

0.04 0.02 0 −0.02

for comparison. To confirm the accuracy of the DDA results, we perform two important tests. First, we show that the chosen DS (5 nm when the diameter d ≤ 100 nm, 25 nm for 100 nm < d ≤ 600 nm, and 50 nm for 600 nm < d ≤ 2000 nm) gives adequate convergence. Second, we compare the results with a known analytical solution, the analytical Lorenz–Mie solution for a sphere. The convergence of DDA with decreasing DS is shown in Figs. 1(a) and 1(b). The chosen DS is marked with red ovals; at these values, the DDA results are accurate enough to test the accuracies of the EBCM and PM method. For cylinders, we cannot compare with an analytical solution, but we can confirm convergence of DDA for both small and large cylinders. The axial optical force calculated by DDA with the chosen DS is compared with that obtained with half of the DS, as shown in Figs. 2(a) and 2(b). To ensure convergence pof both the EBCM and PM, we choose N max  kr  3 3 kr [13], where r is half of the diagonal length of the cylinder (and thus the radius of the circumscribing sphere). This truncation criterion results in a larger error than the Wiscombe criterion [12], but the relative error of approximately 10−6 is adequately small [7]. If we use a larger N max than necessary, the ill conditioning of the resulting linear systems becomes worse, and errors increase. The convergence of DDA with N max is tested as well; the resting error is much smaller than the DDA discretization error. The incident beam is a Gaussian beam propagating along the positive z axis in water, with free-space wavelength λo  1064 nm and beam power P o  100 mw. The refractive index of the particles is chosen to be np  1.59 (that of polystyrene). The numerical aperture of the objective lens forming the trap is NA  1.02. x 10

5

0.06

(a)

Dipole spacing =25nm Dipole spacing =12.5nm

(b)

3

d=2000nm L=200nm

2 1 0

0.03

−1

0

z (m)

1

2

3

4 −6

x 10

Fig. 3. Comparisons of the axial optical force calculated by DDA, EBCM, and PM for cylinders with different ARs: (a) flattened cylinder (AR  0.15), (b) regular cylinder (AR  1), and (c) elongated cylinder (AR  2).

We use all three methods to calculate T-matrices for cylinders with diameter d  1000 nm and lengths L  150, 1000, and 2000 nm, which then were used to calculate axial optical forces. The results are shown in Figs. 3(a), 3(b), and 3(c), respectively. It is found that the optical-force profiles have large difference when the AR is very different from 1; namely, for the flattened (AR  0.15) and elongated cylinders (AR  2). But for the regular cylinder with AR  1, these three methods match very well. In addition, we compared the computation times of calculating T-matrices for cylinders (AR  1) by our MATLAB implementations of PM, EBCM, and DDA in Table 1, which indicates that DDA is significantly slower when diameter increases than the EBCM and PM. The calculations were carried out in MATLAB on a 4.0 GHz PC. We further investigate the different varying trends of accuracy with the size and AR for both Rayleigh and wavelength-scale cylinders. The dependence of the accuracy on size and AR was explored. This revealed that none of the ND calculated by the EBCM and PM linearly increase with the AR for Rayleigh cylinders and wavelength-scale cylinders. Furthermore, the accuracy of both the EBCM and PM for varying AR are very different for Rayleigh cylinders (d < λo ∕10) and wavelength-scale cylinders (λo ∕10 < d < 5λo ). As shown in Fig. 4(a), the EBCM is more accurate than the PM method when

Diameter (nm)

PM

EBCM

DDA

0.02

−1 0.01

−2 −3 −4

−2

Computation Times for Calculating T-Matrices

0.04 F (N)/pN

F (N)/pN

Dipole spacing=50nm Dipole spacing=25nm

0.05

d=200nm L=200nm

4

−3

Table 1. Computation Times of Calculating TMatrices for Cylinders (AR  1) by PM, EBCM, and DDAa

−3

6

−0.04 −4

−3

−2

−1

0 z (m)

1

2

3

4

−6

x 10

0 −4

−3

−2

−1

0 z (m)

1

2

3

4

−6

x 10

Fig. 2. Comparison of axial optical force for cylinders obtained by the DDA method with the current and half of the current DS: (a) for a small cylinder with d  200 nm and L  200 nm, DS is 25 and 12.5 nm; (b) for a large cylinder with d  2000 nm and L  200 nm, DS is 50 and 25 nm.

200 400 600 800 1000 1500 2000 a

0.1078 0.2192 0.4407 0.7874 1.4220 5.3424 7.1898

s s s s s s s

0.7155 2.4949 7.3998 17.7922 35.7306 1.7751 4.2339

s s s s s m m

1.592 s 20.313 s 3.2620 m 18.5143 m 1.04 h 1.532 h 5.4011 h

The calculations were carried out in MATLAB on a 4.0 GHz PC.

August 15, 2014 / Vol. 39, No. 16 / OPTICS LETTERS

calculating the axial optical force for Rayleigh cylinders with diameter d  40 nm. Moreover, the flattened and elongated cylinders with the AR at small and large regions, respectively, have relatively lower accuracy compared with the region in between. Compared with Fig. 4(b) where the ND versus AR for wavelength-scale cylinders are given, the accuracy for Rayleigh cylinders varies more slowly with AR. Furthermore, to directly reflect the differences in the axial optical force obtained by the EBCM, PM, and DDA at various accuracies, the axial optical-force profiles calculated by these three methods are shown as insets in Figs. 4(a) and 4(b). For Rayleigh cylinders with AR of 1.0 (the length is 40 nm), ND are 0.33 and 1.15; namely, the error is 33% and 115%, respectively. It is observed that both the F z profiles from the EBCM and PM are very different from the DDA result. However, when the AR is 2.0, the error for EBCM reduces to 0.0054%, which results in a nearly perfect match to the F z profile. But the PM result still suffers from an error of 51%. For wavelength-scale particles, both the EBCM

Fig. 4. Relationship between the accuracy of the EBCM and PM method and the ARs of cylinders. (a) Rayleigh cylinders. (b) Wavelength-scale cylinders.

4829

and PM have better accuracy in most regions, although flattened and elongated cylinders in this case suffer from large errors. It is also noted that the PM method has very good accuracy in the middle region of the AR, which can be further demonstrated by the F z profiles from the insets in Fig. 4(b). The errors of the EBCM and PM are 0.036% and 0.08%, respectively, when the AR is 1.0, while those for ARs of 2.25 are 0.02% and 11%, respectively. It is observed that the error of 11% is still acceptable. The PM method quickly fails when the AR is larger than 2.75. To extensively present the accuracy of the EBCM and PM method for calculating scattering by cylinders, Figs. 5(a) and 5(b) show the regions where the accuracy is better than 1% and 10% for the EBCM and PM. The black solid curves border the regions with errors smaller than 10%. Within these, the red dotted lines circulate the regions with errors below 1%. It should be noted that the region of error of EBCM > DDA, one can choose the corresponding method for the specific particle size with acceptable accuracy. In conclusion, the accuracy and efficiency of the EBCM, PM, and DDA for calculating the T-matrix in the scattering field by cylinders are compared and analyzed. The relationship between the accuracy and size and AR for Rayleigh and wavelength-scale cylinders using EBCM and PM are investigated. Compared with Rayleigh cylinders, wavelength-scale cylinders suffer from large error when the AR is very different from 1, but

Fig. 5. Accuracy mapping of scattering by cylinders. (a) EBCM method. (b) PM method.

4830

OPTICS LETTERS / Vol. 39, No. 16 / August 15, 2014

has better accuracy at regular AR. Furthermore, we determine the size and AR regimes where EBCM and PM are accurate to better than 1% and 10%. Since the PM method has the highest efficiency, it can be useful to choose the PM method when the particle size is within its region of acceptable error. For larger particles, the EBCM method should be chosen when it provides acceptable error and PM does not. If neither PM or EBCM are sufficiently accurate, DDA can be used. This work was supported by the National Natural Science Foundation of China under grants 61106049 and 61377071 and also by the Australian Research Council grants DP140100753 and DP1095880. References 1. A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Opt. Lett. 11, 288 (1986). 2. F. M. Fazal and S. M. Block, Nat. Photonics 5, 318 (2011). 3. K. Svoboda and S. M. Block, Annu. Rev. Biophys. Biomol. Struct. 23, 247 (1994). 4. T. A. Nieminen, N. du Preez-Wilkinson, A. B. Stilgoe, V. L. Y. Loke, A. A. M. Bui, and H. Rubinsztein-Dunlop, J. Quant. Spectrosc. Radiat. Transfer 146, 59 (2014). 5. T. A. Nieminen, V. L. Y. Loke, A. B. Stilgoe, G. Knoner, A. M. Branczyk, N. R. Heckenberg, and H. Rubinsztein-Dunlop, J. Opt. A 9, S196 (2007). 6. T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Heckenberg, J. Quant. Spectrosc. Radiat. Transfer 79, 1019 (2003). 7. T. A. Nieminen, V. L. Y. Loke, A. B. Stilgoe, N. R. Heckenberg, and H. Rubinsztein-Dunlop, J. Mod. Opt. 58, 528 (2011). 8. S. H. Simpson and S. Hanna, J. Opt. Soc. Am. A 24, 430 (2007).

9. M.I. Mishchenko,N. T.Zakharova, G.Videen, N.G. Khlebtsov, and T. Wriedt, J. Quant. Spectrosc. Radiat. Transfer 111, 650 (2010). 10. L. Ludvig, Kongel. Danske Vidensk. Selsk. Skr., Naturvidensk. Math. Afd. 6, 1 (1890). 11. G. Mie, Ann. Phys. 330, 377 (1908). 12. W. J. Wiscombe, Appl. Opt. 19, 1505 (1980). 13. J. H. Bruning and H. T. Lo, IEEE Trans. Antennas Propag. 19, 391 (1971). 14. D. B. Phillips, S. H. Simpson, J. A. Grieve, G. M. Gibson, R. Bowman, M. J. Padgett, M. J. Miles, and D. M. Carberry, Opt. Express 19, 20622 (2011). 15. F. Borghese, P. Denti, R. Saija, and M. A. Iati, Opt. Express 15, 1418 (2007). 16. P. B. Bareil and Y. L. Sheng, Opt. Express 18, 26388 (2010). 17. D. B. Phillips, M. J. Padgett, S. Hanna, Y. L. D. Ho, D. M. Carberry, M. J. Miles, and S. H. Simpson, Nat. Photonics 8, 400 (2014). 18. F. M. Kahnert, J. Quant. Spectrosc. Radiat. Transfer 79, 775 (2003). 19. T. A. Nieminen, H. Rubinsztein-Dunlop, N. R. Heckenberg, and A. I. Bishop, Comput. Phys. Commun. 142, 468 (2001). 20. T. A. Nieminen, H. Rubinsztein-Dunlop, and N. R. Heckenberg, J. Quant. Spectrosc. Radiat. Transfer 70, 627 (2001). 21. F. M. Kahnert, J. J. Stamnes, and K. Stamnes, Appl. Opt. 40, 3110 (2001). 22. V. L. Y. Loke, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, J. Quant. Spectrosc. Radiat. Transfer 110, 1460 (2009). 23. P. C. Chaumet, A. Rahmani, A. Sentenac, and G. W. Bryant, Phys. Rev. E 72, 046708 (2005). 24. V. Karásek, O. Brzobohatý, and P. Zemánek, J. Opt. 11, 034009 (2009). 25. S. H. Simpson and S. Hanna, Opt. Express 19, 16526 (2011).

Comparison of T-matrix calculation methods for scattering by cylinders in optical tweezers.

The T-matrix method, or the T-matrix formulation of scattering, is a framework for mathematically describing the scattering properties of an object as...
462KB Sizes 0 Downloads 4 Views