Unifying and accelerating level-set and density-based topology optimization by subpixel-smoothed projection

Photonics topology optimization (TopOpt) is a form of “inverse design” [1] in which freeform arrangements of materials are evolved to maximize a given design objective (e.g. transmitted power), with millions of degrees of freedom (DOFs) that can be used to explore nearly arbitrary geometries, including arbitrary topologies (e.g. how many “holes” are present). TopOpt algorithms are divided into an unfortunate dichotomy, however: there are level-set set methods [2] that can exactly describe discontinuous “binary” materials (e.g. silicon in some regions and air in others) but make it more difficult for optimization to change the topology [3]; versus density-based methods [3, 4] that allow the topology to change by smoothly evolving one material to another, but which become more and more difficult to optimize as the structure is “binarized” to be everywhere one material or another. The key difficulty centers around the computation of derivatives: in order to optimize efficiently over many parameters, the design objective must be differentiable with respect to the DOFs, and the derivatives must be well behaved, e.g. lacking sharp “kinks” where the second derivative (Hessian matrix) is suddenly huge in some directions (such an “ill-conditioned” Hessian slows optimization convergence [5, 6]). In density-based methods, binarization involves a “projection” step in which the material is projected towards one extreme or an other, controlled by a steepness hyper-parameter $\beta$, but the projection becomes more ill-conditioned as $\beta$ increases and the structure becomes more binary (more physical). In this paper, we develop a new form of projection that allows us to bridge the gap between density-based and level-set methods: we employ density DOFs $\rho$ that can smoothly describe changes in topology, but our new projection scheme also allows us to increase $\beta$ to $\infty$ (equivalent to a level set) while retaining well-behaved derivatives.

In a level-set method, the DOFs are a level-set function $\phi(\vec{x})$ which denotes one material (e.g. silicon) when $\phi<0$ and another material (e.g. air) when $\phi>0$; unfortunately, when the topology changes (e.g. two boundaries intersect or a new hole appears), differentiation with respect to $\phi$ requires complicated “topological derivative” algorithms [3] that can be onerous to implement, so in practice a level-set implementation may be restricted to the initially chosen topology (which is also known as “shape optimization” with a fixed topology [3]). In density-based TopOpt (reviewed in Sec. 2), the DOFs are a density function $\rho(\vec{x})$ for which $\rho=0$ represents one material, $\rho=1$ represents another material, and intermediate $0<\rho<1$ values represent unphysical “interpolated” materials (e.g. interpolating the refractive index) [7]. This makes the physical solution easy to differentiate with respect to $\rho$ and allows the topology to continuously change (e.g. a new “hole” can appear as $\rho$ goes smoothly from 1 to 0 in some region), but optimization may lead to unphysical designs in which many intermediate materials are present. To combat this, and also to regularize the problem by introducing a minimum lengthscale, density-based methods introduce two additional density fields [8, 9, 10, 11, 12] (Fig. 1a–c): a smoothed $\tilde{\rho}$ given by the convolution of $\rho$ with a minimum-lengthscale low-pass filter, followed by a projection $\hat{\rho}$ that passes $\tilde{\rho}$ through a step function to binarize it to either $0$ or $1$. A step function is not differentiable, however, so in practice the projection step uses a smoothed “sigmoid” function characterized by a steepness parameter $\beta$, as shown in the top row of Fig. 2. As $\beta$ increases towards $\infty$, the structure becomes more binarized, at the cost of slower optimization convergence because the derivative approaches zero everywhere except at the step where it diverges (and its second derivative becomes ill-conditioned). Managing this tradeoff requires cumbersome empirical tuning of the rate at which the hyper-parameter $\beta$ is increased during optimization, in order to obtain a mostly binarized structure with good performance in a reasonable time.

To overcome this challenge, we present a new density projection ($\hat{\tilde{\rho}}$), “subpixel-smoothed projection” (SSP), which supports a seamless transition from density-based to level-set optimization. Our approach offers a drop-in replacement for existing projection functions, and allows designers to leverage all their existing topology-optimization tooling—the only implementation difference is that our $\hat{\tilde{\rho}}(\mathbf{x})$ depends on both $\tilde{\rho}$ and $\|\nabla\tilde{\rho}\|$ at $\mathbf{x}$, rather than just on $\tilde{\rho}$ as for the old $\hat{\rho}$. Let the computational mesh spacing/resolution be of order $\Delta x$. Conceptually, SSP as if the filtered density $\tilde{\rho}$ were first interpolated and projected to yield the conventional $\hat{\rho}$ at “infinite” resolution (regardless of $\Delta x$), which is then convolved with a $\Delta x$-diameter smoothing filter to yield a density $\hat{\tilde{\rho}}$ (depicted in Fig. 2 and Fig. 1d–e) that is binary almost everywhere (for $\beta=\infty$) as the mesh resolution increases ($\Delta x\to 0$), but which changes differentiably as $\rho$ changes. As depicted in Fig. 1, we make this conceptual algorithm practical by instead computing $\hat{\tilde{\rho}}$ directly from the filtered density $\tilde{\rho}$, with the “infinite-resolution smoothing” performed analytically via a planar/linear approximation of the $\tilde{\rho}=0.5$ level set. In particular, Sec. 3 proposes an efficient algorithm, compatible with any dimensionality, that yields a mostly twice-differentiable pointwise function $\hat{\tilde{\rho}}(\tilde{\rho},\|\nabla\tilde{\rho}\|)$, except at changes in topology where it reduces to $\hat{\rho}$ (Sec. 3.1). We demonstrate the flexibility of SSP by applying it in Sec. 4 to various TopOpt problems with three different Maxwell solvers: finite-difference time domain (FDTD) [13], a Fourier modal method (FMM) [14], and a finite-element method (FEM) [15]; the same approach is also directly applicable to non-photonics TopOpt problems (Sec. 5). Several test problems are taken from our recently published photonics-testbed suite [16], and SSP results in final designs comparable to traditional TopOpt, but with faster convergence at large $\beta$ and more rapid $\beta$ increases. We also show how SSP can be applied to both topology optimization (by incrementing the projection steepness $\beta$ in a few steps to $\beta=\infty$) as well as to shape optimization (in which the binary design is optimized for a fixed topology with $\beta=\infty$, which is now differentiable (for a fixed topology). We demonstrate that SSP’s better-conditioned formulation of the projection problem can require fewer optimization iterations at large $\beta$ than the conventional density-projection scheme. By eliminating the density/level-set dichotomy, we believe that SSP should not only improve performance and simplicity—$\beta$ can be increased to $\infty$ as soon as the topology stabilizes, without sacrificing convergence rate—but may additionally lead to improvements in other aspects of topology optimization, such as techniques for imposing lengthscale constraints or more accurate discretization schemes (Sec. 5).

Conceptually, density-based topology optimization maps some array of density values to a real-space grid that is discretized in some fashion (e.g. finite differences, finite elements). Each pixel within this array is then continuously optimized by interpolating between two or more constituent materials. To ensure the final design is binary (and thus, physically realizable), the threshold condition ($\beta$) of a nonlinear projection function is gradually increased. The canonical optimization problem explored throughout this work corresponds to a frequency-domain, differentiable maximization problem over a single figure of merit (FOM) $f(\mathbf{E})$:

where $\rho$ is the spatially varying material “density” (parameterized by the design degrees of freedom or “DOFs” $\boldsymbol{\rho}$), $\mathbf{E}_{n}$ is the time-harmonic electric-field solution at frequency $\omega_{n}$, $\mathbf{\epsilon}_{r}$ is the relative permittivity as a function of the density $\rho$ at each point in space, and $\mathbf{J}_{n}$ is the time-harmonic current density (the source term). While we solve problems constrained by Maxwell’s equations throughout this paper, we emphasize that our new smoothed-projection approach is agnostic to the underlying physics and applicable to any TopOpt problem. Adjoint methods (also known as “backpropagation” or “reverse-mode” differentiation) allow one to compute the gradient of any FOM with respect to all DOFs using just two Maxwell solves (the “forward” and “adjoint” problems) [1, 4].

As discussed earlier, to map the “latent” design variables ($\rho$) to the permittivity profile ($\varepsilon$), one first filters the design variables such that

where $\tilde{\rho}$ is the filtered design field, $\tilde{k}$ is the filter kernel, and $*$ is a problem-dependent multidimensional convolution [8, 12]. The filter radius $\tilde{R}$ (or bandwidth $\sim 1/\tilde{R}$) is typically determined by the minimum desired lengthscale in the design, and is much larger than the underlying grid/mesh spacing $\Delta x$—this regularizes the optimization problem to ensure convergence with resolution and is also related to enforcement of fabrication constraints [11, 17]. Alternatively, one can omit an explicit filtering step by representing the degrees of freedom directly in terms of $\tilde{\rho}$, represented by some smooth/band-limited basis/parameterization [18, 19, 20].

Next, the filtered design variables are projected onto a nearly binary density field, $\hat{\rho}$, using the sigmoid projection function, $P_{\beta,\eta}(\mathbf{\tilde{\rho}})$:

where $\beta$ controls the strength of the binarization ($\beta\to\infty$ corresponds to a Heaviside step function of $\tilde{\rho}-\eta$, while $\beta\to 0^{+}$ gives the identity mapping $\hat{\rho}\to\tilde{\rho}$) and $\eta$ describes the location of the threshold point (the level set $\tilde{\rho}=\eta$). Finally, these projected densities are linearly interpolated between two (possibly anisotropic) material tensors $\mathbf{\varepsilon_{0}}$ and $\mathbf{\varepsilon_{1}}$ at each point $\mathbf{x}$, most commonly via

Such linear interpolation between two materials works well if the materials are both non-metallic (positive $\varepsilon$). Different interpolation techniques are employed when optimizing with metals to avoid nonphysical zero crossings in the interpolated $\varepsilon$ [21], and more generally any differentiable pointwise mapping from $\hat{\rho}$ to the material properties could be employed.

Optimization typically proceeds as follows: (i) the design is initialized with some guess (e.g. a uniform “gray” region $\rho=0.5$) and a relatively low value for $\beta$; (ii) an optimization algorithm iteratively updates this design according to the value of the FOM and its gradient at each step; (iii) once optimization is sufficiently converged, $\beta$ is increased and the optimization algorithm is restarted using the best performing topology from the previous “epoch”. This process is repeated until the final design is acceptably binary and the performance meets the desired specifications. This “$\beta$-scheduling” procedure is highly problem dependent, requiring significant effort and experience-gained heuristics to avoid slow convergence. If one increases $\beta$ too rapidly, then the optimization problem quickly becomes stiffer and stiffer, slowing down the overall convergence. Conversely, if one evolves $\beta$ too slowly, then finding a high-performing binary design is also computationally expensive because of the large number of optimization epochs.

Importantly, when differentiating through the physics of the problem, one first computes the gradient with respect to $\hat{\rho}$ (which directly determines the physical parameters), rather than $\boldsymbol{\rho}$ (the design variables). As such, one must use the chain rule (“backpropagation”) to compute the gradient with respect to $\boldsymbol{\rho}$, which is needed by the optimization algorithm. Consequently, the step-function projection is a problem: the resulting gradient approaches zero almost everywhere as $\beta\to\infty$ (and becomes non-differentiable at the step itself). In other words, standard density-based TopOpt methodologies cannot produce useful gradients for purely binary designs. The amount of binarization is directly related to the value of $\beta$, but the designer must avoid making this value too large (regardless of the underlying schedule); at the final step, if $\beta$ is sufficiently large, the structure can hopefully be “rounded” to the closest binary structure without greatly degrading performance. In contrast, level-set methods avoid this issue by establishing sophisticated methods to differentiate with respect to an interface, which becomes especially cumbersome when the topology changes (e.g. when interfaces cross). By introducing subpixel smoothing into the density-based formulation, we will remove the non-differentiable discontinuity induced by the projection function and allow the designer to arbitrarily specify $\beta$.

As $\beta\to\infty$, $\hat{\rho}$ becomes a discontinuous function of $\tilde{\rho}$, and even for large finite $\beta$ one encounters ill-conditioning problems due to the exponentially large second derivatives. These issues are closely related to the fact that $\hat{\rho}$ is a discontinuous function of space $\mathbf{x}$: as $\tilde{\rho}$ changes, the location of the level set $\tilde{\rho}=\eta$ changes smoothly (except at changes in topology), but the discontinuity of $\hat{\rho}$ in space means that this smooth interface motion translates to a discontinuous dependence $\hat{\rho}(\tilde{\rho}(\mathbf{x}))$ at any given $\mathbf{x}$. However, if we could smooth the discontinuities in space—but only near the interfaces so that the projection remains mostly binary—then $\hat{\rho}$ would become a differentiable function of $\tilde{\rho}(\mathbf{x})$ due to the smooth interface motion (except at topology changes, discussed in Sec. 3.1), even for $\beta=\infty$. This strategy corresponds conceptually to a post-projection smoothing, but in practice must be implemented in conjunction with projection as explained below.

That is, we could conceptually remove the $\hat{\rho}$ discontinuity (in space, and hence in $\tilde{\rho}$) by convolving $\hat{\rho}$ with a smoothing kernel, if this convolution could somehow be evaluated at infinite spatial resolution, as depicted by the dashed “not practical” arrow in Fig. 1. That is, imagine that we could evaluate the exact convolution ($*$) over the computational domain $\Omega$:

where $\hat{k}(\mathbf{x})\geq 0$ is a localized (radius-$\hat{R}$) smoothing kernel, with $\int k\mathop{}\!\mathrm{d}\Omega=1$ and $\hat{k}=0$ for $\|\mathbf{x}\|>\hat{R}$, and where the support diameter $2\hat{R}\gtrsim\Delta x$ is proportional to the spacing/resolution $\Delta x$ of our computational grid or mesh (below, we chose $\hat{R}=0.55\Delta x$ to make it slightly larger than a voxel). (Because $\hat{\tilde{\rho}}$ will be sampled onto a computational grid or mesh, or at quadrature points in a finite-element method [15], we want $2\hat{R}$ to be slightly greater than the maximum sample spacing $\Delta x$ to avoid “frozen” gradients where the derivative is zero at all sampled $\mathbf{x}$.) If we could do this, then $\hat{\tilde{\rho}}$ would have two desirable properties: (i) it would remain $0$ or $1$ as $\beta\to\infty$ except in a radius-$\hat{R}$ (“subpixel”) neighborhood of interfaces (i.e. binary almost everywhere as the computational resolution increases: $\Delta x\to 0$); and (ii) if $\hat{k}$ is sufficiently smooth, then $\hat{\tilde{\rho}}$ will be a differentiable function of $\mathbf{x}$ and also of the underlying $\rho$ even for $\beta\to\infty$ (as long as the position of the $\tilde{\rho}=\eta$ level set changes smoothly with $\rho$). The challenge is to make the computation of such a $\hat{\tilde{\rho}}$ practical and efficient.

The key fact is that the underlying filtered field $\tilde{\rho}$ is, by construction, smooth and slowly varying (on the scale of the filter radius $\tilde{R}\gg\hat{R}\sim\Delta x/2$), which allows us to smoothly interpolate it to any point $\mathbf{x}$ in space—thus implicitly defining $\hat{\rho}(\mathbf{x})=P_{\beta,\eta}(\tilde{\rho}(\mathbf{x}))$ at infinite resolution—and furthermore it implies that $\tilde{\rho}$ can normally be linearized in a small neighborhood $\|\mathbf{x}^{\prime}-\mathbf{x}\|\leq\hat{R}\ll\tilde{R}$:

This allows us to vastly simplify the smoothing computation of $\hat{\tilde{\rho}}$ as explained below, and will ultimately allow us to completely replace the integral (5) with a simple local function of $\tilde{\rho}$ and $\|\nabla\tilde{\rho}\|$. Moreover, in the linearized approximation, the signed distance $d$ from any point $\mathbf{x}$ to the (approximately planar) $\tilde{\rho}=\eta$ level set (i.e., the $\hat{\rho}$ discontinuity at $\beta=\infty$) is simply:

as depicted in Fig. 3. Note that we define $d>0$ or $d<0$ for $\tilde{\rho}$ below or above $\eta$, corresponding to $\hat{\rho}(\mathbf{x})=0$ or $1$ for $\beta\to\infty$, respectively. (For cases where optimization causes two interfaces to meet, changing the topology of the level set, the effect of this linearization is discussed in Sec. 3.1.)

The next step is to choose the smoothing kernel to be rotationally symmetric $\hat{k}(\mathbf{x})=\hat{k}(\|\mathbf{x}\|)$, and to rewrite the convolution (5) in spherical coordinates. (Henceforth, we will always view the convolution as being performed in 3d: even if the underlying problem is 1d or 2d, we will treat it as conceptually “extruded” into 3d, e.g. a 2d $\rho(x_{1},x_{2})$ is treated as a 3d $x_{3}$-invariant function.) In particular, we center the spherical coordinates at $\mathbf{x}$, where the “$z$” axis ($\theta=0$, $z\neq x_{3}$) is oriented in the $\nabla\tilde{\rho}$ direction, and approximate the integral (5) in terms of the linearized $\tilde{\rho}$ from Eq. (6):

In principle, Eq. (8) is already computationally tractable: at each point $\mathbf{x}$ in our grid/mesh where we want to compute $\hat{\tilde{\rho}}$, we evaluate $\tilde{\rho}$ and $\nabla\tilde{\rho}$, and then apply a numerical-quadrature routine to evaluate the integral. This is especially cheap for $\beta\to\infty$ because the integral is simply $1$ or $0$ except for points $\mathbf{x}$ where $|d|<\hat{R}$, so quadrature need only be employed at a small set of points near the interface. However, we can make further efficiency gains by the choice of kernel $\hat{k}$ and other approximations described below.

First, consider the limit $\beta\to\infty$, for which $P_{\infty,\eta}$ is a Heaviside step function. Suppose that the point $\mathbf{x}$ is within $|d|<\hat{R}$ of the $\hat{\rho}$ interface, and is on the $\hat{\rho}=0$ side ($d\geq 0$) as in Fig. 3. In this case, the integral (8) simplifies even further to a function $F(d)$ that depends on $d$ alone:

For any given kernel $\hat{k}(r)$, this function $F(d)$ could be precomputed. If the kernel $\hat{k}$ is simply a “top-hat” filter $\hat{k}(r)=\text{constant}$ for $r\leq\hat{R}$ and $\hat{k}(r)=0$ for $r>\hat{R}$, then $F(d)$ is a “fill factor” of the $\hat{\rho}=1$ region inside the sphere, and can be computed analytically. For $d<0$, we let $F(-d)=1-F(d)$: the fill factor is inverted (equivalent to inverting $\hat{\rho}\to 1-\hat{\rho}$). However, rather than specifying $\hat{k}(r)$ and computing $F(d)$, it is attractive to instead choose $F(d)$ and infer $\hat{k}(r)$, because $F(d)$ directly determines the smoothness of $\hat{\tilde{\rho}}$ as a function of $\rho$ (recalling that $d$ is a smooth function of $\tilde{\rho}$, which in turn is a smooth function of $\rho$). The algebraic details are given in Appendix A, where we construct an inexpensive, twice-differentiable, monotonic, piecewise-polynomial choice of $F$, denoted as $F_{2}(d)$:

which is continuous with a continuous derivative and second derivative, satisfies $F_{2}(-d)=1-F_{2}(d)$, and turns out to correspond to a parabolic kernel $\hat{k}_{2}(r)$ proportional to $1-(r/\hat{R})^{2}$ for $r\leq\hat{R}$. We also show how to construct an infinitely differentiable $F$ if desired: see Eq. (25).

Thus, for $\beta\to\infty$, this $\hat{\tilde{\rho}}$ function is twice differentiable (in $d$ and hence $\rho$) by construction (except at changes in topology as discussed in Sec. 3.1), is almost-everywhere binary (in the limit of $\Delta x\to 0$), and is notably cheap to compute at any desired point $\mathbf{x}$ in space:

1. 1.

   Given the smooth $\tilde{\rho}$ function, interpolate it and its gradient to $\mathbf{x}$ (e.g. via bilinear interpolation in a gridded finite-difference method),

2. 2.

   Compute $d$ by Eq. (7),

3. 3.

   Compute $\hat{\tilde{\rho}}(\mathbf{x})=F(d)$ for the chosen $F$, e.g. Eq. (10) or (25). This is generalized to Eq. (11) for finite $\beta$ below.

Each of these steps is differentiable, so backpropagating derivatives to $\rho$ for adjoint differentiation is straightforward. Much like the computation of $\hat{\rho}$ in standard TopOpt, $\hat{\tilde{\rho}}$ is a purely local function of $\tilde{\rho}$ at each $\mathbf{x}$, with the only difference being that it now depends on both $\tilde{\rho}$ and $\|\nabla\tilde{\rho}\|$. (In a discretized grid/mesh, such local quantities are determined by interpolating adjacent mesh points.) Moreover, even though $\hat{\tilde{\rho}}\to\hat{\rho}$ as $\Delta x\to 0$, it remains a well-conditioned function of the $\tilde{\rho}$ values (it has bounded second derivatives) on any discretized mesh, as discussed in Appendix B.

For finite $\beta$, one could fall back to numerical integration of Eq. (8). However, we can greatly speed up the finite-$\beta$ computation by making two approximations, which still result in a differentiable $\hat{\tilde{\rho}}$ that reproduces the convolution (8) above as $\beta\to\infty$. First, for points $|d|\geq\hat{R}$ far from the level set, we simply omit the smoothing and choose $\hat{\tilde{\rho}}(\mathbf{x})=\hat{\rho}(\mathbf{x})$ as in conventional TopOpt: far from an interface, the projection $\hat{\rho}$ is nearly constant within the smoothing radius $\hat{R}$, so there is not much point in smoothing it further. Second, for points $|d|<\hat{R}$ close to the level set, we employ the $\beta\to\infty$ fill-fraction function $F(d)$ as above, with the modification that we average values of $\hat{\rho}$ evaluated on the two sides of the level set, denoted by $\hat{\rho}_{-}$ and $\hat{\rho}_{+}$ and defined below, instead of averaging $0$ and $1$. Rather than an approximation for Eq. (8), this can be viewed as simply a new definition of $\hat{\tilde{\rho}}$, namely:

where $\hat{\rho}_{\pm}$ are obtained by projecting the linearized approximation of $\tilde{\rho}$ at two points along the $\nabla\tilde{\rho}$ direction (so that they still depend only on $\tilde{\rho}$ and $\|\nabla\tilde{\rho}\|$ at $\mathbf{x}$) within $\hat{R}$ of $\mathbf{x}$:

corresponding to two points on either side of the level set separated by a distance $\hat{R}F(-d)+\hat{R}F(d)=\hat{R}$. This $\hat{\tilde{\rho}}$, demonstrated in Fig. 4, is constructed to make Eq. (11) monotonic and twice differentiable in $d$, and also so that Eq. (11) reduces to $\hat{\tilde{\rho}}\to F(d)$ (as above) in the limit $\beta\to\infty$ (in which case $\hat{\rho}_{-}\to 0$ and $\hat{\rho}_{+}\to 1$) as well as to $\hat{\tilde{\rho}}\to\tilde{\rho}$ as $\beta\to 0^{+}$ (no projection), as can be shown via straightforward algebra. In particular, note that Eq. (12) re-uses our function $F(d)$ purely for convenience, exploiting the fact that $F$ goes smoothly from $1$ at $d=-\hat{R}$ to $0$ at $d=\hat{R}$ to make Eq. (11) transition smoothly to the $\hat{\rho}(\mathbf{x})$ case as $d\to\pm\hat{R}$.

Our new subpixel-smoothed projection $\hat{\tilde{\rho}}$ (11) is computationally inexpensive while offering two fundamental benefits over traditional projection (3): (i) the ability to optimize even as $\beta\to\infty$; and (ii) improved convergence for large finite values of $\beta$ (because it has a bounded second derivatives regardless of $\beta$, from Appendix B, leading to better-conditioned optimization problems [6]). In Sec. 4, we quantify these advantages with numerical experiments, a simple example of which is shown in Fig. 5.

Here, we present numerical experiments using three different Maxwell discretization techniques to illustrate the flexibility of our new SSP method. First, we design a 2d silicon-photonic crossing using a hybrid time-/frequency-domain adjoint solver [7] compatible with Meep, a free and open-source finite-difference time-domain (FDTD) code (FDTD) [23]. The SSP formula (11) is implemented as a vectorized function in Python, automatically differentiable via the autograd software [24]. We demonstrate that SSP is directly compatible with shape optimization (where an initial design topology is already identified and optimized at $\beta=\infty$) and with freeform topology optimization (at finite $\beta$ then $\beta=\infty$ as discussed in Sec. 3.1) on an example 2d waveguide-crossing [25] design problem. We also show greatly accelerated convergence compared to traditional density projection (3) for a large finite $\beta=32$.

Next, we illustrated the versatility of SSP by applying it to other problems taken from our photonics-optimization test suite [16], as well as to other computational-electromagnetics methods. We design a 3d diffractive metasurface using FMMAX [26], a differentiable and GPU-accelerated Fourier Modal Method (FMM) code written in Python. We then design a 2d metallic nanoparticle to enhance near-field concentration using a finite-element electromagnetics solver implemented in the Julia software package Gridap [27], demonstrating that the same method is applicable to discretizations employing unstructured meshes.

We believe that the subpixel-smoothed projection (SSP) scheme introduced in this paper should be widely useful for topology optimization (TO) in many fields—nothing about the construction of our $\hat{\tilde{\rho}}$ is specific to electromagnetism. SSP both simplifies and accelerates density-based TO, allowing $\beta$ to be increased very quickly without concern for ill-conditioning or “frozen” gradients, and at the same time guarantees almost-everywhere binarization of the final structure (without the need for additional penalty terms in problems where the physics might otherwise push $\tilde{\rho}$ to $\approx\eta$ in order to counteract finite-$\beta$ projection [43, 44, 45]). Furthermore, because SSP simply replaces $\hat{\rho}(\tilde{\rho})$ with a new function $\hat{\tilde{\rho}}(\tilde{\rho},\nabla\tilde{\rho})$, our approach should be easy to “drop in” to existing TO software without extensive modifications, and we have shown that it works well in both uniform-grid (e.g. FDTD and FMM) and unstructured-mesh (FEM) discretizations.

Furthermore, we believe that there are many opportunities for future developments building off of the smoothed-projection ideas presented in this paper. We have already begun to revisit methods to impose fabrication-lengthscale constraints [11, 17], which we believe can be further simplified and refined now that one can simply set $\beta=\infty$ to ensure binarization. As mentioned in Sec. 3.1, there is the potential for a more sophisticated smoothed-projection scheme that allows one to differentiate through changes in topology for $\beta=\infty$ (instead of requiring $\beta$ to be increased in a few stages) by relaxing the $\tilde{\rho}$-linearization approximation. (Even with the current scheme, we sometimes observe changes in topology at $\beta=\infty$ where the optimizer “steps over” the non-differentiability, e.g. to delete a small hole, but this cannot be relied on with algorithms whose convergence proofs assume differentiability.) Another exciting opportunity is to exploit the level-set information in $\tilde{\rho}$ to improve accuracy of the PDE discretization of the material interfaces—in other work, we’ve shown that one can improve accuracy in electromagnetism using anisotropic smoothing [46, 47], and there has also been other research on high-order discretization schemes exploiting subpixel interface knowledge [48]. In the finite-element context, we are working to both improve accuracy of the matrix assembly (quadrature over the mesh elements) and allow smaller subpixel smoothing radii $\hat{R}$, by adaptively subdividing the mesh in the vicinity of the level set, analogous to methods developed for immersed-boundary FEM techniques [49].