Open Access
Add to BibSonomyAbstract
An efficient finite-difference frequency-domain method is developed for calculating electromagnetic fields in the neighborhood of subwavelength dielectric and metallic structures. The method is used to investigate two-dimensional near-field and far-field patterns of a focused beam diffracted from an optical disk, specifically from a DVD (Digital Versatile Disk). It is shown that the polarization of illumination has a significant impact on diffraction patterns as expected and that scalar theory does not provide an accurate analysis of diffraction from a DVD.
©1998 Optical Society of America
1. Introduction
Developments in nano-technology and microelectronics fabrication have led to numerous engineered structures with sizes comparable to or smaller than the wavelength of visible light. In the area of optical data storage systems, the DVD, which is becoming increasingly popular, has a minimum pit length that is smaller than the laser wavelength. To study diffraction from such subwavelength structures, vector diffraction theory must be used because scalar theory becomes inaccurate.
The numerical calculation of rigorous electromagnetic diffraction is a broad field with many different approaches and applications. Two general numerical approaches to solving the Maxwell equations are the integral approach and the differential approach. Examples of some differential methods are finite-difference frequency-domain (FDFD), finite-difference time-domain (FDTD), and finite-element. These differential methods are well-suited to problems involving inhomogeneous media with complicated shapes, whereas the integral approaches tend to be better suited to open-space scattering or diffraction problems. Recently, numerical algorithms such as FDTD [1], boundary element method [2], and the integral method [3] have been used to investigate vector diffraction from optical disks.
In this paper, we describe an efficient FDFD method for calculating the electromagnetic fields in the neighborhood of dielectric and metallic structures. We use this method to model a DVD-ROM optical storage system and to calculate near fields and far fields of a two-dimensional DVD.
2. Numerical methods
Since we are interested in steady state phenomena, the electromagnetic field can be assumed to be time-harmonic and the FDFD method can be used to solve the time-harmonic Maxwell equations. In media that are non-magnetic, isotropic, and inhomogeneous, the source-free, time-harmonic Maxwell equations may be written as
where E⃗ is the electric field, H⃗ is the magnetic field, ∊ is the complex electric permittivity, ω is the angular frequency, and the time dependence e -iωt has been assumed. The information about material properties is included in ∊(y,z). Eliminating the magnetic field, the Maxwell equations become
Equation (2) is in the form of coupled Helmholtz equations and is the main equation we work with. With properly imposed boundary conditions and correctly chosen grid geometry, the divergence equation (3) is automatically satisfied [4]. To solve equation (2), we combine the techniques described below to develop a novel FDFD method.
2.1 Concus-Golub iteration
Concus and Golub have proposed a scheme that uses a fast cyclic-reduction solver to iteratively solve general Helmholtz equations [5]. The cyclic-reduction solver [6,7] explores the sparseness of the matrix form of the discrete Helmholtz equation and can handle the separable Helmholtz equation that describes the layered structure. This solver can produce fast direct solutions at O(3n 2 log n) for an n × n matrix and has a small memory requirement. In equation (2), each component has the form of a general Helmholtz equation
which can be solved with the Concus-Golub iteration method. Specifically, for the k-th iteration, we use the expression
It has been demonstrated [5] that the number of necessary iterations can vary dramatically depending on the function ∊(y,z). Usually, the smoother ∊(y,z) is, the faster the rate of convergence. However, the Concus-Golub iteration method does not always converge in all cases. For nonconverging cases, we use the conjugate gradient method.
2.2 Conjugate gradient method
Conjugate gradient methods are probably the most popular iterative techniques for solving systems of linear equations. These methods generate solutions to A ∙ x⃗ = b⃗ by minimizing quadratic functionals in Krylov subspaces, which are spanned by a series of vectors generated by repeated multiplication by A.
The Concus-Golub iteration can be used as a good preconditioner in conjunction with a robust general minimal residual (GMRES) conjugate gradient algorithm [8] to fully solve for electric fields to a very high accuracy [9]. From numerical experiments, it was found that the internal iterations of the Concus-Golub preconditioner do not need to fully converge in order to achieve optimal computing performance [9].
2.3 Radiation boundary condition
To solve the Maxwell equations in an infinite domain numerically, for example in diffraction or scattering problems, it is usually necessary to limit the computation space to a finite domain and to impose a radiation boundary condition on the outer boundary. The purpose of the radiation boundary condition is to absorb scattered or radiated fields so that free space is at least approximately simulated.
The Engquist-Majda [10] and Mur [11] radiation boundary conditions and their variations are often used in FDTD calculations. Although these boundary conditions are derived in the FDTD framework, it is straightforward to adapt them to the frequency domain. In this work, we use Higdon’s variation of the 2nd-order finite-difference radiation boundary condition [12]:
where α 1 and α 2 are the optimal angles for minimal reflection. We choose the ẑ direction to be the primary direction of wave propagation and the angles α1 and α 2 to be zero in order to minimize the reflection from the ẑ direction. This boundary condition is equivalent to the 2nd-order Engquist-Majda radiation boundary condition [12], but it is better suited for the cyclic-reduction method because it only requires the derivative in one direction.
3. Configuration
The surface of a DVD-ROM disk consists of a plastic substrate coated with aluminum with tracks of pressed pits (see Figure 1). We consider a simplified 2D model to simulate vector diffraction from these pits. When both the incident field and the geometry of the structure are two dimensional, the problem can be decomposed into two independent polarizations. We have TE polarization when the electric field is perpendicular to the plane of incidence and parallel to the pits and TM polarization when the electric field is in the plane of incidence.
The complex indices of refraction are n = 1.6 for the polycarbonate substrate and n = 1.5 + 7.8i for the aluminum reflection layer. The wavelength in vacuum of the incident laser light is 650 nm, which corresponds to a reduced wavelength in the polycarbonate substrate of 406 nm. The numerical aperture is NA = 0.6. The incident field is a 2D Gaussian beam that has a full width at half maximum (FWHM) of approximately 600 nm.
In the numerical calculation there are five pits on the disk. The computation box is set to 12 wavelengths (≈ 4.87 μm) in the y direction and 5 wavelengths (≈ 2.03 μm) in the z direction, with 40 grid points in each wavelength. The grid points are located in a Yee cell [13] in two-dimensional space. The unit of length in all figures is the reduced wavelength in the polycarbonate substrate, λs = 406 nm.
4. Near fields
Figure 2 shows movies of the amplitude and phases of the electric fields near the middle pit for the pit heights between 0.1 and 0.8 λs (≈ 40–320 nm). The total field is shown here, which includes both incident and diffracted fields. For TE illumination there is only an x̂ component, whereas for TM illumination there are both ŷ and ẑ components that are coupled together. Here we only show the ŷ component of the TM case, which is generally the dominant component for this polarization.
For TE illumination, there are strong peaks, not only directly in front of the pit, but also adjacent to the pit. With increasing pit height, a visible resonant mode emerges in the space between the pits. The TM near-field profiles are significantly different from the TE ones. Initially, for small pit heights, the amplitude of the primary component of the electric field is nearly the same for both polarizations. However, this similarity diminishes with increasing pit height.
Figure 2. Movies of near-field profiles of the x̂ component of the electric field for TE illumination and the ŷ component of the electric field for TM illumination. The color shows the phase and the amplitude of near-field profiles. White denotes amplitudes larger than 2.0. The white dotted line shows the interface between the substrate and the aluminum. The pit heights are between 0.1–0.8 λs. The unit of length is the wavelength in the substrate, λs = 406 nm. [Media 1] [Media 2]
5. Far fields
To calculate the far fields, we sample the near-field pattern inside the mesh and use a simple transform based on the far-zone approximation of the Green’s function. For example, for the Ex component,
The transform takes the transverse components of the scattered field along a plane above the surface of the disk and projects these components onto a circle in the far field.
Figure 3 shows the far-field intensities in the normal direction for TE illumination, TM illumination and scalar theory as a function of pit height. Both TE and TM cases differ from the predictions of scalar diffraction theory. The TE case exhibits a larger deviation, perhaps related to the emergence of resonant modes in the near-field.
The phase depth in the TE case is less than in scalar diffraction theory. A significant deviation in the TE case has been previously reported in Ref. 3. The TM case is similar to the scalar theory for small pit heights, but begins to deviate for heights larger than 0.3 λs (≈ 120 nm). The position of minima and maxima for TM illumination is roughly the same as for scalar theory.
Figure 3. Far-field intensities I = |E|2 in the normal direction as a function of pit height for the DVD. The intensities are normalized so that unity corresponds to no pit present. The solid line is for TE illumination, the dashed line for TM illumination, and the dotted line for scalar theory. The unit of length is the wavelength in the substrate, λs = 406 nm.
When pit height is around 0.6 λs in the TE case and 0.5 λs in the TM case, the far-field intensities return to close to unity and the near-fields in the substrate look similar to those without a pit at all. For pit heights greater than 0.5 λs, the intensity predicted by scalar theory is considerably smaller than vector theory. Our initial calculations indicate that this difference is due to a large contribution from evanescent waves in scalar theory.
In a usual DVD device, the incident field is circularly polarized. The results for circularly polarized illumination are approximately the average of the TE and TM cases.
6. Conclusion
We have developed a novel FDFD method for calculating the near-field and far-field diffraction of a focused laser beam incident on a two-dimensional DVD surface. We have found that the results of vector diffraction theory are significantly different from those of scalar diffraction theory and that TE and TM illuminations produce different diffraction patterns. The later phenomena has been experimentally demonstrated by Gerber and Mansuripur [14]. The deviation between TE and TM cases is related to the emergence of resonant modes in the space between pits. Although a three-dimensional model is necessary to rigorously analyze DVD and other high-density optical recording formats, our two-dimensional model is valid for longer pits and provides important insight into vector diffraction effects. A three-dimensional version of this algorithm is under development and can serve as a design tool for future optical storage systems. This method can also be applied to near-field calculations of other optical systems with subwavelength structures.
Acknowledgments
We thank J. H. Eberly for invaluable discussions and advice. This research was supported by the National Science Foundation (Grant PHY-94-15583) and Eastman Kodak Company. Computational facilities were supported by Eastman Kodak Company and the University of Rochester.
References
1. J. B. Judkins, C. W. Haggans, and R. W. Ziolkowski, “Two-dimensional finite-difference time-domain simulation for rewritable optical disk surface structure design,” Appl. Opt. 35, 2477–2487 (1996). [CrossRef] [PubMed]
2. M. Ogawa, M. Nakada, R. Katayama, M. Okada, and M. Itoh, “Analysis of scattering light from magnetic material with land/groove by three-dimensional boundary element method,” Jpn. J. Appl. Phys. 35, 336–341 (1996). [CrossRef]
3. D. S. Marx and D. Psaltis, “Optical diffraction of focused spots and subwavelength structures,” J. Opt. Soc. Am. A 14, 1268–1278 (1997). [CrossRef]
4. B.-N. Jiang, J. Wu, and L. A. Povinelli, “The origin of spurious solutions in computational electromagnetics,” J. Comput. Phys. 125, 104–123 (1996). [CrossRef]
5. P. Concus and G. H. Golub, “Use of fast direct methods for the efficient numerical solution of nonseparable elliptic equations,” SIAM J. Numer. Anal. 10, 1103–1120 (1973). [CrossRef]
6. B. L. Buzbee, G. H. Golub, and C. W. Nielson, “On direct methods for solving Poisson’s equations,” SIAM J. Numer. Anal. 7, 627–656 (1970). [CrossRef]
7. R. A. Sweet, “A cyclic reduction algorithm for solving block tridiagonal systems of arbitrary dimension,” SIAM J. Numer. Anal. 14, 706–720 (1977). [CrossRef]
8. Y. Saad and M. H. Schultz, “GMRES: a general minimal residual algorithm for solving nonsym-metric linear systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986). [CrossRef]
9. C. D. Dimitropoulos and A. N. Beris, “An efficient and robust spectral solver for nonseparable elliptic equations,” J. Comput. Phys. 133, 186–191 (1997). [CrossRef]
10. B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation of waves,” Math. Comput. 31, 629–651 (1977). [CrossRef]
11. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. EMC-23, 377–382 (1981). [CrossRef]
12. R. L. Higdon, “Absorbing boundary conditions for difference approximations to the multidimensional wave equation,” Math. Comput. 47, 437–459 (1986).
13. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. AP-14, 302–307 (1966).
14. R. E. Gerber and M. Mansuripur, “Dependence of the tracking performance of an optical disk on the direction of the incident-light polarization,” Appl. Opt. 34, 8192–8200 (1995). [CrossRef] [PubMed]
