Compact 200 line MATLAB code for inverse design in photonics by topology optimization: tutorial

Rasmus E. Christiansen; Rasmus E. Christiansen; Ole Sigmund; Ole Sigmund

doi:10.1364/JOSAB.405955

1. INTRODUCTION

This paper details a 200 line MATLAB code, which demonstrates how density-based topology optimization (TopOpt) can be applied to photonics design. The code is written for scientists and students with a basic knowledge of programming, numerical modeling, and photonics who desire to start using inverse design in their research. We briefly detail the model of the physics, followed by the discretized TopOpt design problem (Section 2). The MATLAB code is then explained in detail (Section 3), followed by two application examples providing the reader with targets for reproduction (Section 4). Then, a number of possible extensions are discussed, and code-snippets for easy implementation are provided along with design examples (Section 5). Finally, we demonstrate the superiority of using gradient-based TopOpt compared to a genetic algorithm (GA) when solving a photonic design problem (Section 6).

TopOpt [1] as an inverse design tool was first developed in the context of solid mechanics in the late 1980s [2]. Since its inception, the method has developed rapidly and expanded across most areas of physics [3–6]. Over the last two decades, the interest in applying TopOpt for photonics has increased rapidly [7] with applications within cavity design [8,9], photonic demultiplexers [10], metasurfaces [11,12] and topological insulators [13] to name a few. While the interest in TopOpt within the photonics community has grown markedly in recent years, significant barriers hinder newcomers to the field from adopting the tool in their work. These are as follows: the required knowledge of numerical modeling, the required knowledge of advanced mathematical concepts, and scientific programming experience. This paper seeks to lower these barriers by providing the reader with a simple two-dimensional (2D) finite-element-based MATLAB implementation of TopOpt for photonics, which is straightforwardly extendable to a range of other design problems. Within the field of structural optimization in mechanics, similar simple MATLAB codes [14–16] have proven themselves highly successful in raising the awareness of TopOpt and serving as a basic platform for further expansion of the method, thus broadening its application as a design tool and successfully driving the field forward.

For readers who are less interested in the programming and method development aspects of TopOpt as an inverse design tool, we have authored a parallel tutorial paper on TopOpt for photonics applications, utilizing the graphical user interface (GUI)-based commercial finite element software COMSOL Multiphysics as the numerical tool to model the physics and solve the optimization problem [17].

Fig. 1. Model domain, $\Omega$, of height ${h_\Omega}$ and width ${w_\Omega}$ with a designable region, ${\Omega _D}$, of height ${h_{{\Omega _D}}}$ and width ${w_{{\Omega _D}}}$ on top of a substrate, ${\Omega _S}$, of height ${h_s}$.

Download Full Size | PDF

2. PHYSICS AND THE DISCRETIZED OPTIMIZATION PROBLEM

We model the physics in the rectangular domain $\Omega$, with the boundary $\Gamma$ (see Fig. 1), using Maxwell’s equations, assuming time-harmonic temporal behavior. We define a subset of $\Omega$ as the design domain and denote this region ${\Omega _D}$. We assume out-of-plane ($z$ direction) material invariance and that all involved materials are linear, static, homogeneous, isotropic, non-dispersive, non-magnetic, and without inherent polarization. Finally, we assume out-of-plane polarization of the electric field (TE polarization). From these assumptions, we derive a 2D Helmholtz-type partial differential equation for the out-of-plane component of the electric field in $\Omega$,

(1)$$\nabla \cdot \left({\nabla {E_z}({\textbf r})} \right) + {k^2}{\varepsilon _r}({\textbf r}){E_z}({\textbf r}) = F({\textbf r}),\quad {\textbf r} \in \Omega \in {{\mathbb R}^2},$$

where ${E_z}$ denotes the $ z $ component of the electric field, $k = \frac{{2\pi}}{\lambda}$ is the wavenumber with $\lambda$ (= lambda in the top200EM interface) being the wavelength, ${\varepsilon _r}$ denotes the relative electric permittivity, and $F$ denotes a forcing term used to introduce an incident plane wave from the bottom boundary of $\Omega$. We apply first-order absorbing boundary conditions on all four exterior boundaries,

(2)$${\textbf n} \cdot \nabla {E_z}({\textbf r}) = - {\rm i}k{E_z}({\textbf r}),\quad {\textbf r} \in \Gamma ,$$

where ${\textbf n}$ denotes the surface normal and $ i $ the imaginary unit. Note that first-order boundary conditions are not as accurate as certain other boundary conditions, e.g., perfectly matched layers [18]; however, they are conceptually simpler and simpler to implement. Next, we introduce a design field $\xi ({\textbf r}) \in [0,1]$ to control the material distribution in $\Omega$ through the interpolation function,

(3)$${\varepsilon _r}(\xi ({\textbf r})) = 1 + \xi ({\textbf r})\left({{\varepsilon _{r,m}} - 1} \right) - {\rm i}\alpha \xi ({\textbf r})\left({1 - \xi ({\textbf r})} \right),\quad {\textbf r} \in \Omega ,$$

where it is assumed that the background has the value ${\varepsilon _r} = 1$ (e.g., air) and where ${\varepsilon _{r,m}}$ (= eps_r) denotes the relative permittivity of the material used for the structure under design and $\alpha$ is a problem dependent scaling factor. The non-physical imaginary term discourages intermediate values of $\xi$ in the design for the focusing problem at hand, by introducing attenuation [19].

As the baseline example, we consider the design of a monochromatic focusing metalens situated in ${\Omega _D}$. To this end, we select the magnitude of $|{E_z}{|^2}$ at a point in space ${{\textbf r}_p}$ (= targetXY) as the figure of merit (FOM) denoted $\Phi$, i.e.,

(4)$$\Phi (\xi ({\textbf r}),{{\textbf r}_p}) = |{E_z}(\xi ({\textbf r}),{{\textbf r}_p}{)|^2} = {E_z}{(\xi ({\textbf r}),{{\textbf r}_p})^*}{E_z}(\xi ({\textbf r}),{{\textbf r}_p}),$$

where ${\bullet ^*}$ denotes the complex conjugate.

The model equation, boundary conditions, material interpolation function, and FOM are all discretized using the finite element method (FEM) [20] using ${{\cal N}_e}(= {\texttt{nElX}} \cdot \texttt{nElY})$ bi-linear quadratic elements. The discretized model uses nodal degrees of freedom (DOFs) for ${E_z}$ and $F$ as well as elementwise constant DOFs for $\xi$, with ${{\cal N}_D}$ of the elements situated in ${\Omega _D}$. The following constrained continuous optimization problem is formulated for the discretized problem:

(5)$$\begin{split}&\mathop {\max}\limits_\xi :\;\Phi = {\textbf E}_z^\dagger {\textbf P}{{\textbf E}_z}, \\[-4pt] &{\rm s.t.}:\;{\textbf S}({\varepsilon _r}){{\textbf E}_z} = \left({\sum\limits_{e = 1}^{{{\cal N}_e}} {{\textbf S}_e}({\varepsilon _{r,e}})} \right){{\textbf E}_z} = {\textbf F}, \\[-4pt] &\qquad :\;{\varepsilon _{r,j}} = 1 + {{\bar {\tilde \xi}} _j}\left({{\varepsilon _{r,m}} - 1} \right) - {\rm i}{{\bar {\tilde \xi}} _j}\left({1 - {{{\bar {\tilde \xi}}}_j}} \right)\quad \forall \;j \in \{1,2, \ldots,{{\cal N}_e}\} , \\[-4pt] &\qquad :\;0 \lt {\xi _j} \lt 1\quad \forall \;j \in \{1,2, \ldots,{{\cal N}_D}\} , \\[-4pt] &\qquad :\;\xi = 0\quad \forall \;{\textbf r} \in \Omega /\{{\Omega _D},{\Omega _S}\} \vee \xi = 1\quad \forall \;{\textbf r} \in {\Omega _S},\\[-1.3pc]\end{split}$$

where ${{\textbf E}_z}$ and ${\textbf F}$ are vectors containing the nodal DOFs for the electric field and forcing term, and $\xi$ (= dVs) and ${\bar {\tilde \xi}}$ are vectors of element DOFs for the design field and the physical filtered and thresholded field [Eqs. (6) and (7)], respectively.

The diagonal selection matrix ${\textbf P}$ weights the ${{\textbf E}_z}$-DOFs that enter $\Phi$ (= FOM). In the baseline example, it has four $\frac{1}{4}$ entries for the selected element’s nodal DOFs, i.e., it selects the field intensity at the focal point, which for simplicity is taken to be at the center of a single finite element. Finally, ${\bullet ^\dagger}$ denotes the conjugate transpose.

To ameliorate numerical issues, such as pixel-by-pixel design variations, and to introduce a weak sense of geometric length scale, a filter and threshold scheme is applied to $\xi$, before using it to interpolate the material parameters [21–23]. First, the following convolution-based filter operation is applied:

(6)$$\begin{split}{\tilde \xi _h} &= \frac{{\mathop {\sum}\limits_{k \in {{\cal B}_{e,h}}} w({{\textbf r}_h} - {{\textbf r}_k}){A_k}{\xi _k}}}{{\mathop {\sum}\limits_{k \in {{\cal B}_{e,h}}} w({{\textbf r}_h} - {{\textbf r}_k}){A_k}}},\\ w({\textbf r}) &= \left\{{\begin{array}{*{20}{l}}{{r_f} - |{\textbf r}|\quad \forall |{\textbf r}| \le {r_f}}\\0\end{array}} \right.,\quad {r_f} \ge 0,\quad {\textbf r} \in \Omega.\end{split}$$

Here ${A_k}$ denotes the area of the $k$th element, ${r_f}$ (= fR) denotes the desired spatial filtering radius, and finally ${{\cal B}_{e,h}}$ denotes the $h$th set of finite elements whose center point is within ${r_f}$ of the $h$th element. Second, a smoothed approximation of the Heaviside function is applied to the filtered design variables as

(7)$$\begin{split}&{{\bar {\tilde \xi}} _h} = H({\tilde \xi _h}) = \frac{{\tanh (\beta \cdot \eta) + \tanh (\beta \cdot ({{\tilde \xi}_h} - \eta))}}{{\tanh (\beta \cdot \eta) + \tanh (\beta \cdot (1 - \eta))}},\\& \beta \in [1,\infty [,\eta \in [0,1].\end{split}$$

Here $\beta$ and $\eta$ control the threshold sharpness and value, respectively.

Adjoint sensitivity analysis [6,24] is carried out to compute the design sensitivities, utilizing the chain rule for the filter and threshold steps as [23,25]

(8)$$\frac{{{d}\Phi}}{{{d}{\xi _h}}} = \mathop {\sum}\limits_{k \in {{\cal B}_{e,h}}} \frac{{\partial {{\tilde \xi}_k}}}{{\partial {\xi _h}}}\frac{{\partial {{{\bar {\tilde \xi}}}_k}}}{{\partial {{\tilde \xi}_k}}}\frac{{{d}\Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}},\quad \frac{{{d}\Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}} = 2\Re \left({{\lambda ^{T}}\frac{{\partial {\textbf S}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}{{\textbf E}_z}} \right),$$

where $\Re$ denotes the real part, ${\bullet ^T}$ the transpose, and $\lambda$ a vector obtained by solving

(9)$${{\textbf S}^{T}}\lambda = - \frac{1}{2}{\left({\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} - {\rm i}\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}}} \right)^{T}}\quad {\rm with}\quad {{\textbf E}_z} = {{\textbf E}_{z,\Re}} + {\rm i}{{\textbf E}_{z,\Im}},$$

where $\Im$ denotes the imaginary part. The $m$th entry of the right-hand side in Eq. (9) is given by

(10)$$\left({\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} - {\rm i}\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}}} \right)_m^{T} = {{\textbf P}_{m,m}}\left({2({{\textbf E}_{z,\Re}}{)_m} - 2{\rm i}{{({{\textbf E}_{z,\Im}})}_m}} \right).$$

The derivations of the expression for $\frac{{{d}\Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}}$ in Eq. (8) and for Eq. (9) are given in Appendix A.

The fundamental advantage of using the adjoint approach to compute the sensitivities of the FOM with respect to the design variables is that only one single additional system of equations [namely Eq. (9)] must be solved. Hence, the sensitivity information is obtained at approximately the same computational cost as the one required to compute the field information itself. In fact, for the examples treated here, it is possible to reuse the LU-factorization used to compute the field information [second line of Eq. (5)], making the computational cost associated with computing the sensitivity of the FOM almost ignorable.

3. MATLAB CODE

The design problem stated in Eq. (5) is implemented in top200EM (see the full code in Code 1, Ref. [26]), which has the following interface:

function [dVs,FOM] = …

top200EM(targetXY,dVElmIdx,nElX,nElY, dVini,eps_r,lambda,fR,maxItr);

The function takes the input parameters:

targetXY: two-values 1D-array with the $ x $ and $ y $ position of the finite element containing the focal point.

dVElmIdx: 1D-array of indices of the finite elements, which are designable, i.e., ${\Omega _D}$.

nElX: Number of finite elements in the $ x $ direction.

nElY: Number of finite elements in the $ y $ direction.

dVini: Initial guess for discretized $\xi$-field. Accepts a scalar for all elements or a 1D-array of identical length to dVElmIdx.

eps_r: Relative permittivity for the material constituting the structure under design.

lambda: The targeted wavelength, $\lambda$, measured in number of finite elements.

fR: Filter radius, ${r_f}$, measured in number of finite elements.

maxItr: Maximum number of iterations allowed by fmincon for solving the optimization problem in Eq. (5).

And returns the following output parameters:

dVs: 1D-array of the optimized discretized $\xi$-field in the design domain, ${\Omega _D}$.

FOM: Value of the figure of merit.

During execution, the data related to the physics, the discretization, the filter, and threshold operations are stored in the structures phy, dis, and filThr, respectively.

The spatial scaling, threshold sharpness $\beta$, and threshold level $\eta$ are hard coded in top200EM as

% SETUP OF PHYSICS PARAMETERS

phy.scale = 1e-9; % Scaling finite element side length to nanometers

% SETUP FILTER AND THRESHOLDING PARAMETERS

filThr.beta = 5; % Thresholding sharpness

filThr.eta = 0.5; % Thresholding level

Note: For simplicity, the code uses the unit of nanometers to measure length, and the finite elements are taken to have a side length of 1 nm. This may be changed by changing the scaling parameter phy.scale.

The algorithm used to solve the design problem is MATLAB’s fmincon,

[dVs,~] = fmincon(FOM,dVs(:),[],[],[], [],LBdVs,UBdVs,[],options);

with the design-variable bounds and options set up as

LBdVs = zeros(length(dVs),1); % Lower bound on design variables

UBdVs = ones(length(dVs),1); % Upper bound on design variables

options = optimoptions(‘fmincon’,‘Algorithm’,‘interior-point’,…

‘SpecifyObjectiveGradient’,true,‘HessianApproximation’,‘lbfgs’,…

‘Display’,‘off’,‘MaxIterations’,maxItr,‘MaxFunctionEvaluations’,maxItr);

The FOM and sensitivities provided to fmincon are computed using the inline function:

FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);

The discretized design field is distributed in the model domain with a hard coded background of air in the top 90% of the domain and solid material in the bottom 10% as

% DISTRIBUTE MATERIAL IN MODEL DOMAIN BASED ON DESIGN FIELD

dFP(1:dis.nElY,1:dis.nElX) = 0; % Design field in physics, 0: air

dFP(dis.nElY:–1:ceil(dis.nElY*9/10),1:dis.nElX) = 1; % 1: material

dFP(dis.dVElmIdx(:)) = dVs; % Design variables inserted in design field

Followed by the application of the filter and threshold operations and the material interpolation,

% COMPUTE MATERIAL FIELD FROM DESIGN FIELD

dFPS = DENSITY_FILTER(filThr.filKer, filThr.filSca,dFP,ones(dis.nElY,dis.nElX));

dFPST = THRESHOLD(dFPS, filThr.beta, filThr.eta);

[A,dAdx] = MATERIAL_INTERPOLATION(phy.eps_r,dFPST,1.0); % Material field

The system matrix for the state equation is constructed,

% CONSTRUCT SYSTEM MATRIX

[dis,F] = BOUNDARY_CONDITIONS_RHS(phy.k,dis,phy.scale);

dis.vS = reshape(dis.LEM(:)-phy.k^2*dis. MEM(:)*(A(:).’),16*dis.nElX*dis.nElY,1);

SysMat = sparse([dis.iS(:);dis.iBC(:)],[dis.jS(:);dis.jBC(:)],[dis.vS(:);dis.vBC(:)]);

The state system is solved using LU-factorization,

% SOLVING STATE SYSTEM: SysMat * Ez = F

[L,U,Q1,Q2] = lu(SysMat); % LU–factorization

Ez = Q2 * (U\(L\(Q1 * F))); Ez = full(Ez); % Solving

The FOM is computed,

% FIGURE OF MERIT

P = sparse(dis.edofMat(dis.tElmIdx,:), dis.edofMat(dis.tElmIdx,:),1/4,…

(dis.nElX+1)*(dis.nElY+1),(dis.nElX+1) *(dis.nElY+1)); % Weighting matrix

FOM = Ez’ * P * Ez; % Solution in target element

The adjoint system is solved by reusing the LU-factorization from the state problem,

% ADJOINT RIGHT HAND SIDE (0th-order quadrature)

AdjRHS = P*(2*real(Ez)–1i*2*imag(Ez));

% SOLVING THE ADJOING SYSTEM: S.’ * AdjLambda = AdjRHS

AdjLambda = (Q1.’) * ((L.’)\((U.’)\ ((Q2.’) * (-1/2*AdjRHS)))); % Solving

The sensitivities in $\Omega$ are computed and filtered, and the values in ${\Omega _D}$ extracted,

% SENSITIVITIES

dis.vDS = reshape(-phy.k^2*dis.MEM(:) *(dAdx(:).’),16*dis.nElX*dis.nElY,1);

DSdx = sparse(dis.iElFull,dis.jElFull, dis.vDS); % Constructing dS/dx

DSdxMulV = DSdx * Ez(dis.idxDSdx); % Computing dS/dx * Field values

DsdxMulV = sparse(dis.iElSens,dis. jElSens,DSdxMulV);

sens = 2*real(AdjF(dis.idxDSdx).’ * DsdxMulV); % Com puting sensitivities

sens = full(reshape(sens,dis.nElY,dis. nElX));

% FILTERING SENSITIVITIES

DdFSTDFS = DERIVATIVE_OF_THRESHOLD(dFPS, filThr.beta, filThr.eta);

sensFOM = DENSITY_FILTER(filThr.filKer, filThr.filSca,sens,DdFSTDFS);

% EXTRACTING SENSITIVITIES FOR DESIGN REGION

sensFOM = sensFOM(dis.dVElmIdx);

Finally, ${E_z}({\textbf r})$ and $\xi ({\textbf r})$ are plotted, and the current FOM-value printed,

% PLOTTING AND PRINTING

figure(1); % Field intensity, |Ez|^2

imagesc((reshape(Ez.*conj(Ez),dis.nElY +1,dis.nElX+1))); colorbar; axis equal;

figure(2); % Physical design field

imagesc(1-dFPST); colormap(gray); axis equal; drawnow;

disp([’FOM:’ num2str(-FOM)]); % Display FOM value

After the TopOpt procedure is finished, a thresholded version of the final design is evaluated, and the resulting $|{E_z}{|^2}$ field and design field are plotted.

% FINAL BINARIZED DESIGN EVALUATION

filThr.beta = 1000;

disp(’Black/white design evaluation:’)

[obj_2,dFPST_2,F_2] = OBJECTIVE_ GRAD(DVini(:),dis,phy,filThr);

Fig. 2. (a) Max-normalized $|{\textbf E}{|^2}$-field in $\Omega$. (b) Metalens design, ${\varepsilon _r} = 3.0$ (black) and ${\varepsilon _r} = 1.0$ (white).

Download Full Size | PDF

Fig. 3. Metalens designs obtained (a) without filtering (${r_f} = 1$), (b) using a filter radius of ${r_f} = 3.0$, (c) using a filter radius of ${r_f} = 6.0$, (d) using a filter radius of ${r_f} = 9.0$. These results illustrate the effect of applying the cone-shaped filter to the design field as part of the optimization process.

Download Full Size | PDF

4. USING THE CODE

Next, we demonstrate how to use top200EM by designing a focusing metalens as follows.

A. Designing a Metalens

First, we define the domain size in terms of the number of finite elements in each spatial direction and the element indices for the design domain.

% DESIGN FIELD INDICES

DomainElementsX = 400;

DomainElementsY = 200;

DesignThicknessElements = 15;

DDIdx = repmat([1:DomainElementsY: DomainElementsX*DomainElementsY],…

DesignThicknessElements,1);

DDIdx = DDIdx+repmat([165:165 +DesignThicknessElements-1]’,1, DomainElementsX);

Second, the optimization problem is solved by executing the command

[DVs,obj]=top200EM([80,200],DDIdx, DomainElementsX,DomainElementsY,…

0.5,3.0,35,6.0,200);

The final binarized design is shown in Fig. 2(b) with black (white) representing solid (air). Figure 2(a) shows the $|{E_z}{|^2}$-field resulting from exciting the metalens in Fig. 2(b) for the targeted incident field, demonstrating the focusing effect at the targeted focal spot. The numerical aperture of the metalens is ${\rm NA} \approx 0.92$, and the transmission efficiency is ${T_A} \approx 0.87$ computed as the power propagating through the lens relative to the power incident on the lens.

B. Solving the Same Design Problem with Different Resolutions

For various reasons, such as performing a mesh-convergence study, it may be necessary to solve the same physical model problem using different mesh resolutions. This may be done with top200EM by multiplying the following inputs by an integer scaling factor: the number of finite elements in each spatial direction, nElX and nElY, the wavelength lambda, and the filter radius fR, and dividing the hard coded value of phys.scale in the code by the same factor.

C. Effect of Filtering

Next, we demonstrate the effect of applying the filtering step [27] by changing the filter radius and designing four metalenses (see Fig. 3) using top200EM. Again, the model domain size and indices for the design domain are defined first,

% DESIGN FIELD INDICES

DomainElementsX = 400;

DomainElementsY = 200;

DesignThicknessElements = 15;

DDIdx = repmat([1:DomainElementsY: DomainElementsX*DomainElementsY],…

DesignThicknessElements,1);

DDIdx = DDIdx+repmat([165:165 +DesignThicknessElements-1]’,1, DomainElementsX);

Then, the optimization problem is solved with the four filtering radii fR$= {r_f} \in \{1.0,3.0,6.0,9.0\}$,

[DVs,obj]=top200EM([80,200],DDIdx,DomainElementsX,DomainElementsY,…

0.5,3.0,35,fR,200);

Looking at the four final binarized designs in Fig. 3, it is clearly observed that as ${r_f}$ is increased, the features in the designs grow. When no filtering is applied [Fig. 3(a)], single-pixel-sized features are observed. Such features may be problematic from a numerical modeling point of view as well as being detrimental to fabrication.

It is noted that without filtering, TopOpt is more prone to identify a local minimum, which performs worse than the ones identified with filtering. In other words, filtering tends to have a convexifying effect, as long as the filter size is not too big.

5. MODIFYING THE CODE

There exists a vast amount of auxiliary tools developed to extend the applicability of density-based TopOpt across a wide range of different problems and physics. The following sections provide examples of how simple some of these tools are to implement in top200EM.

A. Plasmonics

In recent work [28], it was demonstrated that using a refractive index and extinction cross-section-based non-linear material interpolation yielded significantly improved results when designing Au, Ag, and Cu nano-particles for localized field enhancement, with recent applications to enhanced upconversion [29], thermal emission [30], and Raman scattering [31]. This interpolation function avoids artificial resonances in connection with the transition from positive to negative $\varepsilon$ values and reads

\begin{split}\varepsilon (x) &= (n{(x)^2} - \kappa {(x)^2}) - {\rm i}(2n(x)\kappa (x)),\\n(x)& = {n_{{M}_1}} + x({n_{{M}_2}} - {n_{{M}_1}}),\\ \kappa (x) &= {\kappa _{{M}_1}} + x({\kappa _{{M}_2}} - {\kappa _{{M}_1}}).\end{split}

Here $n$ and $\kappa$ denote the refractive index and extinction cross section, respectively. The subscripts ${M_1}$ and ${M_2}$ denote the two materials being interpolated.

The interpolation scheme is straightforward to implement in the code as follows:

First, the following lines of code:

function [A,dAdx] = MATERIAL_INTERPOLATION(eps_r,x,1.0)

A = 1 + x*(eps_r-1)–1i * alpha_i * x.* (1–x); % Interpolation

dAdx = (eps_r-1)*(1+0*x)–1i * alpha_i * (1–2*x); % Derivative of interpolation end

are replaced with

function [A,dAdx] = MATERIAL_INTERPOLATION(n_r,k_r,x)

n_eff = 1 + x*(n_r-1);

k_eff = 0 + x*(k_r-0);

A = (n_eff.^2–k_eff.^2)–1i*(2.*n_eff.*k_eff);

dAdx = 2*n_eff*(n_r-1)-2*k_eff*(k_r-1)-1i*(2*(n_r-1)*k_eff+(2*n_eff*(k_r-1))); end

where for simplicity it is assumed that ${M_1}$ is air, i.e., ${n_{{M}_1}} = 1.0$ and ${\kappa _{{M}_1}} = 0.0$.

Second, the scalar input parameter eps_r is changed to a two-valued 1D-array nk_r.

Third, the line

phy.eps_r=eps_r; % Relative permittivity

is replaced with,

phy.nk_r = nk_r; % Refractive index and extinction coefficient

and fourth, the call to the material interpolation function is changed from

[A,dAdx] = MATERIAL_INTERPOLATION(phy.eps_r,dFPST,1.0);

to,

[A,dAdx] = MATERIAL_INTERPOLATION(phy.nk_r(1),phy.nk_r(2),dFPST);

B. Excitation

The excitation considered in the baseline problem may be changed straightforwardly. As an example, the boundary at which the incident field enters the domain may be changed from the bottom to the top boundary as follows.

First, the index set controlling where the boundary condition is imposed in the vector, ${\textbf F}$ (= F), is changed by moving dis.iRHS = TMP; from line 157 to above line 151.

Second, the values stored in F are changed to account for the propagation direction of the wave by replacing

F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))+1i*waveVector;

F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))+1i*waveVector;

with

F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))-1i*waveVector;

F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))-1i*waveVector;

C. Designing a Metallic Reflector

By introducing the changes presented in Sections 5.A and 5.B, one may design a metallic reflector using topEM200.

First, we set the domain dimensions and design thickness as in the previous examples,

% DESIGN FIELD INDICES

DomainElementsX = 400;

DomainElementsY = 200;

DesignThicknessElements = 15;

DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],…

DesignThicknessElements,1);

DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);

Second, we solve the optimization problem by executing the command,

[DVs,obj]=top200EM([100,200],DDIdx,DomainElementsX,DomainElementsY,…

0.5,[1.9,1.5],35,3.0,200);

Note: For simplicity we selected the following values for the refractive index, $n = 1.9$ (=nk_r(1)), and extinction cross section, $\kappa = 1.5$ (=nk_r(2)) (These values corresponds to values for gold at $\lambda = 350\;{\rm nm}$, i.e., one may think of this choice of material parameters as an unspoken rescaling of space by a factor of 10, i.e., changing the element size to pixels of 10 nm by 10 nm rather than of 1 nm by 1 nm).

Fig. 4. (a) Max-normalized $|{\textbf E}{|^2}$-field in $\Omega$. (b) Metallic reflector design (black) in air background (white).

Download Full Size | PDF

Considering the final binarized reflector design in Fig. 4(b), one can interpret the design as the well-known parabolic reflector broken into pieces to fit the spatially limited design domain. From the max-normalized $|{\textbf E}{|^2}$-field presented in Fig. 4(a), the focusing effect of the reflector is clearly observed.

D. Linking Design Variables

Certain fabrication techniques limit the allowable geometric variations in a design. For example, optical-projection lithography and electron-beam lithography restrict variations in design geometries to 2D patterns, which can then be extruded in the out-of-plane direction. It is straightforward to introduce such a geometric restriction using TopOpt by linking design variables and sensitivities across elements. In top200EM the design field may be restricted to only exhibit in-plane ($ x $ direction) variations as follows.

First, the code-line representing the values of the design variables,

dVs(length(dis.dVElmIdx(:))) = dVini; % Design variables

is replaced with,

dVs(1:nElX) = dVini; % Design variables for 1D design

Second, the code-line transferring the design variables to the elements in the physics model,

dFP(dis.dVElmIdx(:)) = dVs; % Design variables inserted in design field

is replaced with,

nRows=length(dis.dVElmIdx(:))/dis.nElX; % Number of rows in the 1D design

dFP(dis.dVElmIdx) = repmat(dVs,1,nRows)’; % Design variables inserted in design field

and finally the code-line representing individual element sensitivities,

sensFOM = sensFOM(dis.dVElmIdx);

is replaced with,

sensFOM = sensFOM(dis.dVElmIdx);

sensFOM = sum(sensFOM,1); % Correcting sensitivities for 1D design,

which sums the sensitivity contributions from the linked elements.

1. Designing a 1D Metalens

By introducing the changes listed in Section 5.D, topEM200 can be used to design metasurfaces of fixed height with in-plane variations as follows. Again we define the domain dimensions as

Fig. 5. (a) Max-normalized $|{\textbf E}{|^2}$-field in $\Omega$. (b) Metalens design restricted to one-dimensional variations, ${\varepsilon _r} = 3.0$ (black) and ${\varepsilon _r} = 1.0$ (white).

Download Full Size | PDF

% DESIGN FIELD INDICES

DomainElementsX = 400;

DomainElementsY = 200;

DesignThicknessElements = 15;

DDIdx = repmat[1:DomainElementsY:DomainElementsX*DomainElementsY],…

DesignThicknessElements,1);

DDIdx = DDIdx+repmat([165:165+DesignThicknessElements-1]’,1,DomainElementsX);

Followed by the execution of the command,

[DVs,obj]=top200EM([80,200],DDIdx,DomainElementsX,DomainElementsY,…

0.5,3.0,35,3.0,200);

The final binarized design, resulting from solving the design problem, is shown in Fig. 5(b). Here it is clear to see that the design is now restricted to vary only in the $ x $ direction. It is worth noting that the design is still filtered in both spatial directions; hence, the corners of the design appear rounded in Fig. 5(b). The max-normalized $|{\textbf E}{|^2}$-field presented in Fig. 4(a) demonstrates the focusing effect of the lens at the targeted point in the modeling domain, which, due to the reduced design freedom, is not as high as in the unrestricted case (Section 4.A).

E. Continuation of Threshold Sharpness

For some design problems in photonics and plasmonics, e.g., [22,28,31], intermediate values may be present in the final optimized design, i.e., ${{\bar {\tilde \xi}} _{{\rm final}}} \in]0,1[$, despite the applied penalization scheme, as they prove beneficial to optimizing the FOM. However, (in most cases) intermediate values hold no physical meaning, and it is therefore not possible to realize designs containing such intermediate values experimentally. A way to promote that (almost) no design variables take intermediate values in the final design is by using a continuation scheme for the threshold sharpness, gradually increasing it until an (almost) pure 0/1-design is achieved, ensuring that the optimized designs are physically realizable.

To implement the continuation scheme in the code, replace the following lines:

% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON

FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);

$[{\rm dVs},\sim]$ = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);

with,

while filThr.beta<betaMax % Thresholding sharpness bound

% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON

FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);

$[{\rm dVs},\sim]$ = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);

filThr.beta = betaInc * filThr.beta; % Increasing thresholding sharpness

end

and select suitable values for betaMax and betaInc. These values are problem dependent, and some experimentation may be required to identify the best values for a given problem. For the metalens design example in Section 4.A, when considering high material contrast, i.e., large values of eps_r, the values betaMax = 20.0 and betaInc = 1.5 have been found to work well.

6. WHAT ABOUT GENETIC ALGORITHMS?

For a wide range of inverse design problems, where gradient-based optimization methods are applicable [inverse design problems where it is possible to compute the sensitivities of the FOM (and any constraints), e.g., using adjoint sensitivity analysis, a pre-requisite for using gradient-based methods], they have been found to severely outperform non-gradient-based methods, such as GAs [32]. This, both in terms of the computational effort required to identify a local optimum for the FOM and, by direct extension, in terms of the number of design DOFs, it is feasible to consider (a difference of many orders of magnitude [33]). Further, in many cases, gradient-based methods are able to identify better local optima for the FOM [34].

In the following, we provide a demonstration of these claims, by comparing the solution of a metalens design problem obtained using the gradient-based top200EM to the solution obtained using MATLAB’s built-in genetic algorithm ga.

Readers may perform their own comparisons by rewriting top200EM to utilize ga instead of fmincon as follows.

Fig. 6. (a) Design obtained using the GA-based method. (b) Design obtained using the gradient-based method. ${\varepsilon _r} = 3.0$ (black) and ${\varepsilon _r} = 1.0$ (white). (c) Convergence graph. Remark that each design iteration for the GA-based code requires 200 solutions of the physics model equation, while it requires one solution of the physics model equation and one solution of the adjoint equation for the gradient-based code!

Download Full Size | PDF

First, a (near perfect) 0/1-design is ensured by changing the thresholding strength from

filThr.beta = 5; % Thresholding sharpness

to,

filThr.beta = 1e5; % Thresholding sharpness

Second, ga is used instead of fmincon by replacing the following lines of codes:

options = optimoptions(‘fmincon’,‘Algorithm’,‘interior-point’,…

‘SpecifyObjectiveGradient’,true,‘HessianApproximation’,‘lbfgs’,…

‘Display’,‘off’,‘MaxIterations’,maxItr,‘MaxFunctionEvaluations’,maxItr);

% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: FMINCON

FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr);

$[{\rm dVs},\sim]$ = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);

with,

options = optimoptions(’ga’,‘MaxGenerations’,maxItr,‘Display’,‘iter’);

% SOLVE DESIGN PROBLEM USING MATLAB BUILT-IN OPTIMIZER: GA

FOM = @(dVs)OBJECTIVE(dVs,dis,phy,filThr);

rng(1,‘twister’); % SETTING RANDOM SEED TO ENSURE REPRODUCTION OF RESULT

$[{\rm dVs},\sim,\sim,\sim]$ = ga(FOM,length(dVs),[],[],[],[],LBdVs,UBdVs,[],[],options);

where the command rng(1,‘twister’) fixes the random seed for reproducibility.

Third, the computation of the gradient of the FOM is removed by changing

function [FOM,sensFOM] = OBJECTIVE_GRAD(dVs,dis,phy,filThr)

to,

function [FOM] = OBJECTIVE(dVs,dis,phy,filThr)

Fourth, deleting the following lines in OBJECTIVE:

% ADJOINT RIGHT HAND SIDE

AdjRHS = P*(2*real(Ez)–1i*2*imag(Ez));

% SOLVING THE ADJOING SYSTEM: S.’ * AdjLambda = AdjRHS

AdjLambda = (Q1.’) * ((L.’)\((U.’)\((Q2.’) * (-1/2*AdjRHS)))); % Solving

% COMPUTING SENSITIVITIES

dis.vDS = reshape(-phy.k^2*dis.MEM(:)*(dAdx(:).’),16*dis.nElX*dis.nElY,1);

DSdx = sparse(dis.iElFull,dis.jElFull,dis.vDS); % Constructing dS/dx

DSdxMulV = DSdx * Ez(dis.idxDSdx); % Computing dS/dx * Field values

DsdxMulV = sparse(dis.iElSens,dis.jElSens,DSdxMulV);

sens = 2*real(AdjLambda(dis.idxDSdx).’ * DsdxMulV); % Computing sensitivites

sens = full(reshape(sens,dis.nElY,dis.nElX));

% FILTERING SENSITIVITIES

DdFSTDFS = DERIVATIVE_OF_THRESHOLD(dFPS, filThr.beta, filThr.eta);

sensFOM = DENSITY_FILTER(filThr.filKer,filThr.filSca,sens,DdFSTDFS);

% EXTRACTING SENSITIVITIES FOR DESIGNABLE REGION

sensFOM = sensFOM(dis.dVElmIdx);

% FMINCON DOES MINIMIZATION

sensFOM = -sensFOM(:);

and finally changing the code evaluating the final design,

% FINAL BINARIZED DESIGN EVALUATION

filThr.beta = 1000;

disp(’Black/white design evaluation:’);

$[{\rm FOM},\sim]$ = OBJECTIVE(dVs(:),dis,phy,filThr);

to,

% FINAL BINARIZED DESIGN EVALUATION

filThr.beta = 1e5;

disp(’Black/white design evaluation:’);

$[{\rm FOM},\sim]$ = OBJECTIVE(dVs(:),dis,phy,filThr);

Note: The code-lines plotting the design in OBJECTIVE can be removed to reduce wall-clock time.

In order to reduce the computational effort required to reproduce the following example, we consider a problem in a smaller spatial domain, considering a shorter wavelength and fewer DOFs than in the previous examples. This is done by defining a model problem considering 1000 design DOFs as

% DESIGN FIELD INDICES

DomainElementsX = 100;

DomainElementsY = 50;

DesignThicknessElements = 10;

DDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],…

DesignThicknessElements,1);

DDIdx = DDIdx+repmat([35:35+DesignThicknessElements-1]’,1,DomainElementsX);

Followed by the solution of the design problem using the modified GA-based code,

[dVs_GA,FOM_GA]=top200EMGA([10,50],DDIdx,DomainElementsX,DomainElementsY,…

0.5,3.0,20,3.0,500);

and then using the original gradient-based code,

[dVs,FOM]=top200EM([10,50],DDIdx,DomainElementsX,DomainElementsY,…

0.5,3.0,20,3.0,500);

The binarized designs obtained using top200EMGA and top200EM are presented in Figs. 6(a) and 6(b), respectively. The evolution of the FOM value as function of the number of FOM evaluations is plotted in Fig. 6(c) for both methods.

The GA-based method uses 137 design iterations (27,600 FOM evaluations) to identify a solution. In contrast, the gradient-based method uses 208 iterations (208 FOM evaluation) to identify a superior solution. Evaluating the final binarized designs, the GA-based solution has a FOM value of $\Phi \approx 6.85$ while the gradient-based solution has an $\approx 18\%$ better FOM value of $\Phi \approx 8.07$.

Crucially, the GA-based method uses 200 FOM evaluations per design iteration to drive the optimization process. This results in a total of 27,600 FOM evaluations for solving the design problem using the GA-based method, which corresponds to 27,600 solutions of the physics model equation, i.e., Eq. (1). In contrast, when solving the design problem using the gradient-based method, a total of 208 FOM evaluations are performed, corresponding to 208 solutions of the physics model equation plus 208 solutions of the adjoint system [Eq. (9)], which for the baseline example is almost free due to the reuse of the LU-factorization.

Hence, in our example, the identification of a local optimum for the FOM using the gradient-based method required $\approx 1.0\%$ of the computational effort spent by the GA-based method, a ratio that only becomes more pronounced as the size of the design space grows.

This example is based on standard settings for MATLAB’s ga optimizer and can surely be improved and refined. However, this will not change the basic conclusion that non-gradient-based approaches are unsuitable for large-scale TopOpt problems. A more thorough discussion of non-gradient-based methods in TopOpt can be found in Ref. [34]. Readers are invited to form their own conclusions by extending the example and codes provided here.

7. CONCLUSION

We have presented a simple finite element based MATLAB code (downloadable from https://www.topopt.mek.dtu.dk) for TopOpt-based inverse design of photonic structures. We have provided examples of how to use the code, as well as a set of suggestions for code extensions enabling it to handle metallic structures, different model excitations, and linked design variables, and introducing a continuation scheme for the threshold sharpness designed to promote 0/1-designs. Finally, we have demonstrated the superiority of a gradient-based method over a genetic-algorithm-based method when performing inverse design in photonics.

The code can be used for educational purposes as is, and is otherwise meant to serve as a starting point for the reader to develop software to handle their more advanced research applications within photonics. For simplicity and computational speed, the code treats problems in two spatial dimensions; however, it is directly extendable to three spatial dimensions by modifying the finite element matrices, boundary conditions, and index sets appropriately.

APPENDIX A: ADJOINT SENSITIVITY ANALYSIS

The expression for $\frac{{{d}\Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}}$ in Eq. (8) and the expression in Eq. (9) may be derived as follows. First, zero is added to $\Phi$ twice,

(A1)$$\tilde \Phi = \Phi + {\lambda ^{T}}\left({{\textbf S}{{\textbf E}_z} - {\textbf F}} \right) + {\lambda ^\dagger}\left({{{\textbf S}^*}{\textbf E}_z^* - {{\textbf F}^*}} \right),$$

where $({{\textbf S}{{\textbf E}_z} - {\textbf F}}) = 0$ and $\lambda$ is a vector of nodal complex Lagrange multipliers, also called the adjoint variables. Second, one takes the derivative of $\tilde \Phi$ with respect to ${{\bar {\tilde \xi}} _k}$ and exploits that, for the optimization problem in Eq. (5), $\Phi$ does not depend explicitly on ${{\bar {\tilde \xi}} _k}$ and neither $\lambda$ nor ${\textbf F}$ depend on ${{\bar {\tilde \xi}} _k}$ at all, yielding

(A2)$$\begin{split}\frac{{{d}\tilde \Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}} &= \frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}}\frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}} + \frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}}\frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}} \\&\quad+ {\lambda ^{T}}\left({\frac{{\partial {\textbf S}}}{{{d}{{{\bar {\tilde \xi}}}_k}}}{{\textbf E}_z} + {\textbf S}\left({\frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}} + {\rm i}\frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}} \right)} \right) \\ &\quad {+}{\lambda ^\dagger}\left({\frac{{\partial {{\textbf S}^*}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}{\textbf E}_z^* + {{\textbf S}^*}\left({\frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}} - {\rm i}\frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}} \right)} \right).\end{split}$$

Collecting terms including $\frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}$ and $\frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}$ and reducing the remaining terms yields

(A3)$$\begin{split}\frac{{{d}\tilde \Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}} &= \frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}\left({\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} + {\lambda ^{T}}{\textbf S} + {\lambda ^\dagger}{{\textbf S}^*}} \right)\\&\quad + \frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}\left({\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}} + {\rm i}{\lambda ^{T}}{\textbf S} + {\rm i}{\lambda ^\dagger}{{\textbf S}^*}} \right) + 2\Re \left({{\lambda ^{T}}\frac{{\partial {\textbf S}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}{{\textbf E}_z}} \right).\end{split}$$

To eliminate the first two terms in Eq. (A3) containing the derivates $\frac{{\partial {{\textbf E}_{z,\Re}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}$ and $\frac{{\partial {{\textbf E}_{z,\Im}}}}{{\partial {{{\bar {\tilde \xi}}}_k}}}$, the two parenthesis must equal zero,

(A4)$$\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} + {\lambda ^{T}}{\textbf S} + {\lambda ^\dagger}{{\textbf S}^*} = 0,\quad \frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}} + {\rm i}{\lambda ^{T}}{\textbf S} + {\rm i}{\lambda ^\dagger}{{\textbf S}^*} = 0,$$

multiplying the second equation by i, subtracting it from the first, and transposing it yields

(A5)$$\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} - {\rm i}\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}} + 2{\lambda ^{T}}{\textbf S} = 0\Leftrightarrow{{\textbf S}^{T}}\lambda = - \frac{1}{2}{\left({\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Re}}}} - {\rm i}\frac{{\partial \Phi}}{{\partial {{\textbf E}_{z,\Im}}}}} \right)^{T}}.$$

Using Eq. (A5), the expression in Eq. (A3) reduces to the expression for $\frac{{{d}\Phi}}{{{d}{{{\bar {\tilde \xi}}}_k}}}$ in Eq. (8), and the second equation in Eq. (A5) is equal to Eq. (9).

Funding

Danmarks Grundforskningsfond (DNRF147); Villum Fonden (8692).

Disclosures

The authors declare that there are no conflicts of interest related to this paper.

REFERENCES

1. M. P. Bendsøe and O. Sigmund, Topology Optimization (Springer, 2003).

2. M. P. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Comput. Methods Appl. Mech. Eng. 71, 197–224 (1988). [CrossRef]

3. J. Alexandersen and C. S. Andreasen, “A review of topology optimisation for fluid-based problems,” Fluids 5, 29 (2020). [CrossRef]

4. C. B. Dilgen, S. B. Dilgen, N. Aage, and J. S. Jensen, “Topology optimization of acoustic mechanical interaction problems: a comparative review,” Struct. Multidiscip. Optim. 60, 779–801 (2019). [CrossRef]

5. C. Lundgaard, K. Engelbrecht, and O. Sigmund, “A density-based topology optimization methodology for thermal energy storage systems,” Struct. Multidiscip. Optim. 60, 2189–2204 (2019). [CrossRef]

6. J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photon. Rev. 5, 308–321 (2011). [CrossRef]

7. S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics 12, 659–670 (2018). [CrossRef]

8. X. Liang and S. G. Johnson, “Formulation for scalable optimization of microcavities via the frequency-averaged local density of states,” Opt. Express 21, 30812–30841 (2013). [CrossRef]

9. F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, “Maximizing the quality factor to mode volume ratio for ultra-small photonic crystal cavities,” Appl. Phys. Lett. 113, 241101 (2018). [CrossRef]

10. A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics 9, 374–377 (2015). [CrossRef]

11. Z. Lin, V. Liu, R. Pestourie, and S. G. Johnson, “Topology optimization of freeform large-area metasurfaces,” Opt. Express 27, 15765–15775 (2019). [CrossRef]

12. H. Chung and O. D. Miller, “High-NA achromatic metalenses by inverse design,” Opt. Express 28, 6945–6965 (2020). [CrossRef]

13. R. E. Christiansen, F. Wang, O. Sigmund, and S. Stobbe, “Designing photonic topological insulators with quantum-spin-Hall edge states using topology optimization,” Nanophotonics 8, 1363–1369 (2019). [CrossRef]

14. O. Sigmund, “A 99 line topology optimization code written in Matlab,” Struct. Multidiscip. Optim. 21, 120–127 (2001). [CrossRef]

15. E. Andreassen, A. Clausen, M. Schevenels, B. S. Lazarov, and O. Sigmund, “Efficient topology optimization in MATLAB using 88 lines of code,” Struct. Multidiscip. Optim. 43, 1–16 (2011). [CrossRef]

16. F. Ferrari and O. Sigmund, “A new generation 99 line Matlab code for compliance topology optimization and its extension to 3D,” Struct. Multidiscip. Optim. 62, 2211–2228 (2020). [CrossRef]

17. R. E. Christiansen and O. Sigmund, “Inverse design in photonics by topology optimization: tutorial,” J. Opt. Soc. Am. B 38, 496–509 (2021). [CrossRef]

18. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114, 185–200 (1994). [CrossRef]

19. J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: a high-bandwidth low-loss T-junction waveguide,” J. Opt. Soc. Am. B 22, 1191–1198 (2005). [CrossRef]

20. J. Jin, The Finite Element Method in Electromagnetics (Wiley-Interscience, 2013).

21. J. K. Guest, J. H. Prevost, and T. Belytschko, “Achieving minimum length scale in topology optimization using nodal design variables and projection functions,” Int. J. Numer. Methods Eng. 61, 238–254 (2004). [CrossRef]

22. F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidiscip. Optim. 43, 767–784 (2011). [CrossRef]

23. R. E. Christiansen, B. S. Lazarov, J. S. Jensen, and O. Sigmund, “Creating geometrically robust designs for highly sensitive problems using topology optimization—acoustic cavity design,” Struct. Multidiscip. Optim. 52, 737–754 (2015). [CrossRef]

24. D. A. Tortorelli and P. Michaleris, “Design sensitivity analysis: overview and review,” Inverse Prob. Eng. 1, 71–105 (1994). [CrossRef]

25. M. B. Dühring and J. S. Jensen, and O. Sigmund, “Acoustic design by topology optimization,” J. Sound Vib. 317, 557–575 (2008). [CrossRef]

26. R. E. Christiansen and O. Sigmund, “A 200 line topology optimization code for electromagnetism,” Code 1, figshare (2020), https://doi.org/10.6084/m9.figshare.13369355.

27. B. S. Lazarov, F. Wang, and O. Sigmund, “Length scale and manufacturability in density-based topology optimization,” Arch. Appl. Mech. 86, 189–218 (2016). [CrossRef]

28. R. E. Christiansen, J. Vester-Petersen, S. P. Madsen, and O. Sigmund, “A non-linear material interpolation for design of metallic nano-particles using topology optimization,” Comput. Methods Appl. Mech. Eng. 343, 23–39 (2019). [CrossRef]

29. J. Christiansen, J. Vester-Petersen, S. Roesgaard, S. H. Møller, R. E. Christiansen, O. Sigmund, S. P. Madsen, P. Balling, and B. Julsgaard, “Strongly enhanced upconversion in trivalent erbium ions by tailored gold nanostructures: toward high-efficient silicon-based photovoltaics,” Sol. Energy Mater. Sol. Cells 208, 110406 (2020). [CrossRef]

30. Z. A. Kudyshev, A. V. Kildishev, V. M. Shalaev, and A. Boltasseva, “Machine-learning-assisted metasurface design for high-efficiency thermal emitter optimization,” Appl. Phys. Rev. 7, 021407 (2020). [CrossRef]

31. R. E. Christiansen, J. Michon, M. Benzaouia, O. Sigmund, and S. G. Johnson, “Inverse design of nanoparticles for enhanced Raman scattering,” Opt. Express 28, 4444–4462 (2020). [CrossRef]

32. D. E. Goldberg, Genetic Algorithms in Search, Optimization and Learning (Addison, 1989).

33. N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund, “Giga-voxel computational morphogenesis for structural design,” Nature 550, 84–86 (2017). [CrossRef]

34. O. Sigmund, “On the usefulness of non-gradient approaches in topology optimization,” Struct. Multidiscip. Optim. 43, 589–596 (2011). [CrossRef]

Compact 200 line MATLAB code for inverse design in photonics by topology optimization: tutorial

Abstract

Corrections

1. INTRODUCTION

2. PHYSICS AND THE DISCRETIZED OPTIMIZATION PROBLEM

3. MATLAB CODE

4. USING THE CODE

A. Designing a Metalens

B. Solving the Same Design Problem with Different Resolutions

C. Effect of Filtering

5. MODIFYING THE CODE

A. Plasmonics

B. Excitation

C. Designing a Metallic Reflector

D. Linking Design Variables

1. Designing a 1D Metalens

E. Continuation of Threshold Sharpness

6. WHAT ABOUT GENETIC ALGORITHMS?

7. CONCLUSION

APPENDIX A: ADJOINT SENSITIVITY ANALYSIS

Funding

Disclosures

REFERENCES

Supplementary Material (1)

Cited By

Figures (6)

Equations (16)

Journal of the Optical Society of America B