1 1. Introduction
The field of nanophotonics has been the subject of extensive expansion due to the unique capabilities of photonic nanostructures to control the propagation of EM waves. Owing to their constituent nanoscale features, which spectrally, spatially, and even temporally manipulate the optical state of the EM wave, nanophotonic devices extend all the functionalities realized by conventional optical devices in much smaller footprints. Combined with the advances in nanofabrication technologies, these nanostructures have been used to demonstrate devices with enormous potential for groundbreaking technologies addressing major challenges in stateoftheart applications, such as optical communications ^{1}, signal processing ^{2}, biosensing ^{3}, energy harvesting ^{4}, and imaging ^{5}, to name a few. As an example, newlyemerged metasurfaces (MSs) ^{6, 7, 8, 9, 10, 11, 12, 13}, twodimensional planar structures comprising of densely arranged periodic/aperiodic arrays of wellengineered dielectric or plasmonic inclusions, offer profound control of the EM wave dynamics including amplitude, phase, polarization, and frequency in the subwavelength regime^{14, 15, 16, 17}.
Despite extensive achievements in the fabrication and realization of photonic nanostructures, the efforts on the development of accurate and computationally efficient design and optimization approaches for these nanostructures are still at early stages^{18}. With the fast progress in forming more complex nanostructures with several design parameters, the need for new design approaches that can keep pace with the computational requirements for analysis and understanding of all possible design options has become more imminent. In addition, realization of nextgeneration nanodevices with potentially new physics enabled through lightmatter interaction at the nanoscale requires significant knowledge about the role of different design parameters in the functionality of a nanostructure.
Traditional design and optimization approaches for EM nanostructures rely on either using analytical (or semianalytical) modeling ^{19, 20, 21, 22, 23, 24, 18} or bruteforce analysis of the nanostructure through exhaustive search of the design parameter space ^{25}
. The use of these approaches are limited to simple structures that could be either analytically modeled or completely studied by an exhaustive search technique with reasonable computation cost. To improve the computation efficiency of such design and optimization tools, evolutionary approaches (e.g., genetic algorithm
^{26, 27} and particle swarm ^{28}) rely on starting from a random initial guess and converging to the final optimum. While reducing the computation cost compared to bruteforce approaches, such techniques are not guaranteed to converge to the global optimum of a problem (even by allocation of extensive computational resources). They are also limited to a single design problem (i.e., the simulations must be completely repeated when a small change in the nanostructure happens) and are computationally expensive for largescale problems due to the significant amount of iterations to find the optimum design for a given device functionality.More recently, design and optimization approaches based on DL techniques have been proposed and implemented for the design of nanostructures ^{29, 30, 31, 32, 33, 34, 35, 36}. Different reported approaches to date primarily rely on training a neural network (NN)(see Fig. 1(a)) using the response of a set of devices (found by numerical simulations) and using the trained NN to solve the inverse design problem. Despite impressing progress in this area, the reported solutions mostly focus on solving simple problems with reasonably smooth optimization landscapes ^{30}
that have a onetoone mapping of design parameter space to the device response space (i.e., any given response can be obtained by only a single set of design parameters) as shown in Fig. 1(b), where a vector of device response (
) is achieved by a unique vector of design parameters (). Unfortunately, most nanostructures of interest do not have this property. Figure 1(c) shows the optimization landscape of a more general problem in which the onetoone relation between design parameters and output response does not exist. This can result in convergence issues for the NN used for optimization (i.e., finding design parameters for a given output response). Efforts on converting the problem to a onetoone mapping by removing some training data sets (see Fig. 1(c)) ^{37} do not essentially help in solving the problem as most of the design space is not covered by these training datasets. Such approaches at most result in a NN that smooths out the optimization dataset (see Fig. 1(c)) without converging to the global optimum. Other proposed approaches (e.g., the use of tandem networks ^{29}) rely on first training a NN that relates the design space to the response space (i.e., for the forward problem), then cascading it as a pretrained NN with another NN that relates the response space to the design space (i.e., the inverse problem), and finally training the resulting network (from the response space to the design space) to avoid the nononetoone relation. However, such techniques do not solve the main problem; they at best smooth out the optimization landscape as shown in Fig. 1(c). Another notable recent approach is based on using generative adversial networks (GANs) to solve the inverse design problem ^{31}. This technique is built on training a network to solve the forward problem with zero error and use it to generate ground truth data in each iteration. Training such a forwardproblemsolver network with zero error in a general design problem is a major challenge and may require excessive computational resources. In the reported design problem, each desired output needs extensive computation (200,000 iterations to reach the convergence region for each structure)^{31}, which may reduce the value of using GANs if a perfect forward problem solver exists with comparable computation complexity (similar computations can be used to solve the design problem by exhaustive search using the perfect forwardproblem solver). Despite impressing results, the reported GANbased approach will be limited to simple design problems with noncomplex nanostructure. One can also think about breaking the NN for the design problem into two parts where one part is trained by limiting the design space to a smooth (onetoone) region and the second part is trained by the data in the complete design space. The success of such techniques highly depends on the complexity of the problem and the selection of the design parameters in the onetoone region (outside the deadzones in Fig. 1(d)) to converge to acceptable answers. As a result, these approaches can be used to design simple structures, which can also be designed using alternative approaches. Finding a reliable approach to fundamentally address this nonuniqueness issue (without limiting the optimization landscape to the onetoone region (or extrapolating from it, see Fig. 1(d)) is still a major challenge in using DL based approaches for the design of EM nanostructures.Another challenge in using DL techniques to design complex EM nanostructures is the large size of the response and design spaces resulting in the need to train a large NN. As an example, to study the spatial and spectral response of a MS with reasonable accuracy, the response space must constitute the sampled EM intensity in a twodimensional space and in frequency with spatial and spectral resolutions smaller than the smallest spatial and spectral features of the output response, respectively. This typically results in thousands of data points in the response space and quickly rises as the structures with sharper spatial and spectral features are designed. Combined with everincreasing number of design parameters in the nanostructures of recent interest, this results in a very large NN, which is difficult to be trained, even for the problems with onetoone optimization landscapes.
In this paper, we demonstrate a new approach for designing complex EM nanostructures by addressing both the networksize issue and the nonuniqueness issue. Our approach is based on reducing the dimensionality of both the design space and the response space through training multilayer NNs, called autoencoders ^{38}
. Once the dimensionality of the problem is reduced, the problem converts into a onetoone problem in the reduced spaces, which can be solved with considerably less computational complexity. In addition, by reducing the envisioned design parameters to few number of more complex design parameters (e.g., a nonlinear function of the weighted sum of the original design parameters), we can obtain valuable intuitive understanding of the roles of different design parameters in the response of the nanostructure. Such an efficacious approach paves the way for understanding, design, and optimization of complex EM nanostructures with far less computation complexity than the alternative approaches. In addition, a tradeoff between the accepted error and the complexity (and time) of the simulations can be used to solve different problems with desired degrees of computation complexity or to obtain quick (approximate) information about the role of design parameters in the overall device performance. Dimensionality reduction (DR) is a powerful technique in machine learning that has been used to effectively solve problems in a wide range of applications including robotics
^{39}, optical tomography ^{40}^{41}, handwritten digit classification ^{42}, remote sensing^{43}, medical science^{44}, genetics^{45}, and electronics ^{46}.To shows the applicability of our approach, we demonstrate its use for designing novel reconfigurable MS based on PCMs to form wideband amplitude modulation of nearinfrared (nearIR) light. The rest of the paper is organized as follows. Section 2 describes our approach based on DR in detail. Section 3 shows the application of this approach for designing MS for the desired functionality along with the corresponding results. Section 4 describes the advantage of our approach in investigating the underlying physics . Section 5 is dedicated to the discussion of the obtained results and the unique features of the new approach. Final conclusions are made in Section 6.
2 2. Dimensionality reduction of the design and response spaces in designing electromagnetic nanostructures
Figure 2 shows the schematic of the design approach based on DR of the design and response spaces assuming that the optimization landscape is nonunique (or manytoone), i.e., more than one set of design parameters can result in the same response. The original forward problem is shown by path 1 in Fig. 2, where each point in the design space (that includes a vector of dimension D corresponding to a set of design parameters) correspond to a point in the response space (which includes a vector of dimension R) through a manytoone relationship. A NN cannot be trained to inverse this relation as explained above. This is the main complication in the design and optimization problem. In our approach, we first use the DR technique to reduce the dimensionality of the response space as much as possible (i.e., reducing the size of the response vector in Fig.1(b)) while keeping the same number of points in the response space (see path 2 in Fig. 2). This concept is schematically shown in Fig. 3 in which a threedimensional manifold in the response space is reduced to a twodimensional manifold, which includes the same number of points in the response space, but each point is represented by a smaller size vector. Each feature in the reduced response space is related to the features of the original response space through a welldefined nonlinear function. This is a onetoone process.
In the next step, we reduce the dimensionality of the design space as much as possible (see path 4 in Fig. 2). In this process, the redundant nature of the design space is removed resulting in a onetoone relation between the reduced design space and the reduced response space (see path 3 in Fig. 2). After training the relevant DR mechanisms in Fig. 2, the relation between the original response space and the reduced design space (Paths in Fig.2) will be onetoone and thus, it can be simply inverted. Thus, our design problem will relate the desired response to the reduced design parameters (see path 5 in Fig. 2). The reduced design parameters are related to the original design parameters through a onetomany relation that are analytically available through the training process. Thus, we can find several design options by converting the resulting optimum reduced design parameters to several sets of the original design parameters. At this stage, design constraints (e.g., fabrication imperfections, structure robustness, characterization limitations, etc.) can be taken into account to choose the final design parameters.
The heart of our approach is the effective implementation of the DR technique to maximally reduce the dimensionality of both design and response spaces, especially the former. Several DR techniques have been developed in machine learning to facilitate classification, data visualization, reduction of the computation cost, etc. Among different options, principal component analysis (PCA)
^{47}, kernel principal component analysis (KPCA) ^{48}, Laplacian eigen map ^{49}, locally linear embedding ^{50}, and autoencoder ^{38} are the most effective techniques. Considering the features of these techniques, we believe that the autoencoder is the most suitable approach for solving inverse problems in general and designing EM nanostructures in particular.The general schematic of an autoencoder is shown in Fig. 4. Autoencoder is a multilayer NN that can encode the highdimensional data into lowdimensional data (using the encoder part in Fig. 4) and use another NN (see the decoder part in Fig. 4) to decode and recover the highdimensional data. In other words, the autoencoder in Fig. 4 is a feedforward NN where the input layer and the output layer have the same structure and are connected to each other with one or more hidden layers. The number of neurons in the layer with minimum number of neurons represents the dimension of the lowdimensional data (i.e., the dimension of the DR space). This layer is known as the bottleneck of the autoencoder. This way, an autoencoder concentrates the data from a highdimensional manifold in a given space around a lowdimensional manifold or a small set of such manifolds.The goal of an autoencoder is to map an original set of input data
to a lower dimensional set of output data (at the bottleneck) in which and are vectors with size and , respectively(), and contains the essential information of .To find the mapping from highdimensional to lowdimensional data, the autoencoder in Fig. 4 should be trained with a sufficiently large training dataset. The training part of the autoencoder can be considered as an optimization problem where the algorithm minimizes a cost function. The cost function is a measurement of discrepancy between the output of the autoencoder and the input data. The meansquared error (MSE) is used as the cost function of the autoencoder, and the error is minimized using the backpropagation method
^{51}. Assuming the output of the autoencoder structure in Fig. 6 for the input is represented by , the reconstruction MSE of the trained autoencoder is defined as:(1) 
Where represents the number of validation (or test) instances (not used for training) that are used to validate the trained autoencoder (but not used for training) . The number of layers and the topology of the NN is also found using an adhoc method (by trial and error). The training dataset for the design of EM nanostructures is obtained by using numerical simulation of the structure using a random set of input design parameters.
In the approach shown in Fig. 2, we first reduce the dimensionality of the response space by training an autoencoder (see Fig. 5(a)). In the next step, we form a pseudoencoder that relates the original design space to the reduced response space as shown in Fig. 5(b). The reason for naming the structure in Fig. 5(b) a pseudoencoder is the fact that its input and output are from different spaces (in contrast to a conventional autoencoder in Fig. 4). By training the pseudoencoder in Fig. 5(b) to reach the minimum size of the bottleneck layer, we reach the reduced design space. Each parameter in this space is related to the original design parameters through a nonlinear function defined by the NN structure of the pseudoencoder from the original design space to the reduced design space (or the bottleneck) in Fig. 5(b). The training approach is similar to that explained for a general autoencoder in Fig. 4 or Fig. 5(a). The pseudoencoder in Fig. 5(b) corresponds to the paths 3 and 4 in Fig. 2, i.e., these two paths are trained together.
Once the DR of the two spaces are complete, we form a NN by cascading the pseudoencoder in Fig. 5(b) with the pretrained decoder part of Fig. 5(a) to form a completely trained NN for solving the forward problem as shown in Fig. 5(c). The resulting NN in Fig. 5(c) relates the original design parameters to the original response space using a unique set of analytic equations defined by different layers of the NN. While this analytic relation is complicated for a largesize network, it provides extremely valuable information about the roles of different design parameters in the response of the nanostructure with minimal computation complexity (technically by calculating the complex analytic formulas in a conventional environment like MATLAB). However, the goal of this paper is the design of EM nanostructures for which the inverse problem has to be solved. For this purpose, we will use a twostep approach. In the first step, we find the inverse of the part of the NN in Fig. 5(c) that relates the reduced design space to the output space. The resulting inverse network is shown in Fig. 5(d). This is easily achievable as the relation between the reduced design space and the original response space is onetoone (see path 5 in Fig. 2). The NN in Fig. 5(d) allows us to obtain the optimal reduced design parameters for any given desired response. This is the last part in our approach where the DL approaches can be used. The final step is to relate the reduced design parameters to the original design parameters (i.e., the inverse of path 4 in Fig. 2). This is a nonunique relation, i.e., it can provide several sets of design parameters from a given reduced set of design parameters. Fortunately, the encoder part of the pseudoencoder in Fig. 5(b) relates the reduced design parameters to the original design parameters analytically (through the formulation of the underlying NN at different layers). We can use these equations to move layerbylayer from the reduced design parameters to the original design parameters. In this backward process, we can reduce the number of possible solutions by imposing constraints such as fabrication limitations. This approach can provide many possible solutions for a design problem, which is expected due to the nonuniqueness of the problem. Note also that within this design problem, we can use the obtained knowledge about the role of the design parameters (using the forward solver in Fig. 5(c)) and the relation between the reduced design parameters and the original design parameters (using the encoder part of the pseudoencode in Fig. 5(b)) to reduce the complexity in solving the design problem. In this paper, we use the analytic relation between the original and reduced design spaces to completely search the original design space to find the point(s) that correspond to the desired point in the reduced design space.
In addition to solving the nonuniqueness issue, the approach in Fig. 5 considerably reduces the computation cost by reducing the dimensionality of the two spaces. It is clear that the training of the pseudoencoder that relates the design space to the reduced response space (see Fig. 5(b)) requires much less computation compared to training of a NN that relates the design space to the original (nonreduced) response space. Furthermore, the calculation of the inverse NN in Fig. 5(d) does not impose significant computation cost due to its onetoone nature.
3 3. Application to the design of hybrid reconfigurable plasmonicphotonic metasurfaces
To show the applicability of the design approach, we consider a simple design problem for the implementation of a reconfigurable multifunctional MS enabling high performance optical modulation as shown in Fig. 6. The metasurface (MS) in Fig. 6 is composed of a periodic array of gold (Au) nanoribbons fabricated on top of a thin layer of germanium antimony telluride (GeSbTe or in short GST), which is a nonvolatile PCM whose index of refraction can be significantly modified (e.g., from 4.5 to 7 in the nearinfrared region) ^{52}when it undergoes transition from the amorphous to the crystalline state in the near infrared regime or vice versa. In addition, using GST in intermediate states between amorphous and crystalline results in a wide range of tunability for its index of refraction. The GST layer is deposited on a thin film of Au as shown in Fig. 6. By applying electric signals to the Au nanoribbons, the state of the GST underneath that nanoribbon is controlled through resistive heating ^{53}. In addition, by controlling the electric stimulus intermediate states (between amorphous and crystalline) can be obtained for GST ^{54}. We limit the number of GST transition states to 11 (i.e., amorphous, crystalline, and 9 intermediate states)^{55}. The supercell (limited to 3 different building blocks to prevent excitation of high diffraction orders) of the MS in Fig. 6 is composed of three Au nanoribbons with different widths (w, w, and w) and three crystallization levels (l, l, l, corresponding to three indices of refraction, see Methods for more details) of GST underneath with the same height (h). The pitches of the 3 building blocks of the supercell are represented by p, p, and p, respectively. As a result, the MS in this work has 10 design parameters (i.e., dimensionality of the design space is equal to 10) with different units (3 unitless indices of refraction and 7 lengths with units of nanometers).
As a simple functionality, we are interested in amplitude modulation of the incident light at with a considerable bandwidth around the central wavelength. The MS in Fig. 6 is illuminated with a plane wave of light with variable wavelengths in the desired range (from 1250 nm to 1850 nm). The polarization of the incident light is such that the electric field (i.e., E) is perpendicular to the grating direction of the MS. The response of the system is the MS reflectance (calculated as the farfield reflection intensity divided by the intensity of the incident field and integrated over a surface area equal to one supercell in the farfield). The resulting reflectance is sampled at 200 equallyspaced wavelengths in the 12501850 nm range. This results in a response space dimensionality of 200. To obtain the data for training and validation of the DR autoencoders in Fig. 5, we simulate the structure in Fig. 6 with 4000 randomly generated instances (3600 for training and 400 for validation) formed by randomly selecting the design parameters in the acceptable variation ranges shown in the caption of Fig. 6. The simulations were performed using the finite element method (FEM) in the COMSOL Multiphysics environment (see Methods for details).
In the next step, we use the training data to train a series of autoencoders with different numbers of hidden layers (ranging from 3 to 9) to study the DR of the response space from 200 to different values in the range of 1 to 20. For each autoencoder, MSE is calculated by finding the square of the norm of the difference between the reflectance vector obtained from the decoder part of the autoencoder and that obtained from the FEM simulations (also called ground truth data) for each one of the validation data. Note that each vector has 200 elements corresponding to the reflectance at 200 selected wavelengths in the 12501850 nm range.
Figure 7(a) shows the calculated MSE as a function of the dimensionality of the reduced response space. Figure 7(b) shows the comparison of the actual reflectance spectrum for the original data and the reconstructed data using the autoencoder for different dimensionalities of the reduced response space. It is clear from both Figs. 7(a) and 7(b) that the dimension of the response space can be reduced from 200 to 10 with negligible MSE (less than 10). This is a clear advantage of our optimization technique.
Among different autoencoder architectures tested for the response space, the one with 5 layers (with the number of neurons in consecutive layers being 200501050200) is selected based on its low MSE and computation costs to form the pseudoencoder architecture in Fig. 5(b). We also choose 4 layers (102015x) for the encoder part of the pseudoencoder in Fig. 5(b). For each set of values for the dimensions of the reduced response space and the reduced design space, we train the resulting pseudoencoder using the 3600 training data, and we calculate the MSE by comparing the output of the pseudoencoder with the actual output using the 400 validation data. The results for 4 different dimensionalities of the reduced response space are shown in Fig. 8(a). Figure 8(b) shows representative reflectance spectra for 3 different values of the dimensionality of the reduce design space. Using Fig. 8(a), it is clear that the dimension of the design space can be reduced from 10 to 5 without imposing much error.
Using Figs. 7 and 8, we choose the dimensionality of the reduced response space and the reduced design space to be 10 and 5, respectively (10201552030201050200). This considerably reduces the computation time as the dimension of the resulting problem is defined in a 5 10 rather than 10 200. Using these values, the final NN architecture for the analysis (or the solution of the forward problem) of the MS in Fig. 6 is formed according to Fig. 5(c). It is clear that the training of the pseudoencoder that relates the design space to the reduced response space (see Fig. 5(b)) requires much less computation compared to training of a NN that relates the design space to the original (nonreduced) response space.
To form a platform for designing MSs with an arbitrary response, we first find the inverse of the network from the original response space to the reduced design space as shown in Fig. 5(d). This is not computationally extensive due to the onetoone nature of the problem. For this purpose, the pretrained encoder part of the DR algorithm for the response space (left side of Fig. 5(a)) is combined with a NN that connects the reduced response space to the reduced design space. This added NN is trained using the same 3600 training data to form the inverse network that relates the desired response to the reduced deign parameters. The resulting onetoone trained platform (see Fig. 5(d)) results in finding the 5 reduced design parameters for the desired response. To find the 10 original design parameters, we solve the onetomany problem through an analytical search approach using the encoder part of the pseudoencoder for DR of the design space (first part of the platform in Fig. 5(b)). This encoder part relates the original design parameters analytically (through the NN formulation) to the reduced design parameters. Thus, the exhaustive search of the design space is not computationally extensive. We use MATLAB to perform this calculation (sweeping each parameter over 10 possible values) using the minimization of the MSE (defined by the integral of the square of the difference between the desired and the resulting light intensities over the operation bandwidth) as the optimization goal.
Figure 9 shows the results for the design of prefect light absorber for operation in the 15001700 nm wavelength range using the MS in Fig. 6. The desired response is zero reflectivity over the entire operation bandwidth. The overall MSE for the response of the optimal structure in Fig. 9 is 0.0147 MSE. The reflectance for two other (nonoptimal) designs with considerably different design parameters are also shown in Fig. 9. The set of design parameters along with the MSE for these three structures are listed in Table 1.
4 4. Understanding the physics of lightmatter interaction
A main advantage of our approach is the possibility of investigating the underlying physics of the device operation and obtaining intuitive information about the roles of different design parameters on its response. To show this capability, we use our approach with a pseudoencoder (10410502010) to model the MS in Fig. 6. Figure 10(a) shows the resulting pseudoencoder with the dimension of the reduced design space being 4 with green and red arrows representing positive and negative weights, respectively. Note that the DR of the design space is performed with only one encoder layer. Figure 10(b) shows the values of the weights for the monolayer encoder. Each weight is multiplied by its corresponding design parameter to form the inputs to the node of the bottleneck layer. The larger the weight, the stronger the contribution of the corresponding design parameter will be. This strength is also shown in Fig. 10(a) by the thickness of the arrows that connect the nodes of the two layers. As shown in Fig 10(b), the height of the structure (h) plays an important role in changing the response compared to other design parameters as connects to all 4 nodes in the bottleneck layer with reasonably strong weights. Moreover, the crystallization levels l, l, and l can only change one of the reduced response features as they mainly connect to only one node (the purple node) in the bottleneck layer. As a result, as long as the total input to that purple node is fixed, the response will stay the same regardless of how the values of l, l, and l change. This conclusion is reached by assuming a small error in training the pseudoencoder and neglecting the small weights (or arrows in Fig. 10(a)) that connect l, l, and l to the nodes of the bottleneck layer. To test this conclusion, we vary l, l, and l while keeping their weighted sum (according to the trained pseudoencoder) and all other 7 design parameters for the MS in Fig. 6 constant, and we calculate the response of the MS using bruteforce COMSOL simulations (no pseudoencoder intervention). The results for two different weighted sums of l, l, and l are shown in Fig. 10(c). Figure 10(c) clearly confirms our observation from the trained pseudoencoder that l, l, and l effectively act as one design parameter (through their weighted sum). Figure 10(d) shows the results of COMSOL simulations when the structure height () is changed while keeping all other 9 design parameters fixed. The large range of variation of the response in Fig. 10(d) clearly shows the importance of h as a design parameter. It is interesting to see from Fig. 10(d) that different responses for different values of h have low correlation while the responses for different values of the weighted sum of l, l, and l (i.e., blues curves and red curves in Fig. 10(c)) show a similar trend with different locations of peaks and valleys. This suggests that the parameter can be used to obtain different classes of responses while the weighted sum of l, l, and l can be used to finely tune a given class of response. The details of the design parameters for each case are shown in Table S2.
The important observations about the role of different design parameters were obtained from our deeplearning approach without taking any information about the physics of the structure into account. Nevertheless, these observations agree with the physical intuition about the structure in Fig. 6. Each unit cell in this structure is composed of three plasmonic building blocks formed between the Au layer underneath and each Au nanoribbon on the top GST laye (see Fig. 6). Since the supermode of each building block is formed by coupling of the surface plasmon polaritons at the two Au layers, its properties strongly depend on the height of the GST layer (h), which directly controls the coupling strength^{56}. Thus, strong dependence of the MS response on is expected. Figures 11(a) and 11(b) show the electric field patterns for the unit cell structure in Fig. 6 for two different values of h, confirming the strong dependence of the spatial mode profile on h.
To consider the effect of variation of the crystallization fractions (l, l, and l), we note that the reflection response of the overall MS is essentially the sum of three responses defined by the three plasmonic resonators in each unit cell. By combining three wideband resonances with different resonance wavelengths, a wideband reflection response is obtained. Figure 11 shows the variation of the reflectance of the MS with frequency for a given set of l, l, and l values (60%, 70%, and 80%, respectively). The inset shows the field profiles of the three plasmonic resonators within each unit cell at different wavelength regimes.
Due to pronounced lightmatter interaction of the supermode with the GST layer at the higher wavelength (i.e., 16501850 nm), we expect that most of the resistive loss occurs in the building block with high crystallization level (i.e., l) accommodating more free charge carriers. This effect is clarified in the inset of Fig. 11 at higher wavelengths (red border) showing that a good portion of absorption takes place in the rightmost building block (i.e., l). Figure 11 also shows that the absorption loss in the middle wavelength window (i.e., 14501650 nm, shown by green) occurs mostly in the building block with lower crystallization level (i.e., the leftmost building blocks with (i.e., l and l). Finally, Fig. 11 shows similar contributions from the three building blocks at lower wavelengths (e.g., 12501450 nm). This is due to the fact that by increasing the level of crystallization in this regime, the optical constant of GST varies significantly leading to decrease in the lightmatter interaction. This explains the collective role of l, l, and l observed through training the pseudoencoder. Note that obtaining this observation from the basic device properties was not as trivial as that of the role of .
While some of the conclusions about the role of design parameters in Fig. 6 obtained by training the pseudoencoder could also be obtained by the underlying mode properties of the DLbased (e.g., by analyzing the modes of the plasmonic resonators), the ability of our approach in providing useful information about the physics of wavematter interaction in nontrivial structures (e.g., nonlinear and dispersive metamaterials) will be extremely valuable. Indeed, by using this approach to find and understand new phenomena in such nontrivial structures, new ideas for forming novel device can be generated. This is a major advantage of our approach over all existing design approaches, especially those that rely on multiple bruteforce simulations of the structure for different design parameters.
5 5. Discussion
Figure 9 and Table 1 clearly show the ability of our approach in designing MSs with considerably reduced computation complexity. Figure 9 obviously verifies the remarkable modulation depth (between ON and OFFstate) of the optimized structure. They also show the importance of the understanding of the manytoone nature of the design problem. Considering the many options for the original design parameters that correspond to a single set of reduced design parameters, we can easily enforce the fabrication restrictions and other design preferences in the last part of the design approach and find a set of design parameters that results in closetooptimal response. It is important to note that the availability of the analytic relation between the original and reduced design spaces makes the bruteforce optimization (e.g., using analytic search) computationally feasible even for a large number of design parameters. Nevertheless, more sophisticated constrained optimization techniques can be used to solve the last (i.e., manytoone) part of the design problem with explicit inclusion of the fabrication and other designrelated constraints. Such techniques are currently under investigation and will be the subject of future publications.
Design Parameters  
Design  h  l  l  l  p  p  p  w  w  w  MSE 
190  0.5  0.6  0.7  650  650  550  350  500  200  0.0147  
190  0  0.2  0.8  650  650  350  450  250  250  0.0149  
190  0.5  0.1  0.7  650  450  450  200  350  300  0.0152  
190  0.3  0.6  0.8  650  550  550  250  300  450  0.0172 
A unique feature of our DRbased approach is the computation simplicity while appreciating the manytoone nature of the problem. By not considering the latter explicitly, several other existing NNbased techniques are technically limited to only smoothenough problems, or they require apriori assumptions to limit the search for the optimal solution in to a given region in the design space, where the relation to the response space is onetoone (as discussed in Section 1). Nevertheless, by reducing the dimensionality of the problem, our approach requires less computation than any other alternative. For example, in the design problem studied here, we reduced a 10 200 dimensional problem to a 5 10 one.
Compared to bruteforce optimization approaches (e.g., exhaustive search), our technique requires far less computation. For example, by assuming only 10 possible values for the 7 analog variables (i.e., h, w, w, and w, p, p, and p) and 11 values for the discrete ones (i.e., l, l, l) in the design problem in Fig. 6, the exhaustive search algorithm requires the complete EM simulation of the structure for more than 10 times, which is essentially intractable. However, our optimization requires only 4000 EM simulations along with the training process that requires far less computations. Indeed, the entire training of the forward and inverse parts of the platform in Fig. 5 for the design problem in Fig. 6 (results shown in Fig. 9) took less than 3 hours using a simple personal computer with a 3.4 GHz core i76700 CPU and 8 GB of random access memory (RAM).
The main issue that can happen in the nonoptimal use of the DR platform is the nonuniqueness (or nononetoone) relation between the original response space and the reduced design space in the platform shown in Fig. 5(d). This can happen by not carefully considering different numbers for the reduced dimensions in the DR algorithms. Nevertheless, the DR algorithm considerably reduces the possibilities of the manytoone problem by limiting the dimensionalities of the design and response spaces. Thus, the resulting structure in Fig. 5(d) can provide a closetooptimal design even for the cases that the resulting architecture in Fig. 5(d) is not absolutely onetoone. A rigorous mathematical study of the conditions for the dimensionality of the reduced spaces from machinelearning point of view can provide more solid guidelines in selecting the dimensionality of the reduced spaces. However, such a rigorous mathematical study is outside the scope of this paper.
6 6. Conclusion
We demonstrated here a new DLbased approach for the design of EM nanostructures with a wide range of design possibilities. We showed that by reducing the dimensionality of the response and design spaces using an autoencoder and a psuedoencoder, we can convert the initial manytoone problem into a onetoone (or in the worst case, close to onetoone) problem plus a simple onetomany problem that can be solved using bruteforce analytical formulas. The resulting approach considerably reduces the computational complexity of both the forward problem and the inverse (or design) problem. In addition, it allows for the inclusion of the design restrictions (e.g., fabrication limitations) without adding computational complexities. It also provides valuable information about the roles of design parameters in the response of the EM structure, which can potentially enable novel phenomena and devices. Finally, this technique can be extended to solve many different optimization problems in a wide range of disciplines as long as enough data for training the incorporated NNs are provided.
7 7. Methods
The fullwave EM simulations were carried out using the finite element method (FEM) enabled by linking the commercial software package COMSOL Multiphysics 5.3 (wave optics module) with MATLAB to expedite the design, optimization, and analysis processes. Floquet periodic and perfectlymatched layer boundary conditions were used along transverse xaxis and the zaxis in Fig. 6, respectively. The structure was assumed infinite along the ydirection. A linearly polarized planewave of light, excited the MS in the wavelength range of 12501850 nm. The refractive index (n) and absorption coefficient (k) data for amorphous and crystalline GST, Au, and SiO were obtained from the literature ^{57, 58}. The computation domain was meshed using triangular elements with a maximum size of /(10 n) in SiO and GST , and of /50 in Au with being the freespace wavelength. The effective dielectric constant associated with the intermediate states of GST were approximated via the effective medium theory. Among different options, LorentzLorenz formula more accurately describes effective permittivity as^{59}:
(2) 
where and are the permitivitties of the crystalline and amorphous GST, respectively, and , ranging from 0 (amorphous) to 1 (crystalline), is the crystallization fraction of GST.
8 8. Acknowledgements
The authors thank Ali A. Eftekhar for helpful discussion.
9 9. Author contributions
YK and SA contributed equally to this work. The initial idea was developed by YK and SA, and its implementation was discussed by all authors. YK performed the training optimization of the autoencoder and the pseudoencoder, and SA developed the simulation results for training and validation. SA proposed the initial idea for the electromagnetic nanostructure for light absorption, which was discussed in further details by all authors. AA managed the project. All authors participated in writing the manuscript.
10 10. Competing interests
The authors declare no competing financial interest.
References
 Melikyan et al. 2014 Melikyan, A.; Alloatti, L.; Muslija, A.; Hillerkuss, D.; Schindler, P. C.; Li, J.; Palmer, R.; Korn, D.; Muehlbrandt, S.; Van Thourhout, D. Highspeed plasmonic phase modulators. Nature Photonics 2014, 8, 229.
 Zhu et al. 2017 Zhu, T.; Zhou, Y.; Lou, Y.; Ye, H.; Qiu, M.; Ruan, Z.; Fan, S. Plasmonic computing of spatial differentiation. Nature communications 2017, 8, 15391.
 Rodrigo et al. 2015 Rodrigo, D.; Limaj, O.; Janner, D.; Etezadi, D.; De Abajo, F. J. G.; Pruneri, V.; Altug, H. Midinfrared plasmonic biosensing with graphene. Science 2015, 349, 165–168.
 Liu et al. 2011 Liu, X.; Tyler, T.; Starr, T.; Starr, A. F.; Jokerst, N. M.; Padilla, W. J. Taming the blackbody with infrared metamaterials as selective thermal emitters. Physical review letters 2011, 107, 045901.
 5 Huang, L.; Chen, X.; Mühlenbernd, H.; Zhang, H.; Chen, S.; Bai, B.; Tan, Q.; Jin, G.; Cheah, K.W.; Qiu, j. v. p. y. p., ChengWei Threedimensional optical holography using a plasmonic metasurface.
 Yu et al. 2011 Yu, N.; Genevet, P.; Kats, M. A.; Aieta, F.; Tetienne, J.P.; Capasso, F.; Gaburro, Z. Light propagation with phase discontinuities: generalized laws of reflection and refraction. science 2011, 1210713.
 Kildishev et al. 2013 Kildishev, A. V.; Boltasseva, A.; Shalaev, V. M. Planar photonics with metasurfaces. Science 2013, 339, 1232009.
 Arbabi et al. 2015 Arbabi, A.; Horie, Y.; Bagheri, M.; Faraon, A. Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission. Nature nanotechnology 2015, 10, 937.
 Khorasaninejad et al. 2016 Khorasaninejad, M.; Chen, W. T.; Devlin, R. C.; Oh, J.; Zhu, A. Y.; Capasso, F. Metalenses at visible wavelengths: Diffractionlimited focusing and subwavelength resolution imaging. Science 2016, 352, 1190–1194.
 Jahani and Jacob 2016 Jahani, S.; Jacob, Z. Alldielectric metamaterials. Nature nanotechnology 2016, 11, 23.
 Lin et al. 2014 Lin, D.; Fan, P.; Hasman, E.; Brongersma, M. L. Dielectric gradient metasurface optical elements. science 2014, 345, 298–302.
 Taghinejad et al. 2018 Taghinejad, M.; Taghinejad, H.; Xu, Z.; Lee, K.T.; Rodrigues, S. P.; Yan, J.; Adibi, A.; Lian, T.; Cai, W. Ultrafast Control of Phase and Polarization of Light Expedited by HotElectron Transfer. Nano letters 2018, 18, 5544–5551.
 AbdollahRamezani et al. 2015 AbdollahRamezani, S.; Arik, K.; Khavasi, A.; Kavehvash, Z. Analog computing using graphenebased metalines. Optics letters 2015, 40, 5239–5242.
 Sun et al. 2012 Sun, S.; Yang, K.Y.; Wang, C.M.; Juan, T.K.; Chen, W. T.; Liao, C. Y.; He, Q.; Xiao, S.; Kung, W.T.; Guo, G.Y. Highefficiency broadband anomalous reflection by gradient metasurfaces. Nano letters 2012, 12, 6223–6229.
 Arbabi et al. 2017 Arbabi, A.; Arbabi, E.; Horie, Y.; Kamali, S. M.; Faraon, A. Planar metasurface retroreflector. Nature Photonics 2017, 11, 415.
 Decker et al. 2015 Decker, M.; Staude, I.; Falkner, M.; Dominguez, J.; Neshev, D. N.; Brener, I.; Pertsch, T.; Kivshar, Y. S. Highefficiency dielectric Huygens surfaces. Advanced Optical Materials 2015, 3, 813–820.
 Chen et al. 2018 Chen, W. T.; Zhu, A. Y.; Sanjeev, V.; Khorasaninejad, M.; Shi, Z.; Lee, E.; Capasso, F. A broadband achromatic metalens for focusing and imaging in the visible. Nature nanotechnology 2018, 13, 220.
 Molesky et al. 2018 Molesky, S.; Lin, Z.; Piggott, A. Y.; Jin, W.; Vucković, J.; Rodriguez, A. W. Inverse design in nanophotonics. Nature Photonics 2018, 12, 659.
 Piggott et al. 2017 Piggott, A. Y.; Petykiewicz, J.; Su, L.; Vučković, J. Fabricationconstrained nanophotonic inverse design. Scientific Reports 2017, 7, 1786.
 Lu and Vučković 2013 Lu, J.; Vučković, J. Nanophotonic computational design. Optics express 2013, 21, 13351–13367.
 Su et al. 2017 Su, L.; Piggott, A. Y.; Sapra, N. V.; Petykiewicz, J.; Vuckovic, J. Inverse design and demonstration of a compact onchip narrowband threechannel wavelength demultiplexer. ACS Photonics 2017, 5, 301–305.
 Frellsen et al. 2016 Frellsen, L. F.; Ding, Y.; Sigmund, O.; Frandsen, L. H. Topology optimized mode multiplexing in silicononinsulator photonic wire waveguides. Optics express 2016, 24, 16866–16873.
 Piggott et al. 2014 Piggott, A. Y.; Lu, J.; Babinec, T. M.; Lagoudakis, K. G.; Petykiewicz, J.; Vučković, J. Inverse design and implementation of a wavelength demultiplexing grating coupler. Scientific reports 2014, 4, 7210.
 Englund et al. 2005 Englund, D.; Fushman, I.; Vuckovic, J. General recipe for designing photonic crystal cavities. Optics express 2005, 13, 5961–5975.
 Seidel and Rappaport 1994 Seidel, S. Y.; Rappaport, T. S. Sitespecific propagation prediction for wireless inbuilding personal communication system design. IEEE transactions on Vehicular Technology 1994, 43, 879–891.
 Gondarenko and Lipson 2008 Gondarenko, A.; Lipson, M. Low modal volume dipolelike dielectric slab resonator. Optics express 2008, 16, 17689–17694.
 Håkansson and SánchezDehesa 2005 Håkansson, A.; SánchezDehesa, J. Inverse designed photonic crystal demultiplex waveguide coupler. Optics Express 2005, 13, 5440–5449.
 Ma et al. 2013 Ma, Y.; Zhang, Y.; Yang, S.; Novack, A.; Ding, R.; Lim, A. E.J.; Lo, G.Q.; BaehrJones, T.; Hochberg, M. Ultralow loss single layer submicron silicon waveguide crossing for SOI optical interconnect. Optics express 2013, 21, 29374–29382.
 Liu et al. 2018 Liu, D.; Tan, Y.; Khoram, E.; Yu, Z. Training deep neural networks for the inverse design of nanophotonic structures. ACS Photonics 2018, 5, 1365–1369.
 Peurifoy et al. 2018 Peurifoy, J.; Shen, Y.; Jing, L.; Yang, Y.; CanoRenteria, F.; DeLacy, B. G.; Joannopoulos, J. D.; Tegmark, M.; Soljačić, M. Nanophotonic particle simulation and inverse design using artificial neural networks. Science advances 2018, 4, eaar4206.
 Liu et al. 2018 Liu, Z.; Zhu, D.; Rodrigues, S.; Lee, K.T.; Cai, W. A Generative Model for the Inverse Design of Metasurfaces. Nano letters 2018,
 Tahersima et al. 2018 Tahersima, M. H.; Kojima, K.; KoikeAkino, T.; Jha, D.; Wang, B.; Lin, C.; Parsons, K. Deep Neural Network Inverse Design of Integrated Nanophotonic Devices. arXiv preprint arXiv:1809.03555 2018,
 Zhang et al. 2018 Zhang, T.; Liu, Z.; Dai, J.; Han, X.; Li, J.; Zhou, Y.; Xu, K. Spectrum prediction and inverse design for plasmonic waveguide system based on artificial neural networks. arXiv preprint arXiv:1805.06410 2018,
 Ma et al. 2018 Ma, W.; Cheng, F.; Liu, Y. DeepLearning Enabled OnDemand Design of Chiral Metamaterials. ACS nano 2018,
 Qu et al. 2018 Qu, Y.; Jing, L.; Shen, Y.; Qiu, M.; Soljacic, M. Migrating Knowledge between Physical Scenarios based on Artificial Neural Networks. arXiv preprint arXiv:1809.00972 2018,
 Inampudi and Mosallaei 2018 Inampudi, S.; Mosallaei, H. Neural network based design of metagratings. Applied Physics Letters 2018, 112, 241102.
 Kabir et al. 2008 Kabir, H.; Wang, Y.; Yu, M.; Zhang, Q.J. Neural network inverse modeling and applications to microwave filter design. IEEE Transactions on Microwave Theory and Techniques 2008, 56, 867–879.
 Hinton and Salakhutdinov 2006 Hinton, G. E.; Salakhutdinov, R. R. Reducing the dimensionality of data with neural networks. science 2006, 313, 504–507.
 Ciocarlie et al. 2007 Ciocarlie, M.; Goldfeder, C.; Allen, P. K. Dimensionality reduction for handindependent dexterous robotic grasping. 2007,
 Bhowmik et al. 2016 Bhowmik, T.; Liu, H.; Ye, Z.; Oraintara, S. Dimensionality reduction based optimization algorithm for sparse 3D image reconstruction in diffuse optical tomography. Scientific reports 2016, 6, 22242.
 He et al. 2005 He, X.; Yan, S.; Hu, Y.; Niyogi, P.; Zhang, H.J. Face recognition using laplacianfaces. IEEE transactions on pattern analysis and machine intelligence 2005, 27, 328–340.
 Hinton et al. 1997 Hinton, G. E.; Dayan, P.; Revow, M. Modeling the manifolds of images of handwritten digits. IEEE transactions on Neural Networks 1997, 8, 65–74.
 Efremenko et al. 2014 Efremenko, D.; Doicu, A.; Loyola, D.; Trautmann, T. Optical property dimensionality reduction techniques for accelerated radiative transfer performance: Application to remote sensing total ozone retrievals. Journal of Quantitative Spectroscopy and Radiative Transfer 2014, 133, 128–135.

Breger et al. 2017
Breger, A.; Ehler, M.; Bogunovic, H.; Waldstein, S.; Philip, A.; SchmidtErfurth, U.; Gerendas, B. Supervised learning and dimension reduction techniques for quantification of retinal fluid in optical coherence tomography images.
Eye 2017, 31, 1212.  Kim and Tidor 2003 Kim, P. M.; Tidor, B. Subsystem identification through dimensionality reduction of largescale gene expression data. Genome research 2003, 13, 1706–1718.

Choi et al. 2017
Choi, S.; Shin, J. H.; Lee, J.; Sheridan, P.; Lu, W. D. Experimental demonstration of feature extraction and dimensionality reduction using memristor networks.
Nano letters 2017, 17, 3113–3118.  Jolliffe 2002 Jolliffe, I. T. Springer series in statistics; SpringerVerlag New York, 2002; Vol. 29.

Schölkopf et al. 1998
Schölkopf, B.; Smola, A.; Müller, K.R. Nonlinear component analysis as a kernel eigenvalue problem.
Neural computation 1998, 10, 1299–1319.  Belkin and Niyogi 2003 Belkin, M.; Niyogi, P. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation 2003, 15, 1373–1396.
 Roweis and Saul 2000 Roweis, S. T.; Saul, L. K. Nonlinear dimensionality reduction by locally linear embedding. science 2000, 290, 2323–2326.
 Rumelhart et al. 1985 Rumelhart, D. E.; Hinton, G. E.; Williams, R. J. Learning internal representations by error propagation; 1985.
 Wuttig et al. 2017 Wuttig, M.; Bhaskaran, H.; Taubner, T. Phasechange materials for nonvolatile photonic applications. Nature Photonics 2017, 11, 465.
 Tuma et al. 2016 Tuma, T.; Pantazi, A.; Le Gallo, M.; Sebastian, A.; Eleftheriou, E. Stochastic phasechange neurons. Nature nanotechnology 2016, 11, 693.
 Ríos et al. 2015 Ríos, C.; Stegmaier, M.; Hosseini, P.; Wang, D.; Scherer, T.; Wright, C. D.; Bhaskaran, H.; Pernice, W. H. Integrated allphotonic nonvolatile multilevel memory. Nature Photonics 2015, 9, 725.
 Feldmann et al. 2017 Feldmann, J.; Stegmaier, M.; Gruhler, N.; Ríos, C.; Bhaskaran, H.; Wright, C.; Pernice, W. Calculating with light using a chipscale alloptical abacus. Nature Communications 2017, 8, 1256.
 Maier 2007 Maier, S. A. Plasmonics: fundamentals and applications; Springer Science & Business Media, 2007.
 Shportko et al. 2008 Shportko, K.; Kremers, S.; Woda, M.; Lencer, D.; Robertson, J.; Wuttig, M. Resonant bonding in crystalline phasechange materials. Nature materials 2008, 7, 653.
 Huang et al. 2016 Huang, Y.W.; Lee, H. W. H.; Sokhoyan, R.; Pala, R. A.; Thyagarajan, K.; Han, S.; Tsai, D. P.; Atwater, H. A. Gatetunable conducting oxide metasurfaces. Nano letters 2016, 16, 5319–5325.
 Chu et al. 2016 Chu, C. H.; Tseng, M. L.; Chen, J.; Wu, P. C.; Chen, Y.H.; Wang, H.C.; Chen, T.Y.; Hsieh, W. T.; Wu, H. J.; Sun, G. Active dielectric metasurface based on phasechange medium. Laser & Photonics Reviews 2016, 10, 986–994.
Comments
There are no comments yet.