1. Introduction

Electromagnetic metamaterials have been used to realize many exotic scattering responses over the last two decades. Effects including negative refractive index and cloaking have generated significant interest and have served to drive the community [1–5]. A more applied, but still relevant metamaterials achievement, is that of graded designs [6,7]. It was shown early on that unit-cells which form metamaterials may be designed with a spatial dependence across a surface and / or volume, and various lensing effects were shown that utilize spatial degrees of freedom. A principle of gradient metasurfaces is that their scattering properties change slowly as a function of spatial coordinate. Broadband metamaterial absorbers [8,9], and metamaterial spatial light modulators [10], also make use of dissimilar neighboring unit-cells, however there is no such requirement to make neighbors as alike as possible, and designs are simply cobbled together to achieve the desired response. This fact highlights a design feature of sub-wavelength metal-based metamaterials, i.e. their scattering response is primarily connected to their geometry and – due to minimal neighbor interaction – the unit-cell largely governs the electromagnetic properties of the array. Conventional optimization techniques built-in to modern day electromagnetic mode solvers are sufficient to achieve a desired response, and designs of differing units may be assembled with little change to the overall response.

In principle, all-dielectric metasurface (ADM) unit-cells may also be used to tesselate a surface in an arbitrary fashion similar to that achieved with metal-based metamaterials. However, ADMs resonators are only marginally sub-wavelength, and modes utilized are often not confined within their physical bound – with an evanescent tail lying just outside their surface – resulting in significant neighbor interaction. Further, surface modes related to the periodicity are not restricted to nearest neighbors, and may persist over several spatial periods. Despite the success of ADMs [11,12], it is conceivable that still-richer electromagnetic scattering can be achieved if more complex geometries are employed. However, physical understanding for such metasurfaces is poor – simple functional relationships, or even heuristic guidance regarding super unit-cell geometry and final electromagnetic properties is unavailable. The only contemporary means to estimate such metasurface properties, given a candidate geometry, is simulation or fabrication. However, given the vast space of potential designs and the speed of conventional simulation and fabrication, it is completely infeasible to iterate over designs in order to achieve a desired response.

A viable design alternative to numerical simulation for artificial electromagnetic materials (AEMs), including metamaterials [13–17], photonic band gap [18–21], and plasmonics studies [22–24], is deep learning. Deep neural networks (DNNs) have successfully learned a forward mapping $s=f(g)$ between a metasurface geometry $g$ and the resulting frequency dependent electromagnetic scattering $s$, where $f$ is an unknown complex (e.g., highly nonlinear) function [25]. A DNN – once trained – can effectively act as a high-speed simulator that may be used to find the electromagnetic scattering of a candidate geometry substantially faster (e.g., by six orders of magnitude [16]) than conventional simulation.

While deep learning enables substantially faster evaluation of candidates, given the vast number of possibilities for many problems (e.g., approx. $10^{12}$ in our case), it is nonetheless difficult or impossible to iterate over all or most candidate designs. From a design perspective, what would be of greatest utility would be to instead specify a desired frequency dependent electromagnetic response $s$, and have a model or solver compute a specific approximate geometry $\hat {g}$, which yields $s$ (here $\hat {}$ denotes approximate). The solution used to search for such an ideal geometry may be cast as an inverse problem written as $\hat {g}=\hat {f}^{-1}(s)$. Hadamard described a solution method as well-posed if it met three criteria – $H_1$: existence (at least one solution), $H_2$: uniqueness (only one solution), and $H_3$: stability (solution must depend continuously on $s$). A solution approach is then ill-posed if one of the $H_1-H_3$ conditions fails [26]. Past work has shown that finding a specific ADM consisting of distinct neighboring resonators in a super unit-cell geometry for a desired scattering state is an ill-posed non-linear inverse problem and, in particular, conditions $H_1$ and $H_2$ are not met [13].

2. Results

Here we explore a deep inverse modeling approach, termed the neural-adjoint (NA) method, [27] which outperformed many other methods for solving ill-posed inverse problems. We take on the challenging task of finding the ADM geometry which yields a desired spectra. Notably, we chose a design space with 14 free parameters – significantly larger than most recent works, i.e. $5-10$ free parameters in [14,16,18]. Furthermore, existing inverse methods often involve complex training procedures, and ultimately produce sub-optimal solutions. In contrast, the NA only requires that we train a single conventional feed-forward neural network, and – as we show – appears to find close approximations to the globally optimal solution within just one minute, even for our complex ADM design problem. Furthermore, we demonstrate how the NA method can be utilized to expand the design search space, essentially providing a form of active learning that is specifically tailored to solve inverse problems. We explore an example inverse problem, where our desired $s$ is a frequency dependent infrared emissivity $\mathcal {E}(\omega )$ where, from Kirchoff’s law of thermal radiation, we may write $\mathcal {E}(\omega )=A(\omega )=1-R(\omega )-T(\omega )$. Here $R=|r|^2$ and $T=|t|^2$ are the reflectance and transmittance, respectively, and $r$ and $t$ are the reflectivity coefficient and transmissivity coefficient, respectively.

We seek to obtain a match to $\mathcal {E}(\omega )$ using an ADM consisting of a square array of 2$\times$2 resonators, each with an elliptical generatrix lying perpendicular to its directrix (line of length $h$) – depicted in Fig. 1. Although past experience and single unit-cell simulations guide us to choose approximate metasurface dimensions which yield resonances in the desired spectral range, our inverse problem is ill-posed, since we do not know if a solution exists, i.e. condition $H_1$ is not met. Further, we do not meet the condition of uniqueness, $H_2$, since many metasurface solutions (disparate geometries) may yield the same spectra – at least within the accuracy of our chosen evaluation metric, and numerical precision. The proposed ADMs consists of a geometry space of fourteen parameters: [$h$, $p$, $r_{x_1}$–$r_{x_4}$, $r_{y_1}$–$r_{y_4}$, $\theta _1$–$\theta _4$]. As shown in Fig. 1, $h$ is the height of all four elliptical resonators, $p$ is the periodicity of the super cell, $r_{x_1}$ – $r_{x_4}$ and $r_{y_1}$ – $r_{y_4}$ are the x-axis and y-axis radii of four elliptical resonators respectively, and $\theta _1$ – $\theta _4$ are the rotational angles with respect to the center axis of each elliptical resonator. All geometry values are randomly sampled from the data grid shown in Table S1, which is included in Supplement 1. From Table S1 of Supplement 1, $h$, $p$, and $\theta _1$–$\theta _4$ each have five unique values from the range we have chosen. $r_{x_1}$–$r_{x_4}$ and $r_{y_1}$–$r_{y_4}$ have nine unique values for random sampling. The resulting total possible combinations is $8^9\times 6^5=1.04\times 10^{12}$. The metasurface material is simulated as lossy silicon carbide (SiC) with no substrate, and more details can be found in Supplement 1.

 

Fig. 1. (a) 3D illustration of an all-dielectric metasurface absorber. A schematic of the all-dielectric metasurfaces super-cell is shown – top right. (b) Training phase of the neural-adjoint method which uses fixed geometry inputs and corresponding simulated spectrum. (c) The inference phase of the neural-adjoint method, in which the target spectrum drives the search of the optimal geometry. The geometry is shown in a green box here to indicate it is the only trainable parameter.

Download Full Size | PDF

The DNN forward model is trained by a number of pairs consisting of 14 geometrical parameters, and 2000 frequency points from 100-500 THz. The neural network consists of twelve fully connected linear layers, four 1D transpose convolutional layers for upsampling, and one final 1D convolutional layer for spectrum smoothing. The linear fully connected layers have the following structure [14, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 500], and all hidden layers except the last linear layer are batch normalized, activated by leaky rectified linear unit (ReLU) – see Supplement 1 for more information. Results of the forward model trained on 60k $\{g,s\}$ pairs (approximately $5.75\times 10^{-8}$ of the entire geometry space) are shown in Fig. 2(a)-(c), where red curves are numerical simulations, and blue curves are the DNN predictions. The absolute error between simulation and prediction is plotted as the gray curve on the right vertical axis of Fig. 2(a)-(c). In Fig. 2(d) we show the histogram of the mean-squared-error (MSE) for the entire validation set, and find that 95% have an MSE$\leq 2.49\times 10^{-3}$ (dashed gray line), and an average MSE of $1.2\times 10^{-3}$.

 

Fig. 2. Forward model mean-squared-error (MSE) between DNN predicted spectra and numerical simulations. The DNN result with (a) the best, (b) average, and (c) below average performance. (d) Histogram detailing the distribution of MSE for the entire validation set.

Download Full Size | PDF

Having an accurate trained forward model we next turn toward the inverse problem. The NA method [27] finds the optimal inverse solution by fixing all the weights and biases of $\hat {f}$, and computing the forward model’s gradient solely with respect to the input to the network (i.e., the geometry), starting from randomly chosen values, denoted $\hat {g}_0$. It is important to highlight that $\hat {f}$ is a closed-form differentiable expression, and thus calculation of $\partial \hat {f}/\partial g$ is trivial. Further, we can estimate the gradient of the input geometry with respect to a loss function $\mathcal {L}$ that we are free to specify (e.g., mean squared error). Therefore, if $s$ is our desired spectrum, and $\hat {g}_i$ is our current best estimate of the metasurface geometry, we can iteratively move along the loss surface to find a better solution using,

Notably, while the user can specify many different choices for $\mathcal {L}$, the NA method prescribes that a so-called “boundary loss", $\mathcal {L}_{bnd}$ should be added to any user-chosen loss [27] and is given by: $\mathcal {L}_{bnd} = ReLU(|\hat {g}-\mu _g|-\frac {1}{2}R_g)$, where rectified linear unit (ReLU) is the activation function, $\mu _g$ is the mean of the geometry training data, and $R_g$ is its range. This boundary loss punishes the inference process with increasing loss if the geometry search process steps out of the space of the training data, where the forward model may produce inaccurate estimates of scattering parameters. In our experiments we use the following total loss: $\mathcal {L}= (s-\hat {f}(\hat {g}))^2 + \mathcal {L}_{bnd}$.

Well-posed inverse problems. As an initial test of the NA inverse method, we feed in frequency dependent absorptivities $A(\omega )$ where we know apriori that a solution $g$ exists, i.e. $s$ is a numerical simulation from which $\hat {f}$ originates from. In Fig. 3 we show results of the NA inverse method and each sub-plot shows characteristic results of (a) the best results, (b) average results, and (c) and (d) below average results, all based on MSE. In all of the these examples, the NA method identifies a close approximation to the correct solution. This is impressive given the complexity of the spectra present in Fig. 3. We suspect the small remaining errors in the predicted design are due largely to the limited precision of gradient descent as it nears solutions (i.e., minima) in the error space; due to the non-zero learning rate in Eq. (1), it cannot converge to the exact minimum point. We note however that learning rate can be gradually reduced during the search process, at the cost of additional computation time, until a solution of desired precision is obtained.

 

Fig. 3. Neural adjoint inverse results for prediction of target spectra that exist within the geometry space. Example spectra from geometry predictions that have (a) the best, (b) average, and (c)-(d) below average performance compared to the target spectrum.

Download Full Size | PDF

Ill-posed inverse problems. We next turn to the significantly more challenging task of applying the NA method to a spectra where we are unaware if a solution exists within the chosen geometrical parameter space, i.e. criterion $H_1$ for inverse problems is violated. We explore matching the metasurface emission to the external quantum efficiency (EQE) of gallium antimonide (GaSb), shown as the gray curves of Fig. 4, as a potential energy harvesting application using thermophotovoltaics. The metasurface will operate at elevated temperature, and thus we consider the so-called graybody spectra – also termed the spectral exitance $M_{e,\nu }(T)$ – which is given by the blackbody radiation curve times the absorptivity. We keep the top 16,000 neural adjoint solutions (spectra) and determine $M_{e,\nu }(T)$ for each of these at 100 temperatures between 1500 and 2500 k – a total of $1.6\times 10^6$ candidates.

 

Fig. 4. Neural adjoint inverse results for matching the EQE of GaSb. The target spectrum (gray), DNN prediction for $M_{e,\nu }(T=2100K)$ (blue), for the expanded geometry space (a), and the original geometry space (b) explored. Neural adjoint predictions versus $p$ and $h$ for (c) the expanded geometry space, and (d) the original space. Symbol colors indicate MSE, with values given in the colorbar.

Download Full Size | PDF

The shape of the EQE curve differs significantly from typical spectra we see in our geometry space (Fig. 3). None-the-less we find a best solution resulting from the NA method at a temperature of $T=2100k$ that consists of a geometry of [$h=0.566$, $p=1.440$, $r_{x_1}=0.180$, $r_{x_2}=0.155$, $r_{x_3}=0.214$, $r_{x_4}=0.278$, $r_{y_1}=0.285$, $r_{y_2}=0.253$, $r_{y_3}=0.146$, $r_{y_4}=0.256$, $\theta _1=-0.901^{\circ }$, $\theta _2=-20.677^{\circ }$, $\theta _3=-37.982^{\circ }$, and $\theta _4=39.046^{\circ }$]. The spectral exitance resulting from this geometry – calculated from $\hat {f}$ – is shown as the blue curve of Fig. 4(b), and we find an MSE, compared to the EQE of GaSb, of 1.06$\times$10$^{-2}$. We also apply a weighting function $W(\nu )=1\chi _{[100,275]}+0\chi _{(275,300]}$ on the MSE forcing the NA method to focus on the region of interest for energy harvesting purposes. To verify our neural adjoint results, we numerically simulate the predicted geometry and plot the resulting $M_{e,\nu }(T)$ in Fig. 5(a) as the red curve – again compared to the EQE of GaSb (gray curve). As can be seen, the simulated curve has many relatively sharp peaks that are not present in Fig. 4(b). This is because, as noted earlier, $\hat {f}$ does not perfectly match the numerical simulator, and this will introduce errors in the design process. Thus since the NA method relies on $\hat {f}$ to search for designs, it is also limited by the accuracy of the forward model estimate. We also found that since $\hat {f}$ is trained from geometries constrained to a grid, the discrepancy between NA solutions and numerical simulation arises because NA solutions are not confined to the grid, where our model is most accurate. None-the-less we find that our simulated $M_{e,\nu }(T)$ achieves an MSE of 1.65$\times 10^{-2}$, as shown in Fig. 5(a).

 

Fig. 5. Numerically simulated $M_{e,\nu }(T=2100K)$ (red) of the optimal geometry predicted by the neural-adjoint method to match the EQE of GaSb (gray), for the original geometry space (a), and the expanded geometry space (b).

Download Full Size | PDF

A fundamental challenge of ill-posed inverse problems, and a major obstacle here, is that our design search space does not contain a geometry that can realize our target spectrum (i.e., Hadamard’s criteria $H_1$). This is suggested by the fact that our best NA solution, shown in Fig. 4(a), still does not match our target spectrum. However, we can use the NA output to identify where we should expand our search space so that it will include better designs. We can do this by visualizing the error of all inverse solutions returned by NA, and look for trends e.g., we may find that all the best solutions are bunched up against some edge of our initial search space, suggesting that expanding along that dimension may yield better results. However, since we have a 14 dimensional design space, we are unable to easily visualize these data. To address this problem, we use the Uniform Manifold Approximation and Projection (UMAP) [28], which is a type of dimensionality reduction method permitting us to visualize the distribution of our inverse solutions performance in 2D, so that we may more easily identify patterns. From this investigation with UMAP in Fig. S1, shown in Supplement 1, we find that our best NA inverse solutions are limited by height. Shown in Fig. 4(d) are NA solutions as a function of height and periodicity color mapped by corresponding MSE values. It is evident that not only are our best solutions grouped at the maximum height allowed in our geometrical space, $h=0.6$ $\mu$m, but also that the solutions improve as a function of height.

Encouraged by these results, we expanded our original geometry space to include increased height values from 0.6 to 0.75 $\mu$m, by simulating an additional 24k $\{g,s\}$ pairs. The NA model now trained on the expanded geometrical space indeed finds an improved solution, shown as the blue curve in Fig. 4(a), where we realize an MSE that is reduced by a factor of 2.7. The simulated red curve in Fig. 5(b) further validates the result that the MSE of numerical simulations is also reduced – here by a factor of 4.8. A plot of the 1000 best NA solutions in the expanded geometrical space shown in Fig. 4(c), however, indicate that we may be able to make continued improvements, since we still have a gradient pushing for greater heights – although the periodicity seems to be honing in on a value of 1.2 $\mu$m.

3. Discussion and conclusions

Kirchoff’s law of thermal radiation dictates that the emissivity of a surface is equal to its absorptivity for the same frequency, direction, and polarization [29,30]. However, the metasufaces we explore here realize a degree of asymmetry, which may lead to cross-polarization (CP) scattering, thus potentially obfuscating the absorptivity and emissivity [31,32]. We performed 500 simulations (not shown) to explore CP in our metasurfaces, and found that the averaged CP transmittance, and averaged CP reflectance had values of 0.8%, and 1.1%, respectively. It is important to note that the above results represent the average of a distribution of values similar to that shown in Fig. 2(d), (and also averaged from 100-500 THz). However, for the specific metasurface identified from our neural adjoint inverse method, shown in Figs. 4(a) and 5(b), both the CP reflectance and CP transmittance, realize lower values of 0.5% and 0.8%, respectively. Thus the cross polarized transmittance and cross polarized reflectance are negligible.

The periodicity of the 2$\times$2 unit cell ranges in values from 1.0 $\mu$m to 1.5$\mu$m, and so is sub-wavelength for some of the frequency range investigated. We have performed a numerical investigation of representative metasurfaces (not shown) and have found diffraction only becomes measureable above 400 THz, but is small, with peak values of 0.7% at 400 THz and 0.5% at 500 THz. Therefore diffraction found from representative metasurfaces investigated is small, and thus negligible.

We presented the neural-adjoint inverse design method [27] which accurately predicts the high-dimensional all-dielectric metasurface geometry needed to produce a targeted infrared absorptivity spectra. When the geometry needed to produce a desired spectrum lies outside of the bounds, the NA method appears to find the best possible solution within the permissible search space. Unlike other adjoint inverse approaches [33], the NA method does not require any domain knowledge specific to the problem. The neural adjoint inverse method is general, and may be used beyond the two cases demonstrated. For example, if a specific scattering state (transmission, reflection, absorption) is desired, the neural adjoint inverse method only requires a single forward neural network model to find good solutions. In the event that the inverse solution does not exist in the parameter space explored, NA may be used to guide one to a better solution, through expansion of the design parameter search space. This may help to reduce the initial required number of numerical simulations, and to instead use NA guided simulation exploration. By tailoring the absorptivity spectrum of the ADMs, we can achieve highly efficient thermophotovoltaic cells with the application of ADMs to reduce the amount of waste energy. The NA method provides a swift and accurate approach to the inverse design of such ADMs’ emitters. With its exceptional computational speed, high accuracy, and potential use in active learning that is explored here, the neural-adjoint method has an impressive future in not only ADM thermal emitters, but also any ADM inverse problem. The NA method is not restricted to the case presented, but may be applied to many other systems including photonic crystals, artificial dielectrics, and plasmonics.