Vector finite element modeling of optical tweezers

Computer Physics Communications 128 (2000) 558–564 www.elsevier.nl/locate/cpc

Vector finite element modeling of optical tweezers ✩ Daniel A. White 1 Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, L-560. P.O. Box 808, Livermore, CA 94551, USA Received 16 November 1999

Abstract It has been established that under certain conditions microscopic dielectric objects can be trapped by a tightly focused laser beam. This phenomena is commonly referred to as an optical tweezer, and it has important applications in biological instrumentation. In this paper we describe how the vector finite element method can be used to compute optical trapping efficiencies for arbitrarily shaped dielectric objects, and we present new results for trapping efficiencies versus dielectric constant.  2000 Elsevier Science B.V. All rights reserved.

1. Introduction It is well known that light waves have momentum and that this momentum can be transferred to a solid object. It has been shown experimentally that a microscopic dielectric object can be trapped by a tightly focused laser beam. This is referred to as a single-beam optical gradient trap, or “optical tweezers”. The optical trapping force has been measured for dielectric spheres under a variety of conditions [1–4]. The trapping force depends upon many factors including the size and dielectric constant of the object, the wavelength and polarization of the laser beam, and the position of the object with respect to the beam focus spot. The predominate application of optical traps is microbiology, where optical traps have been used for the manipulation of cells [5–7] and the manipulation of viruses and bacteria [8,9]. As another biological application, a properly calibrated optical trap can be used to ✩

This work was supported in part by the U.S. Department of Energy under Grant No. W-7405-Eng-48. 1 E-mail: [email protected].

measure minute biological forces [10–12] within living cells. In this paper, we describe describe how a full-wave vector finite element method can be used to compute optical trapping efficiencies for arbitrary 3D dielectric objects. The method solves Maxwell’s equations on an unstructured three-dimensional grid using edge vector finite elements as a basis for the electric field and face vector finite elements for the magnetic flux density. The fields are integrated in time using the leapfrog method. In order to verify the utility of the vector finite element method we model the optical trapping of dielectric microspheres for which measured data is available. The computed trapping efficiencies are shown to compare favorably with measured values for 1 µm silica spheres, and we also present trapping efficiencies for a range of dielectric constant.

2. Overview of the vector finite element method There are a variety of ways of expressing Maxwell’s equations, in this paper the electric field EE and the magnetic flux density BE are used as the principle vari-

0010-4655/00/$ – see front matter  2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 0 ) 0 0 0 0 7 - 2

D.A. White / Computer Physics Communications 128 (2000) 558–564

ables. Maxwell’s equations for the electromagnetic fields in an inhomogeneous volume Ω are µ−1

ε

∂ E B = −µ−1 ∇ ∂t E × EE − µ−1 σM µ−1 BE − µ−1 M,

∂ E E = ∇ × µ−1 BE − σE EE − JE, ∂t

(1) (2)

∇ · εEE = 0,

(3)

∇ · BE = 0.

(4)

Note that we have assumed zero charge density. The dielectric permittivity ε, the magnetic permeability µ, and the conductivities σE and σM are assumed to be symmetric positive-definite tensors, which are functions of position only. The initial conditions and boundary conditions are given by E = 0) = BEic B(t

E = 0) = EE ic , E(t nˆ × EE = EEBC

on Γ.

in Ω,

(5) (6)

In addition, for the problem to be well-posed we require that the independent magnetic and electric current sources satisfy E = 0, ∇ ·M

∇ · JE = 0.

(7)

Maxwell’s Eqs. (1) and (2) are discretized using the Galerkin method with the vector finite elements defined in [13]. In this approach different finite element basis functions are used for the electric and magnetic fields. The electric field is approximated as a linear combination of edge elements, EE =

Ne X

Ei , αi (t)W

(8)

i=1

where Ne is the number of edges in the mesh. The electric field given by Eq. (8) is a first order approximation to the true field. These functions are called edge functions because they have the property Z E i · tˆj = δij , (9) W where tˆj is the unit tangent along edge j . Hence the degrees-of-freedom αi are given by Z (10) αi = EE · tˆi

559

and can be interpreted as the voltage along edge i of the mesh. The magnetic flux density is approximated as a linear combination of face elements, BE =

Nf X

βi (t)FEi ,

(11)

i=1

where Nf is the number of faces in the mesh. The magnetic flux density given by Eq. (11) is a first order approximation to the true flux density. These functions are called face functions because they have the property Z (12) BEi · nˆ j = δij , where nˆ j is the unit normal to face j . Hence the degrees-of-freedom βi are given by Z (13) βi = BE · nˆ i da and can be interpreted as the magnetic flux through face i of the mesh. Using the approximations (8) and (11) in the Galerkin form of Maxwell’s equations yields a coupled system of ordinary differential equations ∂ (14) C α(t) = KT β(t) − Gα(t) − Rj (t), ∂t ∂ (15) D β(t) = −Kα(t) − Pβ(t) − Dm(t), ∂t where the matrices are given by
Ei , W Ej , (16) Cij = εW
E i , FEj , (17) Kij = µ−1 ∇ × W
Ej , Ei, W (18) Gij = σE W
E i , FEj , (19) Rij = W
−1 E E (20) Dij = µ Fi , Fj ,   σM E E Fi , Fj , Pij = (21) µ2 where (u, v) represents the integral of u · v over the volume Ω. The above system of ODE’s is integrated in time using the second-order central difference “leapfrog” method. In this method the electric degrees-of-freedom are known at whole time steps, the magnetic degreesof-freedom are known at the half time steps. Using superscript n to denote the time level, the ODE is integrated according to

560

D.A. White / Computer Physics Communications 128 (2000) 558–564

  1t G α n+1 C+ 2

  1t G α n − Rj n+1/2 , = 1tKT β n+1/2 + C − 2 (22)   1t P β n+1/2 D+ 2   1t n P β n−1/2 − 1tDmn+1/2 . = −1tKα + D − 2 (23) This integration requires the solution of large, sparse, linear systems at every time step. The accuracy, stability, and efficiency of the vector finite element are described in [14–16].

3. Modeling the optical tweezer A three-dimensional hexahedral grid is used to model the dielectric object and the surrounding medium. In practice the dielectric object can be of any shape; in this paper the dielectric object is a simple sphere and an unstructured hexahedral grid is used that conforms to the surface of the spherical object. Physically there is no outer boundary to the electromagnetic scattering problem but in practice the computational grid must be of finite extent. The computational grid used for all computations in this paper is shown in Fig. 1. The initial electromagnetic fields in the computational domain are zero. A time varying source is used to launch the laser beam into the computational domain. The source is designed such that the laser beam comes to a focus at the center of the spherical grid. The frequency of the laser and the beam cone angle are user-specified parameters. Since the computational domain is finite, the cross section of the incident beam is also finite. Rather than abruptly truncating a Gaussian or polynomial profile, which would introduce artificial high-frequency components, the beam is modeled with a smooth cos2 cross section. The time dependence of the source is given by     t 2 sin(ωt). (24) E(t) = 1 − exp − 2τ The parameter τ determines the rise time of the beam, typically τ is chosen to be four periods. The simulation

Fig. 1. Sample hexahedral grid of a dielectric sphere embedded in a spherical volume.

Table 1 Conductivity values for Maxwellian absorber Layer

1

2

σ

1.25

5.0

3 11.25

4 20.0

5 31.25

6 45.0

is typically run for twelve periods at a sampling rate of 50 samples per period. This is enough time for the simulation to reach steady state. The computational grid is terminated with a Maxwellian absorber in order to simulate an infinite medium [17,18]. In this paper a Maxwellian absorber is defined as several layers of artificial anisotropic media, with the conductivity of the layers graded in such a way as to reduce reflections from the boundary. This can be interpreted as an impedance matching device. In a Maxwellian absorber the electric conductivity σE and the magnetic conductivity σM are equal. In this paper a uniaxial tensor conductivity is employed with the optical axis normal to the boundary of the computational grid. The conductivity in along the optical axis is zero, the conductivity normal to this axis is given in Table 1. The permittivity and permeability are that of free space. The choice of a spherical outer boundary simplifies the implementation of the Maxwellian absorber and results in excellent absorption of the nearly spherical scattered fields. Note that a spherical outer boundary is not required,

D.A. White / Computer Physics Communications 128 (2000) 558–564

561

and that the method used here does not require an orthogonal or structured grid. The time dependent source used to launch the wave into the domain is applied to the grid surface separating the absorber from the interior of the domain, this is the outer surface shown in Fig. 1. Given the electromagnetic fields in the vicinity of the dielectric object the net force on the object is calculated using the Maxwell stress tensor. The Maxwell stress tensor T is given by E · EE + BE · HE )I − (D E EE + BE HE ), T = 1 (D (25) 2

E and BE are the electric where I is the identity tensor, D and magnetic flux densities, and EE and HE are the electric and magnetic fields. The net force on the dielectric object is then I E ˆ (26) F = T · ndΓ, Γ

where Γ is the surface surrounding the object and nˆ is the unit surface normal. The integral is computed using a two-dimensional trapezoidal rule. The Maxwell stress tensor is not continuous across the dielectric discontinuity, in the evaluation of (26) the fields are evaluated on the vacuum side of the interface. The basis for this approach is that it gives the correct answer for the extreme case of scattering by a perfect conductor, which would have zero fields on the inside surface. In order for there to be a significant net force, there must be a significant gradient in the electromagnetic field, hence the terminology optical gradient trap. It is important to note that the force computed via (26) is time varying, but the inertia of the dielectric object is such that it cannot possibly respond to the rapidly fluctuating fields. Thus the optical force that is measured in the laboratory is a time averaged force. The calculation of the stress and the net force is straightforward; any error in the force calculation is do to error in the electromagnetic fields, which is determined by the grid spacing the and the performance of the Maxwellian absorber. Based on previous computational experiments in which computed electromagnetic fields for a dipole antenna radiating into free space are compared to analytical results it is estimated that the electromagnetic fields are correct to within ten percent [14] for a grid density of 6 cells per wavelength. A dielectric object can be trapped below the focus spot on the axis of the laser beam. This is referred to

Fig. 2. Geometry of axial trapping experiments. The laser beam is propagating in the negative x direction.

as axial trapping. The geometry of axial trapping is illustrated in Fig. 2. The incident laser beam is both reflected and refracted by the dielectric object, resulting in a backscattered field and a forward scattered field. The net force on the dielectric object will be towards the focus spot only if there is a significant amount of forward scattering. The amount of forward scattering depends upon the size and dielectric constant of the object, the cone angle of the laser beam, and the position of the object with respect to the focus spot. For axial trapping the force is independent of the polarization of the laser. In the following computational experiments λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ . The laser rise time τ was equal to four periods and the fields were sampled 50 times per period, which is well below the time step required for stability. The small time step is employed for an accurate computation of the time-average force on the object. The fields were updated for 12 periods, or 600 time steps. The simulation was performed using a variety of dielectric constants in order to determine the effect of dielectric constant on the optical trapping efficiency. The optical trapping efficiency defined as cF , (27) nP where c is the speed of light, n is the index of refraction, and P is the power of the laser beam. The timeaverage net force on the dielectric sphere was computed via (26) for dielectric constants ranging from 1.0 to 2.0, and the corresponding trapping efficiency is

Q=

562

D.A. White / Computer Physics Communications 128 (2000) 558–564

Fig. 3. Computed axial trapping efficiency Q vs. relative dielectric constant. Parameters λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ .

shown versus dielectric constant in Fig. 3. The dielectric sphere is trapped (positive Q) for epsilon ranging from 1.12 to 1.55. For ε < 1.1 the dielectric sphere does not refract the fields enough to have a significant effect on the propagation of the beam, for ε > 1.55 there is significant backscattering of the beam and the dielectric sphere is pushed away from the focus. The computed Q values presented here are comparable to the measured results presented in [2]. Note that in [2] the axial force was measured in water, hence the dielectric constant used in this paper should be considered relative to the background dielectric. The amorphous silica spheres used in [2] had a relative dielectric constant of 1.12. For this value of ε, we compute a Q of 0.008, compared to a measured Q of 0.006 ± 0.001. The difference between the measured and computed trapping efficiencies is likely due to approximating the incident laser beam by a cos2 cross section, which may lead to a tighter focus than can be achieved by the actual optical system used in the physical experiment. A dielectric object can also be trapped transversely, with the object adjacent the focus spot rather than below it. The geometry of transverse trapping is illustrated in Fig. 4. The parameters for this experiment are identical to the axial trapping experiment above, except for the orientation of the laser beam. In addition, this computational experiment is performed for two different polarizations of the incident laser beam since the geometry is not symmetric. The force was computed via (26) for dielectric constants ranging from 1.0 to 4.0. For the transverse trapping experiments the computed force has both x and z components, and is

Fig. 4. Geometry of transverse trapping experiments. The laser beam is propagating in the positive z direction.

it difficult to define a single trapping efficiency Q. Instead, we define cFx cFz , Qz = , (28) nP nP where Fx and Fz are the x and z components of the force, respectively. These Q’s are shown in Figs. 5 and 6. A positive Qx means that the x component of the force is towards the focus spot, and Fig. 5 illustrates that the dielectric sphere will be pulled towards the focus for particular values of the dielectric constant. Comparing to the axial trapping experiment, we note that the magnitude of the optical forces for the transverse experiment are significantly greater than that obtained for the axial experiment. This is in qualitative agreement with the measured results in [2], which report transverse trapping efficiencies of 0.15, compared to 0.006 for the axial experiment. In addition, the sphere is pulled towards the focus for a larger range of dielectric constant, approximately 1.1 < ε < 3.7. This phenomena has not yet been verified experimentally. An important difference between the computational experiment and the physical experiment is that in the physical experiment the dielectric sphere is free to move with respect to the focus spot whereas in the computational experiment the sphere is at a fixed location. In the physical experiment the laser drags the sphere, hence it is the maximum trapping force that is measured. Therefore it is expected that the computational results underestimate the trapping efficiency, Qx =

D.A. White / Computer Physics Communications 128 (2000) 558–564

563

Fig. 5. Computed transverse trapping efficiency Qx vs. relative dielectric constant, x polarization. Parameters λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ .

Fig. 7. Computed transverse trapping efficiency Qy vs. relative dielectric constant, y polarization. Parameters λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ .

Fig. 6. Computed transverse trapping efficiency Qz vs. relative dielectric constant, x polarization. Parameters λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ .

Fig. 8. Computed transverse trapping efficiency Qz vs. relative dielectric constant, y polarization. Parameters λ = 1.0 µm, d = 1.0 µm, x = 0.4 µm, and θ = 45◦ .

since the distance d = 1.0 µm used in the computations may not have been optimal. In addition the physical experiment is such that the sphere moves in the z direction until the z component of the trapping force is zero. The above transverse trapping experiment was repeated for the case of electric field polarized in the y direction. The time-average net force on the dielectric sphere is again computed for a range of dielectric constants from 1.0 to 4.0. The x and z directed forces are converted to Q’s via (28) and these Q’s are shown in Figs. 7 and 8. The z component of Q is quite similar to that obtained for the x polarization experiment, it is negative for all values of dielectric constant. The x component of Q differs from that obtained for the x polarization experiment in that it remains positive for

all values of dielectric constant. This means that the sphere is always pulled toward the axis of the beam.

4. Summary We demonstrate that the vector finite element method is useful for numerical modeling of optical tweezers. A 3D unstructured hexahedral grid is used to model the dielectric object and the surrounding space. A laser beam is launched into the computational grid by a time varying source. As the electromagnetic fields are evolved, the net force on the dielectric object is computed by integrating the Maxwell stress tensor over the surface of the object. The net force is a function of time, but since the dielectric object

564

D.A. White / Computer Physics Communications 128 (2000) 558–564

cannot possible respond to to optical frequencies the time-average force is computed. The computed optical trapping efficiencies are computed for both axial and transverse trapping geometries for a 1 µm sphere, using a range of dielectric constants. A simple sphere is used so that the computed results can be compared to previously measure data. For the axial trapping experiment where the computational model best agrees with the experiment set-up the computed trapping efficiency agrees with the measured to within 14%. The discrepancy can be attributed to differences between the modeled laser beam profile (cos2 ) and the actual laser beam profile (unknown). The results for the transverse trapping experiment show that the trapping efficiency Qx is positive for a large range of dielectric constant. In addition, at relatively large dielectric constants the transverse trapping efficiency shows a strong polarization dependence. This polarization dependence is a new phenomena that has not yet been verified experimentally. Given the computed trapping efficiencies for the transverse experiments Qx and Qz it is not clear exactly when the dielectric sphere is trapped since the x component of the force is positive and the z component of the force is negative. In order to determine if the sphere is ultimately trapped it is necessary to move the sphere in both the x and z directions until Qz is zero and Qx is maximized. This was beyond the scope of this initial project. Future work in this area could include improved modeling of incident laser beams (including aberrations present in actual optical systems), moving the object with respect to the focus spot in order to determine the maximum trapping efficiency, and modeling the trapping of actual biological organisms. References [1] A. Ashkin, J.M. Dziedzic, E.J. Bjorkholm, S. Chu, Observation of a single-beam gradient optical trap for dielectric particles, Opt. Lett. 11 (1986) 288–290. [2] W.H. Wright, G.J. Sonek, M.W. Berns, Parametric study of the force on microspheres help by optical tweezers, App. Opt. 33 (9) (1994) 1753–1748.

[3] C. D’Helon, E.W. Dearden, H. Rubenstein-Dunlop, N.R. Heckenberg, Measurement of the optical force and trapping range of a single-beam gradient optical trap for micron-sized latex spheres, J. Mod. Opt. 41 (3) (1994) 595–601. [4] R. Omori, T. Kobayashi, A. Suzuki, Observation of a singlebeam gradient force optical trap for dielectric particles in air, Opt. Lett. 22 (11) (1997) 816–818. [5] T.N. Buican, M.J. Smith, H.A. Crissman, G.C. Salzman, C.C. Stewart, J.C. Martin, Automated single-cell manipulation and sorting by light trapping, Appl. Opt. 26 (1987) 5311–5316. [6] A. Ashkin, J.M. Dziedzic, T.M. Yamane, Optical trapping and manipulation of single cells using infrared laser beams, Nature 330 (1987) 769–771. [7] Y. Tadir, W.H. Wright, O. Vafa, T. Ord, R.H. Asch, M.W. Berns, Micromanipulation of sperm by a laser generated optical trap, Fert. Ster. 57 (1989) 870–873. [8] A. Ashkin, J.M. Dziedzic, Optical trapping and manipulation of viruses abd bacteria, Science 235 (1987) 1517–1520. [9] S.M. Block, D.F. Blair, H.C. Berg, Compliance of bacterial polyhooks measured with optical tweezers, Cytometry 12 (1991) 492–496. [10] A. Ashkin, K. Schultze, J.M. Dziedzic, U. Euteneur, M. Schliwa, Force generation of organelle transport measured in vivo by an infrared laser trap, Nature 348 (1990) 346–348. [11] S.C. Kuo, M.P. Sheetz, Force of single kinesin molecules measured with optical tweezers, Science 260 (1993) 232–234. [12] S. Block, L.S.B. Goldstein, B.J. Schnapp, Using optical tweezers to investigate kinesin-based motility in vitro, J. Cell Biol. 109 (1989) 81a. [13] J.C. Nedelec, Mixed finite elements in R3, Numer. Math. 35 (1980) 315–341. [14] D. White, Discrete time vector finite element methods for solving Maxwell’s equations on 3D unstructured grids, Ph.D. dissertation, University of California at Davis (1997). [15] G. Rodrigue, D. White, A vector finite element time-domain method for solving Maxwell’s equations on unstructured hexahedral grids, SIAM J. Sci. Comput., submitted. [16] D.A. White, Solution of capacitance systems using incomplete Cholesky fixed point iteration, Comm. Numer. Meth. Engng. 15 (1999) 375–380. [17] R.W. Ziolkowski, The design of maxwellian absorbers for numerical boundary conditions and for practical applications using engineered artificial materials, IEEE Trans. Ant. Prop. 45 (4) (1997) 656–671. [18] Z. Sacks, D. Kingsland, R. Lee, J. Lee, A perfectly matched anisotropic absorber for use as an absorbing boundary condition, IEEE Trans. Ant. Prop. 43 (12) (1995) 1460–1463.