Nanophotonic devices can provide solutions to challenges in energy conversion, information technologies, chemical or biological sensing, quantum computing, and secure communications. The realization of practical optical structures and devices is a complex problem due to the multitude of constraints on their optical performance, materials, scalability, and experimental tolerances, all of which are requirements implying large optimization spaces. However, despite the complexity of the process, to date, almost all nanophotonic structures are designed either intuitively or based on a priori selected topologies, and by adjusting a limited number of parameters. These intuition-based models are limited to ad hoc needs and have narrow applicability and predictive power, with the exhaustive parameter searches often performed manually. Since the comprehensive search in hyper-dimensional design space is highly resource-heavy, multi-objective optimization has so far been almost impossible. Humans' restrained capacity to think hyper-dimensionally also limits the perception of multivariate optimization models, and, therefore, advanced machinery is needed to manage the multi-domain, hyper-dimensional design parameter space. In this work, we merge the topology optimization method with deep learning algorithms, such as adversarial autoencoders, and show substantial improvement of the optimization process in terms of computational time (4900 times faster) and final devices efficiencies (∼98%) by providing unparalleled control of the compact design space representations. By enabling efficient, global optimization searches within complex landscapes, the proposed compact hyperparametric representations could become crucial for multi-constrained problems. The proposed approach could enable a much broader scope of the optimal designs and data-driven materials synthesis that goes beyond photonic and optoelectronic applications.

Realization of practical optical structures and devices is an inherently complex problem due to multi-faceted requirements with the manifold of stringent constraints on the optical performance, materials, scalability, and experimental tolerances. These multiple requirements inevitably open up an enormously large optimization space. Despite the complexity of the available parametric space, almost all nanophotonic structures to date are designed either intuitively or based on a priori selected topologies, and by adjusting a very limited number of parameters (e.g., the periodicity, the trivial geometrical shapes, and dimensions of the resonant elements). Such intuition-based models are only useful for ad hoc needs and have limited applicability and predictive power. The exhaustive parameter sweeps are often done “by hand.” Since the comprehensive search in hyper-dimensional design space is highly resource-heavy, multi-objective optimization has so far been almost impossible. Moreover, humans' restrained capacity to think hyper-dimensionally limits our perception of multivariate optimization models. Thus, advanced machinery is needed to manage the multi-domain, hyper-dimensional design parameter space.

The innovatory field of the inverse design has recently been transforming conventional nanophotonics and allowing for the discovery of unorthodox optical structures via computer algorithms rather than engineered “by hand.”1 The realization of “non-intuitive” designs requires truly new approaches combined with already established diverse optimization and sensitivity methods such as genetic algorithms2–4 and different variations of the adjoint method.5–12 Particularly, topology optimization (TO) that previously revolutionized mechanical and aerospace engineering13–15 by providing unexpected solutions to constrained material distribution problems, has recently emerged as a powerful architect for photonic design9,16–20 that offers broader parameter space and flexible incorporation with different computational methods. However, the gradient descent nature of the TO method, in which an initial device design is locally optimized within the parametric space, significantly relies on the initial guess of the material distribution inside the optimization domain. Lack of intuition in choosing the right initial geometries leads to multiple TO runs in order to select the best-performing solution. Since TO is very computationally expensive, this substantially limits its applicability to multi-constrained problems that require a significant expansion of the optimization space to larger parametric domains that include mechanical, chemical, and optical properties.

Recently, different aspects of machine learning (ML) have attracted major interest in the field of nanophotonics.21–25 Various discriminative deep learning models have been adapted to find the solution to direct and inverse electrodynamics problems.26–29 Unlike conventional electrodynamic simulation methods, which require intensive, time-consuming computations, ML algorithms enable almost instantaneous solution searches due to the learning process performed during the training phase. Along with the pure discriminative model, various generative networks, such as generative adversarial networks (GANs)30 and variation autoencoders (VAEs),31 have been used for nanophotonic design optimization. Recently, GANs have been coupled with the TO method for optimizing diffractive dielectric gratings.32,33 It has been shown that adversarial networks could be efficiently trained on topology optimized designs for the rapid generation of large families of highly efficient grating designs in the significantly smaller timescale.

Within this work, we demonstrate that adversarial autoencoders (AAEs) can be efficiently adapted for rapid nanophotonic design optimization. Mainly, we show that an AAE based optimization engine (i) enables >4900 times faster optimization search within the compressed design space (latent space), and (ii) ensures advanced control over the latent space distribution. The latter is essential for multi-constrained optimization problems, where a compact hyper-parametric representation becomes critical for efficient optimization searches within a complex landscape. To showcase our AAE assisted method, we optimize a metasurface thermal emitter for thermophotovoltaic (TPV) applications. Compared to an adjoint-based TO design with 92% efficiency of the thermal emission reshaping, the proposed method provides three times the speedup and gives 98% efficiency. The proposed approach can be adapted to a broader scope of problems in optics, chemistry, and mechanics.

Generative models raised significant interest in the machine learning community due to a fundamentally new approach to data interpretation. Generative networks aim to learn dataset distribution of the training set and generate new data with some additional variations, unlike discriminative models, which learn hard/soft boundaries between classes of the data.

Within this work, we coupled an AAE generative network with the conventional TO method for a highly efficient topology optimization of nanophotonic devices. The chosen AAE consists of three coupler neural networks: encoder, decoder/generator, and discriminator (Fig. 1), following an initial AAE concept.34 Similar to VAEs,35 the main goal of the encoder in AAE is to compress a given input pattern into a compact, continuous design space, a so-called latent space. Then, the decoder learns to reconstruct the real space patterns based on a given latent space coordinate, a latent vector. In contrast to VAEs, where the latent space distribution is assumed to be a standard normal distribution, AAE performs adversarial learning (like in GANs36) by applying the discriminator to force latent space to pre-defined model distribution. AAE can be considered as a combination of VAE and GAN networks.

FIG. 1.

AAE-assisted topology optimization. Starting from a discrete set of the topology optimized designs of nanoantenna that serves as the metasurface building block, the encoder compresses each nanoantenna design into a point in the latent space, a compressed, continuous design space. The decoder reconstructs the design based on the input coordinate in the latent space. The discriminator forces the encoder to construct the latent space with a pre-defined distribution. The trained decoder is then used as a generator, which takes the latent space coordinate as an input and generates a large set of designs. The structure refinement procedure is applied to the generated set to eliminate unstable, low-efficient designs.

FIG. 1.

AAE-assisted topology optimization. Starting from a discrete set of the topology optimized designs of nanoantenna that serves as the metasurface building block, the encoder compresses each nanoantenna design into a point in the latent space, a compressed, continuous design space. The decoder reconstructs the design based on the input coordinate in the latent space. The discriminator forces the encoder to construct the latent space with a pre-defined distribution. The trained decoder is then used as a generator, which takes the latent space coordinate as an input and generates a large set of designs. The structure refinement procedure is applied to the generated set to eliminate unstable, low-efficient designs.

Close modal

Figure 1 shows the proposed flow of the optimization process. AAE assisted optimization consists of the three main steps: (i) data generation using adjoint TO method; (ii) training AAE network on the generated dataset; (iii) structure refinement.

The first step aims at generating efficient designs using the TO and construction of an appropriate training dataset. In the second step, the AAE network is trained on the constructed dataset of TO designs. During the training process, the encoder is forced to produce (i) a latent space that the decoder can use for reconstruction, and (ii) a latent sampling that could pass the discriminator as a sample from the pre-defined model distribution. Once the AAE network is trained, then the decoder can be used as a generator that takes the latent vector as an input and generates a design pattern in the real space. In the third step, a design refinement procedure filters out unstable, low-efficient designs.

The choice of the AAE network over other deep generative network architectures, such as VAEs and GANs, is motivated by its important advantages. First, the AAE method provides the neural network with the dense (continuous) latent space representation34 that becomes critical for interpolating the hyper-dimensional parametric space. Such a continuous representation enables a much broader variety of the generated designs. Second, the AAE approach sets no specific limits on the pre-defined model distributions, hence enabling extremely flexible control over the latent space configuration.34 Finally, the AAE networks have better trainability in comparison with GANs37 because the discriminator in the AAE networks is applied to a compressed continuous latent representation in comparison with GANs, where adversarial learning is applied directly to the patterns/images. Details on the quantitative comparison between GAN and AAE networks can be found in the supplementary material. Within this work, we showcase the proposed approach by optimizing a metasurface thermal emitter design for TPV applications. Section III highlights the main constraints of the problem.

The conventional TPV engine aims at generating electrical power by radiative heat transfer and usually consists of a heater and a photovoltaic (PV) cell array [Fig. 2(a)]. Without losing generality, we consider TPV systems utilizing GaSb PV cells with a working band ranging from λmin=0.5μm to λmax=1.7μm. To ensure efficient electrical power generation, thermal emission of the heater should significantly overlap with the working band of the PV cell. Hence, the temperature of the heater should exceed 1000 °C. Figure 1(b) shows the emission spectrum of the blackbody at 1800 °C. However, even for the appropriately high temperature, only a small portion of the emission overlaps with the PV cell working band [Fig. 2(b), green area], while most of the emission energy remains outside the band [Fig. 2(b), red area]. While in-band radiation produces electron-hole pairs, out-of-band radiation causes undesirable heating of the PV cell, significantly reducing its quantum efficiency and device lifetime. By patterning the surface of the heater with a properly designed thermal emitting metasurface, it is possible to spectrally reshape the emissivity ε(λ) of the heater to maximize the in-band and minimize the out-of-band radiation. The ideal thermal emitter has a step function type profile of the emissivity with ε(λminλλmax)=1 and zero elsewhere [dashed blue contour, Fig. 2(b)].

FIG. 2.

Thermal emitter for TPV applications. (a) Schematic of a TPV engine: a heater patterned with a thermal emitter array and a PV cell. (b) Blackbody radiation of the bare heater (solid black curve) corresponding to emission of blackbody at 1800 °C. The gray rectangular region highlights the GaSb PV cell working band. Only in-band radiation is converted into electrical power (green area), while out-of-band radiation causes heating of the PV cell (red area). Blue dashed contour corresponds to an ideal thermal emitter's emissivity/absorption spectrum. (c) Absorption/emissivity spectrum of the optimized cylindrical gap plasmon thermal emitter (shaded region corresponds to the GaSb PV cell's working band). The cylindrical emitter's unit cell size is 145 nm, radius and height of the cylinder are 50 and 30 nm, respectively, the Si3N4 spacer is 40 nm thick, and the top Si3N4 cover is 90 nm thick. Inset shows the 3D and a side view of the structure.

FIG. 2.

Thermal emitter for TPV applications. (a) Schematic of a TPV engine: a heater patterned with a thermal emitter array and a PV cell. (b) Blackbody radiation of the bare heater (solid black curve) corresponding to emission of blackbody at 1800 °C. The gray rectangular region highlights the GaSb PV cell working band. Only in-band radiation is converted into electrical power (green area), while out-of-band radiation causes heating of the PV cell (red area). Blue dashed contour corresponds to an ideal thermal emitter's emissivity/absorption spectrum. (c) Absorption/emissivity spectrum of the optimized cylindrical gap plasmon thermal emitter (shaded region corresponds to the GaSb PV cell's working band). The cylindrical emitter's unit cell size is 145 nm, radius and height of the cylinder are 50 and 30 nm, respectively, the Si3N4 spacer is 40 nm thick, and the top Si3N4 cover is 90 nm thick. Inset shows the 3D and a side view of the structure.

Close modal

Due to high-temperature operation, thermally emitting metasurfaces should be designed utilizing high-temperature stable material platforms. Recently, it has been demonstrated that transition metal nitrides (TiN, ZrN) exhibit metal-like optical properties and plasmonic attributes on par with noble metals in the visible and near-infrared spectral regions.38–42 In contrast to noble metals conventionally used in plasmonics, transition metal nitrides are stable at very high temperatures.43 In this work, we focus on the TiN/Si3N4 material combination for metasurface thermal emitter designs.44 For more details on dielectric permittivity functions of TiN and Si3N4, see the supplementary material (Sec. S1).

Recently, various selective emitter designs have been investigated including rare-earth oxides,45 photonic crystals,46–48 and metamaterial/metasurface-based emitters.49,50 One of the most commonly used designs of the thermal emitter is a gap plasmon structure,51–54 which consists of a back-reflector, dielectric spacing material, and a top array with plasmonic resonators of simple geometrical shapes.55 This design offers simple fabrication, as well as intuitive design. However, the intuitive simple shapes substantially reduce the degrees of freedom for optimization; and as a result, significantly limit achievable efficiencies. As a reference, we use the parametric optimization of a gap plasmon56 structure that comprises an array of TiN cylindrical resonators. TiN cylinders are deposited on top of a Si3N4 spacer layer that covers an optically thick TiN back-reflector [see the inset in Fig. 2(c)]. Optimization of the array period, dielectric spacer thickness as well as the dimensions of the TiN cylinder (radius and height) is performed with particle swarm optimization method,57 minimizing the norm difference between the emissivity/absorption spectrum of the structure and the ideal emitter emissivity. Figure 2(c) shows the obtained absorption spectrum. While the out-of-band emissivity is substantially suppressed, the mean in-band emissivity/absorption reaches only 84% due to a limited number of resonant in-band modes. This reference case demonstrates that even though the parametric design space is large enough, the trivial initial shape of the resonator fundamentally limits the ultimate achievable efficiency.

As the next step, a material distribution within the simulation domain can be used as a sub-set of the design parameter space. Hence, by applying a gradient-based TO technique to such an optimization domain, it can converge to an optimal, non-intuitive binary material distribution that enables highly efficient device performance. In the next section, we adapt such a density-based TO technique to construct a training dataset for an AAE network to optimize the metasurface thermal emitter design.

To generate a training set with topology optimization, we consider a gap plasmon metasurface configuration. The three-layer structure comprises an optically thick back TiN reflector, a 40-nm-thick Si3N4 spacer, and a 120-nm-thick top layer, with a fixed 280-nm-period of the unit cell in both lateral directions [Fig. 3(a)]. The comparison of the TO of the thermal emitter with different unit cell sizes can be found in the supplementary material. The top layer is defined as the optimization region and discretized with 60 × 60 optimization elements. Initially, the material distribution in the optimization region is set to be a random, smooth dielectric function within a single quadrant of the unit cell. The remaining quadrants are filled using mirror-symmetric reflections along the interior sides. In this case, we enforce a desired mirror-symmetry along x and y directions. Such symmetry choice is allowed to simulate ¼ of the unit cell and use symmetry boundary conditions during the full-wave analysis of TO. In the more general case, setting up a random initial guess for the entire unit cell will lead to topology optimized structures with no symmetry constraints. Interpolation of the material distribution is done using a non-linear interpolation scheme proposed in Ref. 58. TO is realized using an adjoint optimization scheme, which requires two full-wave simulations per iteration (forward and adjoint) and employs a direct commercial full-wave solver built on a finite-difference time-domain (FDTD) approximation of the Maxwell equations (Lumerical FDTD) controlled with a Matlab host script. The spatial distribution of the forward and adjoint fields determines the “heat map” of the dielectric function perturbation inside the optimization domain, which maximizes the figure of merit gradient at a given iteration. After multiple iterations, the material distribution inside the optimization region converges to a binary structure (air/TiN). One of the critical constraints of the optimization procedure is compatibility with available fabrication techniques (i.e., the stability of the final design to fabrication imperfections and achievable fabrication tolerances).6,59 This constraint requires incorporating a two-step robustness algorithm into the optimization procedure: (i) elimination of the sub-precision features by averaging the permittivity of the design over neighboring regions and pushing the design to binary structures with the next optimization evolution; (ii) averaging the figure of merit over the perturbed geometries of the device, hence, reducing the impact of geometric variability on the device efficiency. Thus, we apply spatial filtering, which eliminates the sub-30 nm features from the designs in each 10th iteration of the optimization algorithm, while the total number of TO iterations is set to 50. More details on the implementation of the direct TO can be found in the supplementary material (Sec. S2).

FIG. 3.

Direct Topology Optimization for training dataset generation. (a) The base structure under consideration consists of a 300-nm-thick TiN back-reflector, 40-nm-thick Si3N4 dielectric spacer, and 120-nm-thick top TiN patterned layer in a 280 × 280 nm2 unit cell. The top TiN layer is set to be the optimization region (outlined by red box). (b) Convergence plot of the topology optimization with corresponding evolution of the material distribution. Black and red curves correspond to the in-band and out-of-band average absorption respectively. (c) Designs obtained by the direct topology optimization, including the design with the highest efficiency (92%) framed in the red box (black color corresponds to air, white to TiN). (d) Statistics of efficiency distribution of 200 topology optimized designs. Normalized efficiency is defined as a ratio of in-band emissions of the TO thermal emitter and of an ideal emitter at 1800 °C. (e) The absorption/emissivity spectrum of the best TO (black) and cylindrical (red) thermal emitters. (f) Corresponding emissivity spectra at 1800 °C: the best TO emitter (black), cylindrical emitter (red), and blackbody emission (dashed black curve).

FIG. 3.

Direct Topology Optimization for training dataset generation. (a) The base structure under consideration consists of a 300-nm-thick TiN back-reflector, 40-nm-thick Si3N4 dielectric spacer, and 120-nm-thick top TiN patterned layer in a 280 × 280 nm2 unit cell. The top TiN layer is set to be the optimization region (outlined by red box). (b) Convergence plot of the topology optimization with corresponding evolution of the material distribution. Black and red curves correspond to the in-band and out-of-band average absorption respectively. (c) Designs obtained by the direct topology optimization, including the design with the highest efficiency (92%) framed in the red box (black color corresponds to air, white to TiN). (d) Statistics of efficiency distribution of 200 topology optimized designs. Normalized efficiency is defined as a ratio of in-band emissions of the TO thermal emitter and of an ideal emitter at 1800 °C. (e) The absorption/emissivity spectrum of the best TO (black) and cylindrical (red) thermal emitters. (f) Corresponding emissivity spectra at 1800 °C: the best TO emitter (black), cylindrical emitter (red), and blackbody emission (dashed black curve).

Close modal

For the TO process, the figure of merit (FOM) is defined as the spectrally weighted average of the in-band absorption and the out-of-band reflectivity. The weighting of the FOM spectrum is done with respect to the absorption and reflection spectrum of the ideal emitter (supplementary material, Sec. S2). The FOM corresponds to the case of maximized in-band and minimized out-of-band absorption values, which in the ideal case should converge to the step-function type absorption spectrum [Fig. 2(b)]. Figure 3(b) shows the convergence plot of the TO FOM, illustrating the material evolution inside the unit cell. The spikes of the FOM convergence plot occur due to the applied filtering algorithm. Some of the generated 200 designs are shown in Fig. 3(c), where the implemented filtering algorithm cuts off all sub-30-nm features. Here we only use the designs with: (i) more than 85% of average in-band and (ii) less than 40% of mean out-of-band absorption.

To be able to compare the performance of the generated designs, we have defined the efficiency of the thermal emitter as a product of in-band (effin) and out-of-band (effout) efficiencies. effin is an in-band radiance of the emitter normalized to the in-band radiance of ideal emitter at 1800 °C, while out-of-band efficiency effout is defined as a ratio of the out-of-band radiance of back reflector and radiance of the TO design. The later reflects the fact that the response of the gap plasmon structures in the long-wavelength limit is fully determined by the material properties of the back reflector and the out-of-band emissivity is fundamentally limited by the optical losses of TiN. More details on the efficiency calculation can be found in supplementary material, Sec. S3. Figure 3(d) depicts the statistics of the design efficiencies. The normalized efficiency of the cylindrical emitter is 83%, while the best TO design has 92% efficiency. Figure 3(e) shows the absorption/emissivity spectra of the best TO (black) and cylindrical (red) designs. The obtained TO emitter designs with non-trivial shapes enable more uniform, higher in-band absorption while providing the same rapidly decaying tail in the out-of-band region. We note that the plasmonic back-reflector entirely defines the behavior of the out-of-band tail of the absorption spectrum at longer wavelengths. Figure 3(f) shows the corresponding emission spectra of the best TO thermal emitter (black), cylindrical (red) vs the blackbody emission spectrum at 1800 °C. A dense population of in-band modes in the TO design enables a higher in-band emission that in the ideal case should match the blackbody in-band emission. The out-of-band emission is substantially suppressed in both design cases and is equal to 32% (TO design) and 29% (cylindrical emitter) of the out-of-band blackbody radiation at 1800 °C. The performance of the thermal emitter can be further optimized by adjusting optical properties of the back-reflector material (i.e., increasing the reflectivity and decreasing attenuation at long wavelength range).

Once the AAE network is trained on the obtained TO designs, the encoder takes a 64 × 64 binary, greyscale image/pattern of the TO design as input and compresses it into the 15-dimensional vector, representing a position in the 15-dimensional latent space. The decoder reconstructs the resonant layer design from the latent space by taking the latent coordinates as input. During the training process, the discriminator forces the encoder to form the latent space that would match the pre-defined distribution. In the initial case, we use Gaussian distribution as a pre-defined model.

Common to all deep neural networks, for efficient AAE training, a much larger designs dataset (∼10k) is required, compared with the available set of 200 designs. Since TO is time-consuming, the direct generation of thousands of designs is not practical. To overcome this issue, we used data augmentation, which takes into account the physical and symmetry properties of metasurfaces. Periodicity of the thermal emitter design allows for translational perturbation along a longitudinal direction without affecting the optical response of the structure. The freedom in cross-polarization selection allows for a 90-degree rotation. We use 20 random lateral translations of the original patterns and a 90-degree rotated pattern per design that allowed us to enlarge the training set to 8400 samples. The AAE network is then trained on the expanded dataset. More details on the training process and specifics of the AAE structure are given in supplementary material (Sec. S4).

After the training phase, the decoder is used as a generator of new designs. The generated designs are refined/filtered with two different approaches (Fig. 4): (i) using additional iterations of TO; (ii) using a pre-trained convolutional neural network for predicting the efficiency and robustness of generated designs.

FIG. 4.

Structure refinement and filtering schemes. The trained AAE generator is used for thermal emitter design production. To be able to perform the refinement of the obtained designs, two different approaches are developed: (i) exploiting additional TO based refinement, which includes an additional 15 iterations of the TO process; (ii) the filtering of robust, high efficiency designs using the pre-trained VGGnet type of convolutional neural network (CNN). The VGGnet based filtering included additional Gaussian filtering, for elimination of small and sharp features from the AAE generated designs.

FIG. 4.

Structure refinement and filtering schemes. The trained AAE generator is used for thermal emitter design production. To be able to perform the refinement of the obtained designs, two different approaches are developed: (i) exploiting additional TO based refinement, which includes an additional 15 iterations of the TO process; (ii) the filtering of robust, high efficiency designs using the pre-trained VGGnet type of convolutional neural network (CNN). The VGGnet based filtering included additional Gaussian filtering, for elimination of small and sharp features from the AAE generated designs.

Close modal

The structure refinement with the additional TO process ensures the stability of the final designs, as well as helps with eliminating sub-30 nm features. Using the decoder, we generate 1000 designs with a condition of having at least a 40% TiN filling factor within the domain. This requirement helps to prevent low-efficiency of the designs after the refinement due to a low plasmonic material fraction in the domain. The refinement uses the same constraints on the final design efficiency as the direct TO, that is, (i) at least 85% of mean in-band absorption and (ii) less than 40% of the out-of-band [Fig. 5(a)]. Figure 5(a) indicates that the refinement removes the material “blurring,” transforming the design into a binary structure. The AAE + TO optimization provides a mean efficiency of 90% for the set of 200 designs vs 82% for the set obtained with the direct TO [Fig. 5(b)]. The best AAE + TO designs [red box in Fig. 5(a)] provide the top efficiency of 98%, while the efficiency of the best pattern from the direct TO set is only 92%. The AAE + TO design exhibits almost unit-level in-band absorption while having the same decaying tail of the out-of-band absorption spectrum as the direct TO and trivial emitters [Fig. 5(c)]. The AAE + TO designs enable almost all available in-band emission (98%) while significantly suppressing out-of-band radiation [30% of out-of-band blackbody radiation at 1800 °C; Fig. 5(d)].

FIG. 5.

AAE for design optimization of thermal emitters. (a) Example of designs generated by the trained AAE (left panel); the same designs after the structure refinement process (right panel). (b) Statistics of the efficiency distribution for 200 designs obtained with AAE + TO (blue bars). The same statistics of the 200-design set obtained with the direct TO (gray bars). (c) The absorption/emissivity spectra of the best AAE + TO design in the set (blue curve), the best direct TO design (black), and the optimized cylindrical emitter (red). The inset shows the unit cell configuration of the best design in the set. (d) The corresponding emission spectra of the heaters at 1800 °C temperature. Black dashed curve corresponds to blackbody emission at 1800 °C.

FIG. 5.

AAE for design optimization of thermal emitters. (a) Example of designs generated by the trained AAE (left panel); the same designs after the structure refinement process (right panel). (b) Statistics of the efficiency distribution for 200 designs obtained with AAE + TO (blue bars). The same statistics of the 200-design set obtained with the direct TO (gray bars). (c) The absorption/emissivity spectra of the best AAE + TO design in the set (blue curve), the best direct TO design (black), and the optimized cylindrical emitter (red). The inset shows the unit cell configuration of the best design in the set. (d) The corresponding emission spectra of the heaters at 1800 °C temperature. Black dashed curve corresponds to blackbody emission at 1800 °C.

Close modal

Within the AEE + TO refinement optimization approach, most of the computational time is spent on structure refinement due to additional TO iterations. As an alternative, here we propose to use pre-trained CNN based filtering of highly efficient and robust designs within an AAE generated design set. The CNN based structure filtering process consists of two main steps: (i) applying Gaussian filtering and binarization function to an AAE generated design for the elimination sub-30 nm feature (see supplementary material, Sec. S2) and (ii) robustness and efficiency estimation based on pre-trained CNN. Specifically, the CNN architecture utilized here is a smaller version of the VGGnet network, introduced in Ref. 60 [Fig. 6(a)]. CNN takes a 64 × 64 image of the design as an input and passes it through three hidden layers, which consist of convolutional layers with ReLU activation functions. Each hidden layer is followed by the maximum pooling layer, which ensures the down-sampling of the feature maps. The stack of convolutional layers is followed by one fully connected layer. The base VGGnet architecture is followed by two different activation functions of the final layer: (i) the “soft-max” with “cross-entropy” loss function for the robustness classification (“robust” or “not robust”); (ii) “linear” activation function with “mean squared error” loss function for efficiency prediction (regression). More details on the smaller VGGnet network can be found in the supplementary material.

FIG. 6.

AAE + VGGnet optimization. (a) Schematics of the smaller VGGnet used for design efficiency prediction and robustness classification. (b) Regression map of the VGGnet performance on the training set. (c) Test results of the VGGnet based robustness classification scheme. (d) Statistics of the efficiency distribution for 1000 designs obtained directly from the AAE network (gray bars) and AAE coupled with VGGnet filtering (blue bars). (e) The absorption/emissivity spectra of the best AAE + VGGnet design in the set (blue curve), the best direct TO design (black), and the optimized cylindrical emitter (red). The inset shows the unit cell configuration of the best design in the set. (f) The corresponding emission spectra of the heaters at 1800 °C temperature. Black dashed curve corresponds to blackbody emission at 1800 °C.

FIG. 6.

AAE + VGGnet optimization. (a) Schematics of the smaller VGGnet used for design efficiency prediction and robustness classification. (b) Regression map of the VGGnet performance on the training set. (c) Test results of the VGGnet based robustness classification scheme. (d) Statistics of the efficiency distribution for 1000 designs obtained directly from the AAE network (gray bars) and AAE coupled with VGGnet filtering (blue bars). (e) The absorption/emissivity spectra of the best AAE + VGGnet design in the set (blue curve), the best direct TO design (black), and the optimized cylindrical emitter (red). The inset shows the unit cell configuration of the best design in the set. (f) The corresponding emission spectra of the heaters at 1800 °C temperature. Black dashed curve corresponds to blackbody emission at 1800 °C.

Close modal

The VGGnet has been trained on 5000 designs generated by AAE. The efficiency and robustness ground truth labels have been assessed via FDTD simulation of each AAE generated design. Seventy percent and 30% of the dataset have been used for training and testing, respectively. Figure 6(b) shows the performance of the trained VGGnet regression model. The heatmap shows the statistics of the efficiency value prediction as a function of true labels. Here we can see that the VGGnet is able to predict efficiency with high accuracy. The coefficient of determination (r2 coefficient) is equal to 0.93, which in the ideal case should be equal to 1. To be able to set up the robustness classification problem, we have labeled the design as “robust” if F=max(|effidealefferoded|/effideal,|effidealeffdilated|/effideal)0.95, and as “not robust” if F<0.95, where effj are efficiency values for ideal, eroded, and dilated designs. We have applied ±10nm perturbation of the ideal structure. Figure 6(c) shows the test results for the robustness classification. Here we can see that VGGnet ensures 80% accuracy for “robust” devices classification and 79% for “not robust.”

Once both variations of VGGnet are trained, we have generated 1000 designs using the AAE + VGGnet approach with the two main constraints: high efficiencies (>80%) and robustness. Figure 6(d) shows the statistics of the efficiency distribution of the generated AAE + VGGnet designs (blue bars) and 1000 designs generated directly from AAE (gray bars). For the elimination of the sub-30 nm features, AAE designs have been passed through the Gaussian filter. As we can see, the design set generated directly from an AAE has almost the same distribution as a TO training set, while AAE + VGGnet ensures high efficiency of the generated set due to additional filtering. Almost 88% of the AAE + VGGnet generated designs have efficiency over 80%. The best AAE + VGGnet design has 95.5% efficiency, while the best design within the AAE set has 94.6%. The absorption spectrum, as well as emissivity spectra of the best AAE + VGGnet generated design, are shown in Figs. 6(d)–6(e). For the sake of comparison, spectra of the best TO as well as a cylindrical emitter designs are shown as well.

Along with the higher efficiency, AAE assisted optimization ensures a faster optimization search in comparison with conventional TO. Figure 7(a) shows the comparison between computational costs of direct TO (black), AAE + TO refinement (blue), and AAE + VGGnet (red) approaches. Here we can see that, for the generation of 100 designs, direct TO requires 164 h, while the AAE + TO refinement approach needs 54 h. The direct TO needs 1.64 h per design optimization on average, while the AAE based optimization needs only 31 min. The AAE optimization time consists of the decoder generation time (<1 s per design) and the refinement time (∼31 min). In comparison, the AAE + VGGnet approach requires only 2 min for the generation of 100 highly efficient (>80%) designs, which is 1620 times faster than the AAE + TO approach and more than 4900 times faster than direct TO. The training of the AAE network on topology optimized designs takes 24 min, while the training of VGGnet for robustness classification and efficiency regression takes 45 and 53 min, respectively. All numerical simulations are done on a cluster node with two 12-core Intel Xeon Gold “Sky Lake” processors (24 cores per node) and 96 GB of memory. Direct full-wave simulation at each iteration is done in parallel, while the filtering, calculation of gradients, and material distribution updates are performed in a sequential manner. Figure 7(b) depicts the results of all used optimization methods within this work. Here we would like to outline the main difference between two AAE assisted optimization approaches. The AAE optimization followed by additional TO refinement ensures the best solution and almost maximum possible result of the problem under consideration. However, a TO based refinement procedure substantially limits the proposed approach in terms of optimization speedup, while the AAE + VGGnet approach ensures tremendous speedup of the optimization search and relatively high overall device efficiency.

FIG. 7.

Optimization search efficiencies. (a) Dependence of the computational time of the direct TO (black), AAE + TO based optimization (blue), and AAE + VGGnet (cyan) on the number of the optimized high-efficiency resonant patterns. For the generation of 100 designs, direct TO requires 164 h, the AAE+TO refinement approach needs 54 h, while the AAE + VGGnet approach requires only 2 min for the generation of highly efficient (>80%) 100 designs. The AAE + VGGnet based approach is 1620 times faster than the AAE + TO approach and more than 4900 times faster than direct TO. (b) Comparison of the efficiencies of the obtained best designs for all used optimization methods within this work. The AAE optimization followed by additional TO refinement ensures the best solution and almost maximum possible result of the problem under consideration. While the AAE + VGGnet approach ensures tremendous speedup of the optimization search and relatively high overall device efficiency.

FIG. 7.

Optimization search efficiencies. (a) Dependence of the computational time of the direct TO (black), AAE + TO based optimization (blue), and AAE + VGGnet (cyan) on the number of the optimized high-efficiency resonant patterns. For the generation of 100 designs, direct TO requires 164 h, the AAE+TO refinement approach needs 54 h, while the AAE + VGGnet approach requires only 2 min for the generation of highly efficient (>80%) 100 designs. The AAE + VGGnet based approach is 1620 times faster than the AAE + TO approach and more than 4900 times faster than direct TO. (b) Comparison of the efficiencies of the obtained best designs for all used optimization methods within this work. The AAE optimization followed by additional TO refinement ensures the best solution and almost maximum possible result of the problem under consideration. While the AAE + VGGnet approach ensures tremendous speedup of the optimization search and relatively high overall device efficiency.

Close modal

The structure of the AAE network ensures unparalleled control over latent space distribution by adjusting the pre-defined model during the training phase. Such control can be adapted for the realization of global optimization techniques directly inside the compressed space. One possibility of realizing such global optimization is mapping the latent space distribution into the pre-defined surrogate model and using Bayesian optimization. Additionally, such control over the latent space could be used for the determination of sub-latent space with a lower dimension, corresponding to the best designs in the set by using principal component analysis. The analysis of such lower-dimensional spaces will allow determining the key requirement for achieving the best possible performance of the photonic/plasmonic device within a given design space. Within the aforementioned global optimization frameworks, the control over the configuration of the latent space is critical. Moreover, such latent space control is of significant importance for the multi-constrained problems, which require careful engineering of the compressed hyper-dimensional design spaces for the realization of efficient optimization searches.

To showcase such control, we train the AAE network on the TO design set along with probing it with the two types of pre-defined distributions: “3-Gaussian mixture” [Fig. 8(a)] and “Swiss roll” [Fig. 8(b)]. “3-Gaussian mixture” is a mapping of three 2D Gaussian distributions into a parametric 2D plane (x̃,ỹ):

x̃=xcos(2π3ξ)ysin(2π3ξ),ỹ=xsin(2π3ξ)+ycos(2π3ξ),
(1)

where x,y are random real numbers with the Gaussian distribution, ξ is a random integer with “discrete uniform” distribution from a list ξ{1,2,3}. In contrast the “Swiss roll” mapping projects a single random real number with the Gaussian distribution (x) into parametric 2D plane (x̃,ỹ):

x̃=8x+ξcos(8πx+ξ),ỹ=8x+ξsin(8πx+ξ),
(2)

where the random integer ξ is discussed.

FIG. 8.

Latent space engineering. The structure of the AAE network ensures unparalleled control over latent space distribution by adjusting the pre-defined model during the training phase. Such control can be adapted for the realization of global optimization techniques directly inside the compressed space. To showcase such control, we train the AAE network on the TO design set along with probing it with the two types of pre-defined distributions: (a) 2D latent space distribution with a “3 Gaussian mixture” sampling of the pre-defined model; (b) the same as in (a) but for “swiss roll” distribution. Insets show corresponding data sampling used as a pre-defined model.

FIG. 8.

Latent space engineering. The structure of the AAE network ensures unparalleled control over latent space distribution by adjusting the pre-defined model during the training phase. Such control can be adapted for the realization of global optimization techniques directly inside the compressed space. To showcase such control, we train the AAE network on the TO design set along with probing it with the two types of pre-defined distributions: (a) 2D latent space distribution with a “3 Gaussian mixture” sampling of the pre-defined model; (b) the same as in (a) but for “swiss roll” distribution. Insets show corresponding data sampling used as a pre-defined model.

Close modal

The AAE is trained to compress/reconstruct 64 × 64 input patterns into/from 2D latent space. A “3 Gaussian mixture” sampling indicates the formation of the three-lobe distribution of the design space, enforced by the pre-defined model distribution [Fig.8(a)]. In a more sophisticated case, the pre-defined model is a complex spiral-shape distribution that is also reconstructed by the AAE accordingly [Fig. 8(b)].

The synergy between the inverse design methods and advanced machine learning techniques opens up a new paradigm to address highly complex, multi-constrained problems. Here, we merge the adjoint-based topology optimization with the AAE network and demonstrate faster optimization searches and unparalleled control over the latent space configuration. The latter is crucial for the realization of efficient optimization over high dimensional parametric landscapes, which is required for the design of multiconstrained, multifunctional photonic devices. Specifically, we optimize the design of a thermal emitter metasurface with high-efficiency thermal emission reshaping. We show that AAE + TO optimization-based emitter designs enable thermal reshaping with efficiencies up to 98%. Along with the better efficiency, the proposed AAE + VGGnet approach demonstrates ∼4900 times faster optimization searches in comparison with a conventional direct TO method. The proposed method can transform the area of optical design as well as data-driven materials synthesis for a plethora of applications in photonics, optoelectronics, MEMS, and biomedical synthetics.

See the supplementary material for details on dielectric permittivity functions of TiN and Si3N4, additional information on topology optimization framework, adversarial autoencoder structure, training process, and data augmentation procedure; additional data on AAE based optimized designs; and structure of the VGGnet.

This work was supported in part by the AFOSR (Grant No. FA9550-17-1–0243) and the award on trans-dimensional photonics, the NSF (Grant No. 0939370-CCF), and Defense Advanced Research Projects Agency/Defense Sciences Office Extreme Optics and Imaging (EXTREME) Program (Grant No. DARPA/DSO HR00111720032).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
S.
Molesky
,
Z.
Lin
,
A. Y.
Piggott
,
W.
Jin
,
J.
Vucković
, and
A. W.
Rodriguez
,
Nat. Photonics
12
,
659
(
2018
).
2.
M. M.
Spuhler
,
B. J.
Offrein
,
G. L.
Bona
,
R.
Germann
,
I.
Massarek
, and
D.
Erni
,
J. Light. Technol.
16
,
1680
(
1998
).
3.
A.
Gondarenko
and
M.
Lipson
,
Opt. Express
16
,
17689
(
2008
).
4.
A. V.
Kildishev
,
U. K.
Chettiar
,
Z.
Liu
,
V. M.
Shalaev
,
D.-H.
Kwon
,
Z.
Bayraktar
, and
D. H.
Werner
,
J. Opt. Soc. Am. B
24
,
A34
(
2007
).
5.
F.
Callewaert
,
V.
Velev
,
P.
Kumar
,
A. V.
Sahakian
, and
K.
Aydin
,
Sci. Rep.
8
,
1358
(
2018
).
6.
D.
Sell
,
J.
Yang
,
S.
Doshay
,
R.
Yang
, and
J. A.
Fan
,
Nano Lett.
17
,
3752
(
2017
).
7.
B.
Shen
,
P.
Wang
,
R.
Polson
, and
R.
Menon
,
Optica
1
,
356
(
2014
).
8.
T. P.
Xiao
,
O. S.
Cifci
,
S.
Bhargava
,
H.
Chen
,
T.
Gissibl
,
W.
Zhou
,
H.
Giessen
,
K. C.
Toussaint
,
E.
Yablonovitch
, and
P. V.
Braun
,
ACS Photonics
3
,
886
(
2016
).
9.
J. S.
Jensen
and
O.
Sigmund
,
Laser Photonics Rev.
5
,
308
(
2011
).
10.
P.
Seliger
,
M.
Mahvash
,
C.
Wang
, and
A. F. J.
Levi
,
J. Appl. Phys.
100
,
034310
(
2006
).
11.
C. Y.
Kao
,
S.
Osher
, and
E.
Yablonovitch
,
Appl. Phys. B
81
,
235
(
2005
).
12.
Y.
Elesin
,
B. S.
Lazarov
,
J. S.
Jensen
, and
O.
Sigmund
,
Photonics Nanostructures - Fundam. Appl.
10
,
153
(
2012
).
13.
L.
Martinelli
and
A.
Jameson
,
J. Heat Transfer
135
,
011002
(
2013
).
14.
M. P.
Bendsøe
and
O.
Sigmund
,
Arch. Appl. Mech.
69
,
635
(
1999
).
15.
L.
Yang
,
A. V.
Lavrinenko
,
J. M.
Hvam
, and
O.
Sigmund
,
Appl. Phys. Lett.
95
,
261101
(
2009
).
16.
P. I.
Borel
,
A.
Harpøth
,
L. H.
Frandsen
,
M.
Kristensen
,
P.
Shi
,
J. S.
Jensen
, and
O.
Sigmund
,
Opt. Express
12
,
1996
(
2004
).
17.
J. S.
Jensen
and
O.
Sigmund
,
Appl. Phys. Lett.
84
,
2022
(
2004
).
18.
Y.
Tsuji
and
K.
Hirayama
,
IEEE Photonics Technol. Lett.
20
,
982
(
2008
).
19.
L. F.
Frellsen
,
Y.
Ding
,
O.
Sigmund
, and
L. H.
Frandsen
,
Opt. Express
24
,
16866
(
2016
).
20.
H.
Men
,
K. Y. K.
Lee
,
R. M.
Freund
,
J.
Peraire
, and
S. G.
Johnson
,
Opt. Express
22
,
22632
(
2014
).
21.
Y.
Shen
,
N. C.
Harris
,
S.
Skirlo
,
M.
Prabhu
,
T.
Baehr-Jones
,
M.
Hochberg
,
X.
Sun
,
S.
Zhao
,
H.
Larochelle
,
D.
Englund
, and
M.
Soljačić
,
Nat. Photonics
11
,
441
(
2017
).
22.
R.
Hamerly
,
L.
Bernstein
,
A.
Sludds
,
M.
Soljačić
, and
D.
Englund
,
Phys. Rev. X
9
,
021032
(
2019
).
23.
T. W.
Hughes
,
M.
Minkov
,
Y.
Shi
, and
S.
Fan
,
Optica
5
,
864
(
2018
).
24.
Y.
Zuo
,
B.
Li
,
Y.
Zhao
,
Y.
Jiang
,
Y.-C.
Chen
,
P.
Chen
,
G.-B.
Jo
,
J.
Liu
, and
S.
Du
,
Optica
6
,
1132
(
2019
).
25.
X.
Lin
,
Y.
Rivenson
,
N. T.
Yardimci
,
M.
Veli
,
Y.
Luo
,
M.
Jarrahi
, and
A.
Ozcan
,
Science
361
,
1004
(
2018
).
26.
W.
Ma
,
F.
Cheng
, and
Y.
Liu
,
ACS Nano
12
,
6326
(
2018
).
27.
J.
Peurifoy
,
Y.
Shen
,
L.
Jing
,
Y.
Yang
,
F.
Cano-Renteria
,
B. G.
DeLacy
,
J. D.
Joannopoulos
,
M.
Tegmark
, and
M.
Soljačić
,
Sci. Adv.
4
,
eaar4206
(
2018
).
28.
I.
Malkiel
,
M.
Mrejen
,
A.
Nagler
,
U.
Arieli
,
L.
Wolf
, and
H.
Suchowski
,
Light Sci. Appl.
7
,
60
(
2018
).
29.
K.
Yao
,
R.
Unni
, and
Y.
Zheng
,
Nanophotonics
8
,
339
(
2019
).
30.
Z.
Liu
,
D.
Zhu
,
S. P.
Rodrigues
,
K. T.
Lee
, and
W.
Cai
,
Nano Lett.
18
,
6570
(
2018
).
31.
W.
Ma
,
F.
Cheng
,
Y.
Xu
,
Q.
Wen
, and
Y.
Liu
,
Adv. Mater.
31
,
1901111
(
2019
).
32.
J.
Jiang
,
D.
Sell
,
S.
Hoyer
,
J.
Hickey
,
J.
Yang
, and
J. A.
Fan
,
ACS Nano
13
,
8872
(
2019
).
33.
J.
Jiang
and
J. A.
Fan
,
Nano Lett.
19
,
5366
(
2019
).
34.
A.
Makhzani
,
J.
Shlens
,
N.
Jaitly
,
I.
Goodfellow
, and
B.
Frey
, arXiv:1511.05644 (
2015
).
35.
Y.
Pu
,
Z.
Gan
,
R.
Henao
,
X.
Yuan
,
C.
Li
,
A.
Stevens
, and
L.
Carin
,
Adv. Neural Inf. Process. Syst.
2352
(
2016
).
36.
A.
Creswell
,
T.
White
,
V.
Dumoulin
,
K.
Arulkumaran
,
B.
Sengupta
, and
A. A.
Bharath
,
IEEE Signal Process. Mag.
35
,
53
(
2018
).
37.
I. J.
Goodfellow
,
J.
Pouget-Abadie
,
M.
Mirza
,
Bing
Xu
,
D.
Warde-Farley
,
S.
Ozair
,
A.
Courville
, and
Y.
Bengio
, in
Generative Adversarial Networks, NIPS'14: Proceedings of the 27th International Conference on Neural Information Processing Systems
(
NIPS
,
2014
), Vol. 2, pp.
2672
2680
.
38.
N.
Liu
,
M.
Mesch
,
T.
Weiss
,
M.
Hentschel
, and
H.
Giessen
,
Nano Lett.
10
,
2342
(
2010
).
39.
H. A.
Parsamyan
,
K. V.
Nerkararyan
, and
S. I.
Bozhevolnyi
,
J. Opt. Soc. Am. B
36
,
2643
(
2019
).
40.
G.
Liu
,
X.
Liu
,
J.
Chen
,
Y.
Li
,
L.
Shi
,
G.
Fu
, and
Z.
Liu
,
Sol. Energy Mater. Sol. Cells
190
,
20
(
2019
).
41.
D. G.
Baranov
,
Y.
Xiao
,
I. A.
Nechepurenko
,
A.
Krasnok
,
A.
Alù
, and
M. A.
Kats
,
Nat. Mater.
18
,
920
(
2019
).
42.
A. S.
Roberts
,
M.
Chirumamilla
,
D.
Wang
,
L.
An
,
K.
Pedersen
,
N. A.
Mortensen
, and
S. I.
Bozhevolnyi
,
Opt. Mater. Express
8
,
3717
(
2018
).
43.
W.
Li
,
U.
Guler
,
N.
Kinsey
,
G. V.
Naik
,
A.
Boltasseva
,
J.
Guan
,
V. M.
Shalaev
, and
A. V.
Kildishev
,
Adv. Mater.
26
,
7921
(
2014
).
44.
H.
Reddy
,
U.
Guler
,
Z.
Kudyshev
,
A. V.
Kildishev
,
V. M.
Shalaev
, and
A.
Boltasseva
,
ACS Photonics
4
,
1413
(
2017
).
45.
B.
Bitnar
,
W.
Durisch
,
J. C.
Mayor
,
H.
Sigg
, and
H. R.
Tschudi
,
Sol. Energy Mater. Sol. Cells
73
,
221
(
2002
).
46.
Y. X.
Yeng
,
M.
Ghebrebrhan
,
P.
Bermel
,
W. R.
Chan
,
J. D.
Joannopoulos
,
M.
Soljacic
, and
I.
Celanovic
,
Proc. Natl. Acad. Sci.
109
,
2280
(
2012
).
47.
S. Y.
Lin
,
J. G.
Fleming
,
D. L.
Hetherington
,
B. K.
Smith
,
R.
Biswas
,
K. M.
Ho
,
M. M.
Sigalas
,
W.
Zubrzycki
,
S. R.
Kurtz
, and
J.
Bur
,
Nature
394
,
251
(
1998
).
48.
I.
Celanovic
,
N.
Jovanovic
, and
J.
Kassakian
,
Appl. Phys. Lett.
92
,
193101
(
2008
).
49.
A. V.
Kildishev
,
A.
Boltasseva
, and
V. M.
Shalaev
,
Science
399
,
1232009
(
2013
).
50.
W.
Li
,
U.
Guler
,
N.
Kinsey
,
G. V.
Naik
,
A.
Boltasseva
,
J.
Guan
,
V. M.
Shalaev
, and
A. V.
Kildishev
,
Adv. Mater.
26
,
7959
(
2014
).
51.
M. G.
Nielsen
,
A.
Pors
,
O.
Albrektsen
, and
S. I.
Bozhevolnyi
,
Opt. Express
20
,
13311
(
2012
).
52.
F.
Ding
,
Y.
Yang
,
R. A.
Deshpande
, and
S. I.
Bozhevolnyi
,
Nanophotonics
7
,
1129
(
2018
).
53.
J. T.
Heiden
,
F.
Ding
,
J.
Linnet
,
Y.
Yang
,
J.
Beermann
, and
S. I.
Bozhevolnyi
,
Adv. Opt. Mater.
7
,
1801414
(
2019
).
54.
W.-Y.
Tsai
,
C.-M.
Wang
,
C.-F.
Chen
,
P. C.
Wu
,
Y.-H.
Chen
,
T.-Y.
Chen
,
P. R.
Wu
,
J.-W.
Chen
, and
D. P.
Tsai
,
Sci. Rep.
7
,
42076
(
2017
).
55.
X.
Ni
,
N. K.
Emani
,
A. V.
Kildishev
,
A.
Boltasseva
, and
V. M.
Shalaev
,
Science
335
,
427
(
2012
).
56.
W.
Sabra
,
S. I.
Azzam
,
M.
Song
,
M.
Povolotskyi
,
A. H.
Aly
, and
A. V.
Kildishev
,
Opt. Lett.
43
,
4815
(
2018
).
57.
R.
Poli
,
J.
Kennedy
, and
T.
Blackwell
,
Swarm Intell.
1
,
33
(
2007
).
58.
R. E.
Christiansen
,
J.
Vester-Petersen
,
S. P.
Madsen
, and
O.
Sigmund
,
Comput. Methods Appl. Mech. Eng.
343
,
23
(
2019
).
59.
F.
Wang
,
J. S.
Jensen
, and
O.
Sigmund
,
J. Opt. Soc. Am. B
28
,
387
(
2011
).
60.
K.
Simonyan
and
A.
Zisserman
, in
Proceedings of the Third International Conference on Learning Representations, ICLR 2015, Conference Track Proceedings
(
2015
).

Supplementary Material