## Abstract

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

## 1. 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*.

## 2. 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 { (**r**, *t*), ℋ (**r**,*t*)}. Let 𝒫(**r**,*t*) [ℳ(**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]:

*i*,

*j*, and

*k*stand for either

*x*,

*y*or

*z*and

*∂*denotes the derivative with respect to time. This is a generalized form of the Lorentz force [21, 22]. Using Maxwell’s equations the

_{t}*i*-th component of the electric and magnetic forces can be written as:

*ε*is the Levi-Civita tensor. However, the above expression for the force is incomplete. We need to add the radiative reaction force,

_{ijk}*i.e.*the flux of momentum lost by the particle when it radiates electromagnetic waves [23, 24]

^{he}(

**r**,

*t*)] and third terms [ℱ

^{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 [ℱ

^{pe}(

**r**,

*t*)] and fourth terms [ℱ

^{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 time-dependent electromagnetic wave with envelop ℐ(*t*) (Fourier transform *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 (**r**, *t*) = 𝒢^{−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:

*N*subunits. We can also define the momentum imparted to the object by the electromagnetic field as ${\mathcal{Q}}_{i}\left(t\right)={\int}_{0}^{t}{\mathcal{F}}_{i}\left(t\right)\text{d}t$ where we have assumed that the origin of time is chosen such that the electromagnetic fields at the object are zero for

*t*< 0.

## 3. 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 free-space 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 6*N* × 6*N* 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].

## 4. 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*= cos(

_{y}*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:

^{he}and ℱ

^{hm}but as seen in Eqs. (11)–(12) the spatial derivatives of the electromagnetic field introduce a term proportional to sin(2

*kz*– 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 differentiation 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:

## 5. Time-dependant optical force on a discretized object

*5.1. 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:

*ω*

^{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:

*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.*${\omega}^{\text{pe}}={\omega}^{\text{pm}}={w}_{0}\sqrt{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 time-harmonic 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 ℱ

^{pm}(

*t*) = ℱ

^{pe}(

*t*) and ℱ

^{hm}(

*t*) = ℱ

^{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.

#### 5.2. 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(a) where the total momentum 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.

#### 5.3. Influence of the value of the plasma frequency

As the resonance of the sphere is around *ε* = *μ* = −1 when *ω ^{p}* =

*ω*=

^{pe}*ω*go far $\sqrt{2}{\omega}_{0}$ the total momentum imparted to the sphere decrease particularly when

^{pm}*ω*is shifted toward the low frequency, Fig. 4(a), as for the low frequencies (

^{p}*ε*,

*μ*) 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.

## 6. 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, three-dimensional 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.

## 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). [CrossRef]

**2. **B. T. Draine, “The discrete-dipole approximation and its application to interstellar graphite grains,” Astrophys. J. **333**, 848–872 (1988). [CrossRef]

**3. **P. C. Chaumet, A. Rahmani, and G. W. Bryant, “Generalization of the coupled dipole method to periodic structure,” Phys. Rev. B **67**, 165,404–5 (2003). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**6. **A. Rahmani and G. W. Bryant, “Spontaneous emission in microcavity electrodynamics,” Phys. Rev. A **65**, 033,817–12 (2002). [CrossRef]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**9. **B. T. Draine and J. C. Weingartner, “Radiative Torques on Interstellar Grains: I. Superthermal Spinup,” Astrophys. J. **470**, 551–565 (1996). [CrossRef]

**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. A **18**, 1944–1953 (2001). [CrossRef]

**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). [CrossRef]

**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). [CrossRef]

**13. **A. Rahmani and P. C. Chaumet, “Optical Trapping near a Photonic Crystal,” Opt. Express **14**, 6353–6358 (2006). [CrossRef] [PubMed]

**14. **P. C. Chaumet and C. Billaudeau, “Coupled dipole method to compute optical torque: Application to a micro-propeller,” J. Appl. Phys. **101**, 023,106–6 (2007). [CrossRef]

**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). [CrossRef]

**16. **P. C. Chaumet, K. Belkebir, and A. Rahmani, “Coupled-dipole method in time domain,” Opt. Express **16**, 20,157–20,165 (2008). [CrossRef]

**17. **P. C. Chaumet, K. Belkebir, and A. Rahmani, “Optical forces in time domain on arbitrary objects,” Phys. Rev. A (2010). [CrossRef]

**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). [CrossRef]

**19. **B. D. H. Tellegen, “Magnetic-Dipole models,” Am. J. Phys. **30**, 650–652 (1962). [CrossRef]

**20. **L. Vaidman, “Torque and force on a magnetic dipole,” Am. J. Phys. **58**, 978–983 (1990). [CrossRef]

**21. **M. Mansuripur, “Radiation pressure and the linear momentum of the electromagnetic field in magnetic media,” Opt. Express **15**, 13,502–13,518 (2007). [CrossRef]

**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). [CrossRef]

**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. Express **17**, 2224–2234 (2009). [CrossRef] [PubMed]

**26. **P. C. Chaumet and A. Rahmani, “Efficient iterative solution of the discrete dipole approximation for mag neto-dielectric scatterers,” Opt. Lett. **34**, 917–919 (2009). [CrossRef] [PubMed]

**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). [CrossRef]

**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). [CrossRef]

**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). [CrossRef] [PubMed]

**32. **J. R. Arias-González 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). [CrossRef]

**33. **R. W. Ziolkowski, “Pulsed and CW Gaussian beam interactions with double negative metamaterial slabs,” Opt. Express **11**, 662–681 (2003). [CrossRef] [PubMed]