1. Introduction

Metasurfaces are artificial periodic or non-periodic 2D-structures containing meta atoms on subwavelength scale, enabling us to control the intrinsic properties of light [1,2]. By limiting these nano scale structures to surfaces rather than bulk metamaterials, their characterization becomes significantly easier and their fabrication complexity is reduced. Nevertheless, by using optically thin materials we simultaneously limit their available functionalities [3]. This problem can be overcome by using high index dielectric [4] or by layering metasurfaces [5,6]. The design of the latter is the problem we are addressing here.

Layering metasurfaces adds many degrees of freedom and the resulting stacks show a large range of possibilities in controlling light for various applications, such as asymmetric transmission [7,8], wavefront control [9], or broadband circular polarizers [10]. Due to the nature of layered or stratified media, stacking allows to functionalize Fabry-Perot resonances by tailoring both layer distances and particle resonances. The overall response of a stack is then the combination of several physical effects resulting both from the layering and structural properties. These include but are not limited to Fano resonances or anisotropic field coupling in between layers [11]. However, while simple geometries can be realized experimentally and single metasurfaces are investigated to great extend; efficient multi-objective optimization tools for layered metasurface stacks are still lacking. Among the many optimization techniques machine learning shows great promise for computationally efficient optimization. The data accumulated during design projects in previous years can be used to train neural networks more systematically [12], or it can be used for other optimization techniques like statistical approaches [13]. One problem when inversely designing metamaterials is, that different designs can lead to the same functionality. In particular, in the context of layered metasurfaces, in many cases, the ordering of the surfaces constituting the stack doesn’t matter. Hence, the solution is not unique. This can be problematic for traditional optimization procedures and several attempts in the inverse design of metamaterials, as e.g generative models like variational autoencoders [14], tandem Neural Networks [15,16] or adaptive genetic algorithms [17], have been developed to overcome that problem. The Neural Network based methods often require large sets of training data. However, efficiently constructing stacks of metasurfaces for a desired functionality still poses a challenge, to our knowledge all proposed methods so far treat stacks with a predefined number of layers, merely addressing the geometrical parameters within them.

When treating metasurface stacks of incommensurable periods, i.e. periods of non-rational ratios, rigorous simulation methods become very expensive due to quickly expanding super-cells, necessary to describe whole stack. On top of that the number of possibilities on how to stack these metasurfaces scales exponentially with the number of available metasurfaces Most design approaches are either inefficient as they treat the stack as a whole, making optimization steps computationally expensive [18] or are restricted to a predetermined number of layers [19,20]. Many design algorithms result in complex geometries not suitable for fabrication [21–23]. None of the mentioned methods allow for a variable number of layers while simultaneously optimizing the metasurfaces within the stack. Recently, Li et al. [24] proposed a method allowing a variable number of layers by coupling the scattering matrix theory of multi-layer optics with a genetic algorithm. By using a genetic algorithm during optimization the problem of non-uniqueness was circumvented. Here, we propose a method that efficiently solves the inverse design problem of stacked metasurfaces. We link a metasurface database of Fourier modal method (FMM) simulated data with a modified genetic algorithm (GA) and apply the semi analytical stacking algorithm (SASA) [25]. The algorithm architecture of the modified GA enables us to fully control the geometrical properties of the metasurfaces constituting the stack while tackling the problem of non-uniqueness. By using SASA which operates in fundamental mode approximation (FMA), stacks with an arbitrary number of layers can be found in a matter of seconds. The FMA requires us to set a minimal distance between the metasurfaces which enables us to return results of reliable accuracy even for a high number of layers. The core of our algorithm, is the modified GA featuring two instead of one dimensional chromosomes. While one dimension represents the individual layers and can grow and shrink dynamically during the optimization process, the other dimension optimizes the features within the layers. Therefore, the number of layers and their properties are optimized simultaneously. The metasurface properties include e.g the shape of the metaatoms, their period, width, material, and any other desired parameters. Hence the 2D-encoding enables us to optimize the metasurface parameters individually, resulting in a more structured approach compared to guessing arbitrary geometries.

2. Implementation

The objective of our algorithm is the efficient optimization of layered metasurfaces with simple structures for easy fabrication. Opposed to other optimization approaches [19,20], our algorithm optimizes the number of layers dynamically and, simultaneously, the geometry parameters of the layers within the stack are optimized. The optimization of metasurface stacks with our algorithm is based on three numerical tools. The first is a rigorous FMM [26,27] which provides the scattering matrices of the 2D metasurfaces. This is done in preparation and decoupled from the optimization loop itself. The second tool is a modified GA, which is handling the optimization process. The third numerical tool is a SASA [25] for layered metasurfaces. An overview of the whole procedure and how the numerical tools work together is summarized in Fig. 1.

 

Fig. 1. Algorithm for the optimization of multilayered stacks. First a set of random individuals is created and encoded such, that each metasurface within the stack can get decoded into its corresponding scattering matrix , which is previously calculated via FMM and stored in a look-up table (left). Then, using SASA, the scattering matrix response of the full stack is generated. After that, the typical GA procedure of fitness assessment, mutating, breeding and selecting individuals follows. For the crossover, first the chromosomes (metasurfaces) are exchanged along a random crossover point among the stacks. Then, the arithmetic mean for the continuous parameters is evaluated. Since the individuals can be of different length the mean is only considered for the number of layers of the shorter individual (right).

Download Full Size | PDF

2.1 Semi analytical stacking algorithm (SASA)

SASA treats metasurfaces of identical meta-atoms, with dimensions and distances at subwavelength scale. By using SASA we enable fast simulation of any combination in a set of pre-computed scattering matrices representing metasurface geometries. It handles the underlying physics and is based on the scattering matrix theory of multilayer optics in FMA [28], in the case of periodic metasurfaces meaning that only the zeroth diffraction order is considered. For the FMA to be valid, coupling of the evanescent higher orders with the next metasurface has to be ruled out [29], thus as a limit for the FMA to hold. We require for propagation in $z$ direction that the modulus of the amplitude of the evanescent first diffraction order at distance $z = L$ decreases sufficiently, specified by an attenuation factor $\sigma$. From this limitation the critical distance $d_{\mathrm {crit}}$ can be estimated as the smallest distance for which the FMA holds. This critical distance depends on the wavelength $\lambda$, the period $\Lambda$ of the metaatoms, and the refractive index $n$ of the embedding material,

The attenuation factor $\sigma$ relates to the necessary precision of the semi-analytic model and is contingent upon the accuracy of the measured and simulated data. To access properties of the multilayer stack in FMA for a given wavelength, we need $\lambda <<n\Lambda$. The far-field response of a metasurface is fully contained in the zeroth diffraction order in transmission and reflection [28], and can be expressed via the scattering matrix [30] as

It connects the incoming $\mathbf {E}_{in}$ and outgoing $\mathbf {E}_{out}$ complex two-component electric field vectors in the front (f) and the back (b) of the metasurface (superscripts). That means that e.g in $\hat {\mathbf {T}}^{fb}$, connects the light that comes from the front and is transmitted to the back of the metasurface, while $\hat {\mathbf {R}}^{ff}$ connects the light that comes from the front and is reflected back. In particular $\hat {\mathbf {T}}^{ij}$ and $\hat {\mathbf {R}}^{ii}$ are Jones matrices, i.e 2x2 block matrices. The Stacking procedure within SASA is achieved by applying Redheffer’s start product [31]. In optics we use it to elegantly calculate propagation through to adjacent layers, taking into account inter layer reflections that would be excluded in a conventional transfer matrix approach. The strong simplification we achieve comes from the physical approximation to a 4-port system by use of the fundamental mode approximation. This is the core of our implementation with SASA. This procedure provides a significant speedup compared to rigorous simulations treating the stack [25] as a whole. In particular, rotation of individual metasurfaces with respect to each other are simple rotation matrix operations. Moreover, adjusting the distances in-between metasurfaces using a homogeneous medium of thickness $d$, refractive index $n$, and free space wave-number $k_0 = 2\pi /\lambda$ is represented by the propagator

2.2 Genetic algorithm (GA)

A genetic algorithm (GA) is a metaheuristic [33] optimization approach usually used on discrete data, inspired by the biological process of evolution by mimicking the Darwinian theory of the survival of the fittest. In analogy to genetics the data is usually encoded in a binary string, called the individual. A group of individuals is referred to as population. Each individual has a set of properties, its genotype which is a representing a point in optimization space. Each property is referred to as a gene. The optimization process is motivated by selecting the fittest individuals based on a fitness function determined by the objective of the problem. By randomly mutating a genotype, combining the genotype of different individuals and then selecting the fittest individuals in each generation, the population can be improved iteratively. The strength of a GA lies within its simplicity, flexibility, and being able to solve complex problems rapidly. As a metaheuristic algorithm it doesn’t require tightly sampled data or a smooth optimization space and can optimize many parameters simultaneously, without significant additional computational resources [34]. In the following we highlight how our GA-implementation differs from a common GA.

2.2.1 Encoding of the parameters in a layered metasurface stack

We encode the layered metasurface stacks as displayed in Fig. 1. Instead of encoding the data in binary 1D-arrays, we use 2D-arrays containing integer values. Each column (chromosome) represents the geometric parameters of one layer within the stack, each row (gene) represents one specific geometry parameter, e.g the metasurface parameters as height, width and period of the meta atom inclusions or the height of the dielectric layer. Spacer height and metasurface rotation can be adjusted continuously without significant computational cost, according to Eq. (3), or respectively by adding a rotation matrix to the respective metasurface. The 2D-encoding provides a practical alternative to the standard 1D-encoding in GA’s. Apart from making the result more human readable it also provides the benefit of accessing the geometry parameters independently for individual layers. Enabling the algorithm to adjust each parameter separately i.e. during mutation just one random geometry parameter in a layer is adjusted. Furthermore, this offers the possibility to use unique mutation rates for different geometry parameters, treat some continuously, others discretely, and make use of physical assumptions connected to the geometry parameters e.g during the fitness evaluation or selection process.

2.2.2 Fitness assignment

During the fitness assignment each chromosome is decoded into its corresponding scattering matrix. The scattering matrix of the metasurface is accessed via a look-up table containing the whole set of scattering matrices of possible choices for the metasurface geometries. The look-up table itself is built from rigorous FMM simulations of single metasurfaces. The spacer below the metasurface is computed according to Eq. (3). Then, using SASA, all layers are connected via the Redheffer star product returning the scattering matrix $\mathbf {S}_{\mathrm {full}}$ of the entire stack. Depending on the objective, an appropriately designed fitness function returns a real number proportional to the individual’s fitness. The fitness function is not limited to single objectives but can be designed for multifunctional applications.

2.2.3 Selection

The choice of selection method greatly impacts the convergence behaviour of the algorithm. Which selection method is suited best for the problem depends on the the objective and optimization space [35]. Our GA features a few of the most common selection approaches; namely Roulette Wheel, Elite, Linear Rank, Exponential Rank, and Tournament Selection [36].

2.3 Mutation

During the mutation process the algorithm can drop or spawn layers of random properties within the set of FMM simulations by a small random chance (growth mutation rate). Using 2D-encoding allows the optimization of each geometry parameter with different mutation rates. The mutation rate must be picked with great care. Higher mutation rates tend to converge quicker but to a worse overall fitness, if it is set too high the population might not converge at all. However, the mutation rate should be large enough, such that the individual is not trapped in a local optimum. This issue well known for genetic algorithm and various attempts have been made to solve it, one possible solution to overcome the problem are e.g adaptive mutation rates [37].

2.3.1 Crossover

First, two parents with a potentially different number of layers $I, J$ are selected from within the population by random chance. Then a simple one-point crossover is performed swapping a set of chromosomes (layers) among the parents at one random position. Finally, we perform an arithmetic crossover for the continuous parameters, taking the mean between them as visualized in Fig. 1.

2.3.2 Reduced complexity by dynamic growth

For GA the convergence speed depends mostly on the size of the population and the size of optimization space i.e the number of layers and available metasurfaces. For layered metasurface stacks, the size of the optimization space depends on three factors, the number of layers, the number of available metasurfaces within the lookup table and if the ordering within the stack is of relevance. When constructing the look-up table one has to choose carefully on which metasurfaces to simulate. They define in what range the algorithm can operate, for a discrete solution they can be sampled sparsely but that also means that we might miss well performing metasurfaces, a very high number however decreases the chance to randomly pick a suitable metasurface, if the optimization space is sampled fine enough such that interpolation can yield reasonable results, we would suggest using the continuous optimization for the respective parameter. In general, for a discretely optimized stack where the ordering is relevant, there are

3. Broadband circular polarizer

Chiral metamaterials have attracted increasing attention due to their potential as polarizers, sensors, and detectors. Attempts to achieve chirality include complex three-dimensional metaatoms and bulk materials usually challenging to fabricate and computationally expensive to optimize for [39,40]. Zhao, et al. [10] showed how identical structures of rectangular inclusions can exhibit strong, broadband bianisotropic effects when stacked under rotational twist, constituting a transmissive broadband circular polarizer, where the polarization performance scales with the number of layers. We want to mention that the fabrication effort increases roughly linearly with the number of layers. Using typical e-beam lithography and, say for metallic particles, a lift-off process means one fabrication run per layer. So, fabricating N layers is comparable to fabricating N samples. The great benefit of our approach is an independence on layer alignment. The robustness to misalignment, which may usually be problematic in the fabrication of layered metasurface stacks stems from the fact that the metasurfaces are coupled through the fundamental mode. That means that opposed to effects based on the interaction with light at the inclusion level, the metasurfaces are conceptualized as homogeneous layers. This means that, although each layer requires structuring and so forth, it can be fabricated neglecting the difficulty of alignment, making the process much simpler and reproducible. More so, the idea of combining layers comprised of simple nano-structure geometries allows us to achieve comparable optical functionality as with complex 3D-structures. This is why we chose the circular polarizer as an example. Fabricating chiral 3D-structures like helices poses significantly more fabrication challenges, and parametric dependencies, than multiple layers of twisted nano-wire. The limitations of decreasing the number of layers is always case depended. In our twisted nano-wire stack example, the lower limit of layer numbers is two in order to achieve sufficient asymmetry [7]. Here we optimized a plasmonic layered metasurface stack as circular polarizer in analogy to [10] and then increase the degrees of freedom to see how much the performance can be increased. We use metasurfaces with rectangular gold metaatoms to keep the design simple and easily fabricable. The optimization process is performed on a standard desktop computer (3.00 GHz Intel Core i5-8500 CPU and 16 GB RAM).

We optimized for a transmissive broadband circular polarizer in forward direction. As we are working with circular polarization we need to transform the respective Jones matrix within the scattering matrix, into polar coordinates

 

Fig. 2. a) Structure and reference for geometrical parameters. For the circular polarizer we used periodic rectangular gold meta atoms embedded in $\mathrm {SiO}_2$ separated by $d=100$ nm $\mathrm {SiO}_2$ spacers. The Layers are arranged with a rotational twist of $\theta$ towards each other, b) optimization space for a simplified 4 layer circular polarizer equally spaced and with constant rotational twist in between consecutive layers. The red line indicates the trajectory of a GA optimization run, c) left: transmission spectra of the optimized stack evaluated with SASA, the grey line marks the difference between the spectral components, right: transmission spectra of the optimized stack evaluated with rigorously suing FMM.

Download Full Size | PDF

Here an inter-metasurface distance of $d=100$ nm and rotational twist of $\theta =52^\circ$ leads to the best performing device according to our fitness function Eq. (7). We quickly converge to the optimum in around 9 generations and 9s for a population size of 30 (averaged over 10 iterations of the optimization). Figure 2(b) shows the optimization space and result including the trajectory through optimization space indicated in red for one exemplary run, converging in only 6 populations to the global optimum (the same optimium found in Zhao et al. [10]). Figure 2(c) compares the rigorous (left) with the semi analytical (right) approach, we can see that the general shape of the spectra is similar. Very notable is however, the lowered transmission intensity for the rigorous case, which is most likely due to losses into higher orders neglected in the FMA as we are operating close to SASAs limits. These effects can be diminished by increasing the height of the dielectric spacers in order to decrease near-field coupling or by decreasing the meta atom periods eliminating high-order diffraction. For $d=100$ nm and $P=300$ nm the first diffraction order merely decreases to $\sigma \approx 20\%$ at $\lambda =700$ nm and to $\sigma \approx 14\%$ at $\lambda =1400$ nm as specified by Eq. (1). Next, to increase the degrees of freedom, we allowed our algorithm to pick unique rotation angles in between consecutive layers with a constant distance of $L=100$nm. Simultaneously, we optimize the geometrical properties of the metasurfaces, allowing rectangular metaatoms of different dimensions in each layer. The optimization parameters are summarized in Table 1. Even when excluding the continuous rotation, according to Eq. (4) there are over $P=9.8^{11}$ possible combinations to construct an ordered unique stack out of $N_{\mathrm {min}} = 2$ and $N_{\mathrm {max}} = 5$ layers and $M=250$ different metasurfaces. We ran the algorithm with a population size of 50 individuals for 300 generations. The whole optimization takes around 5 minutes on average, simulating $50\cdot 300=1500$ individuals over the whole wavelength range. To put this into perspective; an FMM run for a just a single 5 layer stack build from metasurfaces of the same period and just single wavelength point takes around 4.5 minutes using 16 Fourier orders on a cluster node. Here we treat stacks with different periods resulting in large unit cells which cannot be any longer easily treated via FMM. The algorithm converges to an optimum around a fitness value of 0.34. Thus the performance of the circular polarizer has improved by a factor of 2 compared to the simplified version as in the reference. The best result out of 30 runs is summarized in Table 2, the spectral response of the stack is illustrated Fig. 3, the inlet visualizes the final structure. We want to note that the SASA simulations are not easily verified by FMM as we treat inconsummerable periods. The algorithm converges to a 5 layer stack, the maximum as defined by the parameter range. This comes to no surprise since the performance of the broadband polarizer scales with the number of layers as discussed in [10]. The geometrical parameters show a symmetry with respect to the center of the stack. However, the overall symmetry in propagation direction is broken by a monotonously increasing twist of the metasurfaces. It is this helix-like asymmetry giving the stack its chiral character and thus its ability to circularly polarize. Interestingly, the stack shows systematic variation of the particle aspect ratio (lower towards the center), height (higher towards the center), inter-layer twist (lower towards the center) and period (lower towards the center). It shows similarities to 3D-helical structures studied intensively in the context of optical chirality [42,43].

 

Fig. 3. Spectra of TLL and TRR components and their difference for the optimized result of the broadband circular polarizer, the inlet visualized the resulting structure to scale. The period, aspect ratio, height and consecutive rotational twist decrease towards the center.

Download Full Size | PDF

Table 1. Optimization Parameters^a

View Table | View all tables in this article

Table 2. Optimization Result

View Table | View all tables in this article

4. Conclusion

We showed that we can use a modified genetic algorithm in combination with the semi-analytical algorithm SASA, to create a functional and efficient tool for the optimization of layered periodic metasurfaces. By employing a genetic algorithm, we overcome the problem of most other optimizers struggling with non-unique solutions. We can optimize any arrangement of periodic metasurfaces for any desired functionality that can be expressed in terms of scattering matrices in FMA for either single or multi objectives. By introducing 2D-encoding of the problem and allowing different mutation rates for each parameter, addressed them individually. During the inverse design process, the algorithm can tailor geometrical parameters of the metasurfaces while it adjusts the number of layers simultaneously. Furthermore, the algorithm can optimize each parameter either continuously or discretely. The approach is especially valuable if inconsummerable periods are treated or rotations of the constituting metasurfaces or are considered as it enables the optimization within seconds avoiding thousands of rigorous simulations in the process. Since we reference the previously in FMM simulated scattering matrices directly in a look up table we have full control of the resulting metaatom shapes, keeping metaatom geometries simple and thus easily fabricable to the current state of art. Using the example of a transmissive plasmonic circular polarizer, we were able to showcase the efficiency, and validity of the algorithm. In future the GA can be applied to a multitude of problems and refined using more advanced mutation or selection methods, enabling scientists to efficiently optimize muti layered metasurface stacks for complex applications.