Using the scattered field finite difference time domain (FDTD) formalism, equations for a plane wave incident from a dense medium onto lossy media are derived. The Richards-Wolf vector field equations are introduced into the scattered field FDTD formalism to model an incident focused beam. The results are compared to Mie theory scattering from spherical lossy dielectric and metallic spheres.
© 2003 Optical Society of America
The finite difference time domain (FDTD) method is popular for analyzing electromagnetic fields in situations that are too complex for analytical methods. Application of FDTD to electromagnetics was first described by Yee . Originally the FDTD method was applied to microwave scattering problems, but it has also become a very useful technique at optical frequencies for calculating electric fields in both the near field and far field which are scattered from nanoparticles, apertures , or antennas . The technique has been shown to correctly model such optical effects as surface plasmon resonances and to provide reliable results at optical frequencies. The FDTD method may be divided into two broad categories generally called the total field approach and the scattered field approach. The total field approach is described in detail in Taflove and Hagness .
The scattered field formalism for FDTD is presented by Kunz and Luebbers . In this formalism the incident electric field is specified throughout the computation space for the time duration of the calculation through a calculation which is separate from FDTD algorithm. This is particularly convenient if the incident field can be determined analytically. Dispersion errors in the incident field are not a factor for the scattered field approach. Another advantage of the scattered field approach is that the absorbing boundary conditions must only handle the scattered field rather than the total field, which is usually easier. Finally, in the case of a focused beam incident upon a small scattering object, the computation space in the scattered field formalism can be limited to just the volume surrounding the object. In the total field formalism the incident field must either be specified over all sides of the computation space, or if the incident field is just specified over one surface of the computation space this surface must be large enough to include the entire region over which there is an appreciable incident field strength, which may be considerably larger than the region in the immediate neighborhood of the scattering object and require correspondingly longer computation times.
The scattered field approach is straightforward for modeling dielectrics and perfect conductors, but is more complicated than the total field approach when modeling materials that exhibit a frequency-dependent permittivity or conductivity. As discussed in Chapter 8 of Ref. , the scattered field formalism can be applied to frequency-dependent materials with a Debye or Lorentz dispersion in a computationally efficient manner. This reference considers the case of an electric field plane wave pulse. In this paper the discussion of Ref.  is extended to include lossy materials that are illuminated with a plane wave and a focused optical beam. The case when the field is incident within an optically dense material is also considered.
2. Plane wave incident on a Debye material
A simple model of a frequency-dependent dielectric function which satisfies the Kramers-Kronig causality relations is the first order Debye equation . The FDTD approach for a plane wave incident upon a material whose dielectric constant is described by the Debye equation is considered in this section. The Debye equation for the dielectric constant is
where ε∞ is the infinite frequency dielectric constant, εs is the zero frequency dielectric constant, ω is the angular frequency, t0 is the frequency-independent relaxation time, and χ(ω) is the frequency-dependent susceptibility. The Fourier transform of the susceptibility gives the corresponding time-dependent susceptibility function ,
where U(t) is the Heavyside unit step function. Following the procedure and notation in Ref. , the recursive equation for updating the scattered electric field at time step n+1 is obtained,
where E⃑ i is the incident field, H⃑ s is the scattered H field, σ is the dc conductivity,
after the substitution ξ=t-Λ.
For nonmagnetic materials with a permeability µ0 of free space, the recursive equation for updating the magnetic field is much simpler,
When modeling absorptive materials such as metals or lossy dielectrics, the time step Δt is typically chosen to be slightly smaller than the Courant time step, tc, to ensure stability of the calculation . The Courant time step is determined by the dimensions of the Yee cell and the speed of light in vacuum. A value of Δt between 0.9·tc and 0.95·tc is generally adequate, while Δt > 0.95tc can often lead to divergences in the calculation. Using shorter time steps unnecessarily increases the computation time and can also lead to larger dispersion errors .
It is convenient to express the time dependence of the electric field as a complex quantity. The amplitude of the plane wave is zero for t≤0,
where k 0 is the wavevector of the plane wave and p is the point of calculation in the computation space. In polar coordinates,
The actual field amplitude can be obtained by choosing either the real or imaginary part of Eq. (9).
In the scattered field formulation of Eq. (3) a convolution integral of the incident field must also be computed. For a sinusoidal field,
In the limit of infinite time, the second term in the square bracket of Eq. (12) goes to zero, leaving only the steady state result from the first term. Therefore, one criterion for reaching the steady state result in the computation is t≫t0. In fact, since the calculation must necessarily be run until the steady state is reached, this second term can be eliminated from the calculation a priori which in turn helps to reduce the computation time.
Plane wave scattering by a metallic sphere can be computed analytically by Mie theory [6,7] (also, see appendix). The scattering from a 100 nm diameter silver sphere was calculated using both FDTD and Mie theory for comparison as shown in Fig. 1. The wavelength was 700 nm corresponding to a frequency of 4.28275×105 GHz. The bulk value for the refractive index of silver, n=0.14+i(4.523), was used . Debye parameters for this refractive index are σ= 8.67651×106(Ω m)-1, t0=6.29065×10-15 s, ε∞=1, and εs=-6163.4. The cell space was 200 by 200 by 200 cells, and each cell was a 2 nm cube. Second order Mur boundary conditions were used on the faces and first order Mur boundary conditions on the edges [5,9], although there are many other boundary condition options available which may further improve the accuracy of the calculation [4,10,11]. However, the scattered field FDTD formulation described in this paper is independent of the particular boundary conditions chosen. The incident field was polarized along the x-axis and propagated in the +z direction. The time step was 0.9 tc, and the calculation was run for 3000 time steps.
The scattering from a 1 µm lossy dielectric sphere with a refractive index n=3.0+i(1.5) in free space at a wavelength of 850 nm was calculated as a second example. The Debye parameters were σ=1.99621×105 (Ωm)-1, t0=1.62952×10-16 s, ε∞=10, and εs=6.3262. The cell space was 200×200×200 cells and each cell was a 20 nm cube. The time step was 0.9 tc for a total of 500 time steps. The result was again compared to Mie theory in Fig. 2.
3. Plane wave incident on a Lorentz material
A very similar formalism can be applied to dispersive materials that obey a Lorentz dispersion equation,
where ωp is the resonant frequency and Δp is the damping coefficient. The corresponding time domain susceptibility function is 
As described in Ref. , a complex susceptibility function can be defined,
If we make the following substitutions into Eq. (17)
then it is in the same form as that derived for Debye materials in Eq. (2). To use the Debye equation format, we calculate the complex values for χ0 and ψ n at each time step using Eqs. (4), (5), and (7), and insert the real part of these numbers into Eq. (3).
This procedure cannot be used for the convolution term in Eq. (3), however. Because the electric field is obtained by taking the real or imaginary part of a complex quantity, the convolution term must be calculated for the real susceptibility function of Eq. (13). The Lorentz susceptibility function can be rewritten as a sum of two terms which are both in the same form as the Debye susceptibility function,
By comparison to the Debye result for the convolution term in Eq. (12), one immediately obtains
Again, assuming that the calculation is properly run to the steady state result, the second term in each of the square brackets can be neglected, so the result simplifies to
A comparison can be made between the Debye and Lorentz FDTD calculations using the previous example of the 1 µm lossy dielectric sphere. In Fig. 3(a) the results are plotted in the xz plane for FDTD calculations using the Lorentz dispersion relation with the parameters ω0=2.27412×1015 Hz, δ=9.20695×1013 Hz, ε∞=1, and εs=2 corresponding to n=3.0+i(1.5) at a wavelength of 850 nm. The cell space is 200×200×200 with a 20 nm cubic cell. The calculation is run for 700 time steps with each step t=0.9 tc. The result is compared along the z axis to the Debye model and Mie theory in Fig. 3(b). In this example, the Debye model at only 500 time steps gives a better match to Mie theory than the Lorentz model at 700 time steps.
4. FDTD in dense media
In optical measurements it is often the case that the light propagates towards a scattering object within a dielectric medium other than air or vacuum. All of the FDTD equations have been formulated in terms of the permittivity ε0 of vacuum. If the computation space is embedded in a dense medium, ε0 is simply replaced in every equation by εm=n2·ε0, the permittivity of the medium. The Courant time step, however, is still calculated for free space for which the speed of light is maximized and the time step minimized to ensure convergence of the calculation even when parts of the computation space still include vacuum.
An example is that of a 100 nm air bubble with n=1.0 embedded in a medium with n=2. A plane wave is incident with a free space wavelength of 700 nm. The computation space is 200×200×200 cells and each cell is a 2.5 nm cube. The plane wave propagates in the +z direction and is polarized along the x axis. Each time step is 0.9 tc and 3000 time steps are calculated. The results are shown in Fig. 4.
A second example is for light incident from a dielectric medium with n=2 onto an embedded 100 nm silver sphere. The incident plane wave has a free space wavelength of 700 nm. The sphere is modeled with a Debye dispersion, and the parameters are σ=1.12658×107(Ωm)-1, t0=7.17096×10-15 s, ε∞=4, and εs=-9120.07. The FDTD computational space is 100×100×100 cells, with a 5 nm cubic cell. The time step is 0.9 tc and 1500 time steps are taken. The FDTD calculation is compared to Mie theory in Fig. 5.
5. Focused field FDTD formulation
Wolf and Richards presented a method for calculating the electric field semi-analytically near the focus of an aplanatic lens [12,13]. The total electric field in the neighborhood of the focus is given by 
where the time dependence eiωt is understood, a(θ,ϕ) is a function describing the incident polarization, p is the point of observation,
λ0 is the free space wavelength, , and ϕ=tan-1(yp/xp ).
Equation (26) for the focused field can be used directly within the scattered field formalism of FDTD. As was done for plane waves, the field is turned on throughout the computation space at time t=0. The incident field is thus determined by its steady state value as calculated by the Richards and Wolf theory multiplied by eiωt. The spatial dependence of the incident field need only be calculated once prior to beginning the FDTD calculation and can be used repeatedly for additional FDTD calculations on other scattering objects illuminated by the same incident focused field. One subtlety involved in calculating the incident field is that the x, y, and z coordinates of the electric field components within the Yee cell are offset in different directions from the origin of the cell . Therefore, the incident field must either be calculated separately for each component Ex, Ey, and Ez in the focused beam at each cell in the computation space, or the incident field can be calculated once for each cell and interpolated to provide the field amplitudes within the cells, which is the approach taken here. The time-dependent part of the convolution term is given by Eqs. (12) and (25).
An example for the focused beam similar to that of Fig. 1 with a 100 nm silver sphere in air is now considered. The Debye parameters of silver are the same as for Fig. 1. The incident beam at a wavelength of 700 nm is linearly polarized along the x-axis prior to focusing and propagates in the -z direction. The incident field amplitude is normalized to unity at the focus. It is focused at the center of the sphere and has a half angle of 60°. The computation space is 200×200×200 cells and each cell is a 2 nm cube. The steady state results were obtained after 3000 time steps where each step was 0.9 tc. The results are compared to Mie theory  in Fig. 6.
As a second example, a focused beam incident upon a 500 nm silver sphere is considered. The wavelength is 700 nm and the incident beam is polarized along the x axis prior to focusing. The same Debye parameters are chosen for silver. The half angle for the focused beam is still 60°. The cell space is 200×200×200 cells and each cell is a 5 nm cube. The steady state results were obtained after 1500 time steps where each step was 0.9 tc. The results are shown in Fig. 7.
A scattered field FDTD formalism has been developed for modeling plane waves or focused beams incident upon Debye or Lorentz scattering objects. The Richards-Wolf vector field formalism is used to compute the incident field for focused beams. The surrounding medium may be free space or a dielectric with n>1. Several examples have been calculated for beams incident upon metallic and lossy dielectric spheres and compared to Mie theory with good agreement.
Mie theory calculations were based on the derivation in Ref. . This reference presents the field equations for the scattered fields outside of the sphere but does not complete the derivation for the fields within the sphere. Continuing with the notation in this reference,
and within the sphere,
Calculations of a focused incident beam using Mie theory were performed by summing the scattered fields from incident plane waves. The incident plane waves were spaced at 2° intervals along the polar direction and 6° intervals along the azimuthal direction. For a specific incident plane wave, the coordinates of the point (x,y,z) were first transformed into a coordinate system for which the plane wave was propagating along the z′-axis and polarized along the x′-axis. The transformation matrix is
After the scattered field was calculated at this point, the incident plus scattered field components were transformed back into the original coordinate system according to
and added to the incident plus scattered fields from the other plane waves.
This work was performed as part of the Information Storage Industry Consortium (INSIC) program in Heat Assisted Magnetic Recording (HAMR), with the support of the U. S. Department of Commerce, National Institute of Standards and Technology, Advanced Technology Program, Cooperative Agreement Number 70NANB1H3056.
References and links
1. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Prop. AP-14, 302–307 (1966).
2. W. A. Challener, T. W. McDaniel, C. D. Mihalcea, K. R. Mountfield, K. Pelhos, and I. K. Sendur, “Light Delivery Techniques for Heat-Assisted Magnetic Recording,” Jpn. J. Appl. Phys. 42, 981–988 (2003). [CrossRef]
4. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method 2nd ed. (Artech House, Boston, 2000).
5. K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics (CRC Press, Boca Raton, 1993).
6. G. Mie, “Beiträge zur optik truber medien, speziell kolloida ler metallösungen,” Ann. d. Physik 25, 377 (1908). [CrossRef]
7. M. Born and E. Wolf, Principles of Optics 5th ed. (Pergamon Press, Oxford, 1975), section 13.5.
8. D. W. Lynch and W. R. Hunter, “Silver (Ag)” in Handbook of the Optical Constants of Solids, E. D. Palik, ed. (Academic, San Diego, 1998).
9. G. Mur, “Absorbing boundary conditions for finite-difference approximation of the time-domain electromagnetic field equations,” IEEE Trans. Electromagn. Compat. 23, 1073–1077 (1981). [CrossRef]
10. Z. P. Liaoet al., “A transmitting boundary for transient wave analyses,” Scientia Sinica 271063–1076 (1984).
11. J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]
12. E. Wolf, “Electromagnetic diffraction in optical systems I. An integral representation of the image field,” Proc. Roy Soc. London Ser. A 253, 349–357 (1959). [CrossRef]
13. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanatic system,” Proc. Roy Soc. London Ser. A 253, 358–379 (1959). [CrossRef]
15. I. K. Sendur and W. A. Challener, “Interaction of focused beams with spherical nanoparticles,” to be published.