I Introduction
This paper details a 200 line MATLAB code, which demonstrates how densitybased Topology Optimization (TopOpt) can be applied to photonics design. The code is written for scientists and students with a basic knowledge of programming, numerical modelling 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 (Sec. II). The MATLAB code is then explained in detail (Sec. III), followed by two application examples providing the reader with targets for reproduction (Sec. IV). Then, a number of possible extensions are discussed and codesnippets for easy implementation are provided along with design examples (Sec. V). Finally, we demonstrate the superiority of using gradientbased TopOpt compared to a genetic algorithm when solving a photonic design problem (Sec. VI).
Topology optimization BOOK_TOPOPT_BENDSOE as an inverse design tool was first developed in the context of solid mechanics in the late 1980s BENDSOE_KIKUCHI_1988. Since its inception the method has developed rapidly and expanded across most areas of physics ALEXANDERSEN_ANDREASEN_2020; DILGEN_ET_AL_2019; LUNDGAARD_2019; JENSEN_SIGMUND_2011. Over the last two decades the interest in applying TopOpt for photonics has increased rapidly MOLESKY_2018 with applications within cavity design LIANG_JOHNSON_2013; WANG_2018, photonic demultiplexers PIGGOTT_2015, metasurfaces ZIN_2019; CHUNG_MILLER_2020 and topological insulators CHRISTIANSEN_NP_2019 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: the required knowledge of numerical modelling, of advanced mathematical concepts and scientific programming experience. This paper seeks to lower these barriers by providing the reader with a simple 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 SIGMUND_2001; ANDREASSEN_2011; FERRARI_2020 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 GUI based commercial finite element software COMSOL Multiphysics as the numerical tool to model the physics and solve the optimization problem CHRISTIANSEN_SIGMUND_COMSOL_2020.
Ii The Physics and the Discretized Optimization Problem
We model the physics in the rectangular domain , with the boundary , (see Fig. 1) using Maxwell’s equations, assuming timeharmonic temporal behaviour. We define a subset of as the design domain and denote this region . We assume outofplane (
direction) material invariance and that all involved materials are linear, static, homogeneous, isotropic, nondispersive, nonmagnetic and without inherent polarization. Finally, we assume outofplane polarization of the electric field (TE polarization). From these assumptions we derive a twodimensional Helmholtztype partial differential equation for the outofplane component of the electric field
in ,(1) 
where denotes the zcomponent of the electric field, is the wavenumber with lambda
in the top200EM
interface being the wavelength, denotes the relative electric permittivity and denotes a forcing term used to introduce an incident plane wave from the bottom boundary of . We apply first order absorbing boundary conditions on all four exterior boundaries,
(2) 
where n denotes the surface normal and the imaginary unit. Note that first order boundary conditions are not as accurate as certain other boundary conditions, e.g. perfectly matched layers BERENGER_1994, however they are conceptually simpler and simpler to implement. Next, we introduce a design field to control the material distribution in
through the interpolation
function,(3) 
where it is assumed that the background has the value (e.g. air) and where eps_r
denotes the relative permittivity of the material used for the structure under design and is a problem dependent scaling factor. The nonphysical imaginary term discourages intermediate values of in the design for the focussing problem at hand, by introducing attenuation JENSEN_SIGMUND_2005.
As the baseline example, we consider the design of a monochromatic focussing metalens situated in . To this end, we select the magnitude of at a point in space targetXY
as the figure of merit (FOM) denoted , i.e.
(4) 
where denotes the complex conjugate.
The model equation, boundary conditions, materialinterpolation function and figure of merit are all discretized using the finite element method (FEM) BOOM_FEM_JIN using nElX
nElY
bilinear quadratic elements. The discretized model uses
nodal degrees of freedom (DOFs) for
and as well as elementwise constant DOFs for , with of the elements situated in . The following constrained continuous optimization problem is formulated for the discretized problem,where and F
are vectors containing the nodal DOFs for the electric field and forcing term and
dVs
and are vectors of element DOFs for the design field and the physical filtered and thresholded field (eqs. (6)(7)), respectively.
The diagonal selection matrix P weights the DOFs that enter FOM
. In the baseline example it has four 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, denotes the conjugate transpose.
To ameliorate numerical issues, such as pixelbypixel design variations, and to introduce a weak sense of geometric length scale, a filter and threshold scheme is applied to , before using it to interpolate the material parameters GUEST_ET_AL_2004; WANG_ET_AL_2011; CHRISTIANSEN_SMO_2015. First, the following convolutionbased filter operation is applied,
(6) 
Here denotes the area of the ’th element and fR
the desired spatial filtering radius, finally denotes the ’th set of finite elements who’s center point is within of the ’th element. Second, a smoothed approximation of the Heaviside function is applied to the filtered design variables as,
(7) 
Here and control the threshold sharpness and value, respectively.
Adjoint sensitivity analysis TORTORELLI_ET_AL_1994; JENSEN_SIGMUND_2011
is carried out to compute the design sensitivities, utilizing the chain rule for the filter and threshold steps as
Duhring_2008; CHRISTIANSEN_SMO_2015,(8) 
where denotes the real part, the transpose and a vector obtained by solving,
(9) 
where denotes the imaginary part. The ’th entry of the righthand side in eq. (9) is given by,
(10) 
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 LUfactorization used to compute the field information (second line of eq. (II)), making the computational cost associated with computing the sensitivity of the FOM almost ignorable.
Iii The MATLAB Code
The design problem stated in eq. (II) is implemented in top200EM
(see the full code in Appendix B), which has the interface:
function [dVs,FOM] = ... top200EM(targetXY,dVElmIdx,nElX,nElY,dVini,eps_r,lambda,fR,maxItr);
The function takes the input parameters:

targetXY
: twovalues 1Darray with the x and yposition of the finite element containing the focal point. 
dVElmIdx
: 1Darray of indices of the finite elements which are designable, i.e. . 
nElX
: Number of finite elements in the xdirection. 
nElY
: Number of finite elements in the ydirection. 
dVini
: Initial guess for discretized field. Accepts a scalar for all elements or a 1Darray of identical length todVElmIdx
. 
eps_r
: Relative permittivity for the material constituting the structure under design. 
lambda
: The targeted wavelength, , measured in number of finite elements. 
fR
: Filter radius, , measured in number of finite elements. 
maxItr
: Maximum number of iterations allowed byfmincon
for solving the optimization problem in eq. (II).
And returns the output parameters:

dVs
: 1Darray of the optimized discretized field in the design domain, . 
FOM
: Value of the figure of merit.
During execution, the data related to the physics, the discretization and the filter and threshold operations are stored in the structures phy
, dis
and filThr
, respectively.
The spatial scaling, threshold sharpness and threshold level are hard coded in top200EM
as,

% SETUP OF PHYSICS PARAMETERS phy.scale = 1e9; % 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 MATLABs fmincon
,

[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options);
with the designvariable 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’,’interiorpoint’,... ’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,... ’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr);
The FOM and sensitivites 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 LUfactorization,

% 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 LUfactorization from the state problem,

% ADJOINT RIGHT HAND SIDE (0thorder 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 are computed and filtered, and the values in 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); % 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 DESIGN REGION sensFOM = sensFOM(dis.dVElmIdx);
Finally, and are plotted and the current FOMvalue 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(1dFPST); 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 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);
Iv Using the Code
Next, we demonstrate how to use top200EM
by designing a focusing metalens as follows.
iv.1 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+DesignThicknessElements1]’,1,DomainElementsX);
Second, the optimization problem is solved by executing the command,

[DVs,obj]=top200EM([200,80],DDIdx,DomainElementsX,DomainElementsY,... 0.5,3.0,35,6.0,200);
The final binarized design is shown in Fig. 2b with black(white) representing solid(air). Figure 2a shows the field resulting from exciting the metalens in Fig. 2b for the targeted incident field, demonstrating the focussing effect at the targeted focal spot. The numerical aperture of the metalens is NA and the transmission efficiency is computed as the power propagating through the lens relative to the power incident on the lens.
iv.2 Solving the same design problem with different resolutions
For various reasons, such as performing a meshconvergence 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 hardcoded value of phys.scale
in the code by the same factor.
iv.3 The Effect of Filtering
Next, we demonstrate the effect of applying the filtering step LAZAROV_ET_AL_2016 by changing the filter radius and designing four metalenses (See Fig. 3) using top200EM
. Again, the modeldomain 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+DesignThicknessElements1]’,1,DomainElementsX);
Then, the optimization problem is solved with the four filtering radii fR
,

[DVs,obj]=top200EM([200,80],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 is increased the features in the designs grow. When no filtering is applied (Fig. 3a) singlepixel sized features are observed. Such features may be problematic from a numericalmodelling 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 ones identified with filtering. In other words, filtering tends to have a convexifying effect, as long as the filter size is not too big.
V Modifying the Code
There exists a vast amount of auxiliary tools developed to extend the applicability of densitybased Topology Optimization 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
.
v.1 Plasmonics
In recent work CHRISTIANSEN_VESTER_2019 it was demonstrated that using a refractive index and extinction cross section based nonlinear material interpolation yielded significantly improved results when designing Au, Ag and Cu nanoparticles for localized field enhancement, with recent applications to enhanced upconversion CHRISTIANSEN_SEMSC_2020. thermal emission KUDYSHEV_ET_AL_2020 and Raman scattering CHRISTIANSEN_OE_2020. This interpolation function avoids artificial resonances in connection with the transition from positive to negative values and reads,
Here and denote the refractive index and extinction cross section, respectively. The subscripts and denote the two materials being interpolated.
The interpolation scheme is straight forward 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_r1)  1i * alpha_i * x .* (1  x); % Interpolation dAdx = (eps_r1)*(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_r1); k_eff = 0 + x*(k_r0); A = (n_eff.^2  k_eff.^2)  1i*(2.*n_eff.*k_eff); dAdx = 2*n_eff*(n_r1)2*k_eff*(k_r1)1i*(2*(n_r1)*k_eff+(2*n_eff*(k_r1))); end
where for simplicity it is assumed that is air, i.e. and .
Second, the scalar input parameter eps_r
is changed to a twovalued 1Darray nk_r
.
Third, the line,

phy.eps_r = eps_r; % Relative permittivity
is replaced with,

phy.nk_r = nk_r; % Refracitve index and extrinction 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); Ψ
v.2 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, 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;
v.3 Designing a Metallic Reflector
By introducing the changes presented in Sec. V.1 and Sec. V.2, 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+DesignThicknessElements1]’,1,DomainElementsX); Ψ
Second, we solve the optimization problem by executing the command,

[DVs,obj]=top200EM([200,100],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, (=nk_r(1)
), and extinction cross section, (=nk_r(2)
).^{1}^{1}1These values corresponds to values for gold at 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. 4b, one can interpret the design as the wellknown parabolic reflector broken into pieces to fit the spatiallylimited design domain. From the maxnormalized field presented in Fig. 4a the focussing effect of the reflector is clearly observed.
v.4 Linking Design Variables
Certain fabrication techniques limit the allowable geometric variations in a design. For example, opticalprojection lithography and electronbeam lithography restricts variations in design geometries to twodimensional patterns, which can then be extruded in the outofplane direction. It is straight forward 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 inplane (xdirection) variations as follows.
First, the codeline 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 codeline 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 codeline 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.
v.5 Designing a 1D metalens
By introducing the changes listed in Sec. V.4 topEM200
can be used to design metasurfaces of fixed height with inplane 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+DesignThicknessElements1]’,1,DomainElementsX); Ψ
Followed by the execution of the command,

[DVs,obj]=top200EM([200,80],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. 5b. Here it is clear to see that the design is now restricted to vary only in the xdirection. It is worth noting that the design is still filtered in both spatial directions, hence the corners of the design appear rounded in Fig. 5b. The maxnormalized field presented in Fig. 4a demonstrates the focussing effect of the lens at the targeted point in the modelling domain, which, due to the reduced design freedom, is not as high as in the unrestricted case (Sec. IV.1).
v.6 Continuation of Threshold Sharpness
For some design problems in photonics and plasmonics, e.g. CHRISTIANSEN_VESTER_2019; CHRISTIANSEN_OE_2020; WANG_ET_AL_2011, intermediate values may be present in the final optimized design, i.e. , 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/1design 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 BUILTIN OPTIMIZER: FMINCON ΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr); Ψ[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options); Ψ
with,

Ψwhile filThr.beta<betaMax % Thresholding sharpness bound ΨΨ% SOLVE DESIGN PROBLEM USING MATLAB BUILTIN OPTIMIZER: FMINCON ΨΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr); ΨΨ[dVs,~] = 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 Sec. IV.1, 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.
Vi What about genetic algorithms?
For a wide range of inverse design problems, where gradientbased optimization methods are applicable^{2}^{2}2Inverse design problems where it is possible to compute the sensitivities of the FOM (and any constraints), e.g. using adjoint sensitivity analysis, a prerequisite for using gradientbased methods., they have been found to severely outperform nongradientbased methods, such as genetic algorithms (GA) GOLDBERG_1989. 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 degrees of freedom it is feasible to consider (a difference of many orders of magnitude AAGE_ET_AL2017). Further, in many cases, gradientbased methods are able to identify better local optima for the FOM SIGMUND_2011.
In the following, we provide a demonstration of these claims, by comparing the solution of a metalens design problem obtained using the gradientbased top200EM
to the solution obtained using MATLABs 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/1design 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’,’interiorpoint’,... Ψ’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,... Ψ’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr); Ψ% SOLVE DESIGN PROBLEM USING MATLAB BUILTIN OPTIMIZER: FMINCON ΨFOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr); Ψ[dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options); Ψ
with,

Ψoptions = optimoptions(’ga’,’MaxGenerations’,maxItr,’Display’,’iter’); Ψ% SOLVE DESIGN PROBLEM USING MATLAB BUILTIN OPTIMIZER: GA ΨFOM = @(dVs)OBJECTIVE(dVs,dis,phy,filThr); Ψrng(1,’twister’); % SETTING RANDOM SEED TO ENSURE REPRODUCTION OF RESULT Ψ[dVs,~,~,~] = 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:’); Ψ[FOM,~] = OBJECTIVE(dVs(:),dis,phy,filThr); Ψ
to,

Ψ% FINAL BINARIZED DESIGN EVALUATION ΨfilThr.beta = 1e5; Ψdisp(’Black/white design evaluation:’); Ψ[FOM,~] = OBJECTIVE(dVs(:),dis,phy,filThr); Ψ
Note: The codelines plotting the design in OBJECTIVE
can be removed to reduce wallclock 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 degrees of freedom than in the previous examples. This is done by defining a model problem considering 1000 design degrees of freedom as,

Ψ% DESIGN FIELD INDICES ΨDomainElementsX = 100; ΨDomainElementsY = 50; ΨDesignThicknessElements = 10; ΨDDIdx = repmat([1:DomainElementsY:DomainElementsX*DomainElementsY],... ΨDesignThicknessElements,1); ΨDDIdx = DDIdx+repmat([35:35+DesignThicknessElements1]’,1,DomainElementsX); Ψ
Followed by the solution of the design problem using the modified GAbased code,

Ψ[dVs_GA,FOM_GA]=top200EMGA([50,10],DDIdx,DomainElementsX,DomainElementsY,... Ψ0.5,3.0,20,3.0,500); Ψ
and then using the original gradientbased code,

Ψ[dVs,FOM]=top200EM([50,10],DDIdx,DomainElementsX,DomainElementsY,... Ψ0.5,3.0,20,3.0,500); Ψ
The binarized designs obtained using top200EMGA
and top200EM
are presented in Fig. 6a and Fig. 6b, respectively. The evolution of the FOM value as function of the number of FOM evaluations is plotted in Fig. 6c for both methods.
The GAbased method uses 137 design iterations (27,600 FOM evaluations) to identify a solution. In contrast the gradientbased method uses 208 iterations (208 FOM evaluation) to identify a superior solution. Evaluating the final binarized designs the GAbased solution has a FOM value of while the gradientbased solution has an better FOM value of .
Crucially, the GAbased 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 GAbased 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 gradientbased method a total of 208 FOM evaulations are performed, corresponding to 208 solution 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 LUfactorization.
Hence, in our example the identification of a local optimum for the FOM using the gradientbased method required of the computational effort spend by the GAbased method. A ratio, which 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 nongradientbased approaches are unsuitable for largescale TopOpt problems. A more thorough discussion of nongradientbased methods in TopOpt can be found in Ref. SIGMUND_2011. Readers are invited to form their own conclusions by extending the example and codes provided here.
Vii Conclusion
We have presented a simple finite element based MATLAB code (downloadable from https://www.topopt.dtu.dk) for TopOptbased 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, linked design variables and introducing a continuation scheme for the threshold sharpness designed to promote 0/1designs. Finally, we have demonstrated the superiority of a gradientbased method over a geneticalgorithmbased 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 in eq. (8) and the expression in eq. (9) may be derived as follows. First, zero is added to twice,
(11) 
where and is a vector of nodal complex Lagrange multipliers, also called the adjoint variables. Second, one takes the derivative of with respect to and exploits that, for the optimization problem in eq. (II), does not depend explicitly on and neither nor F depend on at all, yielding,
(12)  
Collecting terms including and and reducing the remaining terms yields,
(13) 
To eliminate the first two terms in eq. (13) containing the derivates and , the two parenthesis must equal zero,
(14) 
multiplying the second equation by i, subtracting it from the first and transposing it yields,
(15) 
Appendix B MATLAB code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%% A 200 LINE TOPOLOGY OPTIMIZATION CODE FOR ELECTROMAGNETISM %%%%%%%% %  EXAMPLE GOAL  % % Designs a 2D metalens with relative permittivity eps_r capable of % % monocromatic focusing of TEpolarized light at a point in space. % %  FIGURE OF MERIT MAXIMIZED  % % Phi = Ez^2 in a "point" (in the center of a finite element) % %  EQUATION SOLVED  % % \nabla * (\nabla Ez) + k^2 A Ez = F % % With first order absorping boundary condition: % % n * \nabla Ez =  i k Ez on boundaries % % and an incident plane wave propagating from bottom to top % %  DOMAIN AND DISCRETIZATION  % % The equation is solved in a rectangular domain, discretized using % % quadrilateral bilinear finite elements % %  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% Author: Rasmus E. Christansen, June 2020 %%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Disclaimer: % % The authors reserves all rights but does not guaranty that the code is % % free from errors. Furthermore, we shall not be liable in any event % % caused by the use of the program. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [dVs,FOM]=top200EM(targetXY,dVElmIdx,nElX,nElY,dVini,eps_r,lambda,fR,maxItr) % SETUP OF PHYSICS PARAMETERS phy.scale = 1e9; % Scaling finite element side length to nanometers phy.eps_r = eps_r; % Relative permittivity phy.k = 2*pi/(lambda*phy.scale); % Freespace wavenumber % SETUP OF ALL INDEX SETS, ELEMENT MATRICES AND RELATED QUANTITIES dis.nElX = nElX; % number of elements in x direction dis.nElY = nElY; % number of elements in y direction dis.tElmIdx = (targetXY(1)1)*nElY+targetXY(2); % target index dis.dVElmIdx = dVElmIdx; % design field element indices in model of physics [dis.LEM,dis.MEM] = ELEMENT_MATRICES(phy.scale); [dis]=INDEX_SETS_SPARSE(dis); % Index sets for discretized model % SETUP FILTER AND THRESHOLDING PARAMETERS filThr.beta = 5; % Thresholding sharpness filThr.eta = 0.5; % Thresholding level [filThr.filKer, filThr.filSca] = DENSITY_FILTER_SETUP( fR, nElX, nElY); % INITIALIZE DESIGN VARIABLES, BOUNDS AND OPTIMIZER OPTIONS dVs(length(dis.dVElmIdx(:))) = dVini; % Design variables LBdVs = zeros(length(dVs),1); % Lower bound on design variables UBdVs = ones(length(dVs),1); % Upper bound on design variables options = optimoptions(’fmincon’,’Algorithm’,’interiorpoint’,... ’SpecifyObjectiveGradient’,true,’HessianApproximation’,’lbfgs’,... ’Display’,’off’,’MaxIterations’,maxItr,’MaxFunctionEvaluations’,maxItr); % SOLVE DESIGN PROBLEM USING MATLAB BUILTIN OPTIMIZER: FMINCON FOM = @(dVs)OBJECTIVE_GRAD(dVs,dis,phy,filThr); [dVs,~] = fmincon(FOM,dVs(:),[],[],[],[],LBdVs,UBdVs,[],options); % FINAL BINARIZED DESIGN EVALUATION filThr.beta = 1000; disp(’Black/white design evaluation:’) [FOM,~] = OBJECTIVE_GRAD(dVs(:),dis,phy,filThr); end %%%%%%%%%%%%%%% OBJECTIVE FUNCTION AND GRADIENT EVALUATION %%%%%%%%%%%%%%%% function [FOM,sensFOM] = OBJECTIVE_GRAD(dVs,dis,phy,filThr) % 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 % FILTERING THE DESIGN FIELD AND COMPUTE THE MATERIAL 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 % CONSTRUCT THE 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); S = sparse([dis.iS(:);dis.iBC(:)],[dis.jS(:);dis.jBC(:)],[dis.vS(:);dis.vBC(:)]); tic; % SOLVING THE STATE SYSTEM: S * Ez = F [L,U,Q1,Q2] = lu(S); % LU  factorization Ez = Q2 * (U\(L\(Q1 * F))); Ez = full(Ez); % Solving toc; % 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 % 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 FOM = FOM; sensFOM = sensFOM(:); % 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(1dFPST); colormap(gray); axis equal; drawnow; disp([’FOM: ’ num2str(FOM)]); % Display FOM value end %%%%%%%%%%%%%%%%%%%%%%%%%% AUXILIARY FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%% ABSORBING BOUNDARY CONDITIONS AND RIGHT HAND SIDE %%%%%%%%%%%% function [dis,F] = BOUNDARY_CONDITIONS_RHS(waveVector,dis,scaling) AbsBCMatEdgeValues = 1i*waveVector*scaling*[1/6 ; 1/6 ; 1/3 ; 1/3]; % ALL BOUNDARIES HAVE ABSORBING BOUNDARY CONDITIONS dis.iBC = [dis.iB1(:);dis.iB2(:);dis.iB3(:);dis.iB4(:)]; dis.jBC = [dis.jB1(:);dis.jB2(:);dis.jB3(:);dis.jB4(:)]; dis.vBC = repmat(AbsBCMatEdgeValues,2*(dis.nElX+dis.nElY),1); % BOTTOM BOUNDARY HAS INCIDENT PLANE WAVE F = zeros((dis.nElX+1)*(dis.nElY+1),1); % System right hand side F(dis.iRHS(1,:)) = F(dis.iRHS(1,:))+1i*waveVector; F(dis.iRHS(2,:)) = F(dis.iRHS(2,:))+1i*waveVector; F = scaling*F; end %%%%%%%%%%%%%%%%%%%%% CONNECTIVITY AND INDEX SETS %%%%%%%%%%%%%%%%%%%%%%%%% function [dis]=INDEX_SETS_SPARSE(dis) % INDEX SETS FOR SYSTEM MATRIX nEX = dis.nElX; nEY = dis.nElY; % Extracting number of elements nodenrs = reshape(1:(1+nEX)*(1+nEY),1+nEY,1+nEX); % Node numbering edofVec = reshape(nodenrs(1:end1,1:end1)+1,nEX*nEY,1); % First DOF in element dis.edofMat = repmat(edofVec,1,4)+repmat([0 nEY+[1 0] 1],nEX*nEY,1); dis.iS = reshape(kron(dis.edofMat,ones(4,1))’,16*nEX*nEY,1); dis.jS = reshape(kron(dis.edofMat,ones(1,4))’,16*nEX*nEY,1); dis.idxDSdx = reshape(dis.edofMat’,1,4*nEX*nEY); % INDEX SETS FOR BOUNDARY CONDITIONS TMP = repmat([[1:nEY];[2:nEY+1]],2,1); dis.iB1 = reshape(TMP,4*nEY,1); % Row indices dis.jB1 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEY,1); % Column indices TMP = repmat([1:(nEY+1):(nEY+1)*nEX;(nEY+1)+1:(nEY+1):(nEY+1)*nEX+1],2,1); dis.iB2 = reshape(TMP,4*nEX,1); dis.jB2 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEX,1); TMP = repmat([(nEY+1)*(nEX)+1:(nEY+1)*(nEX+1)1;(nEY+1)*(nEX)+2:(nEY+1)*(nEX+1)],2,1); dis.iB3 = reshape(TMP,4*nEY,1); dis.jB3 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEY,1); TMP = repmat([2*(nEY+1):nEY+1:(nEY+1)*(nEX+1);(nEY+1):nEY+1:(nEY+1)*(nEX)],2,1); dis.iB4 = reshape(TMP,4*nEX,1); dis.jB4 = reshape([TMP(2,:);TMP(1,:);TMP(3,:);TMP(4,:)],4*nEX,1); dis.iRHS = TMP; % INDEX SETS FOR INTEGRATION OF ALL ELEMENTS ima0 = repmat([1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4],1,nEX*nEY).’; jma0 = repmat([1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4],1,nEX*nEY).’; addTMP = repmat(4*[0:nEX*nEY1],16,1); addTMP = addTMP(:); dis.iElFull = ima0+addTMP; dis.jElFull = jma0+addTMP; % INDEX SETS FOR SENSITIVITY COMPUTATIONS dis.iElSens = [1:4*nEX*nEY]’; jElSens = repmat([1:nEX*nEY],4,1); dis.jElSens = jElSens(:); end %%%%%%%%%%%%%%%%%% MATERIAL PARAMETER INTERPOLATION %%%%%%%%%%%%%%%%%%%%%%% function [A,dAdx] = MATERIAL_INTERPOLATION(eps_r,x,alpha_i) A = 1 + x*(eps_r1)  1i * alpha_i * x .* (1  x); % Interpolation dAdx = (eps_r1)*(1+0*x)  1i * alpha_i * (1  2*x); % Derivative of interpolation end %%%%%%%%%%%%%%%%%%%%%%%%%%% DENSITY FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [xS]=DENSITY_FILTER(FilterKernel,FilterScaling,x,func) xS = conv2((x .* func)./FilterScaling,FilterKernel,’same’); end function [ Kernel, Scaling ] = DENSITY_FILTER_SETUP( fR, nElX, nElY ) [dy,dx] = meshgrid(ceil(fR)+1:ceil(fR)1,ceil(fR)+1:ceil(fR)1); Kernel = max(0,fRsqrt(dx.^2+dy.^2)); % Cone filter kernel Scaling = conv2(ones(nElY,nElX),Kernel,’same’); % Filter scaling end %%%%%%%%%%%%%%%%%%%%%%%%%%%% THRESHOLDING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ xOut ] = THRESHOLD( xIn, beta, eta) xOut = (tanh(beta*eta)+tanh(beta*(xIneta)))./(tanh(beta*eta)+tanh(beta*(1eta))); end function [ xOut ] = DERIVATIVE_OF_THRESHOLD( xIn, beta, eta) xOut = (1tanh(beta*(xIneta)).^2)*beta./(tanh(beta*eta)+tanh(beta*(1eta))); end %%%%%%%%%%%%%%%%%%%%%%%%%%% ELEMENT MATRICES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [LaplaceElementMatrix,MassElementMatrix] = ELEMENT_MATRICES(scaling) % FIRST ORDER QUADRILATERAL FINITE ELEMENTS aa=scaling/2; bb=scaling/2; % Element size scaling k1=(aa^2+bb^2)/(aa*bb); k2=(aa^22*bb^2)/(aa*bb); k3=(bb^22*aa^2)/(aa*bb); LaplaceElementMatrix = [k1/3 k2/6 k1/6 k3/6 ; k2/6 k1/3 k3/6 k1/6; ... k1/6 k3/6 k1/3 k2/6; k3/6 k1/6 k2/6 k1/3]; MassElementMatrix = aa*bb*[4/9 2/9 1/9 2/9 ; 2/9 4/9 2/9 1/9 ; ... 1/9 2/9 4/9 2/9; 2/9 1/9 2/9 4/9]; end
Funding
This work was supported in by Villum Fonden through the NATEC (NAnophotonics for TErabit Communications) Centre (grant no. 8692) and by the Danish National Research Foundation through NanoPhoton Center for Nanophotonics (grant no. DNRF147).
Disclosures
The authors declare that there are no conflicts of interest related to this article.