## Abstract

The time-average energy density of the optical near-field generated around a metallic sphere is computed using the finite-difference time-domain method. To check the accuracy, the numerical results are compared with the rigorous solutions by Mie theory. The Lorentz-Drude model, which is coupled with Maxwell’s equation via motion equations of an electron, is applied to simulate the dispersion relation of metallic materials. The distributions of the optical near-filed generated around a metallic hemisphere and a metallic spheroid are also computed, and strong optical near-fields are obtained at the rim of them.

© 2007 Optical Society of America

## Corrections

Takashi Yamaguchi and Takashi Hinata, "Optical near-field analysis of spherical metals: Application of the FDTD method combined with the ADE method: errata," Opt. Express**16**, 4375-4375 (2008)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-16-6-4375

## 1. Introduction

Improving optical resolution beyond the diffraction limit is a subject of considerable research effort due to its potential impact upon many diverse areas of science and technology. In near-field optical recording, the data bits are written and read by using optical near-fields [1,2]. When light is introduced into a nanometer-scale object, such as a subwavelength aperture or a scatterer, a localized electromagnetic field (the optical near-field) is generated near the object. The energy intensity and the distribution of the optical near-field depend on the object’s shape, size, and dispersive material.

Matsumoto et al. [3,4] investigated the intensity distribution of the optical near-field generated by a probe with a wedge-shaped metallic plate using the finite-difference time-domain (FDTD) method when a plane wave was incident on the metallic plate. The FDTD method [5,6] is a very useful technique to analyze the electromagnetic fields in both the near field and the far field scattered from a complex geometrical object. When we apply the FDTD method to obtain the scattered fields from a complex geometrical object made of dispersive material, the following considerations are important:

(A-1) Verification of the accuracy of the FDTD results,

(A-2) Coupling method of the dispersion relation with Maxwell’s equation.

Plane-wave scattering by a sphere has been rigorously solved by matching boundary conditions (Mie theory). The solutions of Mie theory can be used for verification of the FDTD results. Challener et al. [7] reported comparisons between the FDTD results and the Mie solutions for a plane-wave incident upon a sphere with the Debye or the Lorentz dispersion. The differences between our paper and Ref. [7] are as follows:

(B-1) The distribution maps of the total electric field computed by the FDTD method were compared with those of the Mie solution. However, the dependence of the accuracy on the cell size was not investigated.

(B-2) The Lorentz-Drude model (LD model), which includes the effects of a free electron and bound electrons to describe the permittivity of metals, was not used.

(B-3) The energy density near the object was not investigated.

Zhu et al. [8] reported comparisons between the FDTD results and the Mie solutions for a gold sphere using a commercial FDTD package (XFDTD). The differences between our paper and Ref. [8] are as follows:

(C-1) The gold sphere was modeled by the modified Debye model.

(C-2) The incident wavelength and the size of the sphere were fixed at 660 nm and 400 nm, respectively.

(C-3) The total electric field was described only on the surface of the sphere.

Tamaru and Miyano [9] investigated the far-field scattering intensity of a silver sphere and a gold spheroid. The differences between our paper and Ref. [9] are as follows:

(D-1) It was described that about one tenth of a wavelength was a sufficient cell size for propagating waves. However, our numerical results show that its result is not valid for the energy density near the object.

(D-2) Silver was characterized by the Drude model (D model), and gold by the D model plus one term of the Lorentz model (L model). In contrast, we use the D model plus five terms of the L model [10] as the LD model to approximate the experimental data in a wide spectral range of light.

(D-3) The dispersion relation of the metals was combined with Maxwell’s equation using the recursive convolution (RC) method in Ref. [9]. We employ an auxiliary differential equation (ADE) method [6,11], which has higher accuracy than the RC method and needs no complex-number arithmetic.

Kashiwa and Fukai [11] analyzed the propagation characteristics in a parallel-plate waveguide filled with an L-model dispersive medium using the FDTD method. In their ADE method, the L model was combined with the constitutive equation by the polarization. In our formulation, the effect of dispersive media is added to Maxwell’s equation through the current to consider the polarization current and the conduction current in the LD-model media.

In this paper, we analyze the time-average energy density of the optical near-field generated by a plane wave incident upon a gold sphere to check the accuracy of the FDTD results [12]. As applications, we also compute the energy density near a gold and a silver hemisphere and spheroid. From the results, we clarify the spatial distributions of the energy density and the wavelength of the incident light providing the maximum energy density of the optical near-field.

The convolutional perfectly matched layer (CPML) [13], which is an applicable PML to anisotropic media and dispersive media [13–17], is used as an absorbing boundary condition to truncate the physical region.

## 2. Formulation

#### 2.1 Optical near-field around a small sphere

Mie theory is applicable to an arbitrarily-sized sphere [18,19]. We used the Mie solution as the exact one to verify the FDTD results. In this paper, the incident electric field *E ^{(i)}* is given by:

where *k*
_{0}=2*π/λ(λ*: incident light wavelength), and *E*
_{0} is the magnitude of the applied electric field. *ux* is the unit vector of the x-direction. When the metallic sphere is small (*k _{0}a*≪1), the electromagnetic fields around the sphere (

*r*>

*a*)

*E*^{(s)},

*Hθ*and

^{(s)}*H*are proportional to the electric dipole moment

_{ϕ}^{(s)}*p0*defined by

The electric dipole moment *p0* is enhanced when the light frequency satisfies the following equation:

where Re[·] denotes the real part of [·]. However, the electromagnetic fields around a nonspherical object are not enhanced when Eq. (3) is satisfied. The enhancement of *p0* depends on not only the frequency and the real part of *ε _{r}(ω)* but also the shape of the object.

#### 2.2 Dispersion relation of metals in the frequency domain

A useful dispersion model for metals such as gold, silver, or aluminum is the LD model. In the time-harmonic case, the relative permittivity of metals is expressed by the LD model as follows [10]:

where *ω _{p}* is the plasma frequency,

*K*in the symbol

*Σ*the number of oscillators with frequency

^{K}_{j=1}*ω*, and collision frequencies

_{j}*ν*and

_{0}*ν*and

_{j}. A_{0}*A*are constants which depend on the material [10].

_{j}When we set *A _{j}*=0 (

*j*=1,2,…,

*K*) in Eq. (5),

*ε*(

_{r}*ω*) becomes the relative permittivity of the D model, and the material is composed of free electrons. When

*A*is zero,

_{0}*ε*becomes the relative permittivity of the L model, and the material consists of bound electrons. Figure 1 shows the relative permittivity of gold parameterized by the LD model (

_{r}(ω)*K*=5) [10,20,21] against wavelengths. The solid curves (——), the dashed curves (–––-), and the dashed-dotted curves (–.–.–) are the numerical results of the LD model, the D model, and the L model, respectively. The circles (○ ○ ○) are the experimental data from Johnson and Christy [22]. As shown in Fig. 1, the numerical results of the LD model in Eq. (4) and the experimental data are in good agreement for the wavelength range from 300 nm to 800 nm.

#### 2.3 FDTD formulation for the LD model

The ADE method is used to couple the LD model with Maxwell’s equation [23]. Discretizing the polarization equation derived from the motion equation of an electron, we have the finite-difference equations for the current and the polarization as follows:

where Δ*t* is the increment of the time domain. The finite difference equations for the electromagnetic field in the metal are expressed as:

The computational order is **H**^{n+1/2}, * E^{n+1}*,

*and*

**J**_{j}^{n+1}*. If the medium is characterized by the D model,*

**P**_{j}^{n+1}*γ*=0 and

_{j}*vanishes from Eqs. (6) and (9). The stability requirement for the difference equations is*

**P**_{j}^{n}*cΔt*<1/{1/(Δ

*x*)

^{2}+1/(Δ

*y*)

^{2}+1/(Δ

*z*)

^{2}}

^{1/2}as if the medium were free space, where

*c*is the velocity of light, Δ

*x*, Δ

*y*, and Δ

*z*are the increments of the spatial coordinates.

Since the scattered electric and magnetic fields cannot be separated in three-dimensional problems, the electromagnetic energy density should be considered rather than only the electric field intensity. In near-field optical recording, the distribution and the intensity of the energy density are more important than the scattering cross-section or the scattered efficiency. We investigate the normalized time-average energy density (the energy density) defined by:

where * E* and

*are the total electromagnetic field outside the object,*

**H***and*

**E**^{(i)}*are the incident electromagnetic field (incident light). Each component of the electromagnetic field is located at different positions in the Yee cell, and there is a half-time-step difference between*

**H**^{(i)}**and**

*E***. We take space and time averages as follows:**

*H*The other components are calculated in the same manner. To obtain 〈*w*〉* _{N}*, the total time-average energy density in a steady state, we examined the following two methods;

(E-1) 〈*w*〉* _{N}* is computed using the steady state of

*and*

**E***.*

**H**(E-2) The total energy density is computed every time step. And the steady state of that is obtained as 〈*w*〉* _{N}* [12].

From our numerical experiments, we found that (E-1) gives us more accurate results than (E-2). Thereby we use (E-1) in this paper.

## 3. Numerical results

Figure 2 shows the three kinds of metallic objects that we analyze.

(a) A sphere with a radius *a*,

(b) A hemisphere with the spherical surface in the region *z* <0 and the flat surface on the *x-y* plane,

(c) A spheroid with a radius *a* in the *x*-direction and a radius *b* in the *y*-direction (*a* > *b*). In Fig. 2, the symbol O is the origin of the coordinate system and the symbol *P* is the observation point.

#### 3.1 Comparisons of the FDTD results with Mie solutions

A gold sphere with a radius of *a*=50.0 nm is analyzed to compare the numerical results by the FDTD method with the Mie solutions.

Figure 3 shows the energy density 〈*w*〉* _{N}* near the sphere as a function of the distance

*r*(>

*a*). The observation points are on the

*x*-axis (

*θ*=90°,

*ϕ*=0°). In Fig. 3, the circles (○ ○ ○) and the rectangles (□ □ □) indicate the FDTD results with the cell size Δ=Δ

*x*=Δ

*y*=Δ

*z*=2.5 nm and Δ=1.25 nm, respectively. The solid line (——) is the Mie solution. The incident light wavelength

*λ*is 550 nm. In the range of 57.5 nm <

*r*<70.0 nm, the relative error defined by |Mie solution — FDTD result|/|Mie solution| is less than 2.5 % when Δ=2.5 nm. However, the accuracy of the FDTD results decreases as the observation point r approaches the surface of the sphere. The reasons for that are caused by:

(F-1) The sphere is modeled by cubic cells,

(F-2) The average procedures in Eq. (13) are used to compute the energy density. When Δ=1.25 nm, the relative error of the FDTD result is 2.0 % at *r*=55.0 nm.

Table 1 shows the cell-size dependence of the FDTD results. The cell size Δ is varied from 1.25 nm to 5.0 nm. When Δ=1.25 nm, the relative error of 〈*w*〉* _{N}* is less than 2.1 % at both

*r*=55.0 nm and

*r*=60.0 nm. We used a multiprocessor (Opteron 270, 2.0 GHz) workstation to solve the problem. When Δ=1.25 nm, the computation time was about 150 hours and the required memory size was about 1.6 GB. When Δ=2.5 nm, they were 20 hours and 310 MB, respectively.

Gold, *a*=50 nm, *λ*=550 nm, *θ*=90°, *ϕ*=0°.

In Fig. 4, we show the relation between the energy density 〈*w*〉* _{N}* near the gold sphere and the incident light wavelength

*λ*. The solid line is the Mie solution. The circles and the rectangles indicate the FDTD results when Δ=2.5 nm and 5.0 nm, respectively. These results tell us the following:

(G-1) The energy density 〈*w*〉* _{N}* of both the FDTD results and the Mie solution reach the maximum at

*λ*=550 nm.

(G-2) Equation (3) is not exactly valid when *a*=50.0 nm, because Re[*ε _{r}*(

*ω*)]=-5.37 at

*λ*=550 nm.

(G-3) The relative errors of the FDTD results are less than 3.8 % when Δ=2.5.

Figure 5 illustrates the energy density 〈*w*〉* _{N}* near the gold sphere against the angle

*θ*. The observation points are fixed at 10 nm from the surface of the sphere and laid on the

*x-z*plane (

*ϕ*=0°). The incident light wavelength is 550 nm. The rectangles, the circles, and the filled circles show the FDTD results with the radius

*a*=20.0 nm and the cell size Δ=1.0 nm,

*a*=50.0 nm and Δ=2.0 nm, and

*a*=100 nm and Δ=2.0 nm, respectively. The dotted lines are the Mie solutions. The energy density 〈

*w*〉

*reaches its maximum about at*

_{N}*θ*=90° when

*a*=20.0 nm. If

*k*≪1 (

_{0}r*r*≥

*a*), 〈

*w*〉

*is given by:*

_{N}Equation (14) shows that the energy density 〈*w*〉* _{N}* attains the maximum at

*θ*=90° and

*ϕ*=0°. The maximum of 〈

*w*〉

*generated by the 50-nm sphere is larger than that by the 20-nm and the 100-nm ones. Since*

_{N}*k*=2

_{0}a*π a/λ≈*0.57 is not much smaller than unity when

*a*=50.0 nm and

*λ*=550 nm, Eq. (14) is inapplicable and the effect of size must be considered. The FDTD results are in good agreement with the Mie solutions. The relative errors at the maximum of 〈

*w*〉

*are 3.5 % for*

_{N}*a*=20.0 nm when Δ/

*λ*=1.818×10

^{-3}, 0.86 % for

*a*=50.0 nm and 0.58 % for

*a*=100.0 nm when Δ/

*λ*=3.64×10

^{-3}.

#### 3.2 Optical near-field around a hemisphere and a spheroid

Figure 6 shows the spatial distribution of the energy density 〈*w*〉* _{N}* generated around a silver hemisphere with

*a*radius a=50.0 nm (see Fig. 2 (b)). The observation plane corresponds to the

*x-z*plane at

*y*=0. The incident light wavelength

*λ*is 500 nm. As shown in Fig. 6, a strong optical near-field is generated at the rim of the hemisphere. The maximum of 〈

*w*〉

*appears around at*

_{N}*θ*=90°,

*ϕ*=0° and

*θ*=90°,

*〈*=180°which correspond to the direction of the polarization of the incident electric filed. Silver, λ=500 nm, a=50 nm, Δ=2.0 nm, y=0 nm.

Figure 7 shows the relation between the energy density 〈*w*〉* _{N}* and the incident light wavelength

*λ*. The observation point is placed at 10 nm from the surface of the hemisphere and on the

*x*-axis (

*θ*=90°,

*ϕ*=0°). The solid curve (——), the dashed curve (–––-), and the circles (○ ○ ○) represent the energy density 〈

*w*〉

*generated by the metallic hemisphere made of gold, silver, and perfect conductor, respectively. The following features of those curves will be noted:*

_{N}(H-1) The differences from the energy density 〈*w*〉* _{N}* for the gold sphere:

(1) The maximum of 〈*w*〉* _{N}* by the hemisphere is about 2.0 times as large as that by the sphere (see Fig. 4).

(2) The light wavelength providing the maximum of 〈*w*〉* _{N}* for the hemisphere is

*λ*=580 nm, which is shifted toward longer than that for the sphere.

(H-2) The maximum of 〈*w*〉* _{N}* for the silver hemisphere is about 2.1 times as large as that for the gold one. The energy loss in silver becomes smaller than that in gold because the imaginary part of relative permittivity of silver is Im[

*ε*(

_{r}*ω*)]=0.93 at

*λ*=580 nm, whereas that of gold is Im[

*ε*(

_{r}*ω*)]=2.15.

(H-3) Any strong energy density cannot be obtained around the perfect conductor, because the dipole moment is not enhanced by the incident light.

Figure 8 shows the spatial distribution of 〈*w*〉* _{N}* around a silver spheroid with radii

*a*=50.0 nm and

*b*=40.0 nm (see Fig. 2 (c)). The observation plane corresponds to the

*x-z*plane at

*y*=0. The incident light wavelength

*λ*is 500 nm. The intensity of 〈

*w*〉

*increases at the projections of the spheroid as well as the hemisphere.*

_{N}Figure 9 shows the relation between the energy density 〈*w*〉* _{N}* and the incident light wavelength

*λ*. The observation point is placed at 10 nm from the surface of the spheroid and on the

*x*-axis (

*θ*=90°,

*ϕ*=0°). The solid curve, the dashed curve, and the circles represent the energy density 〈

*w*〉

*generated by the metallic spheroid made of gold, silver, and perfect conductor, respectively. The distribution of 〈*

_{N}*w*〉

*by the spheroid attains the maximum at*

_{N}*θ*=90° and

*ϕ*=0°, 180°, because the major axis of the spheroid corresponds with the polarized direction of the incident light. The results obtained from Fig. 9 are as follows:

(I-1) The differences from the energy density 〈*w*〉* _{N}* for the gold sphere:

(1) The maximum of 〈*w*〉* _{N}* by the spheroid is about 1.5 times as large as that by the sphere (see Fig. 4).

(2) The light wavelength providing the maximum of 〈*w*〉* _{N}* for the spheroid is

*λ*=570 nm, which is shifted toward longer than that wavelength for the sphere.

(I-2) The maximum of 〈*w*〉* _{N}* for the silver spheroid is about 1.6 times as large as that for the gold one. The reason is the same as (H-2).

(I-3) The maximum of 〈*w*〉* _{N}* for the silver spheroid is smaller than that for the silver hemisphere.

In Fig. 10, we show the relation between the energy density 〈*w*〉* _{N}* and the incident light wavelength

*λ*for the gold spheroids with

*b/a*=0.4, 0.6, 0.8, and 1.0. The observation point and the cell size are the same as Fig. 9 except the cell size is 2.0 nm when the radius

*b*is 20.0 nm. From these results, the maximum of 〈

*w*〉

*increases and appears at a longer light wavelength as*

_{N}*b/a*becomes smaller.

## 4 Conclusions

In this paper, we computed the normalized time-average energy density (energy density) of the optical near-field generated by a metallic sphere, a metallic hemisphere, and a metallic spheroid using the finite difference time domain (FDTD) method. An auxiliary differential equation method was applied to couple the Lorentz-Drude model with Maxwell’s equation. To verify the accuracy of the energy density near a gold sphere computed by the FDTD method, we compared the FDTD results with the Mie solutions. The concluding remark for the comparison is as follows:

(1) The energy density computed by the FDTD method was in good agreement with that of the Mie solution.

(2) When the cell size was 2.5 nm, the relative errors between the FDTD results and the Mie solutions were within 3.8 %.

(3) As the observation point was closer to the surface of the sphere, the accuracy of the FDTD results decreased. In that case, we must make the cell sizes smaller.

As applications, we investigated the energy density around a hemisphere and a spheroid, which were made of gold and silver. It was found from our numerical analysis that:

(1) The maximums of the energy density near the silver hemisphere and the silver spheroid were larger than those near the gold ones, respectively.

(2) The maximum of the energy density near the silver hemisphere was larger than that near the silver spheroid.

(3) The light wavelengths providing the maximum of the energy density near the hemisphere and the spheroid were shifted toward longer than those near the sphere.

(4) The maximum of the energy density near the spheroid increased and appeared at a longer light wavelength as the ellipticity of the spheroid (=*b/a*) becomes smaller.

We are now investigating shapes which produce the large energy density near the object.

## Acknowledgments

The authors would like to thank Professor Tsuneki Yamasaki of Nihon University for his helpful suggestions.

## References and links

**1. **M Ohtsu and K. Kobayashi, *Optical Near Field* (Springer-Verlag, Berlin, 2004).

**2. **M. Ohtsu (Ed.), *Progress in Nano-Electro-Optics III-Industrial Applications and Dynamics of the Nano-Optical System* (Springer, 2004). [PubMed]

**3. **T. Matsumoto, T. Shimano, H. Saga, H. Sukeda, and M. Kiguchi, “Highly efficient probe with a wedge-shaped metallic plate for high density near-field optical recording,” J. Appl. Phys. **95**, 3901–3906 (2004). [CrossRef]

**4. **T. MatsumotoM. Ohtsu, ed. (Springer, 2004).

**5. **T. Uno, *Finite Difference Time Domain Method for Electromagnetic Field and Antennas*, in Japanese (Corona Publishing Co., Ltd, 1998).

**6. ** A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method, Third Edition* (Artech house, 2005).

**7. **W. A. Challener, I. K. Sendur, and C. Peng, “Scattered field formulation of finite difference time domain for a focused light beam in dense media with lossy materials,” Opt. Express **11**, 3160–3170 (2003). http://www.opticsinfobase.org/abstract.cfm?URI=oe-11-23-3160 [CrossRef] [PubMed]

**8. **R. J. Zhu, J. Wang, and G. F. Jin, “Mie scattering calculation by FDTD employing a modified Debye model for Gold material,” Optik **116**, 419–422 (2005). [CrossRef]

**9. **H. Tamaru and K. Miyano, “Localized Surface Plasmon Resonances of Metal Nanoparticles: Numerical Simulations and Their Experimental Verification,” Kogaku Japanese Journal of Optics **33**, 165–170 (2004).

**10. **A. D. Rakić, A. B. Djurišić, J. M. Elazar, and M. L. Majewski, “Optical Properties of Metallic Films for Vertical-Cavity Optoelectronic devices,” Appl. Opt. **37**, 5271–5283 (1998). [CrossRef]

**11. ** T. Kashiwa and I. Fukai, “A treatment by FDTD method of dispersive characteristics associated with electronic polarization,” Microwave Optics Tech. Lett. **3**, 203–205 (1990). [CrossRef]

**12. **T. Yamaguchi, T. Yamasaki, and T. Hinata, “FDTD Analysis of Optical Near-field Around a Metallic Sphere,” in *The Papers of Technical Meeting on Electromagnetic Theory* (IEE, Japan, 2006), pp. 143–148.

**13. **J. A. Roden and S. D. Gedney, “Convolutional PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media,” Microwave Optical Tech. Lett. **27**, 334–339 (2000). [CrossRef]

**14. **J. De Moerloose and M. A. Stuchly, “Behavior of Berenger’s ABC for evanescent waves,” IEEE Microwave Guided Wave Lett. **5**, 344–346 (1995). [CrossRef]

**15. ** J. P. Bérenger, “Improved PML for the FDTD Solution of Wave-Structure Interaction Problems,” IEEE Trans. Antennas Propag. **45**, 466–473 (1997). [CrossRef]

**16. **J. P. Bérenger, “Evanescent Waves in PML’s: Origin of the Numerical Reflection in Wave-Structure Interaction Problems,” IEEE Trans. Antennas Propag. **47**, 1497–1503 (1999). [CrossRef]

**17. ** T. Yamaguchi, “Numerical Analysis of Pulse Reflection from Anisotropic Dielectric Layer,” Special Issue of Nihon Univ. CST 2006 Annual Conf. Short Note **1**,113–116 (2007).

**18. **J. A. Stratton, *Electromagnetic Theory* (McGraw-Hill, 1941).

**19. **T. Hosono, *The Foundation of Electromagnetic Wave Theory*, in Japanese (Shoko-do, 1973).

**20. **M. A. Ordal, L. L. Long, R. J. Bell, S. E. Bell, R. R. Bell, R. W. Alexander. Jr., and C. A. Ward, “Optical properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Pt, Ag, Ti, and W in the infrared and far infrared,” Appl. Opt. **22**, 1099–1119 (1983). [CrossRef] [PubMed]

**21. **J. B. Judkins and R. W. Ziolkowski, “Finite-difference time-domain modeling of nonperfectly conducting metallic thin-film gratings,” J. Opt. Soc. Am. A **12**, 1974–1983 (1995). [CrossRef]

**22. **P. B. Johnson and R.W. Christy, “Optical Constants of the Noble Metals,” Phys. Rev. B **6**, 4370–4379 (1972). [CrossRef]

**23. **M. Bordovsky, P. Catrysse, S. Dods, M. Freitas, J. Klein, L. Kotacka, V. Tzolov, I. Uzunov, and J. Zhang, “Waveguide design, modeling, and optimization - from photonic nano-devices to integrated photonic circuits,” Proc. SPIE **5355**, 65–80 (2004). [CrossRef]