## Abstract

We provide a compact 200 line MATLAB code demonstrating how topology optimization (TopOpt) as an inverse design tool may be used in photonics, targeting the design of two-dimensional dielectric metalenses and a metallic reflector as examples. The physics model is solved using the finite element method, and the code utilizes MATLAB’s *fmincon* algorithm to solve the optimization problem. In addition to presenting the code itself, we briefly discuss a number of extensions and provide the code required to implement some of these. Finally, we demonstrate the superiority of using a gradient-based method compared to a genetic-algorithm-based method (using MATLAB’s *ga* algorithm) for solving inverse design problems in photonics. The MATLAB software is freely available in the paper and may be downloaded from https://www.topopt.mek.dtu.dk.

© 2021 Optical Society of America

## 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].

## 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$,

`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,

`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.,

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:

`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:

`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

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]

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);`

## 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

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).

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

`% 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.

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,

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

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,

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]