Inverse Design of Plasmonic Structures with FDTD
11 Topology Optimization of Plasmonics with the FDTD method in Frequency Domain
Zhou Zeng and Xianfan Xu* School of Mechanical Engineering and Birck Nanotechnology Center Purdue University, West Lafayette, IN 47906 * Corresponding author, [email protected]
Abstract
We present a new method for calculating sensitivity in plasmonics with the FDTD method by direct differentiation of the time-marching system in frequency domain. This new method supports general frequency-domain objective functions, does not relay on implementation details of the FDTD method, works with general isotropic materials and can be incorporated into density-based topology optimization. We demonstrate the accuracy of this method for plasmonics. We also present a framework to carry out density-based topology optimization using the new sensitivity formula. We use non-linear material interpolation to counter the non-physical field amplification, adopt filtering-and-projection regularization to ensure manufacturability of the design. As examples of this approach, successful reconstruction of electric fields of a plasmonic bowtie aperture and optimization of absorption within a lossy medium produced by a plasmonic antenna with the framework are presented.
Keywords : topology optimization; plasmonics; discrete adjoint sensitivity analysis; density-based INTRODUCTION
Topology optimization of plasmonics has only been reported in recent years and is almost exclusively used with frequency-domain solvers such as the Finite Element Method (FEM). A nano-photonic grating coupler was reported to improve coupling efficiency from 50% to 68% and was designed using a density-based topology optimization with the FEM [1]. The same method was also applied to design a Ag antenna-strip yielding a more than 1200 fold field enhancement [2], structures with resonant peak at the specified incident wavelength [3], as well as a wide variety of two-dimensional plasmonic structures [4]. In addition, the level-set based approach with the FEM has shown promising results in design of metal dipole antennas with optimized surface current [5]. Topology optimization with the Finite-Difference Time-Domain (FDTD) [6] method in frequency domain based on topological derivatives was used to design a near-field transducer for heat-assisted magnetic recording [7]. The adjoint variable method (AVM) with FDTD in time domain [8] was demonstrated to accurately calculate sensitivity for two-dimensional plasmonic structures, but was not incorporated into a topology optimization scheme. The FDTD method is not a popular choice in topology optimization of plasmonics due to the problem rooted in plasmonics where numerical error in electromagnetic simulations of plasmonic structures is highest at the metal-dielectric interfaces. These locations feature strong field enhancement and contribute strongly to design sensitivities [9]. By using FEM calculations with conformal boundary representation one can obtain highly accurate design sensitivities. However, this is not applicable to the FDTD method. The problem forces many to use the inherently consistent discrete adjoint sensitivity [10] which gives an accurate evaluation of sensitivity in the presence of local field concentrations. It is compatible with solvers with a matrix representation to enable easy implementation. Here the sensitivity is defined as derivatives of the objective function to design variables for a discretized system of the original PDE. On the other hand, the variational adjoint sensitivity is derived from discretization, using numerical integration/interpolation, of derivatives for the original PDE. Figure 1 illustrates the differences between the continuous adjoint sensitivity analysis and the discrete adjoint sensitivity analysis. Fig. 1 Discrete adjoint sensitivity (left) compared with variational adjoint sensitivity (right). Differentiation of the discretized problem also involves a numerical model. Due to popularity of using the FDTD method in simulating plasmonics, we propose a new method of calculating frequency-domain sensitivity with the FDTD method as a frequency-domain solver and incorporate this method into the density-based topology optimization framework [11].To use the FDTD method as a frequency-domain solver, time-domain results of simulations are transformed into frequency domain by Fourier transform. The method accepts any objective function that is only dependent of frequency-domain variables, i.e., solutions of Maxwell’s equations at a designated frequency. As a result, the method solves similar problems as the wide-spread FEM-based approach which also gives solutions in frequency domain. The paper is structured as follows: We first derive the formula of frequency-domain sensitivity with the FDTD method and validate its accuracy with a numerical example. Subsequently, we introduce a density-based topology optimization method with the filtering-and-projection regularization as well as a non-linear material interpolation scheme to jointly improve robustness. The method is illustrated through reconstruction of electric fields of a plasmonic bowtie aperture and optimization of optical absorption in a lossy medium. DISCRETE ADJOINT SENSITIVITY ANALYSIS
Derivations
We start with the basic formulation of three dimensional FDTD method with lossy isotropic materials. For convenience, we index each field component by the position of the corresponding Yee cell, which is illustrated below. Fig. 2 Component arrangement of the Yee cell. Each E component is surrounded by four H components and each H component is surrounded by four E components. With the leapfrog time-stepping scheme, we can write down evaluations of all field components (cid:1831) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:3435)(cid:4666)(cid:1861) (cid:3397) 0.5(cid:4667)Δ(cid:1876), (cid:1862)Δ(cid:1877), (cid:1863)Δ(cid:1878), (cid:1866)Δ(cid:1872)(cid:3439) (cid:1831) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:4666)(cid:1861)Δ(cid:1876), (cid:4666)(cid:1862) (cid:3397) 0.5(cid:4667)Δ(cid:1877), (cid:1863)Δ(cid:1878), (cid:1866)Δ(cid:1872)(cid:4667) (cid:1831) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:4666)(cid:1861)Δ(cid:1876), (cid:1862)Δ(cid:1877), (cid:4666)(cid:1863) (cid:3397) 0.5(cid:4667)Δ(cid:1878), (cid:1866)Δ(cid:1872)(cid:4667) (cid:1834) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:4666)(cid:1861)Δ(cid:1876), (cid:4666)(cid:1862) (cid:3397) 0.5(cid:4667)Δ(cid:1877), (cid:4666)(cid:1863) (cid:3397) 0.5(cid:4667)Δ(cid:1878), (cid:4666)(cid:1866) (cid:3397) 0.5(cid:4667)Δ(cid:1872)(cid:4667) (cid:1834) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:3435)(cid:4666)(cid:1861) (cid:3397) 0.5(cid:4667)Δ(cid:1876), (cid:1862)Δ(cid:1877), (cid:4666)(cid:1863) (cid:3397) 0.5(cid:4667)Δ(cid:1878), (cid:4666)(cid:1866) (cid:3397) 0.5(cid:4667)Δ(cid:1872)(cid:3439) (cid:1834) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:3404) (cid:1831)(cid:3435)(cid:4666)(cid:1861) (cid:3397) 0.5(cid:4667)Δ(cid:1876), (cid:4666)(cid:1862) (cid:3397) 0.5(cid:4667)Δ(cid:1877), (cid:1863)Δ(cid:1878), (cid:4666)(cid:1866) (cid:3397) 0.5(cid:4667)Δ(cid:1872)(cid:3439) (1a) (1b) (1c) (1d) (1e) (1f) where Δ(cid:1876), Δ(cid:1877), Δ(cid:1878) are space steps and
Δ(cid:1872) is the time step. The central-difference expressions for the space and time derivatives convert the two differential curl equations (cid:2034)(cid:1776)(cid:2034)(cid:1872) (cid:3404) (cid:3398)(cid:1487) (cid:3400) (cid:1779) (cid:3398) (cid:1787) (cid:2034)(cid:1778)(cid:2034)(cid:1872) (cid:3404) (cid:1487) (cid:3400) (cid:1782) (cid:3398) (cid:1784) (2) (3) into two finite difference equations (cid:1854) (cid:3041) (cid:3398) (cid:1854) (cid:3041)(cid:2879)(cid:2869)
Δ(cid:1872) (cid:3404) (cid:3398)(cid:1777) (cid:3032) (cid:1857) (cid:3041) (cid:3398) (cid:1865) (cid:3041) (cid:1856) (cid:3041)(cid:2878)(cid:2869) (cid:3398) (cid:1856) (cid:3041)
Δ(cid:1872) (cid:3404) (cid:1777) (cid:3035) (cid:1860) (cid:3041) (cid:3398) (cid:1862) (cid:3041) (4) (5) where (cid:1856) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1830) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1830) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1830) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1860) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1834) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1834) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1834) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1862) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1836) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1836) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1836) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) (6) and (cid:1854) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1828) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1828) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1828) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1857) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1831) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1831) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1831) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1865) (cid:3041) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:1839) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1839) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:1839) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) (7) are column vectors representing the relevant fields; (cid:1777) (cid:3035) , whose nonzero elements are (cid:3399) 1 Δ(cid:1876)⁄ , (cid:3399) 1 Δ(cid:1877)⁄ and (cid:3399) 1 Δ(cid:1878)⁄ , is a matrix for the curl operator on the H components while (cid:1777) (cid:3032) is a matrix for the curl operator on the E components. The z-transform [12] converts a discrete-time signal, which is a sequence of real or complex numbers, into a complex frequency-domain representation given below: (cid:1876)(cid:3556)(cid:4666)(cid:1878)(cid:4667) (cid:3404) (cid:3533) (cid:1876)(cid:4670)(cid:1866)(cid:4671)(cid:1878) (cid:2879)(cid:3041)(cid:2998)(cid:3041)(cid:2880)(cid:2868) (8) where tilde denotes z-transform of the corresponding variable. Applying z-transform to the two finite difference equations (4), (5) yields (cid:1878)(cid:1856)(cid:4634)(cid:4666)(cid:1878)(cid:4667) (cid:3398) (cid:1856)(cid:4634)(cid:4666)(cid:1878)(cid:4667)Δ(cid:1872) (cid:3404) (cid:1777) (cid:3035) (cid:1860)(cid:3560)(cid:4666)z(cid:4667) (cid:3398) (cid:2836)̃(cid:4666)z(cid:4667) (cid:1854)(cid:3560)(cid:4666)z(cid:4667) (cid:3398) (cid:1878) (cid:2879)(cid:2869) (cid:1854)(cid:3560)(cid:4666)(cid:1878)(cid:4667)Δ(cid:1872) (cid:3404) (cid:3398)(cid:1777) (cid:3032) (cid:1857)̃(cid:4666)z(cid:4667) (cid:3398) (cid:1865)(cid:3557) (cid:4666)(cid:1878)(cid:4667) (9) (10) Popular implementations of linear dispersions [6], including the recursive-convolution method and the auxiliary differential equation method, can also be cast into frequency-domain expressions using z-transform (cid:1856)(cid:4634)(cid:4666)z(cid:4667) (cid:3404) (cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667) (cid:1854)(cid:3560)(cid:4666)(cid:1878)(cid:4667) (cid:3404) (cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667) (11) (12) Note that (cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) and (cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667) are diagonal matrix for orthotropic materials given below (cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) (cid:3404) diag (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:2013) (cid:3051)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2013) (cid:3052)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2013) (cid:3053)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667) (cid:3404) diag (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:2020) (cid:3051)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2020) (cid:3051)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2020) (cid:3051)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) (13) Letting (cid:1878) (cid:3404) (cid:1857) (cid:3036)(cid:3104)(cid:2940)(cid:3047) yields solutions at frequency (cid:2033) . (cid:2033)Δ(cid:1872) is often very small to resolve aliasing in discrete-time fourier transform (DTFT) such that the following approximations can be made (cid:1878) ≅ 1 (cid:3397) (cid:1861)(cid:2033)Δ(cid:1872) (14a) (cid:1878) (cid:2879)(cid:2869) ≅ 1 (cid:3398) (cid:1861)(cid:2033)Δ(cid:1872) (14b) Combining the two difference equations (9), (10), the z-transform of constitutive equations (11), (12) and the two approximations (14a), (14b), we have (cid:3436)(cid:3398)(cid:1861)(cid:2033)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) (cid:1777) (cid:3035) (cid:1777) (cid:3032) (cid:1861)(cid:2033)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3436)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3404) (cid:3436) (cid:2836)̃(cid:4666)(cid:1878)(cid:4667)(cid:3398)(cid:1865)(cid:3557) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (15) which is equivalent to the Finite-Difference Frequency-Domain (FDFD) formulation [13]. Through some lengthy manipulations of index in Equation (1), it can be shown that the system matrix in (15) is symmetric in problems where no Perfectly Matched Layers (PMLs) are present and only lossless boundary conditions, including Perfect Electric Conductor (PEC), Perfect Magnetic Conductor (PMC) and periodic boundary conditions, are used. According to the adjoint sensitivity analysis [14], given an objective function (cid:1832)(cid:3435)(cid:1779)(cid:4666)(cid:2033)(cid:4667), (cid:1782)(cid:4666)(cid:2033)(cid:4667)(cid:3439) ∈(cid:1337) which is numerically evaluated as (cid:1858)(cid:3435)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667), (cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667) (cid:3439) ∈ (cid:1337) , variation of the numerical evaluation (cid:1858) produced by variation of material properties (cid:2012)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667), (cid:2012)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667) is given by (cid:2012)(cid:1858) (cid:3404) Re (cid:4684)(cid:3398) (cid:3436)(cid:1857)̃ (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560) (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3021) (cid:2012) (cid:3436)(cid:3398)(cid:1861)(cid:2033)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) (cid:1777) (cid:3035) (cid:1777) (cid:3032) (cid:1861)(cid:2033)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3436)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667)(cid:3440)(cid:4685) (cid:3404) Re (cid:4684)(cid:3398) (cid:3436)(cid:1857)̃ (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560) (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3021) (cid:3436)(cid:3398)(cid:1861)(cid:2033)(cid:2012)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) (cid:1861)(cid:2033)(cid:2012)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3436)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667)(cid:3440)(cid:4685) (cid:3404) Re (cid:4672)(cid:1861)(cid:2033)(cid:1857)̃ (cid:3028)(cid:3021) (cid:4666)(cid:1878)(cid:4667)(cid:2012)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667)(cid:1857)̃(cid:4666)(cid:1878)(cid:4667) (cid:3398) (cid:1861)(cid:2033)(cid:1860)(cid:3560) (cid:3028)(cid:3021) (cid:4666)(cid:1878)(cid:4667)(cid:2012)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667)(cid:4673) (cid:3404) Re (cid:4684) (cid:3533) (cid:1861)(cid:2033)(cid:1831)(cid:3560) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:1831)(cid:3560) (cid:3028)(cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2012)(cid:2013) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:3398) (cid:1861)(cid:2033)(cid:1834)(cid:3561) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:1834)(cid:3561) (cid:3028)(cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:2012)(cid:2020) (cid:3050)(cid:3036),(cid:3037),(cid:3038)(cid:3050),(cid:3036),(cid:3037),(cid:3038) (cid:4685) (16) where (cid:3435)(cid:1857)̃ (cid:3028) (cid:4666)(cid:1878)(cid:4667), (cid:1860)(cid:3560) (cid:3028) (cid:4666)(cid:1878)(cid:4667) (cid:3439) are the adjoint solution satisfying (cid:3436)(cid:3398)(cid:1861)(cid:2033)(cid:1830) (cid:3032) (cid:4666)(cid:1878)(cid:4667) (cid:1777) (cid:3035) (cid:1777) (cid:3032) (cid:1861)(cid:2033)(cid:1830) (cid:3035) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3436)(cid:1857)̃ (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:1860)(cid:3560) (cid:3028) (cid:4666)(cid:1878)(cid:4667)(cid:3440) (cid:3404) (cid:3436) (cid:1487) (cid:3032) (cid:1858) (cid:1487) (cid:3035) (cid:1858)(cid:3440)(cid:3364)(cid:3364)(cid:3364)(cid:3364)(cid:3364)(cid:3364)(cid:3364)(cid:3364)(cid:3364) (17) It should be noted that transpose is removed in the adjoint system because of symmetry of the matrix. The adjoint currents (cid:1487) (cid:3032) (cid:1858), (cid:1487) (cid:3035) (cid:1858) are first derivatives of the objective function with respect to field variables (cid:1487) (cid:3032) (cid:1858) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1831)(cid:3560) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1831)(cid:3560) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1831)(cid:3560) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) , (cid:1487) (cid:3035) (cid:1858) (cid:3404) (cid:1737)(cid:1736)(cid:1736)(cid:1735) ⋮(cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1834)(cid:3561) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1834)(cid:3561) (cid:3052)(cid:3036),(cid:3037),(cid:3038),(cid:3041) (cid:2034)(cid:1858) (cid:3403) (cid:2034)(cid:1834)(cid:3561) (cid:3053)(cid:3036),(cid:3037),(cid:3038),(cid:3041) ⋮ (cid:1740)(cid:1739)(cid:1739)(cid:1738) (18) 0 The symmetry of the system matrix suggests the same setup can be used for both forward simulation (15) and adjoint simulation (17) and the only difference between them is the source currents. Properly configured PMLs absorb most outgoing waves and have the same effect as an unbounded region, leading to the same sensitivity formula if no material perturbation can occur in PMLs. It should be noted that: 1) the summation in (16) is the FDTD equivalence of integration; 2) the Re operator comes from the fact that the objective function is real; 3) the analysis is done in the frequency domain such that the system is put together as a single matrix equation. In short, we convert a FDTD problem to a FDFD problem which we solve gradients for. Implementation
To solve the adjoint simulation, the objective function needs to be written explicitly with arguments from vectors (cid:1857)̃(cid:4666)(cid:1878)(cid:4667) and (cid:1860)(cid:3560)(cid:4666)(cid:1878)(cid:4667) with respect to which the derivatives can be calculated. In practice, the objective function is always calculated based on fields interpolated from the original electromagnetic fields at the Yee cell. To simplify notation, we assume that the objective function only depends on electric field (cid:1857)̃(cid:4666)(cid:1878)(cid:4667) , which is the case for most problems. Denote the interpolated electric fields by (cid:1857)̃(cid:4632)(cid:4666)z(cid:4667) (cid:3404) (cid:1737)(cid:1736)(cid:1735)(cid:1831)(cid:3560)(cid:3552) (cid:2869) (cid:4666)(cid:1878)(cid:4667)(cid:1831)(cid:3560)(cid:3552) (cid:2870) (cid:4666)(cid:1878)(cid:4667)⋮(cid:1831)(cid:3560)(cid:3552) (cid:3041) (cid:4666)(cid:1878)(cid:4667)(cid:1740)(cid:1739)(cid:1738)
19 Assuming linear interpolation, each interpolated electric field is given by 1 (cid:1831)(cid:3560)(cid:3552) (cid:3036) (cid:3404) (cid:3533) (cid:1875) (cid:3036)(cid:3037) (cid:1831)(cid:3560) (cid:3037)(cid:3037)∈(cid:3015) (cid:3284)
20 where (cid:1840) (cid:3036) is the neighboring element for each interpolated electric field and (cid:1875) (cid:3036)(cid:3037) is the linear interpolation weight for the ith point from its j th neighbor [15]. We discard the original index system (cid:1831) (cid:3051)(cid:3036),(cid:3037),(cid:3038),(cid:3041) to avoid cluttered scripts. By chain rule, we have derivatives of objective function with respect to original electric fields at the Yee cell given by (cid:2034)(cid:1858)(cid:2034)(cid:1831)(cid:3560) (cid:3036) (cid:4666)(cid:1878)(cid:4667) (cid:3404) (cid:3533) (cid:2034)(cid:1858)(cid:2034)(cid:1831)(cid:3560)(cid:3552) (cid:3037) (cid:4666)(cid:1878)(cid:4667) (cid:3036)∈(cid:3015) (cid:3285) (cid:1875) (cid:3037)(cid:3036)
21 where (cid:3105)(cid:3033)(cid:3105)(cid:3006)(cid:3560)(cid:3552) (cid:3285) (cid:4666)(cid:3053)(cid:4667) on the right-hand side can be given an analytic form depending on the user-defined formula of objective function. Mathematically, the above equation exhibits a well-known conecpt originating in multi-grid methods – the transpose of interpolation [16]. The adjoint currents are obtained by differentiation of the objective function to the interpolated fields followed by transpose of interpolation given by chain rule (21). Then the adjoint problem (17) can solved. With solutions from the forward and backward problems, sensitivity is evaluated according to (16). Similar to what is implemented in the MEEP FDTD solver [17], we stop simulation once lingering fields in the volume where sensitivity is evaluated are smaller than 0.5% of their maximum. 2
Example: Absorption of a Block with Varying Permittivity
We demonstrate the accuracy of the sensitivity formula for a function measuring the absorption of a plasmonic scatterer. The measure of absorption is the integral of Poynting flux over a specified region (cid:1832) (cid:3404) (cid:3505)Re(cid:4666)(cid:1779) (cid:3400) (cid:1782)(cid:3365) (cid:4667) ⋅ (cid:1866)(cid:3548)(cid:1856)(cid:1845) (cid:3020) (22) We consider a simple geometry, a block with a size of
30 (cid:3400) 30 (cid:3400) 30 nm . The block has spatially uniform permittivity controlled by a density parameter using linear interpolation (cid:2013)(cid:4666)(cid:2025)(cid:4667) (cid:3404) (cid:2013) (cid:2869) (cid:2025) (cid:3397) (cid:2013) (cid:2870) (cid:4666)1 (cid:3398) (cid:2025)(cid:4667) (23) where (cid:2013) (cid:2869) is permittivity of gold and (cid:2013) (cid:2870) is permittivity of vacuum. According to (16), the adjoint sensitivity is given by ∂(cid:1858)(cid:2034)(cid:2025) (cid:4666)(cid:2025)(cid:4667) (cid:3404) Re (cid:4684) (cid:3533) (cid:1861)(cid:2033)(cid:1831)(cid:3560) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:1831)(cid:3560) (cid:3028)(cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:4666)(cid:2013) (cid:2869) (cid:3398) (cid:2013) (cid:2870) (cid:4667) (cid:3050),(cid:3036),(cid:3037),(cid:3038) (cid:4685) (24) The absorption is calculated by integrating the Poynting flux over a box enclosing the object. A plane wave of wavelength 800 nm is incident upon the block. The layout for the example is shown in Fig. 3. 3 Fig. 3 Layout for the example. A box of size 3
30 (cid:3400) 30 (cid:3400) 30 nm is seated in an open region with its axis aligned with coordinates. An incident plane wave travels in the (cid:1878) (cid:2878) direction. PMLs of 8-cell thickness is used to absorb any outgoing waves. Absorption is calculated on the box surrounding the block. We take the LD4 model [18] as dispersion model of gold. We linearly scale the strength of each susceptibility to implement the linearly interpolated permittivity (23). At 800 nm, the permittivity of gold is (cid:2013) (cid:2869) (cid:3404) (cid:2013) (cid:2868) (cid:4666)(cid:3398)22.3 (cid:3398) 2.03(cid:1862)(cid:4667) calculated from the LD4 model. We use a 2 nm Yee cell length to run forward simulations to calculate the absorption versus the density parameter (cid:2025) and run adjoint simulations to get the sensitivity with respect to (cid:2025) . Result of the absorption is shown in Fig. 4. 4 -1 10
01 10 A b s o r p t i on Fig. 4 The absorption of the block versus the density parameter (cid:2025) . It should be noted that the absolute value of absorption does not reflect actual physical quantity. We then compare the adjoint sensitivity (24) with sensitivity calculated from central difference scheme (cid:2034)(cid:1858)(cid:2034)(cid:2025) (cid:4666)(cid:2025)(cid:4667) ≅ (cid:1858)(cid:4666)(cid:2025) (cid:3397) Δ(cid:2025)(cid:4667) (cid:3398) (cid:1858)(cid:4666)(cid:2025) (cid:3398) Δ(cid:2025)(cid:4667)2Δ(cid:2025) (25) to verify the accuracy of the discrete adjoint sensitivity formula. Fig. 5 compares the two. The obvious agreement between the two confirms the validity of the discrete adjoint sensitivity formula. 5 -3 10 -2 10 -1 10
01 10 (i) -2 10 -1 10
01 10 (ii) -1.6 10 -8 10
08 10 (iii) Fig. 5. (i)-(iii): Comparison of the central difference sensitivity and the discrete adjoint sensitivity. (i): (cid:2025) ∈ (cid:4670)0, 0.1(cid:4671) . (ii): (cid:2025) ∈ (cid:4670)0.1, 0.5(cid:4671) . (iii): (cid:2025) ∈ (cid:4670)0.5, 1(cid:4671) . 6 Comparison with Previous Work
To our knowledge, there are only two original methods of sensitivity analysis of plasmonic structures with the FDTD method to date. The adjoint variable method (AVM) proposed in [19] considers the time marching scheme in the FDTD method as a dynamic system and uses time traversal to simulate the adjoint system. It was shown that this method works reasonably well in both dielectrics [20] and plasmonics [21]. However, even though the adjoint equation can be solved without modifying the existing FDTD solver, this approach only works with permittivity while ours works with both permittivity and permeability. The approach depends on implementation of dispersion while ours transforms the dispersion to a unified z-domain expression. The biggest drawback of the approach is that the objective function is restricted to a time integral function. It is not applicable to functions that has a non-linear dependence on frequency-domain variables. A method based on topological derivatives [22] was applied in design of a near-field transducer for heat-assisted magnetic recording [7]. It evolves a bit-map representation of structure by metalizing bits at the boundary and in the void. Each bit is either metalized or nonmetallized based on changes of the objective function due to perforation of the material domain by an infinitesimal hole. However, it is still a continuum method which is prone to errors in plasmonic structures. 7 DENSITY-BASED TOPOLOGY OPTIMIZATION
Density-based topology optimization uses a material interpolation scheme to relate material property to density parameters (cid:2025) ’s. For example, we can assign permittivity (cid:2013) using linear interpolation (cid:2013)(cid:4666)(cid:1876), (cid:2025)(cid:4667) (cid:3404) (cid:2013) (cid:2869) (cid:2025)(cid:4666)(cid:1876)(cid:4667) (cid:3397) (cid:2013) (cid:2870) (cid:3435)1 (cid:3398) (cid:2025)(cid:4666)(cid:1876)(cid:4667)(cid:3439), (cid:2025) ∈ (cid:4670)0,1(cid:4671) (26) where (cid:2013) (cid:2869) is the permittivity for the material domain and (cid:2013) (cid:2870) is the permittivity for the void domain (usually (cid:2013) (cid:2870) (cid:3404) (cid:2013) (cid:2868) ). Under a specific numerical model of the Maxwell’s equations, each node (line segment, pixel or voxel) is assigned a density parameter and the problem of finding the optimal design amounts to determining the value of (cid:2025) for each node. To utilize gradient-based optimization techniques, density parameters can take intermediate values, resulting in gray transition regions. Such gray transition regions do not represent any realistic material and might lead to unphysical device design. Hence, regularizations must be applied to restrict the design space and to ensure manufacturability of the design.
Non-linear Material Interpolation
It is suggested [2] that in topology optimization of plasmonics the non-linear material interpolation is favored over the linear interpolation in (26) which can cause non-physical field amplification. In the non-linear method, the refractive index is linearly interpolated (cid:2013)(cid:4666)(cid:2025)(cid:4667) (cid:3404) (cid:3435)(cid:3493)(cid:2013) (cid:2870) (cid:3397) (cid:4666)1 (cid:3398) (cid:2025)(cid:4667)(cid:3493)(cid:2013) (cid:2869) (cid:3439) (cid:2870) , (cid:2025) ∈ (cid:4670)0,1(cid:4671) (27) 8 To implement the non-linear material interpolation in the FDTD method, we need to match the permittivity at specified (cid:2025) . Take gold at wavelength 800 nm for example. Since a complex number is in fact a two-dimensional vector, we pick two susceptibility out of the LD4 model [18] and modify their intensity to yield an accurate value for permittivity. Assuming the two susceptibility are given by (cid:2031) (cid:2869)(cid:4593) (cid:3398) (cid:1861)(cid:2031) (cid:2870)(cid:4593)(cid:4593) , (cid:2031) (cid:2870)(cid:4593) (cid:3398) (cid:1861)(cid:2031) (cid:2870)(cid:4593)(cid:4593) ∈ (cid:1319) and their intensities are (cid:2026) (cid:2869) , (cid:2026) (cid:2870) ∈ (cid:1337) respectively, we can calculate the intensities by solving the linear equations below (cid:2013) (cid:4593) (cid:4666)(cid:2025)(cid:4667) (cid:3404) (cid:2013) (cid:2998) (cid:3397) (cid:2026) (cid:2869) (cid:4666)(cid:2025)(cid:4667)(cid:2031) (cid:2869)(cid:4593) (cid:3397) (cid:2026) (cid:2870) (cid:4666)(cid:2025)(cid:4667)(cid:2031) (cid:2870)(cid:4593) (cid:2013) (cid:4593)(cid:4593) (cid:4666)(cid:2025)(cid:4667) (cid:3404) (cid:2026) (cid:2869) (cid:4666)(cid:2025)(cid:4667)(cid:2031) (cid:2870)(cid:4593)(cid:4593) (cid:3397) (cid:2026) (cid:2870) (cid:4666)(cid:2025)(cid:4667)(cid:2031) (cid:2870)(cid:4593)(cid:4593) (28a) (28b) where (cid:2013) (cid:2998) is the instantaneous dielectric constant (must be positive) and (cid:2013)(cid:4666)(cid:2025)(cid:4667) (cid:3404) (cid:2013) (cid:4593) (cid:4666)(cid:2025)(cid:4667) (cid:3398) (cid:1861)(cid:2013) (cid:4593)(cid:4593) (cid:4666)(cid:2025)(cid:4667) . We fix (cid:2013) (cid:2998) to be 1 and calculate (cid:2026) (cid:2869) , (cid:2026) (cid:2870) versus (cid:2025) using the non-linear interpolation scheme for gold at 800 nm (Fig. 6). Notice that this method only works for single frequency and does not generally guarantee to work for every material. It works for gold at 800 nm, which is sufficient for our purpose. Each point in the Yee lattice is assigned a density parameter to facilitate density-based topology optimization. 9 Fig. 6. Susceptibility intensities for non-linear material interpolation scheme of gold at 800 nm. Non-linear dependence of the intensity can be observed.
Density Derivatives
A gradient-based optimization algorithm uses derivatives of design variables to evolve the solution at each iteration. The derivatives of objective function with respect to density parameters are density derivatives. In our approach, a rectangular design domain is constructed with a regular grid (cartesian lattice). Each point in the grid is assigned a density parameter (cid:2025) and linear interpolation is used to map grid density parameters into Yee lattice density parameters. In the 2D example illustrated in Fig. 7, (cid:2025) (cid:3036),(cid:3037) denotes density parameters on the regular grid while (cid:1868) (cid:3051)(cid:3036),(cid:3037) denotes density parameters on the Yee lattice. -20246810 0 0.2 0.4 0.6 0.8 1
0 Fig. 7 Linearly interpolated density parameters in a 2D Yee lattice. Each Yee lattice density parameter (in green) is linearly interpolated from the grid density parameter (in blue). The Yee cell size does not have to be equal to the grid size. The interpolation effectively decouples the geometry from the Yee lattice such that the design region can be arbitrarily sized and positioned. However, the grid size for the design domain cannot be arbitrarily small. It cannot be 2 times smaller than the grid size of the Yee lattice, otherwise it leaves some grid points unused during the interpolation. The sensitivities with respect to 3D Yee lattice density parameters is given by (cid:2034)(cid:1858)(cid:2034)(cid:1868) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:3404) Re (cid:4684) (cid:3533) (cid:1861)(cid:2033)(cid:1831)(cid:3560) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:1831)(cid:3560) (cid:3028)(cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667) (cid:2034)(cid:2013) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:2034)(cid:1868) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:3398) (cid:1861)(cid:2033)(cid:1834)(cid:3561) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667)(cid:1834)(cid:3561) (cid:3028)(cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:4666)(cid:1878)(cid:4667) (cid:2034)(cid:2020) (cid:3050)(cid:3036),(cid:3037),(cid:3038) (cid:2034)(cid:1868) (cid:3050)(cid:3036),(cid:3037),(cid:3038)(cid:3050),(cid:3036),(cid:3037),(cid:3038) (cid:4685) (29) 1 according to the sensitivity formula (16) where (cid:3105)(cid:3084) (cid:3298)(cid:3284),(cid:3285),(cid:3286) (cid:3105)(cid:3043) (cid:3298)(cid:3284),(cid:3285),(cid:3286) , (cid:3105)(cid:3091) (cid:3298)(cid:3284),(cid:3285),(cid:3286) (cid:3105)(cid:3043) (cid:3298)(cid:3284),(cid:3285),(cid:3286) are evaluated based on material interpolation scheme at choice. Following the chain rule in calculus, we propagate sensitivities back to the regular grid of the design domain to be used in optimization algorithm. Filtering and Projection Scheme
The three-field density representation [23] is adopted as the regularization technique in our method, which uses a density filter followed by a projection. Viewing discretized density parameters (cid:2025) as a 3D image, we can reassign each density value (cid:2025) (cid:3036) as a weighted sum of its neighbors (cid:2025)(cid:3548) (cid:3036) (cid:3404) ∑ (cid:1875) (cid:3036)(cid:3037) (cid:2025) (cid:3037)(cid:3037)∈(cid:3015) (cid:3284) ∑ (cid:1875) (cid:3036)(cid:3037)(cid:3037)∈(cid:3015) (cid:3284) (30) We use a Gaussian distribution function for the filter function (cid:1875) (cid:3036)(cid:3037) (cid:1875) (cid:3036)(cid:3037) (cid:3404) (cid:1857) (cid:2879)(cid:2868).(cid:2873)(cid:4678)(cid:3630)(cid:3051) (cid:3284) (cid:2879)(cid:3051) (cid:3285) (cid:3630)(cid:3019) (cid:2870)⁄ (cid:4679) (cid:3118) (31) as suggested in [24]. (cid:1844) denotes the filter radius which is used to truncate the filter. Density filters can effectively remove small features during the optimization process but end up creating more gray transition regions. Therefore, a projection operator is used to project the filtered density to 0/1. Suggested in [25], the following projection operator is used in our method 2 (cid:2025)(cid:3548)̅ (cid:3036) (cid:3404) tanh(cid:4666)(cid:2010)(cid:2015)(cid:4667) (cid:3397) tanh(cid:3435)(cid:2010)(cid:4666)(cid:2025)(cid:3548) (cid:3036) (cid:3398) (cid:2015)(cid:4667)(cid:3439)tanh(cid:4666)(cid:2010)(cid:2015)(cid:4667) (cid:3397) tanh(cid:3435)(cid:2010)(cid:4666)1 (cid:3398) (cid:2015)(cid:4667)(cid:3439) (32) where (cid:2010) is a parameter used to control the strength of the projection and (cid:2015) ∈ (cid:4670)0,1(cid:4671) defines the projection level. Eventually, projected density parameters (cid:2025)(cid:3548)̅ (cid:3036) are linearly interpolated into Yee lattice density parameters (cid:1868) (cid:3050)(cid:3036)(cid:3037)(cid:3038) . To prevent the optimization from getting stuck in local minimum earlier in the process due to the density being projected to 0/1 immediately, projection strength (cid:2010) is increased after the objective function has not changed much for (cid:1866) (cid:3046)(cid:3030) iterations, mimicking the (cid:2010) -continuation scheme suggested in [25]. At the same time, we also allow the filter radius to change. Our continuation scheme can be written in algorithmic form shown below. The method of moving asymptotes (MMA) [26] is applied to solve the optimization in each round.
Initialize filter radius = (cid:2200)
Initialize projection level = (cid:2236)
Run the following for k times: (cid:2191) (cid:3404) (cid:2778)
Initialize design parameters (cid:2251) (cid:2778) while True: Evaluate objective function (cid:2188) (cid:2191) if |(cid:2188) (cid:2191) (cid:3398) (cid:2188) (cid:2191)(cid:2879)(cid:2778) | (cid:3407) (cid:2235)(cid:2188) (cid:2191) for consecutively (cid:2196) (cid:2201)(cid:2185) times then: (cid:2200) (cid:3404) (cid:2200) ⁄ (cid:2185) (cid:2200) (cid:2236) (cid:3404) (cid:2236) (cid:3400) (cid:2185) (cid:2236) break out of the loop run forward and adjoint simulation evaluate sensitivity (cid:2260)(cid:2188) (cid:3403) (cid:2260)(cid:2251)(cid:4666)(cid:2251) (cid:2191) (cid:4667) evolve the solution from (cid:2251) (cid:2191) to (cid:2251) (cid:2191)(cid:2878)(cid:2778) using MMA (cid:2191) (cid:3404) (cid:2191) (cid:3397) (cid:2778) return the solution We programmed a FDTD solver in C++ for handling solutions of forward and adjoint problems. The solver allows specifying material property of each Yee cell point so that we are able implement the non-linear material interpolation. A Python interface was designed to talk to the solver, calculate sensitivities and evolve solutions. This is very similar to COMSOL’s built-in graphical interface of adjoint solver which the user writes a script to incorporate into a topology optimization scheme.
Comparison between FDFD and FDTD in Frequency Domain
Surprisingly, there is no known use of the FDFD method in topology optimization. From an analytical point of view, there is no difference between using the FDFD method and using the FDTD method as a frequency domain solver in topology optimization of plasmonics since the former is the frequency-domain form of the latter. In terms of execution, the straightforward material model of FDFD enables easier implementation of non-linear material interpolation compared to FDTD. However, the ill-conditioned system matrix due to an inappropriate choice of the PML and the high multiplicity of near-zero eigenvalues of system matrix for plasmonic and 4 nanophotnic systems that stagnates iterative methods make it much more difficult to implement an efficient 3D FDFD solver [13]. CASE STUDY 1: RECONSTRUCTION OF ELECTRIC FIELDS OF A BOWTIE APERTURE
Bowtie apertures are known to produce highly localized fields and have potential applications in optical lithography [27], high density data storage [28], etc. The enhanced electric fields are confined within only a tiny region of the nanometer length scale near the surface of the nanostructures and decay significantly thereafter [4]. In this example, reconstruction of electric fields of a plasmonic bowtie aperture is carried out using the proposed sensitivity analysis in Section 2 and density-based topology optimization in Section 3. Fig. 8 Dimensions of the bowtie aperture for generating the objective electric field. 5 A bowtie aperture carved out of a gold layer of thickness 60 nm is placed in an open region where an incident wave of wavelength 800 nm is traveling in direction perpendicular to the layer. Fig. 8 shows the dimensions of the bowtie aperture. The objective electric field (cid:1779) (cid:2868) generated by the aperture is calculated using the FDTD method with a Yee cell length of 4 nm. Although in practice a smaller mesh size like 2 nm is needed to accurately capture the near field [29], we did not use a finer size because we are not interested in the actual physics of the system. Layout for the optimization is shown below Fig. 9 Layout for the optimization (case study 1). All rectangles are squares with their centers alighned. To measure the difference between the design electric field (cid:1779) and the objective electric field (cid:1779) (cid:2868) , the objective function is defined as 6 (cid:1832)(cid:4666)(cid:1779)(cid:4667) (cid:3404) (cid:3505)(cid:4666)‖(cid:1779)‖ (cid:3398) ‖(cid:1779) (cid:2868) ‖(cid:4667) (cid:2870) (cid:1856)(cid:1845) (cid:3020) (33) where (cid:1845) is a rectangular surface at the exit plane of the bowtie aperture (Fig. 9). It is obvious that a smaller value of the objective function implies a better agreement with the objective field. The surface (cid:1845) is chosen at the exit plane, which increases the possibility of the algorithm to recover the original geometry because different structures can produce the same field patterns at a large distance. The design region is a
200 (cid:3400) 200 (cid:3400) 60 nm block centered in the gold layer (Fig. 9). The region is big enough so that the original bowtie geometry can be fitted in. The design region is discretized into a (cid:1840) (cid:3051) (cid:3400) (cid:1840) (cid:3052) (cid:3400) (cid:1840) (cid:3053) rectangular grid. Below show parameters used in the optimization. Note that (cid:1840) (cid:3053) discretized points are equivalent to (cid:1840) (cid:3053) (cid:3398) 1 discretized intervals. Table 1 Parameters used in optimization (case study 1). Parameter [unit] Value (cid:1844) [Yee cell length] 3 (cid:2010) (cid:1863) (cid:2009) (cid:1855) (cid:3045) (cid:1855) (cid:3081) (cid:1840) (cid:3051) , (cid:1840) (cid:3052) , (cid:1840) (cid:3053)
51, 51, 16 (cid:1845) (cid:3051) , (cid:1845) (cid:3052)
41, 41 (cid:1866) (cid:3046)(cid:3030)
3 7 The objective field is obtained by setting up the objective density distribution (cid:2025) (cid:2868) (Fig. 10) in the design region and running the simulation once. Viewing the bowtie geometry as a planar polygon, the objective density is to set to 1 if it is at a point inside the polygon. The optimization algorithm starts with (cid:2025) (cid:3404) 1 , representing pure gold layer, and runs for 5 rounds of optimization with different filter radius and projection strength, which constitutes a total of 105 iterations and 185 function evaluations. We enforce the density to be constant along z direction so that the optimization produces a planar structure. Fig. 10 shows distributions of density in the final iteration, where the projected density (cid:2025)(cid:3548)̅ represents the optimized structure. Although symmetry is not enforced during the optimization, the density distributions are symmetric with a smaller bowtie aperture and four holes in the corner. The smaller bowtie aperture generates the same magnitude of localized electric field in the center as the objective bowtie aperture, suggesting that much of the empty space is redundant in the objective bowtie aperture. Four holes emerge because the objective electric field has four smaller spikes in the four corners of (cid:1845) whose values are about 6 times that of the incident field. 8 Fig. 10 Distributions of the design density (cid:2025) , the filtered density (cid:2025)(cid:3548) , the projected density (cid:2025)(cid:3548)̅ in the final iteration and the objective density (cid:2025) (cid:2868) (case study 1). Fig. 11 compares the final electric field magnitude distribution with the objective electric field magnitude distribution on the surface (cid:1845) . The (cid:1838) (cid:2870) norm difference (cid:1516) (cid:4666)‖(cid:1779)‖(cid:2879)‖(cid:1779) (cid:3116) ‖(cid:4667) (cid:3118) (cid:3031)(cid:3020) (cid:3268) (cid:1516) ‖(cid:1779) (cid:3116) ‖ (cid:3118) (cid:3031)(cid:3020) (cid:3268) bewteen them is about 3%, suggesting a successful recovery of the electric field on the surface (cid:1845) . 9 Fig. 11 Distributions of the final electric field magnitude and the objective electric field magnitude on the surface (cid:1845) (case study 1). The magnitudes are normalized by the incident magnitude. The maximum magnitude at center is 42.6 for both distributions. A closer look at intermediate steps (Fig. 12) of the optimization reveals that only boundary perturbation is taking place in later iterations. The sharp decrease in the objective function occurring between iteration 49 and iteration 50 is induced by the change of projection strength (cid:2010) from 5.8 to 10.5. Higher value of (cid:2010) brings in higher contrast of air and gold to the boundary, which contributes to the generation of highly localized fields. 0 Fig. 12 Distributions of the projected density (cid:2025)(cid:3548)̅ (left) and the field magnitude (right) in some intermediate steps showing evolution of the density (case study 1). Each iteration is marked with diamond in the objective function graph at the bottom. The objective function is normalized so that it is evaluated as 1 in the 1 st iteration. 1 Fig. 13 Distributions of the projected density (cid:2025)(cid:3548)̅ (left) and the design field magnitude (right) in the first 8 iterations (case study 1). During the first 8 iterations (Fig. 13), the general shape of the final design already emerges. Early emergence of a clearly defined shape is also reported in [4][2] using density-based topology optimization. We believe that this phenomenon is related to plasmonic structures being objects of study. Slow convergence in later iterations is due to the objective function highly sensitive to boundary perturbations. For the same reason, the four holes persist through out the optimization 2 and trap the optimization into a local optimum. Better results than 3% can be obtained if these holes are manually removed from the design earlier in the optimization. We did another optimization where we manually removed the four holes at the 19 th iteration as shown in Fig. 14. Fig. 14 The design density (cid:2025) before (left) and after (right) the four holes are removed at the 19 th iteartion (case study 1). The four holes emerge again in the final design (Fig. 15) hence the final structure remains very similar to the previous result in Fig. 12. However, the (cid:1838) (cid:2870) norm difference achieves 0.5% after only 48 iterations, indicating a better result and faster convergence. 3 Fig. 15 The projected density (cid:2025)(cid:3548)̅ at the 48 th iteration in the optimization with four holes removed at the 19 th iteration (case study 1). CASE STUDY 2: MAXIMIZING OPTICAL ABSORPTION OF AN FEPT SLAB
In this example, topology optimization is performed to design a gold antenna to maximize the optical absorption within a
50 (cid:3400) 50 (cid:3400) 20 nm block of a 20 nm thick FePt slab. FePt is used in heat assisted magnetic recording as recording media where optical absorption within a small thermal spot is a very important figure of merit [30]. The antenna sits in parallel to the FePt slab with a 10 nm gap in-between and is placed in an open region where an incident wave of wavelength 800 nm is traveling in the direction perpendicular to the layer. Refractive index for FePt at 800 nm is (cid:4666)n (cid:3404) 3.2, k (cid:3404) 2.6(cid:4667) according to [31]. Again, we use a Yee cell length of 4 nm. Layout for the system is shown below. 4 Fig. 16 Layout for the optimization (case study 2). All rectangles are squares with their centers alighned. We define the objective function as the ratio of the absorption within surface (cid:1845) produced by the design structure to that produced without any structure (pure air): (cid:1832)(cid:4666)(cid:1779), (cid:1782)(cid:4667) (cid:3404) (cid:3505)Re(cid:4666)(cid:1779) (cid:3400) (cid:1782)(cid:3365) (cid:4667) ⋅ (cid:1866)(cid:3548)(cid:1856)(cid:1845) (cid:3020) (cid:3403) (cid:3505)Re(cid:4666)(cid:1779) (cid:3049)(cid:3042)(cid:3036)(cid:3031) (cid:3400) (cid:1782)(cid:3365) (cid:3049)(cid:3042)(cid:3036)(cid:3031) (cid:4667) ⋅ (cid:1866)(cid:3548)(cid:1856)(cid:1845) (cid:3020) (34) Fields (cid:4666)(cid:1779) (cid:3049)(cid:3042)(cid:3036)(cid:3031) , (cid:1782) (cid:3049)(cid:3042)(cid:3036)(cid:3031) (cid:4667) are calculated by setting design density (cid:2025) to be 0. The design region is a
200 (cid:3400) 200 (cid:3400) 60 nm block 10 nm away from the FePt slab (Fig. 16) and is discretized into a (cid:1840) (cid:3051) (cid:3400)(cid:1840) (cid:3052) (cid:3400) (cid:1840) (cid:3053) rectangular grid. Table 2 shows the parameters used in the optimization, which, compared to the previous case, does not change radius after each round. 5 Table 2 Parameters used in optimization (case study 2). Parameter [unit] Value (cid:1844) [Yee cell length] 6 (cid:2010) (cid:1863) (cid:2009) (cid:1855) (cid:3045) (cid:1855) (cid:3081) (cid:1840) (cid:3051) , (cid:1840) (cid:3052) , (cid:1840) (cid:3053)
51, 51, 16 (cid:1866) (cid:3046)(cid:3030)
3 The optimization algorithm starts with (cid:2025) (cid:3404) 0 , representing pure air, and runs for 8 rounds of optimization adding up to a total of 53 iterations and 91 function evaluations. Fig. 17 reveals the final design as a bowtie-like structure with an approximately 50 nm gap between the two tips. The optimized antenna provides an absorption ratio (34) of 4.9. We ran two more optimization with different initial density, (cid:2025) (cid:3404) 0.5 and (cid:2025) (cid:3404) 1.0 , and ended up with almost the same design and absorption in the end, which implies a certain degree of convexity in the optimization problem. In addition, the absorption of a bowtie antenna with the same dimensions to the bowtie aperture in previous section is calculated and its density (cid:2025) (cid:3029) is shown in Fig. 17. The absorption ratio for the bowtie antenna is evaluated to be 2.3 and is half of what the optimized design achieves. We believe gap and tip sizes ought to be close to the target area size to achieve optimal absorption. Fig. 18 shows some of the intermediate steps. 6 Fig. 17 Distributions of the design density (cid:2025) , the filtered density (cid:2025)(cid:3548) , the projected density (cid:2025)(cid:3548)̅ in the final iteration (case study 2), and the bowtie antenna density (cid:2025) (cid:3029) for comparison. 7 Fig. 18 Distributions of the projected density (cid:2025)(cid:3548)̅ in some intermediate steps showing evolution of the density (case study 2). Each iteration is marked with diamond in the objective function graph at the bottom. 8 CONCLUSION
We present a framework for applying topology optimization in plasmonics using the FDTD method as a frequency-domain solver. An original method of calculating discrete adjoint sensitivities of a general frequency-domain function is presented with detailed derivations and validated with a numerical example involving a plasmonic scatterer. Results match very well with sensitivities obtained with the central-difference scheme. A density-based topology optimization method with the FDTD method in frequency domain and the filtering-and-projection regularization is carried out to successfully recover near field patterns of a plasmonic bowtie aperture. We also proposed a procedure to implement the non-linear material interpolation in the FDTD method. This optimization method can be applied to problems with various lossless boundary conditions, including PEC, PMC and periodic boundary conditions.
ACKNOWLEDGEMENTS
Support to this work by ASTC – the Advanced Storage Technology Consortium and the National Science Foundation is acknowledged.
REFERENCE [1] J.A. Andkjaer, O.; Sigmund, S.; Nishiwaki, T. Nomura, Topology optimization of grating couplers for the efficient excitation of surface plasmons, J. B Opt. Phys. 27 (2010) 1828–1832. doi:10.1364/JOSAB.27.001828. [2] R.E. Christiansen, J. Vester-Petersen, S.P. Madsen, O. Sigmund, A non-linear material 9 interpolation for design of metallic nano-particles using topology optimization, Comput. Methods Appl. Mech. Eng. 343 (2019) 23–39. doi:10.1016/j.cma.2018.08.034. [3] Y. Deng, Z. Liu, C. Song, P. Hao, Y. Wu, Y. Liu, J.G. Korvink, Topology optimization of metal nanostructures for localized surface plasmon resonances, Struct. Multidiscip. Optim. 53 (2016) 967–972. doi:10.1007/s00158-015-1388-8. [4] Y. Deng, Z. Liu, C. Song, J. Wu, Y. Liu, Y. Wu, Topology optimization-based computational design methodology for surface plasmon polaritons, Plasmonics. 10 (2015) 569–583. doi:10.1007/s11468-014-9842-9. [5] S. Zhou, W. Li, Q. Li, Level-set based topology optimization for electromagnetic dipole antenna design, J. Comput. Phys. 229 (2010) 6915–6930. doi:10.1016/j.jcp.2010.05.030. [6] A. Taflove, S.C. Hagness, Computational Electrodynamics[1] J.A. Andkjaer, O.; Sigmund, S.; Nishiwaki, T. Nomura, Topology optimization of grating couplers for the efficient excitation of surface plasmons, J. B Opt. Phys. 27 (2010) 1828–1832. doi:10.1364/JOSAB.27.001828. [2] R.E. Christiansen, J. Vester-Petersen, S.P. Madsen, O. Sigmund, A non-linear material 9 interpolation for design of metallic nano-particles using topology optimization, Comput. Methods Appl. Mech. Eng. 343 (2019) 23–39. doi:10.1016/j.cma.2018.08.034. [3] Y. Deng, Z. Liu, C. Song, P. Hao, Y. Wu, Y. Liu, J.G. Korvink, Topology optimization of metal nanostructures for localized surface plasmon resonances, Struct. Multidiscip. Optim. 53 (2016) 967–972. doi:10.1007/s00158-015-1388-8. [4] Y. Deng, Z. Liu, C. Song, J. Wu, Y. Liu, Y. Wu, Topology optimization-based computational design methodology for surface plasmon polaritons, Plasmonics. 10 (2015) 569–583. doi:10.1007/s11468-014-9842-9. [5] S. Zhou, W. Li, Q. Li, Level-set based topology optimization for electromagnetic dipole antenna design, J. Comput. Phys. 229 (2010) 6915–6930. doi:10.1016/j.jcp.2010.05.030. [6] A. Taflove, S.C. Hagness, Computational Electrodynamics