Accurate adjoint design sensitivities for nano metal optics

Paul Hansen; Lambertus Hesselink

doi:10.1364/OE.23.023899

1. Introduction

Metal optical devices, ever more prominent in nano-photonics in the last decade, usually cannot be designed with pencil and paper as can simple lens-based optical systems. Instead, the strong evanescent fields, surface plasmons, intricate metal-dielectric boundary shapes and dispersive materials must be modeled using fully-vectorial FDTD and FEM codes, and these have been a mainstay of research and engineering for decades. Numerical modeling is used to investigate the physics of light-matter interaction at the nano-scale, but also to support the process of device design. An engineer with a CAD representation of a device can use FDTD or FEM to predict how it will perform under certain illumination conditions, and use the insight from simulation to attempt to improve the design. But to gauge the effect of a possible perturbation to the design one must know the sensitivity of the device’s behavior. This design sensitivity can be obtained by brute force, perturbing and re-simulating the design, which is computationally costly when the number of degrees of freedom is large. Adjoint design sensitivity analysis is a method to obtain the sensitivity to all perturbations of interest at once, for the computational cost of a single electromagnetic simulation. This method, adopted from its earlier uses in aeronautical design and control theory, has been making inroads in photonics for a decade, but while it has been straightforwardly implemented for dielectric structures it has seen limited use in metal optics due to the particular numerical challenges of modeling light-metal interaction: both light intensity and numerical error are highest at metal corners, jointly increasing the error of design sensitivities. Consequently design sensitivities are highly susceptible to the numerical error in metal optics.

We present a rigorous derivation of adjoint design sensitivities for optics using the mathematics of partial differential equations (PDEs), used in the preponderance of PDE-constrained optimization papers in other fields such as fluid mechanics. This approach is more extensible and powerful than existing approaches based explicitly on Maxwell’s equations because it does not rely on properties such as reciprocity, so it is applicable directly to other physical systems; and it provides practitioners of optics an introduction to the language of the field.

Furthermore we report the first highly-accurate adjoint design sensitivity analysis for metal optics. By using FEM instead of FDTD, the error of fields at interfaces is much reduced. To reduce the localized field enhancement that leads to errors we round all sharp corners to a small radius of curvature. We demonstrate the accuracy of this method by calculating sensitivities of a nano-structure called a C-shaped engraving (CSE). The CSE, recently employed for optical trapping [1], is a resonant nano-aperture with a reflecting back plane and is used for concentrating light into a sub-wavelength spot. After finding the sensitivities of the structure we optimize it to increase its field enhancement.

2. Related work

Although adjoint sensitivity analysis is not yet widespread in nanophotonic design, it has a long history in other engineering disciplines, having been used in control theory as early as 1971 [2]. It was adopted for fluid mechanics [3], aeronautical design [4–7] and structural design [8] through the 1980s. By the early 2000s it was applied to acoustics [9], to photonic band gap materials [10–12], and to microwave problems using FDTD and FEM [13–15]. Variants were proposed using FDTD and TLM [16–18]. It was applied to microlens design in 2009 [19]. In 2011, the present authors first presented an application of adjoint sensitivity analysis to two-dimensional dielectric and TE-mode metallic structures operating in the near-infrared regime, using FDTD [20, 21]. It was applied as well to dielectric light-trapping structures [22], mode splitters in two dimensions [23], and dielectric grating couplers [24]. However, its following extensions to plasmonic (TM-mode) metallic structures have been limited in scope and much less accurate. TLM has been applied for two-dimensional problems with strictly axis-aligned boundaries [25], but with reported error reaching 100%. Work on two-dimensional plasmonic structures [26] using commercial FDTD solvers have obtained pixel-by-pixel gradients accurate enough for optimization by a pixel-wise iterative metallization technique [17, 27]. Amorphous metal shapes have been optimized in 2D using a continuum design variable [28] and in 3D using a level-set approach [29], but there has been no report of accurate design sensitivities for the sharp, field-enhancing structures that make metal optics so appealing for nano-optics applications.

The obstacle to applying the adjoint method to metal optics lies in the difficulty of obtaining accurate field solutions near sharp edges and corners. FDTD, which has no explicit representation of material boundaries, is known to be inaccurate along interfaces (this has been remedied in the case of smooth dielectric boundaries [30] but not for corners or metal-dielectric interfaces). Furthermore the electromagnetic field at a sharp metal corner diverges, and divergent fields cannot be accurately represented by FDTD or even standard FEM using polynomial basis functions. The order of the divergence directly limits the rate of convergence of the field calculation [30], and the error is largest along interfaces between metals and dielectrics. Unfortunately, the adjoint sensitivity calculation for perturbing interfaces depends solely upon the interface electromagnetic fields (see Eq. (40)). We remedy the metal interface problem in two ways. First of all we completely eschew the low-accuracy FDTD and TLM approaches used for optics in almost all of the works cited above, and calculate the fields using FEM with cubic basis elements (COMSOL Multiphysics). The FEM basis functions can be evaluated at element edges to yield accurate fields on either side of a boundary even when field discontinuities are expected there. Secondly by rounding off sharp corners to a finite radius of curvature we eliminate the field divergences that spoil the accuracy of FEM and FDTD. We apply these techniques to the optimization of a highly field-enhancing metal nanostructure.

In Section 3 we present the basics of adjoint design sensitivity analysis and a rigorous derivation of the sensitivity to boundary perturbations for Maxwell’s equations. In Section 4 we discuss its implementation in FEM. In Section 5 the effects of several corner-rounding schemes on are presented. Finally in Section 6 we show how this method can be used to find design sensitivities for a metal nanostructure used for optical trapping, demonstrate < 5% error in sensitivities, and optimize the nanostructure to increase the energy density in a measurement volume by 70%.

2.1. Alternative optimization schemes

Sensitivity analysis is applicable to local optimization algorithms such as gradient descent and Newton’s method which iteratively improve a design until no further improvements can be found through perturbation. Unless the optimization problem is convex [31], the final design is not guaranteed to be the best one possible. Global optimization methods such as simulated annealing [32] and particle swarm optimization [33] can find globally optimal designs but converge very slowly. Adjoint design sensitivity analysis greatly reduces the runtime of local optimization algorithms but has no direct application to global optimization algorithms that make no use of gradients.

3. Adjoint sensitivity

When an optical system is perturbed by changing its shape, size, materials or illumination, its measurable behaviors change as well. The response of an optical device’s behavior to such a perturbation is called its design sensitivity. When a key behavior of the device varies smoothly with design parameters, the design sensitivity can be expressed as a derivative and solved for using calculus, and applied to device optimization by methods such as gradient descent. Here we consider the change of a behavior F to perturbations of several parameters p = (p₁, p₂,···) (where some parameters may represent shape, others material properties, etc.). In optimization contexts, F is an objective function representing the goodness of a device, and the best design is represented by the parameters p^* that minimize (or maximize) F. Naive calculation of all design sensitivities ∂F/∂p is expensive, as perturbations to each p_i need to be considered in turn, with one numerical simulation per parameter. Adjoint design sensitivity analysis, built on a linear algebra identity called duality [34], is a shortcut for calculating the sensitivity of a device to all perturbations of interest at very low computational cost.

The following derivation of the system sensitivity applies duality to Maxwell’s equations in a solver-agnostic way. Duality can also be applied directly to numerical algorithms such as FDTD and FEM. We extend the scheme of [35] to systems with internal material boundaries, working within the general framework of linear flux-conserving PDEs. (See as well a derivation for dispersive FDTD equations [36], a circuit analysis approach [16], and a derivation from Lorentz reciprocity [37].) Because adjoint methods have been very extensively explored in fluid mechanics, we hew to the language of flux-conserving PDE systems and show how Maxwell’s equations fit into that framework. This should facilitate the transfer of further adjoint-based methods from fluid mechanics into the field of optics (such as correction of functional estimates [38] and a posteriori error analysis [39]).

3.1. Duality in linear algebra

In linear algebra, duality is a relationship between solutions of two linear systems, called the primal system and the dual system. Given four vectors u, v, f and g, and a matrix A, the following identity holds,

g^{†} u = v^{†} f,

given that the vectors are related through

{\begin{array}{l} A u = f & (Primal system) \\ A^{†} v = g & (Dual system) . \end{array}

The primal and dual systems have the same form and same dimensions; the dual system is obtained from the primal system by taking the conjugate transpose, or adjoint, of A (and giving adjoint design sensitivity analysis its name). Duality is useful for sensitivities in the following way. Let A be a function of a vector of parameters p, and consider a parameterized linear system, called the forward system,

A (p) x = b (Forward system) .

We wish to know the sensitivity of some function F(x) to perturbations of each element p_i of p. The sensitivities follow by the chain rule, taking the form of computationally inexpensive inner products, analogous to Eq. (1),

\frac{\partial F}{\partial p_{i}} = \frac{\partial F}{\partial x} \frac{\partial x}{\partial p_{i}} for i \in (1, 2, \dots) .

Of the right-hand terms, ∂F/∂x can be obtained analytically, but ∂x/∂p_i must be solved for separately. Direct differentiation of the forward system (Eq. (3)) yields a set of additional, computationally expensive “primal” linear systems to solve for ∂x/∂p_i,

A \frac{\partial x}{\partial p_{i}} = \frac{\partial b}{\partial p_{i}} - \frac{\partial A}{\partial p_{i}} x for i \in (1, 2, \dots) .

Duality provides an alternative to laboriously solving all of these primal systems one at a time. Making the following identifications for each system in turn,

\begin{array}{l} u_{i} \equiv \frac{\partial x}{\partial p_{i}} \\ f_{i} \equiv \frac{\partial b}{\partial p_{i}} - \frac{\partial A}{\partial p_{i}} x \\ g_{i} \equiv {(\frac{\partial F}{\partial x})}^{†} \end{array}

it turns out that all of the primal systems, differing only in their inhomogeneous terms f_i, share the same dual system,

A {(p)}^{†} v = {(\frac{\partial F}{\partial x})}^{†}

and once v is calculated the sensitivities follow by computationally inexpensive “dual” inner products as in Eq. (1),

\frac{\partial F}{\partial p_{i}} = v^{†} (\frac{\partial b}{\partial p_{i}} - \frac{\partial A}{\partial p_{i}} x) for i \in (1, 2, \dots) .

Adjoint design sensitivity analysis refers to calculating system sensitivities like Eq. (4) more efficiently by replacing many expensive primal system solutions with a single dual system solution. The total computation required for calculation of all elements of ∂F/∂p is one forward calculation (Eq. (3)), one dual calculation (Eq. (7)), and one matrix multiplication and inner product per element of p (Eq. (8)).

3.2. Duality in partial differential equations

Adjoint design sensitivity analysis can also be applied to systems of partial differential equations such as Maxwell’s equations. Some of the analysis follows by analogy with linear algebra: the matrix A of the forward and primal problems is replaced by a partial differential operator ℒ, and the adjoint matrix A^† of the dual problem is replaced by the adjoint partial differential operator ℒ^†, obtained from ℒ not by transposition but by integrating by parts; the inner products g^†u = v^†f (Eq. (1)) become integrals. But PDE systems also feature boundary conditions, and the boundary conditions of the dual system are not necessarily just transposed versions of the primal boundary conditions but depend as well on ℒ and the form of the functional that gives ∂F/∂p_i. In the event that material boundaries and measurement domains can change shape, these too must be incorporated into the sensitivity calculation.

Before stating the dual expressions for PDEs, we introduce the forward problem (the physics of the device, e.g. Maxwell’s equations) and its mathematical representation. We represent a device design as a partition of all of space Ω into two domains Ω₁ and Ω₂ of distinct materials (Fig. 1), with the understanding that more material domains could be handled as well. Material properties of the device may vary smoothly inside each domain, but discontinuities of material properties are located only on the domain boundary ∂Ω₁₂. The domain boundaries and material properties in space are given as functions of parameters p which vary during the design process. The performance of a device is given by a merit functional F which we seek to minimize during the design process. It may depend directly on p but also on physical fields U, as F(U(p), p). Material discontinuities along ∂Ω₁₂ lead to field discontinuities as well, with U = U₁ on the Ω₁ side and U = U₂ on the Ω₂ side. We frequently make use of the jump in fields U²¹ ≡ U₂ − U₁ and the average field $\bar{U} \equiv \frac{1}{2} (U_{1} + U_{2})$ .

Fig. 1 Device design schematic. The entire domain Ω is partitioned into subdomains Ω₁ and Ω₂, with a material discontinuity on the interior interface. During the design process, the material parameters in Ω and the shape of the interface ∂Ω₁₂ and outer boundary ∂Ω may change. The normal vector n points out of outer ∂Ω and from Ω₁ to Ω₂ along ∂Ω₁₂.

Download Full Size | PDF

The fields U satisfy a system of partial differential equations, the forward system:

\begin{array}{l} ℒ (p) U & = S & in Ω \\ B (p) U & = S^{\partial Ω} & on outer \partial Ω \\ B_{12} (p) U^{21} & = S^{\partial Ω_{12}} & on interior \partial Ω_{12} . \end{array}

The physics in the device is represented by ℒ(p), a differential operator depending on design parameters p; and source currents S. The fields satisfy boundary conditions at the outer boundary ∂Ω of all space, and jump conditions along the inner domain boundary ∂Ω₁₂. Boundary conditions are encoded by matrices B which may vary with position. Maxwell’s equations and many other physical systems can be expressed in this framework: the conditions on ∂Ω are often periodic, perfect electric conductor (PEC) or perfect magnetic conductor (PMC); and the conditions on ∂Ω₁₂ are continuity of the tangential E and H fields. The boundary matrices B and B₁₂ are responsible for extracting the tangential field components.

The device performance is quantified by an objective functional F(U) such as a mode overlap integral or the total energy in a volume. Let F be a sum of such integrals over Ω, ∂Ω and ∂Ω₁₂,

F = \int_{Ω} F_{Ω} (U) + \oint_{\partial Ω} F_{\partial Ω} (U) + \int_{\partial Ω_{12}} F_{\partial Ω_{12}} (\bar{U}) .

The sensitivity of F to p_i is the primal functional of the PDE system, an integral of the form

\frac{d F}{d p_{i}} = {(g, u)}_{Ω} + {(C^{†} h, u)}_{\partial Ω} + {(C_{12}^{†} h_{12}, \bar{u})}_{\partial Ω_{12}} .

The measurement matrices C and C₁₂, not present in the discrete duality formulation, restrict the field components measurable on ∂Ω and ∂Ω₁₂ to those contributing to the flux normal to the boundary; this limitation on F is discussed in the appendix.

In the case that F does not directly depend on changes to the shape of Ω, then similarly to the discrete system

\begin{array}{l} g & = {(\frac{d F_{Ω}}{d U})}^{†} \\ C^{†} h & = {(\frac{d F_{\partial Ω}}{d U})}^{†} \\ C_{12}^{†} h_{12} & = {(\frac{d F_{\partial Ω_{12}}}{d \bar{U}})}^{†} . \end{array}

Like for the discrete system, ∂F/∂p can most straightforwardly be obtained by differentiating Eq. (9) with respect to each design parameter p_i in turn. This yields a partial differential equation not for U but for its sensitivity u ≡ dU/dp_i: the primal system

\begin{array}{l} ℒ (p) u & = f & in Ω \\ B (p) u & = e & on outer \partial Ω \\ B_{12} (p) u^{21} & = e_{12} & on interior \partial Ω_{12} . \end{array}

where simplifying variables have been introduced,

\begin{array}{l} f & = - \frac{\partial ℒ}{\partial p_{i}} U + \frac{\partial S}{\partial p_{i}} \\ e & = - \frac{\partial B}{\partial p_{i}} U + \frac{\partial S^{\partial Ω}}{\partial p_{i}} \\ e_{12} & = - \frac{\partial B_{12}}{\partial p_{i}} U^{21} + \frac{\partial S^{\partial Ω_{12}}}{\partial p_{i}} . \end{array}

The adjoint approach is to assemble a dual system,

\begin{array}{l} ℒ {(p)}^{†} v & = g & in Ω \\ β (p) v & = h & on \partial Ω \\ β_{12} (p) v^{21} & = h_{12} & on interior \partial Ω_{12} \end{array}

which is solved for the adjoint field v, and to evaluate the dual functional,

\frac{d F}{d p} = {(v, f)}_{Ω} + {(v, γ^{†} e)}_{\partial Ω} + {(\bar{v}, γ_{12}^{†} e_{12})}_{\partial Ω_{12}} .

The dual system has its own operator ℒ^† and its own boundary matrices β and β₁₂; the dual functional has its own measurement matrices γ and γ₁₂, mimicking the form of the primal functional. These components are determined from the primal system and primal functional. The most straightforward step is obtaining ℒ^† from ℒ through integrating by parts, analogous to transposing A for a linear algebra system. For fields u and v, the inner product (v, ℒu)_Ω is integrated by parts,

{(v, ℒ u)}_{Ω} = {(ℒ^{†} v, u)}_{Ω} + {(v, (n \cdot G) u)}_{\partial Ω} + {(v, (n \cdot G) u)}_{\partial Ω_{1}} + {(v, (n \cdot G) u)}_{\partial Ω_{2}} .

Integration by parts produces ℒ^† and a “flux matrix” n · G that acts along the domain boundaries ∂Ω and ∂Ω₁₂. The remaining dual system matrices β, β₁₂, γ and γ₁₂ are collectively determined from n · G, B, B₁₂, C and C₁₂ to satisfy

\begin{array}{l} n \cdot G = β^{†} C - γ^{†} B & on \partial Ω \\ n \cdot G = - β_{12}^{†} C_{12} = γ_{12}^{†} B_{12} & on \partial Ω_{12} . \end{array}

The full derivation is given in the appendix, and carried out explicitly for Maxwell’s equations in Section 3.4.

3.3. Sensitivity to boundary perturbations

The primal system Eqs. (13)–(14) accounts for perturbations inside Ω through the presence of ∂ℒ/∂p_i, which is nonzero at points in space where physics changes when p_i changes. It also accounts for perturbations on ∂Ω through the presence of ∂B/∂p_i, which is nonzero at points on the edge of Ω where boundary conditions change. To find the sensitivity of the system to boundary deformations, smooth perturbation of ∂Ω to ∂Ω′ can be equivalently cast as a change to the boundary operator B while ∂Ω is treated as though it is fixed.

Imagine a small displacement of the boundary from ∂Ω to ∂Ω′, wherein each point x on Ω shifts to x′ = x + η(x)n ∈ ∂Ω′. The original boundary condition is replaced by a new boundary condition, as follows (Fig. 2):

\begin{array}{l} B U (x) = S^{\partial Ω} (x) & \forall x \in \partial Ω & (original) \\ B^{'} U (x^{'}) = S^{\partial Ω^{'}} (x^{'}) & \forall x^{'} \in \partial Ω^{'} & (perturbed) . \end{array}

Each point x′ on ∂Ω′ “came from” a point x on ∂Ω. Because to each x′ ∈ ∂Ω′ there corresponds an x ∈ ∂Ω (given by x = x′ − η(x)n), we may write the boundary condition on the perturbed boundary as if it were a condition on the original unperturbed boundary,

B^{'} U (x + η (x) n) = S^{\partial Ω'} (x + η (x) n) \forall x \in \partial Ω .

Because ηn is the displacement of ∂Ω, differentiating Eq. (20) with respect to p_i yields the correct boundary conditions for Eq. (13) when boundaries may move. At η = 0, where also B′ = B and S^∂Ω′ = S^∂Ω, we have

\frac{\partial B}{\partial p_{i}} U (x) + B (\frac{\partial η}{\partial p_{i}} \partial_{n} U + \frac{\partial U}{\partial p_{i}}) = \frac{\partial S^{\partial Ω}}{\partial p_{i}} (x)

This is the boundary condition seen above in Eq. (13), but with an extra term arising from perturbation of ∂Ω by η(x)n. The boundary terms for the differentiated system with a movable boundary now read

\begin{array}{l} e & = - \frac{\partial B}{\partial p_{i}} U - B \frac{\partial η}{\partial p_{i}} \partial_{n} U + \frac{\partial S^{\partial Ω}}{\partial p_{i}} & on \partial Ω \\ e_{12} & = - \frac{\partial B_{12}}{\partial p_{i}} U^{21} - B_{12} \frac{\partial η}{\partial p_{i}} \partial_{n} U^{21} + \frac{\partial S^{\partial Ω}}{\partial p_{i}} & on \partial Ω_{12} . \end{array}

Fig. 2 Boundary condition on perturbed boundary. Each x ∈ ∂Ω becomes x′ ∈ ∂Ω′, and U must satisfy a new boundary condition on ∂Ω′ (because, for instance, the matrix to extract tangential components of U depends on the slope of ∂Ω).

Download Full Size | PDF

3.4. Maxwell’s equations and adjoint maxwell’s equations

Maxwell’s equations are the governing physics in optical design. Applying the framework of the previous section, we arrive at the operators ℒ and ℒ^† and suitable boundary conditions for the adjoint system. Due to Lorentz reciprocity, the adjoint PDE can be treated as a time-reversed version of Maxwell’s equations. This means that software for simulating electromagnetics can be applied as well to the adjoint problem without change, obviating the need for a special-purpose PDE solver. However, in nonreciprocal systems ℒ^† will need to be treated separately.

We shall begin by using the operator form of Maxwell’s equations to derive the adjoint operator and the flux matrix. The flux matrix and the primal boundary conditions suggest the form of the measurement matrix (and somewhat constrain the form of F), and the dual boundary conditions follow. The objective functional F can be evaluated once the forward system is solved for U = (E, H). The primal functional ∂F/∂p_i depends on the field sensitivities u = (E′, H′). It is more computationally efficient to obtain system sensitivities from the dual functional, which depends on dual fields v = (ℰ, ℋ). The electric displacement D = εE and magnetic induction B = μH will be used eventually for notational convenience, along with their adjoint counterparts 𝒟 = ε^†ℰ and ℬ = μ^†ℋ.

Orthogonal coordinates are defined along ∂Ω and ∂Ω₁₂, given by the normal direction n and two surface-tangential directions t₁ and t₂, illustrated in Fig. 3. Subscripts n, t₁ and t₂ will refer to field and gradient components along the n, t₁ and t₂ directions in the following.

Fig. 3 Coordinates on ∂Ω and ∂Ω₁₂. Fields and spatial derivatives will be represented in these coordinates.

Download Full Size | PDF

In matrix notation, Maxwell’s equations in the frequency domain read

\begin{array}{l} i ω [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] & = [\begin{matrix} - J^{Ω} \\ M^{Ω} \end{matrix}] & in Ω \\ [\begin{matrix} - t_{2}^{T} & 0 \\ t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] & = M^{\partial Ω} & on \partial Ω \\ [\begin{matrix} 0 & t_{2}^{T} \\ 0 & - t_{1}^{T} \\ t_{2}^{T} & 0 \\ - t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} E^{21} \\ H^{21} \end{matrix}] & = [\begin{matrix} - J^{\partial Ω_{12}} \\ M^{\partial Ω_{12}} \end{matrix}] & on \partial Ω_{12} . \end{array}

We identify the boundary condition matrices B and B₁₂,

B = [\begin{matrix} - t_{2}^{T} & 0 \\ t_{1}^{T} & 0 \end{matrix}], B_{12} = [\begin{matrix} 0 & t_{2}^{T} \\ 0 & - t_{1}^{T} \\ t_{2}^{T} & 0 \\ - t_{1}^{T} & 0 \end{matrix}] .

On the outer boundary ∂Ω, B enforces a PEC boundary condition. We define an integral objective functional F that depends on (E, H) in Ω, surface-tangential H on ∂Ω, and mean surface-tangential (Ē, H̄) on ∂Ω₁₂,

\begin{array}{l} F & = \int_{Ω} F_{Ω} ([\begin{matrix} E \\ H \end{matrix}]) + \oint_{\partial Ω} F_{\partial Ω} (C [\begin{matrix} E \\ H \end{matrix}]) + \oint_{\partial Ω_{12}} F_{\partial Ω_{12}} (C_{12} [\begin{matrix} \bar{E} \\ \bar{H} \end{matrix}]) \\ C & = [\begin{matrix} 0 & t_{1}^{T} \\ 0 & t_{2}^{T} \end{matrix}] \\ C_{12} & = [\begin{matrix} t_{1}^{T} & 0 \\ t_{2}^{T} & 0 \\ 0 & t_{1}^{T} \\ 0 & t_{2}^{T} \end{matrix}] . \end{array}

The measurement matrices C and C₁₂, required by duality [35] (see the appendix), project (E, H) down to only tangential fields. The primal functional comes from directly differentiating F with respect to p_i,

\frac{\partial F}{\partial p_{i}} = {(g, [\begin{matrix} E^{'} \\ H^{'} \end{matrix}])}_{Ω} + {(C^{†} h, [\begin{matrix} E^{'} \\ H^{'} \end{matrix}])}_{\partial Ω} + {(C_{12}^{†} h_{12}, [\begin{matrix} {\bar{E}}^{'} \\ {\bar{H}}^{'} \end{matrix}])}_{\partial Ω_{12}}

with constant fields

\begin{matrix} g = [\frac{{\frac{\partial F_{Ω}}{\partial E}}^{†}}{{\frac{\partial F_{Ω}}{\partial H}}^{†}}], & C^{†} h = [\frac{{\frac{\partial F_{\partial Ω}}{\partial E}}^{†}}{{\frac{\partial F_{\partial Ω}}{\partial H}}^{†}}], & C_{12}^{†} h_{12} = [\frac{{\frac{\partial F_{\partial Ω_{12}}}{\partial \bar{E}}}^{†}}{{\frac{\partial F_{\partial Ω_{12}}}{\partial \bar{H}}}^{†}}] \end{matrix} .

Directly differentiating the forward system with respect to a parameter p_i gives the primal system to solve for E′ = ∂E/∂p_i and H′ = ∂H/∂p_i,

\begin{array}{l} i ω [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] & = f & in Ω \\ [\begin{matrix} - t_{2}^{T} & 0 \\ t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] & = e & on \partial Ω \\ [\begin{matrix} 0 & t_{2}^{T} \\ 0 & - t_{1}^{T} \\ t_{2}^{T} & 0 \\ - t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} E^{21'} \\ H^{21'} \end{matrix}] & = e_{12} & on \partial Ω_{12} . \end{array}

To evaluate e and e₁₂ we use the equations for perturbed boundaries, Eq. (22). Perturbations to the boundary change the surface coordinates such that

\frac{\partial t_{1}}{\partial p_{i}} = n \frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}}

and

\frac{\partial t_{2}}{\partial p_{i}} = n \frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}}

. The perturbations to B and B₁₂ are then

\begin{array}{l} \frac{\partial B}{\partial p_{i}} & = [\begin{matrix} - n^{T} \frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}} & 0 \\ n^{T} \frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}} & 0 \end{matrix}], \\ \frac{\partial B_{12}}{\partial p_{i}} & = [\begin{matrix} 0 & n^{T} \frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}} \\ 0 & - n^{T} \frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}} \\ n^{T} \frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}} & 0 \\ - n^{T} \frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}} & 0 \end{matrix}] . \end{array}

The primal source currents f, e and e₁₂ are thus

\begin{array}{l} f & = - i ω [\begin{matrix} \frac{\partial ε}{\partial p_{i}} & 0 \\ 0 & - \frac{\partial μ}{\partial p_{i}} \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] + [\begin{matrix} - \frac{\partial J^{Ω}}{\partial p_{i}} \\ \frac{\partial M^{Ω}}{\partial p_{i}} \end{matrix}] m \\ e & = [\begin{matrix} (\frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}}) E_{n} + \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} E_{t_{2}} + \frac{\partial M_{t_{1}}^{\partial Ω}}{\partial p_{i}} \\ - (\frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}}) E_{n} - \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} E_{t_{1}} + \frac{\partial M_{t_{2}}^{\partial Ω}}{\partial p_{i}} \end{matrix}], \\ e_{12} & = [\begin{matrix} - (\frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}}) H_{n}^{21} - \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} H_{t_{2}}^{21} - \frac{\partial J_{t_{1}}^{\partial Ω 12}}{\partial p_{i}} \\ (\frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}}) H_{n}^{21} + \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} H_{t_{1}}^{21} - \frac{\partial t_{t_{2}}^{\partial Ω 12}}{\partial p_{i}} \\ - (\frac{\partial}{\partial t_{2}} \frac{\partial η}{\partial p_{i}}) E_{n}^{21} - \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} E_{t_{2}}^{21} + \frac{\partial M_{t_{1}}^{\partial Ω 12}}{\partial p_{i}} \\ (\frac{\partial}{\partial t_{1}} \frac{\partial η}{\partial p_{i}}) E_{n}^{21} + \frac{\partial η}{\partial p_{i}} \frac{\partial}{\partial n} E_{t_{1}}^{21} + \frac{\partial M_{t_{2}}^{\partial Ω 12}}{\partial p_{i}} \end{matrix}] . \end{array}

The differential operator of the primal system is

ℒ = i ω [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] .

Left-multiplying ℒu by v^† and integrating by parts produces the adjoint operator and a surface flux term,

\begin{array}{l} \int_{Ω} {[\begin{matrix} ℰ \\ ℋ \end{matrix}]}^{†} (i ω [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}]) \\ = \int_{Ω} {(- i ω [\begin{matrix} ε^{†} & 0 \\ 0 & - μ^{†} \end{matrix}] [\begin{matrix} ℰ \\ ℋ \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] [\begin{matrix} ℰ \\ ℋ \end{matrix}])}^{†} [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] \\ + \oint_{\partial Ω} {[\begin{matrix} ℰ \\ ℋ \end{matrix}]}^{†} [\begin{matrix} 0 & - n \times \\ - n \times & 0 \end{matrix}] [\begin{matrix} E^{'} \\ H^{'} \end{matrix}] \end{array}

where the adjoint operator is identified as

ℒ^{†} = - i ω [\begin{matrix} ε^{†} & 0 \\ 0 & - μ^{†} \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}]

and the flux matrix n · G ∈ ℝ^6×6 is

n \cdot G = [\begin{matrix} 0 & - n \times \\ - n \times & 0 \end{matrix}] = [\begin{matrix} 0 & t_{1} t_{2}^{T} - t_{2} t_{1}^{T} \\ t_{1} t_{2}^{T} - t_{2} t_{1}^{T} & 0 \end{matrix}] .

This flux matrix is antisymmetric, (n · G)^T = −n · G. As shown in the appendix, under these circumstances β = B, β₁₂ = B₁₂, γ = C and γ₁₂ = C₁₂, providing the boundary components of the dual system.

The dual system is constructed according to Eq. (15),

\begin{array}{l} i ω [\begin{matrix} ε^{†} & 0 \\ 0 & - μ^{†} \end{matrix}] [\begin{matrix} ℰ \\ ℋ \end{matrix}] + [\begin{matrix} 0 & - \nabla \times \\ - \nabla \times & 0 \end{matrix}] [\begin{matrix} ℰ \\ ℋ \end{matrix}] & = g & in Ω \\ [\begin{matrix} - t_{2}^{T} & 0 \\ t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} ℰ \\ ℋ \end{matrix}] & = h & on \partial Ω \\ [\begin{matrix} 0 & t_{2}^{T} \\ 0 & - t_{1}^{T} \\ t_{2}^{T} & 0 \\ - t_{1}^{T} & 0 \end{matrix}] [\begin{matrix} ℰ^{21} \\ ℋ^{21} \end{matrix}] & = h_{12} & on \partial Ω_{12} \end{array}

where g, h and h₁₂ are given by the primal functional.

The dual functional follows Eq. (16),

\begin{array}{l} \frac{\partial F}{\partial p_{i}} & = {([\begin{matrix} ℰ \\ ℋ \end{matrix}], f)}_{Ω} + {([\begin{matrix} ℰ \\ ℋ \end{matrix}], γ^{†} e)}_{\partial Ω} + {([\begin{matrix} ℰ \\ ℋ \end{matrix}], γ_{12}^{†} e_{12})}_{\partial Ω_{12}} \\ = \int_{Ω} [ℰ^{†} ℋ^{†}] (- i ω \frac{\partial}{\partial p_{i}} [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] + \frac{\partial}{\partial p_{i}} [\begin{matrix} - J^{Ω} \\ M^{Ω} \end{matrix}]) \\ + \oint_{\partial Ω} [ℋ_{t_{1}}^{†} ℋ_{t_{2}}^{†}] (E_{n} [\begin{matrix} \partial_{t_{2}} \\ - \partial_{t_{1}} \end{matrix}] \frac{\partial η}{\partial p_{i}} + \frac{\partial η}{\partial p_{i}} \partial_{n} [\begin{matrix} E_{t_{2}} \\ - E_{t_{1}} \end{matrix}] + \frac{\partial}{\partial p_{i}} [\begin{matrix} M_{t_{1}}^{\partial Ω} \\ M_{t_{2}}^{\partial Ω} \end{matrix}]) \\ + \oint_{\partial Ω_{12}} [{\bar{ℰ}}_{t_{1}}^{†} {\bar{ℰ}}_{t_{2}}^{†}] (H_{n}^{21} [\begin{matrix} - \partial_{t_{2}} \\ \partial_{t_{1}} \end{matrix}] \frac{\partial η}{\partial p_{i}} + \frac{\partial η}{\partial p_{i}} \partial_{n} [\begin{matrix} - H_{t_{2}}^{21} \\ H_{t_{1}}^{21} \end{matrix}] - \frac{\partial}{\partial p_{i}} [\begin{matrix} J_{t_{1}}^{\partial Ω 12} \\ J_{t_{2}}^{\partial Ω 12} \end{matrix}]) \\ + \oint_{\partial Ω_{12}} [{\bar{ℋ}}_{t_{1}}^{†} {\bar{ℋ}}_{t_{2}}^{†}] (E_{n}^{21} [\begin{matrix} - \partial_{t_{2}} \\ \partial_{t_{1}} \end{matrix}] \frac{\partial η}{\partial p_{i}} + \frac{\partial η}{\partial p_{i}} \partial_{n} [\begin{matrix} - E_{t_{2}}^{21} \\ E_{t_{1}}^{21} \end{matrix}] + \frac{\partial}{\partial p_{i}} [\begin{matrix} M_{t_{1}}^{\partial Ω 12} \\ M_{t_{2}}^{\partial Ω 12} \end{matrix}]) . \end{array}

By integration by parts, the spatial derivatives can be stripped from ∂η/∂p_i and the surface integrals can be reduced to products of forward and adjoint fields. Consider the integral over ∂Ω:

\begin{array}{l} \oint_{\partial Ω} [ℋ_{t_{1}}^{†} ℋ_{t_{2}}^{†}] (E_{n} [\begin{matrix} \partial t_{2} \\ - \partial_{t_{1}} \end{matrix}] \frac{\partial η}{\partial p_{i}} + \frac{\partial η}{\partial p_{i}} \partial_{n} [\begin{matrix} E_{t_{2}} \\ - E_{t_{1}} \end{matrix}] + \frac{\partial}{\partial p_{i}} [\begin{matrix} M_{t_{1}}^{\partial Ω} \\ M_{t_{2}}^{\partial Ω} \end{matrix}]) \\ = & \oint_{\partial Ω} \partial_{t_{2}} (ℋ_{t_{1}}^{†} E_{n} \frac{\partial η}{\partial p_{i}}) - \partial_{t_{1}} (ℋ_{t_{2}}^{†} E_{n} \frac{\partial η}{\partial p_{i}}) \\ + & \oint_{\partial Ω} \frac{\partial η}{\partial p_{i}} [- (\partial_{t_{2}} ℋ_{t_{1}}^{†}) E_{n} - ℋ_{t_{1}}^{†} \partial_{t_{2}} E_{n} + (\partial_{t_{1}} ℋ_{t_{2}}^{†}) E_{n} + ℋ_{t_{2}}^{†} \partial_{t_{1}} E_{n} + ℋ_{t_{1}}^{†} \partial_{n} E_{t_{2}} - ℋ_{t_{2}}^{†} \partial_{n} E_{t_{1}}] \\ + & \oint_{\partial Ω} ℋ_{t_{1}}^{†} \frac{\partial M_{t_{1}}^{\partial Ω}}{\partial p_{i}} + ℋ_{t_{2}}^{†} \frac{\partial M_{t_{2}}^{\partial Ω}}{\partial p_{i}} . \end{array}

Integrating over a closed surface, ∮ ∂_t₁() = ∮ ∂_t₂() = 0, which eliminates the first surface integral above. The remaining two surface integrals simplify under application of Maxwell’s equations and adjoint Maxwell’s equations in the surface coordinate system:

\begin{array}{l} i ω B_{t_{1}} = - \partial_{t_{2}} E_{n} + \partial_{n} E_{t_{2}} \\ i ω B_{t_{2}} = - \partial_{n} E_{t_{1}} + \partial_{t_{1}} E_{n} \\ i ω 𝒟_{n} = - \partial_{t_{1}} ℋ_{t_{2}} + \partial_{t_{2}} ℋ_{t_{1}} \end{array}

yielding

\oint_{\partial Ω} i ω \frac{\partial η}{\partial p_{i}} (ℋ_{t_{1}}^{†} B_{t_{1}} + ℋ_{t_{2}}^{†} B_{t_{2}} + 𝒟_{n}^{†} E_{n}) + ℋ_{t_{1}}^{†} \frac{\partial M_{t_{1}}^{\partial Ω}}{\partial p_{i}} + ℋ_{t_{2}}^{†} \frac{\partial M_{t_{2}}^{\partial Ω}}{\partial p_{i}} .

The same steps apply to all the surface integrals of Eq. (36), removing all the spatial derivatives so the dual sensitivity functional reduces to

\begin{array}{l} \frac{\partial F}{\partial p_{i}} & = \int_{Ω} [ℰ^{†} ℋ^{†}] (- i ω \frac{\partial}{\partial p_{i}} [\begin{matrix} ε & 0 \\ 0 & - μ \end{matrix}] [\begin{matrix} E \\ H \end{matrix}] + \frac{\partial}{\partial p_{i}} [\begin{matrix} - J^{Ω} \\ M^{Ω} \end{matrix}]) \\ + \oint_{\partial Ω} i ω \frac{\partial η}{\partial p_{i}} (𝒟_{n}^{†} E_{n} + ℋ_{t_{1}}^{†} B_{t_{1}} + ℋ_{t_{2}}^{†} B_{t_{2}}) \\ + \oint_{\partial Ω} ℋ_{t_{1}}^{†} \frac{\partial M_{t_{1}}^{\partial Ω}}{\partial p_{i}} + ℋ_{t_{2}}^{†} \frac{\partial M_{t_{2}}^{\partial Ω}}{\partial p_{i}} \\ + \oint_{\partial Ω_{12}} i ω \frac{\partial η}{\partial p_{i}} (- {\bar{𝒟}}_{n}^{†} E_{n}^{21} + {\bar{ℰ}}_{t_{1}}^{†} D_{t_{1}}^{21} + {\bar{ℰ}}_{t_{2}}^{†} D_{t_{2}}^{21} + {\bar{ℬ}}_{n}^{†} H_{n}^{21} - {\bar{ℋ}}_{t_{1}}^{†} B_{t_{1}}^{21} - {\bar{ℋ}}_{t_{2}}^{†} B_{t_{2}}^{21}) \\ + \oint_{\partial Ω_{12}} - {\bar{ℰ}}_{t_{1}}^{†} \frac{\partial J_{t_{1}}^{\partial Ω 12}}{\partial p_{i}} - {\bar{ℰ}}_{t_{2}}^{†} \frac{\partial J_{t_{2}}^{\partial Ω 12}}{\partial p_{i}} + {\bar{ℋ}}_{t_{1}}^{†} \frac{\partial M_{t_{1}}^{\partial Ω 12}}{\partial p_{i}} + {\bar{ℋ}}_{t_{2}}^{†} \frac{\partial M_{t_{2}}^{\partial Ω 12}}{\partial p_{i}} . \end{array}

Because objective functions are real-valued, only the real part of ∂F/∂p_i need be calculated.

4. Solving the adjoint system

Given a differential operator ℒ and a means to solve it, one may not be equipped yet to solve its adjoint ℒ^†. Because the form of ℒ^† is different than the form of ℒ, a solver written for ℒ (for example a finite-difference algorithm) will not solve ℒ^† and another code must be written for it. However in the case of Maxwell’s equations the form of ℒ^† is quite similar to that of ℒ,

\begin{array}{l} ℒ U = [\begin{matrix} i ω ε & - \nabla \times \\ - \nabla \times & - i ω μ \end{matrix}] U, & ℒ^{†} v = [\begin{matrix} - i ω ε^{†} & - \nabla \times \\ - \nabla \times & i ω μ^{†} \end{matrix}] v . \end{array}

The adjoint operator in this formulation of Maxwell’s equations is thus simply the complex conjugate of ℒ with transposed constitutive tensors. Often the constitutive tensors are diagonal, in which case to solve an adjoint system ℒ^†v = g in the frequency domain one can use two steps: solve ℒw = g^* using an ordinary Maxwell’s equations solver, then find v = w^*. If ε or μ is not a symmetric matrix then an operator ℒ′ equal to ℒ with transposed materials ε^T and μ^T can be used to find w.

Evaluation of all the elements of ∂F/∂p requires one solution of the forward system (Eq. (23)) for (E, H), one solution of the dual system (Eq. (35)) for (ℰ, ℋ), and one set of integrals for each parameter (Eq. (40)).

5. Accuracy of gradients

Use of adjoint sensitivities in photonics to date has mostly been limited to dielectric structures with limited field enhancement along boundaries. Metal-optics applications have mostly been limited to two-dimensions until recently [25, 26, 28, 40] due to the difficulty in arriving at sufficiently accurate field solutions along metal boundaries, edges and corners.

Sharp corners have been long understood to reduce the accuracy of electromagnetic field simulators [30], because the electric field E in the vicinity of a sharp point is weakly singular—it may grow to infinity, albeit remaining square integrable (so the total energy is finite). The approximation E_h found by FDTD or FEM is piecewise polynomial of degree p (for FDTD, linear E_h ∈ P₁(𝒞); for FEM, often p = 2 or p = 3, E_h ∈ P_p(𝒞)). As long as E is smooth, the error E − E_h converges as a power of h [41], with ‖E − E_h‖₂ ∝ h^2p.

However, the electric field at a distance ρ from a conductive point of angle β [42] grows without bound as ρ → 0,

u \propto ρ^{(π / β) - 1} .

Because this function is not smooth on an interval ρ ∈ [0, R], its approximation by a polynomial is accordingly poor in a neighborhood of the left endpoint. The maximum pointwise error ‖E − E_h‖_∞ is infinite and the 2-norm ‖E − E_h‖₂—although it may grow smaller with h—is no longer bounded by h^2p.

When E is inaccurate due to sharp corners, the sensitivity ∂E/∂p_i is inaccurate as well. Rounding off sharp features to a small but finite radius of curvature should thus improve gradient accuracy while reducing computational resources. This is illustrated below with a simple 2D TM simulation, in which the gradients are made more accurate by first chamfering the corners and then rounding the corners with polynomial splines. We follow by applying sensitivity analysis to a metallic nano-resonator with very high field enhancement, the C-shaped engraving (CSE) [1]. By rounding its sharp features with subdivision surfaces, accurate design sensitivities (error < 5%) are obtained, and the design sensitivities are employed optimize the CSE for higher field enhancement. Due to the finite point spread function of nano-fabrication methods such as lithography, practical realizations of nano-structures are inevitably rounded so this artificial corner rounding actually brings the simulations closer to agreement with reality.

5.1. Two-dimensional metal dimer

We use FEM to calculate one design sensitivity of a nanoparticle dimer in two dimensions (Fig. 4). Two 120 nm gold squares (ε = −44.7234 − 3.3160i) are placed corner-to-corner with a 50 nm gap. The squares are illuminated from above (+y), through free space, by a uniform x-polarized plane wave of amplitude E₀ with λ = 1000 nm. The electric field in the gap is strongly enhanced due to both the sharpness of the corners and their close proximity. Because the field enhancement in the gap can be very useful for applications, a useful figure of merit for the structure is the mean enhancement over a measurement box of area A in the gap:

F = \frac{1}{A} \int_{box} \frac{{| E_{x} |}^{2} + {| E_{y} |}^{2}}{E_{0}^{2}} .

We introduce a single design parameter, the displacement of the left particle along the x-axis from its initial position. It is displaced in 2.5 nm steps from −50 nm to 0 nm. The intensity is highest when the halves of the dimer are at their closest spacing (Fig. 5a). At each position a forward FEM calculation obtains F, and an adjoint FEM calculation obtains the adjoint fields needed to evaluate the sensitivity ∂F/∂p by an integral (Eq. (40)) over all sides of the moving particle. The “true” sensitivity of F is calculated numerically from the forward calculations (Fig. 5b).

Fig. 4 Electric field enhancement for translating metal box test. The left particle of the dimer will be translated to the left starting from this initial position. The gold permittivity is ε = −44.7234 − 3.3160i.

Download Full Size | PDF

Fig. 5 Left: F (a) and ∂F/∂x (b) as the leftmost particle of the dimer is swept along the x axis. The derivative ∂F/∂x is calculated numerically from the samples of F. By coincidence the “Chamfer 1” curves are closely overlapped by the “Fillet” curves. Right: Coarsest meshes (c) and finest meshes (d) used to test convergence of F and DF.

Download Full Size | PDF

To illustrate the effects of sharp corners on sensitivity calculations, four dimers are simulated: one with sharp corners, one with singly-chamfered corners (the corners are chopped off, leaving a 5 nm edge and two vertices), one with doubly-chamfered corners (the corners are replaced with three edges and four vertices, roughly a 10 nm radius of curvature), and one has filleted (rounded) corners with a radius of curvature of 5 nm. The filleted corner is represented in FEM by a spline and the local field calculation is accurately adjusted for the curvature. In general rounding the corners reduces F, here by as much as 15%. Reducing the radius of curvature will systematically reduce its impact on F, although for many nanoscale structures round corners are inevitable fabrication side-effects and should be incorporated into modeling from the beginning.

Corner rounding especially affects accuracy of design sensitivities. At each x displacement the adjoint design sensitivity ∂F/∂p is calculated and compared to the numerical derivative ΔF/Δp. The relative error in the sensitivity,

relative error = \frac{{‖ \partial F / \partial p - Δ F / Δ p ‖}_{2}}{{‖ Δ F / Δ p ‖}_{2}}

is expected to converge as a power of the element size h, and its exponent can be obtained by calculating the relative error at progressively higher mesh resolutions. Because fields vary most rapidly near material boundaries, here the mesh is refined only along the perimeter of the dimer and in the measurement volume. The selected mesh elements are refined by powers of two from a coarse mesh (characteristic size h ≈ 5 nm) to a fine mesh (h ≈ 5 nm/32). The meshes for each dimer variant are shown in Fig. 5c,d, and the relative errors in the gradient in Fig. 6.

Fig. 6 Relative error in ∂F/∂p as the mesh is refined (relative error = ‖∂F/∂p − ΔF/Δp‖₂/‖ΔF/Δp‖₂). Slopes are 0.32 (Sharp), 0.59 (Chamfer 1), 0.75 (Chamfer 2), and 2.1 (Fillet).

Download Full Size | PDF

The unrounded dimer is most problematic for sensitivity calculations. Even with sub-angstrom mesh elements the design sensitivity is wrong by 10–50%. Not only is the relative error large, but the order of convergence is deeply sublinear, with relative error ∝ h^0.32. Successive chamferings soften the corner and improve the error and its convergence rate; the doubly-chamfered dimer’s gradients are accurate to about 6.2% at the coarsest mesh, which is very useable for optimization. But rounding the corner with a spline is far and away superior to resolving the corner with straight segments. The relative error for the filleted dimer is about 0.26% on the coarse mesh and it scales as h^2.1.

Two conclusions can be drawn. Firstly, while FEM is generally expected to converge very rapidly with mesh refinement, reducing h is an impractical way to produce an accurate adjoint sensitivity for sharp metal structures. A small change to the geometry is a more practical approach. Secondly, handling round corners with conformal elements such as splines is far more effective than resolving round corners into flat segments, even when the tesselated corner looks very round. Additionally we note that the fields near sharp points in three dimensions grow faster than those around corners in two dimensions [42], and while h-refinement ultimately yielded marginally useable gradients without rounding in this 2D case, it might not be possible at all for some 3D geometries.

With accurate design sensitivities, three-dimensional metal-optics optimization becomes possible. We next turn to a nano-structure of practical interest: the C-shaped engraving (CSE).

6. C-shaped engraving

The CSE is a resonant metal-optical cavity formed by extruding a C-shaped hole partway into a metal layer. When illuminated from the top by x-polarized light at a resonant wavelength (980 nm for the design shown in Fig. 7), charge gathers at the sharp, central ridge on the metal surface and gives rise to strong optical near-fields. Much like for the C-shaped aperture [43], the CSE’s resonant wavelength depends on depth and perimeter; longer perimeters lead to longer resonant λ and Fabry-Perot-like resonances occur for some depths. Arrayed CSEs were recently used as optical traps in a nano-optical conveyor belt [1]. The design used in that work (Fig. 7) began with a canonical C-shaped aperture perimeter with a characteristic size of 60 nm (henceforth the “CSE-60” design), and its depth (150 nm) was determined subsequently by brute force, simulating CSE depths from 100 to 200 nm and selecting the design that maximized field enhancement. A CSE with higher field enhancement will produce a stiffer, more power-efficient optical trap. The CSE-60 was chosen by fixing a perimeter and varying the depth. However, a better design is expected to come of optimizing the shape of the entire inner surface. Gradient descent using adjoint sensitivities is an excellent approach, but as for 2D metal optics the gradients are inaccurate for a sharp-cornered structure.

Fig. 7 Electric field intensity enhancement in the C-shaped engraving (CSE) under 980 nm x-polarized illumination from above (+z). Its initial dimensions are a = 60 nm and d = 150 nm. The CSE has eight designable parameters, corresponding to displacements of its interior faces. The eight transverse displacement parameters F1 through F8 are illustrated.

Download Full Size | PDF

For illustrative purposes the CSE will first be constrained to remain C-shaped with flat, axis-aligned interior faces. The gradient accuracy will be calculated for perturbations of each facet in turn, for sharp-cornered and rounded structures. Then the design constraints will be relaxed and the surface optimized freely.

6.1. Sensitivity of sharp-cornered structures

The CSE’s performance will change as its interior facets F1–F8 are allowed to move in and out. Calculating the design sensitivities for each interior face of the CSE gives the direction to perturb the design to improve its behavior. But the calculated sensitivities of this sharp-cornered structure suffer too much numerical error to be useful, even at a very fine mesh resolution (Table 1). As in the 2D dimer test, the inaccuracy of the field solutions is associated with the inability of FEM basis elements, which are polynomial in nature, to resolve the weakly singular electric fields near sharp conductive corners [42]. Rounding off all the sharp features of the CSE to a finite radius of curvature dramatically improves the accuracy of the design sensitivities.

Table 1. Relative error of ∂F/∂p_i for three CSE variants. *Absolute sensitivity to p₂ is very low, as F varies about 1% from the smallest to largest value of p₂; the relative error in these near-zero values is accordingly very high.

View Table | View all tables in this article

Depending on the structure and the merit function, rounding will have a greater or lesser effect on F. The radius of curvature should be chosen to maintain sufficient accuracy of F while not unduly increasing computational burden. During meshing, the small mesh elements resolving corners will smoothly join the larger mesh elements of the model through a volume of higher-resolution mesh. Smaller corner elements will force the addition of more mesh elements and increase the time and space requirements for the calculation. For greatest computational efficiency the rounding radius should be the largest possible that keeps the change in F below a desired threshold. Nano-fabrication techniques inevitably introduce rounding, so numerically beneficial corner rounding may be desirable as well from an accuracy standpoint.

The merit function is the mean normalized field intensity in a measurement volume V,

F = \frac{1}{V} \int_{V} {‖ E ‖}^{2} / E_{0}^{2},

where E₀ is the electric field magnitude in the incident wave. To calculate the accuracy of the design sensitivities we select four internal faces of the CSE in turn (labeled 1, 2, 4 and 6 in Fig. 7) and translate them in and out by up to 10 nm. At 980 nm, the gradients obtained by the adjoint method are compared with the numerical derivative of the measured value. The sensitivities are grossly incorrect for this sharp-cornered structure (Figs. 8(a) and 8(d)). The source of the problem is apparent when looking at the field inside the CSE (Fig. 7). The maximum field intensity is found along the corners of the ridge of the CSE. If the field singularities were removed by gently rounding off the sharp features of the CSE, the simulation could converge at a lower mesh resolution and the sensitivity integrals could be accurate.

Fig. 8 Objective function, sensitivities and adjoint sensitivities of the CSE. Top: objective functions. Bottom: adjoint (dashed line) and measured (solid line) sensitivities. Left-to-right: sharp, singly-rounded and doubly-rounded structures.

Download Full Size | PDF

We quantify the rounding of the CSE by the radius of curvature of the corners, r. The sensitivity test is repeated now with all corners rounded to r = 10 nm, realistic for nano-fabrication today. The rounded meshes are generated from unrounded meshes by one or two iterations of Loop subdivision [44], corresponding to the single and double chamfers used in the dimer test. The values of F are very similar, but the sensitivities ∂F/∂p_i are only accurate for the rounded structures (Fig. 8).

While the sensitivity for the sharp-cornered structure is extremely inaccurate, the sensitivities for singly- and doubly-rounded CSEs are accurate enough to be applied for optimization. The strongest sensitivity ∂F/∂p₁, due to changing the length of the ridge, is an integral over enhanced fields near trihedral convex corners where the fields are least accurate. The doubly-rounded structure calculates this sensitivity to 3% error. Sensitivity to p₂ is low, and the absolute error in ∂F/∂p₂ is small despite its relatively large relative error. The remaining parameters p₄ and p₆ are significant contributors to changes in F but both correspond to concave interior faces of the CSE, where field enhancement is low and the accuracy is correspondingly high (error < 1.3% for the doubly-rounded structure).

7. Optimization of the CSE

The performance of the CSE will be improved by relaxing its shape constraints, allowing its perimeter to vary freely and even allowing its cross-section to vary with z. A new unconstrained CSE is given by a coarse mesh with 96 vertices (the “control vertices”) which are permitted to move in the x and y directions, so the design is not constrained to have axis-normal faces or even remain C-shaped. The mesh is smoothed everywhere with one iteration of Loop subdivision, and smoothed near the ridge with a second iteration of Loop subdivision to improve accuracy where the field intensity is highest.

The control vertex positions are optimized by steepest descent over 42 iterations (Fig. 9). The resulting design is resonant at 980 nm, and mean enhancement F has increased from 41 to 135. This outperforms the original CSE-60 design (F = 79) by 70%, in spite of its rounded features (Table 2).

Fig. 9 Top: Initial (left) and final optimized (right) unconstrained CSE designs. The optimized design is over 3× better than the initial design and 70% better than the original CSE-60 design. Bottom: progress of optimization over 42 iterations. Arrows mark the designs shown above. The horizontal dashed line at F = 79 marks the performance of the hand-tuned CSE-60, surpassed by the sixth iteration.

Download Full Size | PDF

Table 2. Performance of resonant CSE-60 compared to initial and optimized unconstrained CSE designs.)

View Table | View all tables in this article

Optimization significantly changes the shape of the aperture. The optimized design is undercut so its cross-section is smaller at the entrance (z = 0) than in the center of the aperture (around z = −75 nm). The top of the ridge has extended closer to the flat back of the aperture, strengthening the field in the gap, and the ridge has also become narrower in the y direction and sharpened to a point in the z direction. The sharpening enhances the “lightning rod” behavior of the ridge, gathering electric charge tightly underneath the measurement volume. The concentration of field towards the top of the aperture is visible in the cross-section plots in Fig. 10. The undercut of the ridge also appears to reduce Ohmic losses in the deeper portion of the CSE. As the ridge recedes, the field enhancement in the aperture lowers and the electric current in the aperture walls decreases. Because the electric field is nearly zero at the metal floor of the aperture (z = −150 nm), it is not necessary to pull the ridge back to reduce losses and the aperture perimeter remains essentially unchanged.

Fig. 10 Top: field enhancement in standard CSE-60. Bottom: field enhancement in optimized undercut CSE. The xy cross-section (left) is taken at z = 20 nm. The yz (center) and xz (right) cross-sections are taken through x = 0 and y = 0, respectively. Dimensions are in nm.

Download Full Size | PDF

We note as well that the optimized structure is slightly asymmetrical in the y direction (Fig. 10, bottom-center). The symmetry breaking between the initial and final designs is a result of asymmetry in the computational mesh.

It is exciting to note that even though the C-aperture has been studied since 2000, elongated and squeezed [45], rounded [46], tapered [47], bent [48], and in other variations (and furthermore its kin the bowtie aperture [49], H-aperture [50], L-aperture [51], and so on as well), the undercut design here is a new idea, produced only with the aid of sensitivity analysis. Furthermore, this CSE is only the first metal structure we have sought to improve with adjoint design sensitivity analysis. The potential for new nano-structure designs is evidently vast.

8. Conclusion

Design of metal nano-structures for optics presents an array of challenges to the engineer, such as the lack of analytically solvable models, complicated material properties, and the need to solve Maxwell’s equations without approximations. Computation is a mainstay of nano optical design today, and stands to benefit from the mathematics of adjoint design sensitivity analysis developed over several decades in other fields of engineering. Starting from a general linear PDE framework we have derived the equations of sensitivity analysis with moving internal material boundaries and applied them to Maxwell’s equations. The electromagnetics solvers commonly applied to nano-optics problems cannot accurately resolve the fields near sharp metal corners, making it difficult to use adjoint sensitivities to optimize metal optics devices. Rounding the sharp corners of these structures to a small but finite radius of curvature reduces or eliminates the field singularities and restores the accuracy of the field solvers, drastically improving the accuracy of adjoint sensitivities; rounded structures are also more realizable with practical tools. We use adjoint sensitivities to tune a resonant metal nanostructure, including the sharp ridge that contributes the most field enhancement and the most numerical error. This technique is broadly applicable for electromagnetics, including devices for communications, data storage, biosensing, and optical trapping.

A. Conditions for Duality with Jump Conditions

Let the primal system take the form

\begin{array}{l} ℒ u & = f & in Ω \\ B u & = e & on outer \partial Ω \\ B_{12} (u_{2} - u_{1}) & = e_{12} & on interior \partial Ω_{12} \end{array}

and the dual system

\begin{array}{l} ℒ^{†} v & = g & in Ω \\ β v & = h & on \partial Ω \\ β_{12} (v_{2} - v_{1}) & = h_{12} & on interior \partial Ω_{12} . \end{array}

The primal functional is an inner product with u,

I_{primal} = {(g, u)}_{Ω} + {(h, C u)}_{\partial Ω} + {(h_{12}, C_{12} \bar{u})}_{\partial Ω_{12}}

and the dual functional is an inner product with v,

I_{dual} = {(v, f)}_{Ω} + {(γ v, e)}_{\partial Ω} + {(γ_{12} \bar{v}, e_{12})}_{\partial Ω_{12}}

where the average fields on the interface are

\bar{u} = \frac{1}{2} (u_{1} + u_{2})

and

\bar{v} = \frac{1}{2} (v_{1} + v_{2})

.

Under conditions on B, B₁₂, C, C₁₂, β, β₁₂, γ, and γ₁₂, the primal and dual functionals have the same value. The primal system and primal functional specify ℒ, B, B₁₂, C and C₁₂. Then the adjoint operator ℒ^† and dual system matrices β, β₁₂, γ and γ₁₂ can be found from the primal system. The duality conditions relating the primal and dual systems are found in two steps. First ℒ^† is found from ℒ by integration by parts, which produces surface integrals involving a flux matrix n · G. Second, equating the primal and dual functionals and substituting in terms with u and v from the primal and dual systems yields alternative expressions for the same surface integrals. Matching the integrands relates the primal and dual matrices through the flux matrix n · G, giving the duality conditions. These steps were carried out in [35] for exterior flow problems, and are extended here to the material interface ∂Ω₁₂ and the matrices B₁₂, C₁₂, β₁₂ and γ₁₂.

ℒ^† is obtained by taking the inner product of Eq. (46) with another field v and integrating by parts. A first-order flux-conserving system such as Maxwell’s equations will produce a surface term as a byproduct as well, with a surface normal-dependent flux matrix n · G.

{(v, ℒ u)}_{Ω} = {(ℒ^{†} v, u)}_{Ω} + {(v, (n \cdot G) u)}_{\partial Ω} + {(v, (n \cdot G) u)}_{\partial Ω_{1}} + {(v, (n \cdot G) u)}_{\partial Ω_{2}} .

Because the integrals over ∂Ω₁ and ∂Ω₂ are calculated on opposite sides of ∂Ω₁₂, they can be written as such, denoting the fields on the two sides of the interface by u₁, v₁, u₂ and v₂ and taking the surface normal to be n = n₁ = −n₂.

{(v, ℒ u)}_{Ω} = {(ℒ^{†} v, u)}_{Ω} + {(v, (n \cdot G) u)}_{\partial Ω} + {(v_{1}, (n \cdot G) u_{1})}_{\partial Ω_{12}} - {(v_{2}, (n \cdot G) u_{2})}_{\partial Ω_{12}} .

Next, the difference between the primal and dual functionals is taken to be zero by hypothesis:

0 = I_{primal} - I_{dual} = {(g, u)}_{Ω} - {(v, f)}_{Ω} + {(h, C u)}_{\partial Ω} + {(h_{12}, C_{12} \bar{u})}_{\partial Ω_{12}} - {(γ v, e)}_{\partial Ω} - {(γ_{12} \bar{v}, e_{12})}_{\partial Ω_{12}} .

The constant terms f, g, e, e₁₂, h, and h₁₂ can be substituted in from Eqs. (46) and (47), and one term drawn to the left hand side to match the form of Eq. (51):

{(v, ℒ u)}_{Ω} = {(ℒ^{†} v, u)}_{Ω} - {(γ v, B u)}_{\partial Ω} + {(β v, C u)}_{\partial Ω} - {(γ_{12} \bar{v}, B_{12} (u_{2} - u_{1}))}_{\partial Ω_{12}} + {(β_{12} (v_{2} - v_{1}), C_{12} \bar{u})}_{\partial Ω_{12}} .

The surface integral terms in Eqs. (51) and (53) must match at all points. This produces a duality condition on the exterior boundary [35],

n \cdot G = β^{†} C - γ^{†} B on Ω_{12},

and an equality that must hold over the inner boundaries,

\begin{array}{l} {(γ_{12} \frac{v_{1} + v_{2}}{2}, B_{12} (u_{2} - u_{1}))}_{\partial Ω_{12}} - {(β_{12} (v_{2} + v_{1}), C_{12} \frac{u_{1} - u_{2}}{2})}_{\partial Ω_{12}} \\ = - {(v_{1}, (n \cdot G) u_{1})}_{\partial Ω_{12}} + {(v_{2}, (n \cdot G) u_{2})}_{\partial Ω_{12}} . \end{array}

When expressed as a matrix there are four conditions to solve together,

\int_{\partial Ω_{12}} {[\begin{matrix} v_{1} \\ v_{2} \end{matrix}]}^{T} [\begin{matrix} - n \cdot G \\ - n \cdot G \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \end{matrix}] = \int_{\partial Ω_{12}} {[\begin{matrix} v_{1} \\ v_{2} \end{matrix}]}^{T} \frac{1}{2} [\begin{matrix} β_{12}^{T} C_{12} - γ_{12}^{T} B_{12} & β_{12}^{T} C_{12} + γ_{12}^{T} B_{12} \\ - β_{12}^{T} C_{12} - γ_{12}^{T} B_{12} & γ_{12}^{T} B_{12} - β_{12}^{T} C_{12} \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \end{matrix}] .

Equating the central matrices element by element we obtain the conditions for duality to hold with jump conditions:

n \cdot G = - β_{12}^{T} C_{12} = γ_{12}^{T} B_{12} .

In the case of Maxwell’s equations, the flux matrix is asymmetric, (n · G)^T = −n · G. In such a case the duality conditions may be satisfied by β = B, β₁₂ = B₁₂, γ = C and γ₁₂ = C₁₂.

Not all the components of u contribute to the flux through an interface or boundary, (n · G)u. In the case of Maxwell’s equations, only the four tangential fields E_‖ and H_‖ satisfy a continuity equation at interfaces and contribute to the Poynting vector, and only parallel fields are present in the PEC, PMC and radiation boundary conditions; this is because only these fields are in the row space of n · G. Indeed if necessary the normally-directed fields E_⊥ and H_⊥ may be solved for in terms of E_‖ and H_‖. One of the functions of B and C both is to project the full field u on ∂Ω or ∂Ω₁₂ onto the space of fields that contribute to the flux (the transverse fields in the case of Maxwell’s equations). On exterior boundaries, Eq. (54) implies that B and C together must project out all fields that do not contribute to the flux through ∂Ω. On interior boundaries ∂Ω₁₂, Eq. (55) implies that B₁₂ and C₁₂ must both project u down to its flux-contributing components. Thus B₁₂ must measure all four components of E_‖ and H_‖ in Maxwell’s equations, and C₁₂ likewise.

References and links

1. Y. Zheng, J. Ryan, P. Hansen, Y.-T. Cheng, T.-J. Lu, and L. Hesselink, “Nano-optical conveyor belt, part ii: Demonstration of handoff between near-field optical traps,” Nano letters 14, 2971–2976 (2014). [CrossRef] [PubMed]

2. J. L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, vol. 170 (Springer Verlag, 1971). [CrossRef]

3. O. Pironneau, “On optimum design in fluid mechanics,” J. Fluid Mech. 64, 97–110 (1974). [CrossRef]

4. A. Jameson, “Aerodynamic design via control theory,” J. Sci. Comput. 3, 233–260 (1988). [CrossRef]

5. A. Jameson, “Optimum aerodynamic design using cfd and control theory,” AIAA paper 1729, 124–131 (1995).

6. A. Jameson, L. Martinelli, and N. Pierce, “Optimum aerodynamic design using the navier–stokes equations,” Theor. Comp. Fluid Dyn. 10, 213–237 (1998). [CrossRef]

7. A. Jameson, “Re-engineering the design process through computation,” J. Aircraft 36, 36–50 (1999). [CrossRef]

8. V. Komkov, K. K. Choi, and E. J. Haug, Design Sensitivity Analysis of Structural Systems, vol. 177 (Academic, 1986).

9. O. Sigmund and J. S. Jensen, “Systematic design of phononic band–gap materials and structures by topology optimization,” Philos. Trans. R. Soc. London, Ser. A 361, 1001–1019 (2003). [CrossRef]

10. J. S. Jensen and O. Sigmund, “Systematic design of photonic crystal structures using topology optimization: Low-loss waveguide bends,” Appl. Phys. Lett. 84, 2022–2024 (2004). [CrossRef]

11. P. Borel, A. Harpøth, L. Frandsen, M. Kristensen, P. Shi, J. Jensen, and O. Sigmund, “Topology optimization and fabrication of photonic crystal structures,” Opt. Express 12, 1996–2001 (2004). [CrossRef] [PubMed]

12. L. Frandsen, A. Harpøth, P. Borel, M. Kristensen, J. Jensen, and O. Sigmund, “Broadband photonic crystal waveguide 60 bend obtained utilizing topology optimization,” Opt. Express 12, 5916–5921 (2004). [CrossRef] [PubMed]

13. Y.-S. Chung, C. Cheon, I.-H. Park, and S.-Y. Hahn, “Optimal shape design of microwave device using fdtd and design sensitivity analysis,” IEEE Trans. Microwave Theory Tech. 48, 2289–2296 (2000). [CrossRef]

14. Y.-S. Chung, C. Cheon, I.-H. Park, and S.-Y. Hahn, “Optimal design method for microwave device using time domain method and design sensitivity analysis. ii. fdtd case,” IEEE Trans. Magn. 37, 3255–3259 (2001). [CrossRef]

15. Y.-S. Chung, J. Ryu, C. Cheon, I.-H. Park, and S.-Y. Hahn, “Optimal design method for microwave device using time domain method and design sensitivity analysis. i. fetd case,” IEEE Trans. Magn. 37, 3289–3293 (2001). [CrossRef]

16. N. K. Nikolova, J. W. Bandler, and M. H. Bakr, “Adjoint techniques for sensitivity analysis in high-frequency structure cad,” IEEE Trans. Microwave Theory Tech. 52, 403–419 (2004). [CrossRef]

17. N. K. Nikolova, H. W. Tam, and M. H. Bakr, “Sensitivity analysis with the fdtd method on structured grids,” IEEE Trans. Microwave Theory Tech. 52, 1207–1216 (2004). [CrossRef]

18. N. K. Nikolova, Y. Li, Y. Li, and M. H. Bakr, “Sensitivity analysis of scattering parameters with electromagnetic time-domain simulators,” IEEE Trans. Microwave Theory Tech. 54, 1598–1610 (2006). [CrossRef]

19. Y.-S. Chung, B.-J. Lee, and S.-C. Kim, “Optimal shape design of dielectric micro lens using fdtd and topology optimization,” J. Opt. Soc. Korea 13, 286–293 (2009). [CrossRef]

20. P. Hansen, Y. Zheng, E. Perederey, and L. Hesselink, “Nanophotonic device optimization with adjoint fdtd,” in “CLEO: Applications and Technology,” (Optical Society of America, 2011), p. JTuI61. [CrossRef]

21. P. Hansen, Y. Zheng, E. Perederey, and L. Hesselink, “Adjoint fdtd for nanophotonic device optimization,” in “Joint International Symposium on Optical Memory and Optical Data Storage,” (Optical Society of America, 2011), p. OTuE2.

22. O. D. Miller, V. Ganapati, and E. Yablonovitch, “Inverse design of a nano-scale surface texture for light trapping,” in “CLEO: Science and Innovations,” (Optical Society of America, 2012), pp. CF2J–2.

23. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express 21, 21693–21701 (2013). [CrossRef] [PubMed]

24. A. C. Niederberger, D. A. Fattal, N. R. Gauger, S. Fan, and R. G. Beausoleil, “Sensitivity analysis and optimization of sub-wavelength optical gratings using adjoints,” Opt. Express 22, 12971–12981 (2014). [CrossRef] [PubMed]

25. O. Ahmed, M. Bakr, X. Li, and T. Nomura, “Adjoint variable method for two-dimensional plasmonic structures,” Opt. Lett. 37, 3453–3455 (2012). [CrossRef]

26. S. Bhargava, O. Miller, V. Ganapati, and E. Yablonovitch, “Inverse design of optical antennas for sub-wavelength energy delivery,” in “CLEO: Science and Innovations,” (Optical Society of America, 2013), pp. CM2F–2.

27. M. S. Dadash, “Simulator independent exact adjoint sensitivity analysis of self-adjoint microwave structures,” Master’s thesis, McMaster University (2011).

28. Y. Deng, Z. Liu, C. Song, J. Wu, Y. Liu, and Y. Wu, “Topology optimization-based computational design methodology for surface plasmon polaritons,” Plasmonics 10, 1–15 (2014).

29. S. Bhargava and E. Yablonovitch, “Multi-objective inverse design of sub-wavelength optical focusing structures for heat assisted magnetic recording,” in “SPIE Optical Engineering+ Applications,” (International Society for Optics and Photonics, 2014), 92010M.

30. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. Joannopoulos, S. G. Johnson, and G. Burr, “Improving accuracy by subpixel smoothing in the finite-difference time domain,” Opt. Lett. 31, 2972–2974 (2006). [CrossRef] [PubMed]

31. S. P. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University, 2004). [CrossRef]

32. L. Ingber, et al., “Adaptive simulated annealing (asa): Lessons learned,” Control and cybernetics 25, 33–54 (1996).

33. R. Poli, J. Kennedy, and T. Blackwell, “Particle swarm optimization,” Swarm Intell. 1, 33–57 (2007). [CrossRef]

34. M. B. Giles and N. A. Pierce, “An introduction to the adjoint approach to design,” Flow Turbul. Combust. 65, 393–415 (2000). [CrossRef]

35. M. B. Giles and N. A. Pierce, “Adjoint equations in cfd: duality, boundary conditions and solution behaviour,” AIAA paper 97, 1850 (1997).

36. E. Abenius and B. Strand, “Solving inverse electromagnetic problems using fdtd and gradient-based minimization,” International journal for numerical methods in engineering 68, 650–673 (2006). [CrossRef]

37. O. D. Miller, “Photonic design: from fundamental solar cell physics to computational inverse design,” arXiv preprint arXiv:1308.0212 (2013).

38. N. A. Pierce and M. B. Giles, “Adjoint and defect error bounding and correction for functional estimates,” J. Comput. Phys. 200, 769–794 (2004). [CrossRef]

39. M. B. Giles and E. Süli, “Adjoint methods for pdes: a posteriori error analysis and postprocessing by duality,” Acta Numer. 11, 145–236 (2002). [CrossRef]

40. Y. Zhang, O. S. Ahmed, and M. H. Bakr, “Adjoint sensitivity analysis of plasmonic structures using the fdtd method,” Optics letters 39, 3002–3005 (2014). [CrossRef] [PubMed]

41. J. Jin, The Finite Element Method in Electromagnetics (John Wiley & Sons, 2014).

42. J. D. Jackson, “Classical electrodynamics,” Classical Electrodynamics 1, 832 (1998).

43. X. Shi, L. Hesselink, and R. L. Thornton, “Ultrahigh light transmission through a C-shaped nanoaperture,” Opt. Lett. 28, 1320–1322 (2003). [CrossRef] [PubMed]

44. C. Loop, “Smooth subdivision surfaces based on triangles,” Master’s thesis, University of Utah (1987).

45. X. Shi and L. Hesselink, “Design of a C aperture to achieve λ/10 resolution and resonant transmission,” J. Opt. Soc. Am. B 21, 1305–1317 (2004). [CrossRef]

46. J. B. Leen, P. Hansen, Y.-T. Cheng, and L. Hesselink, “Improved focused ion beam fabrication of near-field apertures using a silicon nitride membrane,” Opt. Lett. 33, 2827–2829 (2008). [CrossRef] [PubMed]

47. J. B. Leen, E. Cho, S.-D. Suh, P. C. Hansen, J.-S. Sohn, S.-H. Choa, and L. Hesselink, “90 bent metallic waveguide with a tapered c-shaped aperture for use in hamr,” in “Optical Data Storage 2007,” (International Society for Optics and Photonics, 2007), pp. 66200R. [CrossRef]

48. P. Hansen, L. Hesselink, and B. Leen, “Design of a subwavelength bent c-aperture waveguide,” Opt. Lett. 32, 1737–1739 (2007). [CrossRef] [PubMed]

49. E. X. Jin and X. Xu, “Obtaining super resolution light spot using surface plasmon assisted sharp ridge nanoaperture,” Appl. Phys. Lett. 86, 111106 (2005). [CrossRef]

50. X. Xu, E. X. Jin, and S. M. Uppuluri, “Enhancement of optical transmission through planar nanoapertures in a metal film,” in “Optical Science and Technology, the SPIE 49th Annual Meeting,” (International Society for Optics and Photonics, 2004), pp. 230–243.

51. J. Xu, T. Xu, J. Wang, and Q. Tian, “Design tips of nanoapertures with strong field enhancement and proposal of novel l-shaped aperture,” Opt. Eng. 44, 018001 (2005). [CrossRef]

Structure	p ₁	p ₂	p ₄	p ₆
Sharp	475	224	62.9	96.8
Singly-rounded	12.5	164^*	0.845	2.84
Doubly-rounded	3.07	126^*	0.69	1.28

Structure	p ₁	p ₂	p ₄	p ₆
Sharp	475	224	62.9	96.8
Singly-rounded	12.5	164^*	0.845	2.84
Doubly-rounded	3.07	126^*	0.69	1.28

Accurate adjoint design sensitivities for nano metal optics

Abstract

1. Introduction

2. Related work

2.1. Alternative optimization schemes

3. Adjoint sensitivity

3.1. Duality in linear algebra

3.2. Duality in partial differential equations

3.3. Sensitivity to boundary perturbations

3.4. Maxwell’s equations and adjoint maxwell’s equations

4. Solving the adjoint system

5. Accuracy of gradients

5.1. Two-dimensional metal dimer

6. C-shaped engraving

6.1. Sensitivity of sharp-cornered structures

7. Optimization of the CSE

8. Conclusion

A. Conditions for Duality with Jump Conditions

References and links

Cited By

Figures (10)

Tables (2)

Equations (57)

Optics Express