Adaptive Genetic Algorithm for Optical Metasurfaces Design

As optical metasurfaces become progressively ubiquitous, the expectations from them are becoming increasingly complex. The limited number of structural parameters in the conventional metasurface building blocks, and existing phase engineering rules do not completely support the growth rate of metasurface applications. In this paper, we present digitized-binary elements, as alternative high-dimensional building blocks, to accommodate the needs of complex-tailorable-multifunctional applications. To design these complicated platforms, we demonstrate adaptive genetic algorithm (AGA), as a powerful evolutionary optimizer, capable of handling such demanding design expectations. We solve four complex problems of high current interest to the optics community, namely, a binary-pattern plasmonic reflectarray with high tolerance to fabrication imperfections and high reflection efficiency for beam-steering purposes, a dual-beam aperiodic leaky-wave antenna, which diffracts TE and TM excitation waveguides modes to arbitrarily chosen directions, a compact birefringent all-dielectric metasurface with finer pixel resolution compared to canonical nano-antennas, and a visible-transparent infrared emitting/absorbing metasurface that shows high promise for solar-cell cooling applications, to showcase the advantages of the combination of binary-pattern metasurfaces and the AGA technique. Each of these novel applications encounters computational and fabrication challenges under conventional design methods, and is chosen carefully to highlight one of the unique advantages of the AGA technique. Finally, we show that large surplus datasets produced as by-products of the evolutionary optimizers can be employed as ingredients of the new-age computational algorithms, such as, machine learning and deep leaning. In doing so, we open a new gateway of predicting the solution to a problem in the fastest possible way based on statistical analysis of the datasets rather than researching the whole solution space.

problem is provided by fitness functions, F(p). A fitness function assigns a figure of merit or a measure of goodness [8] to a specific individual, which is usually a positive number and called fitness. The goal is to find an individual with maximum possible fitness, or the optimal solution. It should be noted that, since GAs work with chromosomes, there should be a decoding step (p → q) before evaluation of the fitness function.
Genetic algorithms consist of three main phases [8], namely, i. initiation, ii. reproduction, and iii. generation replacement. In the initiation phase, a population of randomly generated chromosomes, an initial generation, is formed. This initial generation is then evolved based on the fitness function to the proceeding generations, including more fit individuals, through reproduction and generation replacement processes. Reproduction process, takes advantage of GA operators, including, selection, cross-over, and mutation. One of the most successful selection strategies is tournament selection [5,8], in which a sub-population of individuals are randomly selected from the current GA generation. These selected individuals then compete based on their fitness and the one with the highest fitness is selected. These selected parents (members of the current generation) then go through reproduction process to generate members of the next generation, or children. The most important GA operators in the reproduction phase are cross-over and mutation, depicted in fig.S-1(b). The standard cross-over (recombination) operator takes two parent chromosomes as input and generates two child chromosomes. In order to do so, a random position is selected along the chromosome string, and the segments before and after the selected positions are exchange between the parents (see upper panel of fig.S-1(b)), with a probability of p cross . On the other hand, mutation S-2 operator selects a random element in the chromosome string, with a probability of p mutation , and inverts it. The purpose of cross-over and mutation operators are recombining the genes between members of a population, and searching the parts of the solution domain that are not represented by the current members of the population, respectively, with the hope of generating more fit children. Finally, in the generation replacement of GA, the current generation of parents is replaced with the new generation of children, generated as explained above. It should be mentioned that, in this paper, population size, N pop , has been kept constant for all the generations in an optimization. Number of generations is represented by N gen . Also, we have used an elitist strategy [8] to ensure monotonic increase in the highest fitness values of consecutive generations. In the generation replacement phase, the individual with highest fitness value, in the current generation, is transferred to the next generation to realize such elitism. A schematic of three main GA phases is presented in fig.S-1(c). We have used unit-cell patterns from one of the applications presented in the companion paper, the binary plasmonic metasurface for beam-steering, to illustrate initiation, reproduction, and generation replacement phases of GA.
Flowchart of a conventional genetic algorithm is depicted in fig.S-2. The algorithm starts by initiation of the first generation with N pop random chromosomes. Then, it proceed by assigning fitness values to the members of the current population, by decoding the individuals from chromosome space to parameter space [8], and evaluating the fitness function. If the stop criteria (for instance, specific number of generations or a threshold for the highest fitness value) is met, the algorithm will terminated and the near-optimal solution is obtained. Otherwise, the selection strategy, tournament selection is this paper, is utilized to select eligible parents form current generation for recombination and mutation, with probabilities of p cross and p mutation , respectively. The current generation is then replace with the next generation in the generation replacement phase, taking advantage of elitism, as explained above.

S-2 Design Procedure for the Binary Pattern Reflectarray Application
It is a well-known concept that, for a 2-D beam-steering reflectarray, the distribution of the reflection phase-delay over the aperture of the array should have the following form [4]: where k 0 is the free-space wavenumber and (θ 0 , ϕ 0 ) is the targeted direction of the reflected beam. It is worth mentioning that, the phase distribution of eq.S-1 is periodic both in x and y with periodicities of Λ x = 2π/k x and Λ y = 2π/k y , respectively. The common practice in designing reflectarrays is to quantize the phase-delay distribution is eq.S-1 to a few levels and realize those phase levels by the metasurface elements. As an example, in section 3.1 the companion paper, we have targeted (30 • ,45 • ) as the direction of the reflected beam, and also considered 8 levels of quantization. The resulting quantized phase-delay distribution is plotted in fig.S-3(a), for 3 periods of the phase-delay distribution in both x and y directions. In order to realized the above-mentioned quantized phase distribution by a reflectarray, one needs to take advantage of metasuface elements with retardation phases equal to the quantization levels. In this work, we have employed the AGA technique to design 8 binary plasmonic unit-cells to realized the phase distribution of fig

Conventionl Leaky-Wave Antennas
Leaky wave antennas (LWAs) operate on the principle of diffraction. Conventional LWAs use periodic ridges (gratings) inscribed on top of slab waveguides to couple the guided modes of the slab to radiative modes of free-space. The periodicity (Λ) of the grating is chosen such that one of the diffraction orders (usually -1 indexed diffraction order) of the waveguide mode is directed towards a desired angle (θ 0 ). Assuming the modal index of the waveguide mode is given by (β), the periodicity can be expressed as Λ = λ 0 /(β/k 0 − sinθ 0 ). Note that β depends on the polarization of the waveguide mode. Accordingly, for a given Λ, waveguide modes of each polarization get radiated to different angles. Figure S-6 shows the possible (fixed) beaming angles for TE and TM polarized waveguides modes in dielectric slab of the LWA shown above as a function of periodicity of the grating. The angles are computed using the formula θ = arcsin(β/k 0 − λ 0 /p). As seen, the angular difference between the two polarizations is fixed by the periodicity of an LWA (and almost constant). The GA optimization methodology employed in this paper allows us to choose unrestricted beaming angles.

S-3.2 Another example of an optimized Dual-Beam LWA
In order to further investigate the capability of GAs to design an aperiodic dual-beam LWA introduced in section 3.2 of the companion paper, we have solved another set of targeted radiation directions, namely, +45

S-4 Material Models for Visible-Transparent Infrared-Emitting Metasurface
We have used the FDTD method to accurately characterize the metasurface unit-cell, corresponding to visible-transparent IR-emitting application in section 3.4 of the companion paper, in the IR regime. In order to incorporate material dispersion in FDTD, a dispersion model should be fit to the wavelength dependent permittivity data. Drude dispersion model is employed for silver (Ag) and Lorentz model for polyimide (PI), fused silica(SiO 2 ), and indium tin oxide(ITO). Drude dispersion model can be written as follows [12]: where ε ∞ is the relative permittivity at infinite frequency, ω p is the frequency of the Drude pole, and γ p is the inverse of the pole relaxation time [12]. Also, Lorentz model has the following form: where N poles is the number of Lorentz pole pairs, ∆ε p is the change in the relative permittivity due to the p'th pole pair, ω p is the undamped frequency of the p'th pole pair, and δ p is the damping coefficient [12].
We have employed ε ∞ = 5.0, ω p = 1.3666E+16, and γ p = 2.7332E+13 as Drude model parameters of Ag [11]. Moreover, parameters represented in table S-1 are used in the Lorentz models fit to PI, ITO, and SiO 2 [9] in IR. Relative permittivities of the materials of table S-1 are also plotted in fig.S-8.

S-5 More Details Regarding the Performed Optimizations in the Companion Paper
In this section, we present GA parameters and statistics corresponding to the applications presented in section 3 (Results and Discussion) of the companion paper. An introduction to GAs along with a brief discussion of its main parameters is given in section S-1 of this supporting document. Table  S-3 summarize the parameters that will be presented in this section with short descriptions. There are some well known rules regarding the choice of GA parameters, as follows. A rule of thumb is suggested in ref. [3] for the choice of the population size, N pop as: where l is the chromosome length, k, effectively, is the average number of bits per parameter, and χ = 2 for binary encoding. In our case l = N par , and k = 1. Therefore, N pop = O(2N par ). However, we have been able to achieve good convergence in most cases with N pop ≈ N par . Higher value for N pop leads to more genetic diversity [8], and therefore, faster convergence. it also implies, however, larger fitness evaluation time for each generation. From our experience, number of generations should be chosen as N gen = O(N pop ). For multi-objective optimization, we have split N gen equally among the AGA stages. Probability of crossover, p cross , is typically chosen in 0.6-0.9 range, with an optimal value of 0.7 for most of the problems [8]. Higher values of p cross results in faster search in the chromosome domain. Probability of mutation, p mut , should be chosen small, usually in 0.01-0.1 range [8]. The fitness function for this problem is defined as First and second exponential terms in eq.S-5 represent phase and amplitude objective sub-functions, respectively, and are plotted in fig.S-10. The constant factors in the exponents of the Gaussian sub-functions (-6.0 and -4.0 for phase and amplitude, respectively) are chosen carefully, based on S-13 the physical considerations and goals of the problem. The phase objective function has been chosen to have a sharp peak, giving rise to accurate optimized phases and fast GA convergence. Since one cannot expect perfect reflection amplitudes close to 1.0 for many retardation phases in this application, as will be discussed in section S-5.1.3, we have chosen a smoother Gaussian function as the amplitude objective. Using this amplitude objective function, one can expect high amplitude fitnesses even for amplitudes around 0.7-0.8.

S-5.1.2 Optimization Parameters and Statistics
GA parameters and statistics corresponding to the binary pattern reflectarray metasurface application in section 3.1 of the companion paper are presented in table S-4. Taking advantage of a GPU-enabled parallel FDTD (PFDTD) solver [6], we have been able to solve each GA individual (see section S-1) is 16 s. We have also taken advantage of a computer cluster equipped with GPU nodes to optimize all the 8 distinct phase-delays simultaneously, resulting in a speed-up factor of 8. Although it takes about 45 h to design this metasurface, with a specific set of goals, we have solved 80,000 distinct unit-cell patterns in the course of the optimization. It has been discussed in section 4 of the paper, that this available surplus data can be used to design metasurfaces with alternative goals, without further computational cost.

S-5.1.3 GA Convergence Plots
It has been mentioned in section 3.1 of the paper, that we have employed a two-stage adaptive procedure to optimize the unit-cells of the binary pattern metasurface. In the first stage, we just optimize phase-delay of the unit-cells, and in the second stage we optimize both phase-delay and magnitude of the reflection coefficient. Results of the adaptive optimization are presented in fig.S-11. As it can be seen in fig.S-11 (a), at the end of the first stage, optimal values for all of the targeted phase-delays have been achieved. In the second stage, we consider a small weight for the magnitude objective sub-function, in eq.S-5. Consequently, as it can be observed in fig.S-11 (b), we have been able to improve the reflection magnitudes for almost all of the targeted phase-delays, and also the optimized phases do not deviate too much from the optimal values achieved in the first phase of the optimization, and remain in a acceptable range from the targeted values.
Furthermore, GA convergence plots, showing maximum fitness values in each GA generation versus generation index are plotted in fig.S-12. A discontinuity can be observed in all the convergence plots at the 51st generation, where the fitness function is changed based on the adaptive optimization. It can also be seen in fig.S-12 that, for almost all of the targeted phase-delays, an improvement in the overall fitness is achieved at the end of the second phase of the adaptive optimization (except for φ = -45 • ). S-16

S-5.2.1 Fitness Function
Cross-correlation of two arbitrary functions v(t) and w(t), that can in general be complex-valued, is defined as [2] R The cross-correlation operator can be normalized as follows: so that it returns 1.0 if both function are the same and values less than 1.0 if they are different. For this application, the fitness function is defined as where f y and f x have the same form. For instance, f y can be written as In other words, f y (p) is the zero-lag normalized cross-correlation of |E y | and |E y | target .

S-5.2.2 Optimization Parameters and Statistics
GA parameters and statistics corresponding to the dual-beam aperiodic LWA, presented in section 3.2 of the companion paper, are presented in table S-5.

S-5.3.1 Fitness Function
We have employed the following fitness function for the compact birefringent all-dielectric metasurface in section 3.3 of the companion paper. Identical objective sub-functions f x and f y correspond to x and y components of the reflected field, respectively, and are plotted in fig.S-14. We have chosen the constant factor -6.0 in the Gaussian sub-functions carefully to assure sharp peaks at the targeted phase-delays, and therefore accurate optimized phase-delay for each component of the reflected field. S-18

S-5.3.2 Optimization Parameters and Statistics
Table S-6 presents the GA parameters and statistics corresponding to the compact birefringent alldielectric metasurface in section 3.3 of the companion paper. It is worth mentioning that, taking advantage of a computer cluster equipped with GPU-enabled nodes, we have been able to optimize all of the 26 distinct phase-delay combinations, corresponding to perpendicular components of the reflected field, simultaneously. This has provided us with a significant advantage, as, we have been able to design the whole birefringent metasurface in an optimization time equal to that of a single phase-delay combination (a speed-up factor of 26). Moreover, 32,500 GA individuals have been solved in about 27 h of optimization time, giving rise to a pool of surplus data that can be employed for other design goals, without further computational cost (see section 4 of the companion paper).     fig.S-16. As it has been discussed in the companion paper, we have employed a two-stage adaptive optimization procedure for this problem. We just optimize the metasurface based on the IR objective sub-function for the first 15 generations, and then consider a small weight for the visible objective in eq.(8) of the companion paper. The discontinuity in the maximum fitness plot of fig.S-16 at the 15th generation is due to the change in the fitness function based on the adaptive optimization. As it can be seen, an improvement in the overall fitness value in the second phase of the optimization has been achieved. Maximum values of the IR and visible sub-functions in each generation are also plotted S-23 in the figure. In the first phase of the adaptive optimization, maximum overall fitness is equal to the maximum IR fitness. However, in the second stage, maximum values of f IR and f vis do not necessarily correspond to the maximum overall fitness. As, an individual might possess the maximum f IR value, but not necessarily the highest overall fitness.