Inverse design of optical mode converters by topology optimization: tutorial

This tutorial details the use of topology optimization (TopOpt) for the inverse design of electromagnetic mode-converters. First, the design problem under consideration is stated. Second, suitable models for the geometry and physics are formulated and third the TopOpt method is outlined. Then follows three increasingly advanced design examples. In the first, the mode converter is allowed to consist of a non-physically-realizable material distribution, leading to a design exhibiting near perfect conversion from the input mode i to the output mode o in terms of power conversion , providing a performance benchmark. Then follows two examples demonstrating the imposition of relevant restrictions on the design, first ensuring a physically realizable device blueprint, and second introducing feature-size control and ensuring device connectivity. These examples demonstrate how TopOpt can be used to design device blueprints that only require a minimum of post-processing prior to fabrication, which only incur a minor reduction of performance compared to the initial unconstrained design. A software tool is provided for reproducing the first design example. This tool may be extended to implement the other design examples in the paper, to explore other device configurations or, given sufficient computational resources, to design 3D devices.


Introduction
This paper provides a tutorial for the application of densitybased topology optimization (TopOpt) [1][2][3] to the design of optical mode converters [4,5]. Three design examples of increasing complexity are provided, demonstrating that device blueprints supporting near perfect mode conversion, designed to adhere to fabrication limitations, may be created using TopOpt. Although only considering examples of conversion of optical modes between waveguides, the method is general and may be applied to a wide range of mode-conversion problems by adjusting relevant details in the tutorial appropriately.
While having been proven capable of solving a wide range of electromagnetic design problems [6,7], TopOpt in undeniably a non-trivial method to apply, potentially barring researchers, scientists and engineers alike from utilizing the tool for their applications. This work serves to lower any such barriers, and to this end a software tool, based on COMSOL Multiphysics [8], is provided along with this text, allowing the reader to reproduce the first design example without any implementation work required.
TopOpt, and other related inverse design methods, have experienced rapidly growing interest in recent years for a variety of electromagnetics applications, ranging from its early adaptation for design of photonic-crystal-based devices [9] over photonic cavity design [10][11][12], the design of optical lenses [13] and concentrators [14], through more exotic applications such as designing topological insulators [15,16] to the design of optical multiplexer and mode-converters [17,18] to name but a few. The latter being most directly relevant to this paper.
While the basics of the TopOpt method are presented in the following, not all aspects are explained in detail for brevity. The interested reader is instead recommended to explore these details by consulting prior works, such as tutorials on the basics of the TopOpt method in the context of photonics [3,19] and/or the review paper by Jensen and Sigmund [6] and references therein.
In brief, the density-based TopOpt method is a widely applicable, large-scale (in terms of design degrees of freedom) [20] inverse design method utilizing gradient-based optimization and adjoint sensitivity analysis [21] for efficient gradient calculations. It relies on a differentiable model for the design problem (geometry and physics) in question, most often a set of partial differential equations along with appropriate boundary conditions defined and solved on a particular modeling domain; on a mathematical field, discretized using a (high) number of design variables, that controls the geometry of the device under design; on a set of auxiliary tools used to constrain and modify the mathematical field prior to interpolating device geometry, using said field to ensure physical realizability, device fabricability etc for the optimized design. Different realizations of the method, with more-or-less functionality, have been implemented across a variety of research, opensource, and commercial codebases.
As outlined above, TopOpt may be used for a wide class of photonics device design problems. That being said, the goal of this tutorial is to demonstrate how TopOpt may be used specifically to design optical mode converters. The procedure explained in this work may however be adapted to other photonic device design problems by adjusting the individual steps appropriately.
To limit the computational costs associated with the design examples, thus allowing readers to reproduce results on a standard laptop/desktop computer, the examples are restricted to two spatial dimensions. The procedure detailed in the following is however directly applicable for design problems in both one, two and three spatial dimensions. Thus an interested reader, with access to sufficient computational resources, will be able to extend the examples to three spatial dimensions by making the appropriate minor modifications.

Design problem example
The mode-converter design problem (MCDP) considered as the example in the following is, MCDP: Given an optical waveguide A supporting the propagating mode i at the angular frequency ω, and an optical waveguide B supporting the propagating mode o at ω, sketched in figures 1(C) and (D) respectively: Design a device which, when inserted between A and B as sketched in figure 1 A point-by-point procedure for solving MCDP using simulation-based design tools may be stated as,  figure 1). 4. Excite the input waveguide in the mode-converter design model problem (at Γ P,in in panel E of figure 1) using the field profile of mode i calculated in step 2, and compute the resulting electromagnetic field, E C , in the model domain Ω C . 5. Using a suitable measure, calculate the overlap between the electromagnetic field propagating along the output waveguide B found in step 4, and the field profile for mode o calculated in step 3. 6. Systematically adjust the mode-converter design (situated in Ω d in figure 1(E)), to improve the overlap calculated in step 5. 7. Repeat steps 4 through 7 until a stopping criterion is reached.
To solve MCDP using TopOpt, the design problem must be recast as a continuous optimization problem based on the models for the geometry and physics as detailed in the following sections.

Model geometries
The model geometries employed in the solution of MCDP are sketched in figure 1 panels A-E. These are used in the computation and analysis of the input mode i (panels A-B), the output mode o (panels C-D) and the field resulting from the input mode propagating from waveguide A, through the mode converter under design, to waveguide B (panel E  (h B,wg ) and extends along the x-direction through Ω A (Ω B ), centered along the y-direction. A port, at which mode i(mode o) is coupled into the waveguide, is located at the center of the domain, oriented normal to the waveguide. The model domain in panel E, labeled Ω C , consists of a central region of width w MCDP and height h MCDP bounded to the left(right) by a light green(blue) region of width w PML , labeled Ω PML . The domain contains the input waveguide A (green) of height h A,wg , connected to a design domain Ω d where the mode converter is situated, which in turn is connected to the output waveguide B (blue) of height h B,wg . The design domain is surrounded on all sides by a narrow region (transparent light gray) employed in the inverse design procedure as outlined in section 2.6. The port at which the input(output) mode is excited(recorded) is placed δ Port from the left(right) edge of the center domain. The parameters required to setup the five model geometries, and the values assumed in the design examples, are given in table 1

Physics models
The electromagnetic field is modeled using Maxwell's equations [22], assuming linear, static, homogeneous, isotropic, non-dispersive, non-magnetic materials as well as timeharmonic electromagnetic field behavior, i.e. E = Ee −iωt and H = He −iωt , where E(H) is the spatial distribution of the electric(magnetic) field, i is the imaginary unit, ω is angular frequency and t is time.
First the input(output) mode i(mode o) in waveguide A(B) is computed along with the associated propagation constant. This is done by solving an eigenvalue problem, with appropriate boundary conditions, on a waveguide cross-section embedded in a sufficiently large region of surrounding background material to avoid significant boundary effects, as sketched in panel A(B) of figure 1. For the two-dimensional problem considered in the following, the waveguide crosssection is one dimensional and is infinitely extruded in the zdirection by employing periodic boundary conditions, labeled Γ periodic,1 (2) , and truncated using perfect electric conductors in the y-direction, labeled Γ PEC . The solution to the eigenvalue problem will be a set of modes and propagation constants, from which the appropriate mode and constant is selected and stored.
For the geometries in panels C and D of figure 1, periodic boundary conditions are imposed between the left and right boundaries, labeled Γ periodic,1 (2) . First-order scattering boundary conditions are imposed on the top and bottom boundaries, labeled Γ Scatt . A port boundary condition is imposed along the green(blue) line at the center of the domain, labeled Γ P,A (Γ P,B ), at which the mode i(mode o) computed in the previous step is introduced. Note that it is not necessary to compute the electromagnetic field in Ω A (Ω B ), as only the boundary mode profile along Γ P,A (Γ P,B ) for mode i(mode o) is required in the design of the mode converter. However, calculating the electromagnetic field in Ω A (Ω B ) serves as a tool for visualizing and investigating the propagating mode thus ensuring that the model is correctly implemented.
For the model geometry in panel E, a first order scattering boundary condition is imposed along all outer boundaries, labeled Γ Scatt . A perfectly matched layer [23] is introduced in the light green/blue regions, labeled Ω PML , to minimize reflections from the termination of the waveguides. The input mode i is excited at the port boundary condition (dark green line) labeled Γ P,in , and the output mode o is evaluated at the port (dark blue line), labeled Γ P,out .
The physics model equations solved to obtain the profiles for mode i and mode o, and to compute the electromagnetic field distribution in the domains Ω A , Ω B and Ω C , are, Eigenvalue problem solved to compute the cross-section of mode i in waveguide A yielding E A,i (r) =Ẽ A,i (y)e −βi x , r ∈ Γ P,A and β i (ω): Eigenvalue problem solved to compute the cross- Transmission problem modeling mode i propagating along waveguide A yielding E A (r), r ∈ Ω A at β i (ω), ω, ε r (r): Transmission problem modeling mode i propagating along waveguide B yielding E B (r), r ∈ Ω B at β o (ω), ω, ε r (r): Transmission problem modeling the field propagation in Here, n is the surface normal, r is the spatial coordinate, β * are modal propagation constants, c is the speed of light in vacuum and ε r ∈ {ε r,wg , ε r,bg } is the relative electric permittivity for the waveguides and mode converter and for the background medium, respectively. Note that when solving the design problem using TopOpt the relative permittivity is actually computed from the refractive index, n, and extinction cross-section, κ, using equation (17). The magnetic field, H, is computed from the electric field as, where µ 0 is the vacuum permeability. For the design examples, it is assumed that the waveguides and mode converter consist of silicon and that the background material is air (vacuum). The list of parameter values used when solving the physics model problems is given in table 2

The waveguide (Eigen)modes
To solve MCDP the electric-and magnetic-field profiles of mode i and mode o and the associated propagation constants must be known. As outlined in section 2.3, these are calculated by solving the eigenvalue problems in equations (1) and (2). In all design examples, the lowest order TE mode in waveguide A and the first higher order TE mode in waveguide B are selected as mode i and mode o. These modes are visualized propagating along waveguide A and B by solving the time-harmonic model problems stated in equations (3) and (4), using the previously computed mode profiles as the excitation at the port boundaries. The resulting electric field magnitude for E A (r) and E B (r) in Ω A and Ω B are presented in panels A and B of figure 2 on a max-normalized colormap.

The objective function-mode conversion
With the design problem stated (section 2.1) and the geometry and physics models defined (sections 2.2 and 2.3), the next step is to recast the design problem as a mathematical optimization problem, which in turn may be solved using a gradient-based optimization algorithm. To this end, the objective function to be optimized is selected as follows.
Given the model problem in equation (5), an optimal solution to MCDP consists of a device, which when introduced into Ω d in panel E of figure 1, losslessly converts mode i propagating in waveguide A to mode o propagating in waveguide B. Another way of stating this is, that all power propagating in mode i in waveguide A (introduced in the model through the port condition at Γ P,in ), is transferred through the device in Ω d to mode o in waveguide B. Labeling the time-averaged power flow coupled in to the domain through the port at Γ P,in as P in and the resulting timeaveraged power flow in mode o in waveguide B as P o,B one can define the figure of merit (FOM), which equals unity for the perfect transmission of power from the external input to the desired output mode o in waveguide B. 3 If power is lost to scattering or absorption, or if power is coupled to another mode in waveguide B, Φ will take a value between zero and one. The FOM in equation (7) can now be used to formulate the objective function for the TopOpt problem through the following steps. First, recall that the time-averaged power flow through a surface Γ, is computed from the time-averaged Poynting vector as, where ℜ denotes the real part and • * the complex conjugate. Thus, to compute the power flow through a given surface, one needs to know the electric and magnetic fields at said surface. Second, under normal circumstances an electromagnetic field propagating along a waveguide can be expanded in an infinite series of orthogonal modes, where the modal coefficients e k and h k may be computed as, with Γ being a plane intersecting the waveguide. Note, in practice a truncation using a finite number of terms is sufficient to expand the field to sufficient accuracy. Utilizing the modal decomposition 4 , and exploiting the orthogonality of the modes, the time-averaged power flow of the electromagnetic field may be written as a sum of the individual modal contributions, with the power flow in mode k calculated as, 4 Assuming that the electromagnetic field is sufficiently accurately captured by a finite number of terms in the expansion ensuring interchangeability of the modal sum and power flow integral.
Inserting equation (10) into equation (12) and assuming that the power coupled into waveguide A is coupled perfectly into mode i, P in = P i,A , one may write the FOM in equation (7) as, Here E C (H C ) is the electric(magnetic) field obtained by solving the model problem stated in equation (5).
at Γ P,in (only the desired mode is excited in waveguide B) the expression in equation (13) reduces to the ratio of the power in mode o at Γ P,out to the power in mode i at Γ P,in . Hence if these powers are equal, then no scattering or absorption occur and perfect mode conversion is achieved.
The expression in equation (13) can be evaluated based on the solution to the physics model problems (section 2.3). It is taken as the objective function for the TopOpt problem.

The basic TopOpt problem
To solve MCDP using TopOpt [3], the problem is recast and solved as a mathematical optimization problems of the form, Here, the scalar field ξ(r) constitutes the optimizable quantity (the design field), that is used to control the material layout in Ω d . That is ξ controls the geometry of the mode converter, with ξ = 1 corresponding to the device material and ξ = 0 to the background material. The functions c e,k and c i,l are used to impose a set of, problem dependent, equality and in-equality constraints on the optimization problem. In short, the problem of determining the geometry of the mode converter is recast as the problem of maximizing a function Φ (equation (13)) through iterative and systematic modification of ξ, while respecting all constraints imposed on the problem. Note that, crucially to the efficiency of the TopOpt method, ξ is allowed to vary continuously between zero and one. This choice enables the use of efficient gradientbased optimization algorithms in the solution of equation (14). The main challenge that this choice introduces, which will be demonstrated in the first design example, is that the optimized values of ξ(r) likely consists of large areas of intermediate, and thus non-physically meaningful values. A suite of tools has been developed to eliminate the problem of intermediate design field values, as will be demonstrated in the second design example.
A standard filtering procedure [25] is applied to the design field ξ to limit rapid spatial oscillations using the equation, where r f is the filter radius. Note that the filtering procedure is carried out over Ω d,e to enable feature-size control along the edges of the design domain Ω d (see section 2.11. The filter step is followed by the application of a smoothed approximation of a threshold operator [26], one of the tools used to eliminate Here β is the threshold strength and η the threshold level. By applying the threshold operation, along with a gradual increase of β during the iterative solution of the optimization problem, it is possible to recover a design field consisting solely of device and background materials The filtered and thresholded design fieldξ is coupled to the physics model via the material interpolation [27], n(ξ) = n bg +ξ(n wg − n bg ), Through the steps outlined above ξ(r) controls the geometry of the mode converter, and thus any change to the design field will cause a change in the mode-converter layout.
The parameters related to solving the optimization problem and the values used in the examples are listed in table 3 Note that TopOpt is almost always employed to solve design problems that, when recast as optimization problems, are non-convex [28]. This has, among others, the two implications that there is no guarantee that the solution identified by TopOpt will constitute the globally optimal solution, and that the final design geometry is often found to be sensitive to the initial guess. In fact, in practice one will almost never discover the global optimum for the design problem. For sufficiently geometrically sensitive photonic design problems, exhibiting multiple local extrema with significantly varying FOM-values, examples being grating [29] or cavity design [10], a carefully chosen starting guess, or a large number of starting guesses may be required to attain satisfactory device performance when applying TopOpt. That being said, in practice it is authors experience that even for such difficult problems, most often only a few different initial guesses are needed to achieve satisfactory performance. For less sensitive photonic design problems, such as coupler design [30], metalens design [31] beam-splitter design [32] and the present case of modeconverter design, a simple uniform initial guess is often sufficient to achieve (near) optimal device performance. In cases where one already has a working device geometry for the problem at hand, this can also be used as the initial guess for the optimization, however this risk getting the inverse design process stuck in a local optimum as the working device geometry might already be (nearly) 'locally' optimal. The fact that the design problems are generally non-convex means that solving them using different initial guesses are likely to result in different optimized geometries. That being said, for the three design examples considered in the following, only a single initial guess was employed and near perfect mode conversion achieved.

Practical implementation.
The numerical implementation of the physics models (sections 2.2 and 2.3) was achieved using the finite element method [33]. The optimization problem (section 2.6) was implemented based on the discretized physics model and was solved using the gradient-based optimization algorithm, the method of moving asymptotes [34]. The gradients of the objective function and constraints were computed using adjoint sensitivity analysis [21] carried out based on the discretized physics model (see appendix for an outline of the approach.). For the examples in the following the mode-converter design domain was discretized using square elements with a side length of 20 nm, first order basis functions and nodal design variables.
The tutorial software supplied with this work is implemented in COMSOL Multiphysics [8] in a similar style to the tutorial software provided for the more basic TopOpt tutorial in [3]. Readers interested in the details of the underlying finiteelement implementation of the model and TopOpt problems, are recommended to consult [19] and associated software.
Here an electromagnetic metalens design problem is solved using a freely available 200 line MATLAB code.

Example 1 -the naive approach
As the first design example, MCDP is solved using the models and procedure presented in sections 2.2-2.6, with the parameters in tables 1-3. This results in the design presented in panel A of figure 3 showing the input waveguide (black) to the left, the mode converter (gray-scale) in the middle, and the output waveguide (black) to the right. The optimized material distribution in Ω d clearly consists of a non-physical mixture (gray) of air (white) and silicon (black). The presence of large regions of gray-scale means that the device blueprint is not physically realizable, outside perhaps further development and application of non-standard fabrication techniques, such as gradient index direct laser writing [35].
Setting a side for a moment the issue of realizing the design. Visual inspection of panels B-C, showing |E C (r)| and the ycomponent of E C (r) respectively, suggests near perfect mode conversion from the input to the output waveguide. Computing the reflectance and transmittance for the device, one obtains R ≈ 0.000 025 and T ≈ 0.999 with a modal overlap between the targeted mode o and the input mode i, measured in terms of the power flow, of P o,B /P i,A ≈ 0.996. That is, nearly all power flowing through the mode converter to waveguide B is transmitted into mode o as requested in MCDP.
Returning to the realizability of the design, naively binarizing the design blueprint in panel A of figure 3 by thresholding it around the mean value ofξ = 0.5, result in a poorly performing design as illustrated in panels D-F of figure 3.
Panel D shows the binarized design, which is now seen to consist solely of air (white) and silicon (black). Panels E and F show |E(r)| and the E y (r)-component of the electric field, respectively. It is immediately obvious from these panels that the device no longer functions as an efficient mode converter. Indeed, computing reflectance, transmittance and mode-conversion efficiency one obtains, R ≈ 0.23, T ≈ 0.58 and P o,B /P i,A ≈ 0.017. Clearly, a naive binarization does not result in the design of high-performance devices for MCDP. Fortunately a suite of TopOpt tools, which ensure physically realizable (pure black and white) designs with high performance, have been developed. The application of a subset of these tools resolve the binarization issue, as will be demonstrated next.

Continuation and pamping
When solving transmission-dependent design problems, like MCDP, a combination of penalization through damping (or pamping) [36] and continuation of the threshold strength [26] result in physically realizable designs with high performance.
In brief, pamping consists of introducing artificial attenuation in the physics model for intermediate value ofξ. For the design problem at hand, this is done by modifying the material interpolation scheme (equation (17)) as, n(ξ) = n bg +ξ(n wg − n bg ), where α i (=0.01 in the following) is a coefficient controlling the magnitude of the attenuation introduced whenξ takes values other than 0 or 1. Continuation of the threshold strength in equation (16) is implemented by increasing β every n β design iterations from an initial value β ini to a final value β final . The ideal rate at which β is increased along with the values of β ini and β final are design problem dependent. In the following β is increased every n β = 50 design iterations in five increments with an increase of β by 50% per increment starting at β ini = 5.

Example 2 -obtaining a physically realizable design
For the second design example, the TopOpt tool is modified by replacing the material interpolation in equation (17) by equation (18) and by using the continuation procedure for the threshold operation. Otherwise the solution of MCDP is unchanged from the previous example. This results in the optimized design presented in panel A of figure 4, which is seen to consist solely of air (white) and silicon (black).
Panel B shows |E C (r)| and panel C the y-component of E C (r) for the design in panel A. Studying these fields, the design appears to achieve near-perfect mode conversion. Corroborating this qualitative observation are the reflectance, transmittance and mode-conversion efficiency, which are calculated to R ≈ 0.000 025, T ≈ 0.999 and P o,B /P i,A ≈ 0.999, respectively. Thus, employing pamping and continuation in the TopOpt procedure had near-zero influence on the device performance, while their use resulted in a design blueprint that is now physically realizable. That being said, closer inspection of the design reveals other potential challenges to the fabrication of the device. Firstly, a number of fine details are observed in the design blueprint, which might not be amenable to reliable and accurate fabrication. Secondly, the device contains a disconnected island of silicon towards its lower left corner. If the mode converter is to be realized as a silicon-on-insulator device, 'free floating' islands of material are permissible as they will rest on a substrate. However, if the device is to be membranized, such islands cannot be realized.
These observations are manifestations of issues inherent to the TopOpt tool presented thus far. That is, there is nothing in the current formulation of the method that prohibits disconnected islands of material, nor that ensures a minimum size of individual features. Both of which may prohibit accurate fabrication. Next, it is explained and demonstrated how these issues are resolved.

Feature-size control and solid-feature connectivity
Different fabrication tools have different limits to the design features they can accurately manufacture. Therefore, the ability to specify and control the length-scales of a device directly in the design process is of great value. A useful measure of feature-size is the maximum radius of the brush (ball), with which it is possible to accurately 'paint' the feature. To demonstrate feature-size control using TopOpt, the method proposed by Zhou et al [37] is employed in the third design example. This method has, among others, proven useful for fabrication of optimized designs by electron beam lithography [12]. The method uses the indicator functions I s and I v to define two integrals, which measure if the solid(void) parts of a design contain features smaller than the specified brush (ball) radius. The indicator functions and integrals may be written as, , , 0} ] 2 dr, Here η e and η d , together with the radius, r f , of the filter operation (equation (15)), are used to determine the minimum feature-sizes on the solid(void) features as detailed in [38]. The constant c LS is a tuning parameter, the value of which is selected to make the indicator functions numerically well-behaved, as detailed in [37]. The integrands in equation (20) are strictly non-negative and hence the value of the integrals is zero if and only if no feature in a given design is below the specified minimum feature-size. Hence, one may specify the constraints, where ϵ s > 0 and ϵ v > 0 are introduced to relax the constraints, which enable gradually imposing the constraints during the iterative design process allowing for the design to initially develop without having to adhere to the constraints. Further, the relaxation allow for unavoidable numerical errors in the evaluation of the integrals. In the third example ϵ s and ϵ v are gradually, monotonically decreased from 1 (inactive constraints) to 10 −5 over the course of 10 continuation steps with 50 design iterations per step. The minimum feature-size imposed on the final design for the third example is 50 nm. This is achieved by changing the filter-radius to r f = 100 nm and selecting η e = 0.75 and η d = 0.25. Note that in the example the integrals in equation (20) are evaluated over the extended design domain, Ω d,e , while the design is only free to change inside the design domain, Ω d . This ensures that the specified minimum feature-sizes are also respected along the edges of Ω d , which would otherwise not be the case.
When designing suspended [12] multilayered [14] or fully three-dimensional devices, physics dictates that no solid features are allowed to be disconnected from the rest of the device, since free-floating members are not possible to realize. The way of prohibit free-floating islands of material, when designing a device using TopOpt, is to impose a connectivity constraint as part of the optimization problem. In this work the connectivity constraint is formulated using a heat-transfer problem. Conceptually the constraint may be understood as follows. Consider any (sufficiently) solid material in Ω d,e as a heat source which is also highly conductive, and consider any background material as being insulating. Next, define boundaries to which the solid material must be connected as perfect heat sinks and all other boundaries as perfect insulators. Now, if a solid feature is connected to a heat sink, (nearly) all heat generated by said feature will be conducted to the sink. Thus, if all solid features are connected to heat sinks the temperature everywhere in the device will be low. On the contrary, if any solid feature is disconnected from the heat sinks, it will generate heat that cannot be conducted away, as the surrounding background material insulate it, hereby creating a region of high temperature. Thus, by integrating the temperature field over Ω d,e one obtains a measure of connectivity. If the integrated temperature is below a certain threshold value all solid features will be connected, while if one or more solid features are completely disconnected from the rest of the design, the integrated temperature will exceed the threshold.
In practice the constraint is implemented using the following system of equations defined on Ω d,e , η)) , Here C F denotes the artificial temperature (or connectivity) field. The constants c 0 (=10 −6 ), c 1 (=10 10 ) denote the artificial conductivity of the background and design material, respectively. The constants S 0 (=0), S 1 (=10 20 ) denote the artificial heat generated by the background and design material, respectively. β CF = β S = 50 and η CF = η S = 0.55 are the threshold strength and threshold level used to determine what parts of the developing design,ξ, that is counted as solid material for the evaluation of the constraint. Imposing the constraint that, where ϵ C is a sufficiently small number, ensures that all solid material (silicon) is connected to the boundaries Γ D where the heat sinks are situated. A detailed explanation of the connectivity constraint and its implementation may be found in [39]. In the following design example, the constraint ensure that all solid features are connected to waveguide A or waveguide B. To achieve this the heat sink (zero Dirichlet) and insulator (zero Neumann) boundary conditions, are imposed as sketched in panel H of figure 1.

Example 3 -ensuring a fabricable design
Imposing the feature-size constraints (equation (21) and the connectivity constraint (equation (22)) along with using pamping and threshold continuation as in the second design example, the design problem MCDP is solved for a third time. This results in the optimized design presented in panel A of figure 5. The design again consists solely of air (white) and silicon (black). Further, it is now observed that all silicon features are connected and that the specified feature-sizes are respected (illustrated using the red and orange discs).
Panels B and C of figure 5 show |E C (r)| and the y-component of E C (r), respectively. Again (near) perfect mode conversion seemingly occur from waveguide A to waveguide B. The reflectance, transmittance and mode-conversion efficiency are computed to R ≈ 0.0001, T ≈ 0.995 and P o,B /P i,A ≈ 0.995, respectively. Thus, the performance dropped by less than half a percent in order to ensure that the design blueprint is physically realizable and adheres to the specified feature-size and connectivity limitations. Notably without changing any other parameters in the model, such as the size of the design domain or the initial guess.

Conclusion
A step-by-step tutorial for how TopOpt can be applied as a tool for the inverse design of optical mode converters was provided. It was demonstrated that TopOpt is capable of designing high-performance mode-converters, which are physically realizable, respect specified feature-sizes and ensuring that the device geometry is connected. All examples are kept two-dimensional to reduce the required computational effort, allowing readers to reproduce the examples in a few hours on a standard laptop. However, given sufficient computational resources it is straight-forward to extend all examples, and the provided tutorial software, to full three-dimensional problems by modifying the model geometries and physics problems appropriately. Treating the full 3D problem it may not be possible to achieve as extreme performance as the 2D examples in this work, since light can then be lost through out-of-plane scattering. Further, it is straight-forward to adapt the approach outlined here to other mode-conversion systems, such as metasurface-based mode-converters, also by adjusting the geometry and physics models accordingly. Finally, additional constraints can be imposed on the design problem, such as requiring a specific reflectance back into the input waveguide, or designing mode-converters that operate across specified bandwidths.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Funding
This work was supported by the Danish National Research Foundation through the NanoPhoton-Center for Nanophotonics, Grant No. DNRF147.

Conflict of interest
The author declare no conflicts of interest.

Appendix. Discrete adjoint sensitivity analysis
Being able to efficiently compute the change in the FOM and constraints with respect to design perturbations, i.e. the gradient with respect to the design variables, is essential for solving the inverse design problem efficiently using a gradientbased optimization algorithm. To this end, the TopOpt method employs adjoint sensitivity analysis [21]. Adjoint sensitivity analysis may be performed either directly on the model equations before numerical discretization (differentiate then discretize), or on the discretized model equations (discretize then differentiate). In this work the latter approach is taken, as illustrated with an example in the following.
The design field, ξ, is discretized using a set of coefficients (design variables) along with a set of basis functions as, The system of model equations in equation (5) is discretized using the finite element method [33]. The resulting linear system of equations may be written as, where S is the design dependent system matrix and E(F) is a vector of degrees of freedom for the electric field(forcing). In the discretized model, the FOM may be written as a function of the electric field vector, Φ(E). Using the technique for integration of composite functions (the chain rule), and for the sake of simplicity assuming that no operations (smoothing, thresholding etc) are applied to the design field, the gradient of the FOM with respect to design variable j is computed as, where E ℜ (E ℑ ) denotes the real(imaginary) part of the electric field vector. The gradient in equation (25) is rewritten as follows. First, zero is added to the FOM twice, where λ is a vector of complex values unknowns (Lagrange multipliers), also called adjoint variables. Then, one takes the derivative ofΦ with respect to ξ j , exploiting that neither λ nor F depend on ξ j . This yields, ) .
Next, collecting the terms including ∂E ℜ ∂ξj and ∂E ℑ ∂ξj and reducing the expression yields, The first two terms in equation (27), containing the expensive to compute derivates ∂E ℜ ∂ξj and ∂E ℑ ∂ξj , may be eliminated by requiring that, ∂Φ ∂E ℜ + λ T S + λ † S * = 0, ∂Φ ∂E ℑ + iλ T S − iλ † S * = 0, (28) multiplying the second equation by i, subtracting it from the first equation and transposing the result yields, Requiring that equation (29) is satisfied, the expression in equation (27) reduces to, Thus, all that is needed to compute the gradient of the FOM with respect to the design field is to calculate the right hand side in equation (29) and solve this equation system once, independent on the number of design variables.
An example of the derivation of the right hand side in equation (29) follows here. For simplicity in this example, a transverse electric polarization is assumed for the 2D physics model, resulting in the electric and magnetic fields being given as,