Discrete dipole approximation for time-domain computation of optical forces on magnetodielectric scatterers

We present a general approach, based on the discrete dipole approximation (DDA), for the computation of the exchange of momentum between light and a magnetodielectric, three-dimensional object with arbitrary geometry and linear permittivity and permeability tensors in time domain. The method can handle objects with an arbitrary shape, including objects with dispersive dielectric and/or magnetic material responses. © 2011 Optical Society of America OCIS codes: (290.5850) Scattering, particles; (350.4855) Optical tweezers or optical manipulation. References and links 1. E. M. Purcell and C. R. Pennypacker, “Scattering and absorption of light by nonspherical dielectric grains,” Astrophys. J.186, 705–714 (1973). 2. B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J.333, 848–872 (1988). 3. P. C. Chaumet, A. Rahmani, and G. W. Bryant, “Generalization of the coupled dipole method to periodic structure,” Phys. Rev. B67, 165,404–5 (2003). 4. P. C. Chaumet and A. Sentenac, “Numerical simulations of the electromagnetic field scattered by defects in a double-periodic structure,” Phys. Rev. B 72, 205,437–8 (2005). 5. A. Rahmani, P. C. Chaumet, and F. de Fornel, “Enrironment-induced modification of spontaneous emission: Single-molecule near-field probe,” Phys. Rev. A 63, 023,819–11 (2001). 6. A. Rahmani and G. W. Bryant, “Spontaneous emission in microcavity electrodynamics,” Phys. Rev. A 65, 033,817–12 (2002). 7. F. Bordas, N. Louvion, S. Callard, P. C. Chaumet, and A. Rahmani, “Coupled dipole method for radiation dynamics in finite photonic crystal structures,” Phys. Rev. E 73, 056,601 (2006). 8. A. Rahmani, P. C. Chaumet, and G. W. Bryant, “Discrete dipole approximation for the study of radiation dynamics in a magnetodielectric environment,” Opt. Express 18, 8499–8504 (2010). 9. B. T. Draine and J. C. Weingartner, “Radiative Torques on Interstellar Grains: I. Superthermal Spinup,” Astrophys. J.470, 551–565 (1996). 10. A. G. Hoekstra, M. Frijlink, L. B. F. M. Waters, and P. M. A. Sloot, “Radiation forces in the discrete-dipole approximation,” J. Opt. Soc. Am. A18, 1944–1953 (2001). 11. P. C. Chaumet, A. Rahmani, and M. Nieto-Vesperinas, “Photonic force spectroscopy on metallic and absorbing nanoparticles,” Phys. Rev. B 71, 045,425 (2005). 12. P. C. Chaumet, A. Rahmani, A. Sentenac, and G. W. Bryant, “Efficient computation of optical forces with the coupled dipole method,” Phys. Rev. E 72, 046,708–6 (2005). 13. A. Rahmani and P. C. Chaumet, “Optical Trapping near a Photonic Crystal,” Opt. Express 14, 6353–6358 (2006). 14. P. C. Chaumet and C. Billaudeau, “Coupled dipole method to compute optical torque: Application to a micropropeller,” J. Appl. Phys. 101, 023,106–6 (2007). 15. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: An overview and recent developments,” J. Quant. Spect. Rad. Transf. 106, 558–589 (2007). 16. P. C. Chaumet, K. Belkebir, and A. Rahmani, “Coupled-dipole method in time domain,” Opt. Express 16, 20,157– 20,165 (2008). #129020 $15.00 USD Received 3 Jun 2010; revised 13 Sep 2010; accepted 4 Oct 2010; published 25 Jan 2011 (C) 2011 OSA 31 January 2011 / Vol. 19, No. 3 / OPTICS EXPRESS 2466 17. P. C. Chaumet, K. Belkebir, and A. Rahmani, “Optical forces in time domain on arbitrary objects,” Phys. Rev. A (2010). 18. P. C. Chaumet and M. Nieto-Vesperinas, “Time-averaged total force on a dipolar sphere in an electromagnetic field,” Opt. Lett.25, 1065–1067 (2000). 19. B. D. H. Tellegen, “Magnetic-Dipole models,” Am. J. Phys. 30, 650–652 (1962). 20. L. Vaidman, “Torque and force on a magnetic dipole,” Am. J. Phys. 58, 978–983 (1990). 21. M. Mansuripur, “Radiation pressure and the linear momentum of the electromagnetic field in magnetic media,” Opt. Express15, 13,502–13,518 (2007). 22. M. Mansuripur and A. R. Zakharian, “Maxwell’s macroscopic equations, the energy-momentum postulates, and the Lorentz law of force,” Phys. Rev. E 79, 026,608–10 (2009). 23. E. E. Radescu and G. Vaman, “Exact calculation of the angular momentum loss, recoil force, and radiation intensity for an arbitrary source in terms of electric, magnetic, and toroid multipoles,” Phys. Rev. A 65, 046,609– 47 (2002). 24. E. E. Radescu and G. Vaman, “Toroid moments in the momentum and angular momentum loss by a radiating arbitrary source,” Phys. Rev. A 65, 035,601–3 (2002). 25. P. C. Chaumet and A. Rahmani, “Electromagnetic force and torque on magnetic and negative-index scatterers,” Opt. Express17, 2224–2234 (2009). 26. P. C. Chaumet and A. Rahmani, “Efficient iterative solution of the discrete dipole approximation for mag netodielectric scatterers ,” Opt. Lett. 34, 917–919 (2009). 27. P. C. Chaumet and A. Rahmani, “Coupled-dipole method for magnetic and negative refraction materials,” J. Quant. Spect. Rad. Transf. 110, 22–29 (2009). 28. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302–307 (1969). 29. A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Transactions on Microwave Theory and Techniques 23, 623–630 (1975). 30. A. Taflove, “Application of the finite-difference time-domain method to sinusoidal steady state electromagnetic penetration problems,” IEEE Trans. Antennas Propagat. 22, 191–202 (1975). 31. P. C. Chaumet, “Comment on “Trapping force, force constant, and potential depths for dielectric spheres in the presence of spherical aberrations”,” Appl. Opt. 43, 1825–1826 (2004). 32. J. R. Arias-Gonź alez and M. Nieto-Vesperinas, “Optical forces on small particles: attractive and repulsive nature and plasmon-resonance conditions,” J. Opt. Soc. Am. A 20, 1201–1209 (2003). 33. R. W. Ziolkowski, “Pulsed and CW Gaussian beam interactions with double negative metamaterial slabs,” Opt. Express11, 662–681 (2003).


Introduction
The discrete dipole approximation (DDA, also called coupled dipole method) is a general computational method that is widely used to study scattering problems in electrodynamics and photonics.The method was introduced by Purcell and Pennypacker [1,2] to address the scattering of light by interstellar grains with arbitrary shapes.Since its initial formulation, the DDA been generalized to study the scattering of light by periodic structures [3,4], spontaneous emission [5,6,7,8], optical forces [9,10,11,12,13] and optical torques [14].A recent overview of the DDA can be found in Ref. [15].
In the standard formulation of the DDA the electromagnetic fields are computed in the frequency domain.While this is most appropriate for finding the steady state solution of an electromagnetic problem, this restriction means that one cannot exploit the advantages of the DDA (ability to consider arbitrary scatterers, spatial discretization limited to the scatterer and its immediate neighborhood) to explore transient phenomena.This is a particularly limiting factor for optomechanical studies where one might be interested in modelling the optical force experienced by a scatterer illuminated by a time-varying electromagnetic (EM) field.
Previously we developed a formulation of the DDA to study light scattering, in time domain, by a linear, dispersive, lossy objects [16].More recently, we generalized the time-domain DDA to address optical forces in non-magnetic systems [17].In this article we describe how these two aspects (time varying fields and optical forces on magnetodielectric scatterers) can be combined in a more general formulation of the DDA, that can be used to study optical forces on arbitrary three-dimensional magnetodielectric object with material dispersion and/or losses, in time domain.

DDA for optical force in time domain
Consider an object (scatterer) in vacuum.We start by discretizating the object, over a cubic lattice, into a collection of N polarizable subunits.If the lattice spacing is small enough compared to the spatial variation of the electromagnetic fields, one can treat each subunit within the dipole approximation.Hence, to compute the optical force on an arbitrary object we first need to derive the optical force on a subunit, in the dipole approximation.Consider a small (compared to the spatial variations of the fields) particle located at position r with permittivity ε and permeability µ (we assume that the background medium is vacuum).The particle is illuminated by an arbitrary, time-dependent incident EM field {E (r,t), H (r,t)}.Let P(r,t) [M (r,t)] be the time-dependent electric (magnetic) dipole moment induced in the particle by the electric (magnetic) field of the incident EM wave, then the i-th component of the force experienced by the particle is the sum of the force on the electric dipole [18] plus the force on the magnetic dipole [19,20]: where i, j, and k stand for either x, y or z and ∂ t denotes the derivative with respect to time.This is a generalized form of the Lorentz force [21,22].Using Maxwell's equations the i-th component of the electric and magnetic forces can be written as: where ε i jk is the Levi-Civita tensor.However, the above expression for the force is incomplete.We need to add the radiative reaction force, i.e. the flux of momentum lost by the particle when it radiates electromagnetic waves [23,24] Therefore, the optical force experienced by the particle is the sum of five contributions.The first [F he (r,t)] and third terms [F hm (r,t)] on the right-hand-side of Eq. ( 2) correspond to the "standard" terms that appears in the time harmonic case for the electric dipole and magnetic dipole, respectively.The second [F pe (r,t)] and fourth terms [F pm (r,t)] of Eq. ( 2) whose time average vanish if the incident wave is time harmonic are related to the Poynting vector.In the absence of material dispersion, these two terms are proportional to the time derivative of the Poynting vector.The last term, Eq. (3) which does not vanish for a time harmonic wave [25], is a self-interaction term.Within the framework of the DDA let us consider an arbitrary object illuminated by a timedependent electromagnetic wave with envelop I (t) (Fourier transform I(ω) = G [I (t)]), and a spectrum centered at frequency f 0 .The general aspects of the time-domain formulation of the DDA have been detailed in a previous article, therefore only a brief description will be given here [16].We first solve the linear system associated to the time-harmonic DDA, for M values of the frequency with the incident field I(ω)E(ω), to find the value the electric dipole and magnetic dipole at each subunits location for the time harmonic problems [26,27].Then we compute the temporal inverse Fourier transform of the fields in accordance with the Nyquist-Shannon sampling theorem such that E (r,t) = G −1 [I(ω)E(r, ω)] [16].
We now focus on how to compute the optical force in time domain on a subunit located at r.The terms related to the harmonic force can obviously be computed through: The derivation of the terms related to the time derivative of the Poynting vector is slightly more involved.We first compute these terms without the time derivative: The time derivative is calculated in the frequency domain using the inverse Fourier transform: The recoil term can also be obtained by performing the differentiation in the frequency domain: Using the DDA, once the contributions to the optical force are known for each subunit, the optical force on the scatterer is obtained by adding the optical force on the N subunits.We can also define the momentum imparted to the object by the electromagnetic field as Q i (t) = t 0 F i (t)dt where we have assumed that the origin of time is chosen such that the electromagnetic fields at the object are zero for t < 0.

Advantages of the DDA for the computation of optical forces in time domain
Traditionally, the finite difference time domain (FDTD) method has been used to study electromagnetic scattering problems in time domain [28,29,30].In the FDTD, the differential forms of the time-dependent Maxwell equations are discretized in both space and time.As a result, the entire computational window must be discretized and not just the scatterer.Moreover, in the FDTD, boundary conditions at the edges of the computational window must be handled carefully (usually using perfectly matched layer techniques).Furthermore, numerical dispersion and/or instabilities may occur when the fields are propagated over large distances (large object compared to the wavelength).
By contrast, the DDA circumvent many of these issues.First, there is no computational window per se since only the scatterer is discretized.Once the fields inside the scatterer are computed, the fields anywhere outside can be computed in a straightforward way using a freespace susceptibility tensor that is known analytically.In the FDTD, the fields are computed either within the computational window or in the far-field limit using near-to-far-field transformations.Furthermore, since unlike the FDTD, the DDA is based on the integral form of Maxwell's equations and, therefore, is "built" around the concept of field susceptibility tensor, certain type of geometries such as a semi-infinite substrate, or a multilayer stack can be handled analytically which is not possible with the FDTD.
Another advantage of the DDA is how easy it is to account for material dispersion as showed in the next section which deals with plasmon resonance.This is because each frequency component is dealt with separately, and therefore there is no risk of error propagation in time.Anisotropic scatterers are also easily accounted for as they only amount to the introduction of a dyadic tensor for the electric and/or magnetic polarizabilities.While the computation of optical forces in time domain has been done using FDTD for a 2D structure with a permeability equal to one (non-magnetic case) [22], we are not aware, at the present time, of a general FDTD treatment of optical forces in time domain for a 3D magneto-dielectric scatterer with material dispersion and losses.
Finally, let us mention that in the DDA one needs to solve linear systems of size 6N × 6N where N is the number of discretization subunits for the scatterers under study.Except for the smallest values of N these systems cannot be solved by direct matrix inversion, rather, one should use an efficient iterative scheme [26] with an appropriate initial estimate of the solution [16].

How many contributions to the time-dependant optical force?
Before considering the case of a discretized scatterer with material dispersion and losses in Sec. 5, to help us understand why we have introduced five force terms, in this section we take a closer look at the different contributions to the time-dependent optical force in the simpler case of a dipole scatterer in the presence of a time-harmonic incident plane wave.Consider a spherical scatterer with radius a (in the dipole approximation) illuminated by a time-harmonic plane wave with normalized amplitude [E x = H y = cos(kz − ωt)].For the sake of simplicity we assume that the sphere is homogeneous with real permittivity and permeability.We also neglect material dispersion.Taking into account the radiative reaction term necessary to satisfy the optical theorem [2,31], we can write the five terms of the time-dependent force introduced previously as: where α m 0 and α e 0 are the quasistatic polarizabilities.Notice that in the time harmonic regime, it is common to split, artificially, the optical force into two contributions: the gradient force and the radiation pressure (or scattering force) [32].In the time domain, however, such a decomposition of the optical force into two terms only is no longer adequate as other contributions emerge.In the time domain picture, the counterparts to the gradient force and radiation pressure are included in the time-dependent terms F he and F hm but as seen in Eqs. ( 11)-( 12) the spatial derivatives of the electromagnetic field introduce a term proportional to sin(2kz − 2ωt) which vanishes on average (there is no gradient force with a plane wave illumination), and a term proportional to sin 2 (kz − ωt) which gives the radiation pressure when the time average is performed.The three other terms [Eqs.( 8)-( 10)], on the other hand, are obtained by differ- entiation of the fields with respect to time and are not related to the concept of gradient force or radiation pressure of the time-harmonic regime.Notice that in Eq. ( 15) we have neglected terms in (ka) 7 and (ka) 10 .Therefore, the total force, in time domain, on a dipole illuminated by a time-harmonic plane wave is: where the second term cannot be neglected since its time-average is not zero (contrary to the first term).In other words, when dealing with a particle in the dipole approximation (or an actual collection of dipoles) it is essential to take into account terms beside the usual gradient force and radiation pressure in calculating the total momentum imparted by the electromagnetic wave to the dipole.

Contributions to the force: single dipole versus discretized scatterer.
As we discussed, the foundation of the DDA is the description of the various scattering processes at the level of a single dipole (subunit).However, most of the time the underlying structure of the scatterer as a collection of dipoles must be viewed as a convenient numerical description as opposed to a true representation of the internal geometry of the scatterer.Because of this, the self-term in Eq. (3) will lose its significance when a finite object is discretized into a large collection of dipoles.To illustrate this point, consider a spherical scatterer made of a double-negative medium for which the material responses are described by lossy Drude models [33].In the frequency domain, the permittivity and permeability are given by: where ω pe , ω pm , Γ e and Γ m denote the corresponding plasma and damping frequencies respectively.Let us start by assuming that the incident field is a plane wave with a Gaussian time envelop of the form: where f 0 = ω 0 /2π = c/λ 0 is the central frequency of the pulse τ = 8/ f 0 is the duration of the pulse.The spectral and time profiles of the incident pulse are plotted in Figs.1(a) and 1(b).Let us first consider a sphere with radius a = λ 0 discretized into N = 33552 subunits.The parameters for this computation are f 0 chosen such that Re[ε(ω 0 )]=Re[µ(ω 0 )]=-1, i.e. ω pe = ω pm = w 0 √ 2, and damping terms Γ e = Γ m = ω pe /10.This corresponds to a sphere made of a lossy, negative-index material and exhibiting a resonance at frequency ω 0 , i.e. in the timeharmonic case the time average of the optical force would be maximum for ω = ω 0 .Notice that due to the large size of the sphere the resonance does not occur for Re[ε(ω 0 )]=Re[µ(ω 0 )]=-2 (plasmon resonance of a sphere much smaller than the wavelength) but is shifted toward Re[ε(ω 0 )]=Re[µ(ω 0 )]-1 (surface plasmon resonance).We plot in Fig. 2(a) the total optical force and its different contribution, as a function of time.Because ε(ω) = µ(ω), and the equivalence between the plasma frequencies and damping terms for the electric and magnetic parts of the material response, we have F pm (t) = F pe (t) and F hm (t) = F he (t).Note that, in the present configuration, the oscillations of the total force are mainly due to the term associated with the Poynting vector.In Fig. 2(b) we can see that the contribution of this term to the transfer of momentum from the EM field to the object vanishes (as expected) at the end of the pulse and only the harmonic contributions to the force remain.The convergence of the method is illustrated in Fig. 2(c) which shows the momentum imparted by the EM wave to the object, versus the number of subunits.The computed value of the momentum only changes by 1% when the number of subunits is increased from N = 4224 to N = 113104.
The main difference between the optical force on a single dipole and on one subunit of a discretized object comes from the recoil (self-interaction) of Eq. ( 3).The recoil (radiation reaction) force showed Fig. 2(d) is very weak compare to the other terms.In fact the magnitude of this term decreases as the number of subunit N increases as illustrated by the two curves plotted for different values of N, Figs.2(d) and 2(e) for the force and the momentum respectively.This is due to the fact that the radiation reaction force for each subunit scales as the volume of the subunit squared (cross product of the electric and magnetic dipole moments), hence when N, the number of subunits, increases the contribution to the total force on the object of this recoil term (summed over all the subunits) decreases like 1/N as illustrated in Fig. 2(f).Note that this term would vanish in the limit where N tends to infinity, however, its contribution to the force would then be taken into account through the other contributions of the force and the multiple scattering between the subunits.In other words, the separation of the total optical force into 5 terms, while helpful in understanding the origin of the force experienced by a small particle, is somewhat artificial in the case of a discretized object.As a result, as the number of discretized subunits tend to infinity, the various contributions of the total forces can be grouped into four terms instead of the five terms, which is consistent with the expression for the generalized Lorentz force derived by Mansirupur [22].
However, we emphasize again that if one is interested in calculating the optical forces on a single, or a collection of particles treated in the dipole approximation it is essential to take into account the recoil term explicitly.This was illustrated in Sect.4, where in Eq. ( 16) the recoil term gives the term involving the product of the electric and magnetic polarizabilities which is of the same magnitude as the other two terms in the square brackets.

Influence of losses and plasmon resonances
If a small sphere is illuminated with a plane wave at frequency f 0 , the spectrum of the optical force would exhibit peaks at two frequencies: the zero frequency and 2 f 0 .This can be seen in Eq. ( 16) where there is a term cos 2 (kz − ωt) in the expression of the total force.Accordingly, for an homogeneous sphere with no dispersion and an illumination given by Eq. ( 18) we observe two peaks: one at zero frequency and the second at 2 f 0 .As showed in Fig. 3(b), an increase in absorption (damping term) produces a slight redshift of the maximum around 2 f 0 , and a decrease of the magnitude of the two peaks, confirmed by Fig. 3 mentum imparted to the object is weaker for higher absorption.This decrease of the momentum transfer with material losses is due to the fact that an increase of the damping term weakens the resonance.

Influence of the value of the plasma frequency
As the resonance of the sphere is around ε = µ = −1 when ω p = ω pe = ω pm go far √ 2ω 0 the total momentum imparted to the sphere decrease particularly when ω p is shifted toward the low frequency, Fig. 4(a), as for the low frequencies (ε, µ) are close to one.If we compare the ration between the maximum at the frequency 2 f 0 and the maximum at the null frequency is higher when the central frequency of the pulse correspond to the resonance of the sphere showing that the observation of oscillation of the force with the frequency of the pulse should be done at a resonance.

Conclusion
In conclusion we have presented a general framework based on the discrete dipole approximation (DDA) for the computation of optical force on arbitrary magnetodielectric, threedimensional objects in time domain.The principal advantage of the method is that only the scatterer and its immediate neighborhood need to be discretized, allowing analytic expressions of the incident fields to be used.From its DDA foundation our approach inherits the ability to handle any material with a linear response, including dispersive, anisotropic and/or lossy magnetodielectric scatterers.

Fig. 1 .
Fig. 1.Spectral (a) and time (b) profiles of the incident field.

Fig. 2 .
Fig. 2. (a) Total force (solid line) versus time and its different contribution, i.e.F pm = F pe (dashed line), F hm = F he (dot-dashed line).(b) Total momentum imparted to the object and its contribution associated.(c) Total momentum imparted to the object by the pulse versus the numbers of subunits to represent the sphere.(d) In solid line (dashed line) F recoil with N = 33552 (N = 4224) subunits to represent the object.(e) Momentum imparted to the object due to F recoil versus the time.(f) Momentum imparted to the object due to F recoil versus the numbers of subunits to represent the sphere.

Fig. 3 .Fig. 4 .
Fig. 3. (a) Total momentum imparted to the object and (b) spectrum of the force, for different values of the damping term Γ.