Open Access
Add to BibSonomy
Abstract
Many single-layer dielectric metasurfaces studied today offer optimal performance by maximizing the refractive index contrast between a meta-atom and a background material, thereby increasing the field confinement and enhancing the optical resonance. Such architectures typically utilize high-index dielectric meta-atoms contrasted against an air background. Patterning can be done through deposition and growth of high-index materials or etching designs into a high-index layer. Unfortunately, such structures may be easily damaged, are susceptible to moisture accumulation, and cannot be readily incorporated into multilayer designs. To improve upon the mechanical and environmental stability of conventional dielectric-in-air metasurfaces, we show that dielectric metasurfaces in a solid host matrix can still achieve large enough resonances to exhibit high reflectance over a broad spectral band. Particle swarm optimization is used to discover different broadband reflector designs, each with different performance advantages including ultra-wide broadband reflectance and polarization independence. All designs exhibit a reflectance of R ≥ ~99% over the range of at least Δλ/λ0 ≥ ~20% centered at the wavelength λ0 = 1.55 µm. Finally, a decoupled method approximation is employed to predict the spectral response for a multilayer stack of metasurfaces, but with the added benefit of a lower computational cost.
© 2021 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
For centuries the bulk properties of a material, such as refractive index, dictated the spectral responses of optical devices. However, with the advancement of optical metamaterials, researchers can now fabricate devices with structures on size scales less than the wavelength of incident light, giving rise to unique optical properties determined by the geometric structuring. By varying the size, shape, material, and periodicity of these "meta-atoms", researchers can selectively fine-tune different resonances (Mie, Fano, Fabry-Perot, Bloch, guided, leaky, etc.) within these materials for specific applications [1–5]. By arranging meta-atoms in a 2D planar fashion, one can create a metasurface with exotic properties not seen in traditional thin films, while also keeping the thickness of the device small enough to exhibit lower optical losses than traditional multilayer Bragg stacks.
The ability to tailor the optical responses achieved by metasurface manipulation is nearly limitless. In practice, however, exploiting this tunability remains a challenge because identifying a specific design with the target optical response is a non-trivial task. Metasurfaces suffer from the curse of dimensionality, insofar that the multidimensional design space cannot be exhaustively searched by brute force due to the sheer number of possible designs [6–8]. For this reason, many researchers have leveraged optimization and machine learning techniques to find optimal metasurface designs, all with different advantages and disadvantages [9–13]. For example, artificial neural networks (ANNs) have been used in various works to link metasurface designs to optical responses [6,8,10,14,15]. ANNs can be used to create a library of designs [16,17] or act as surrogate models in other optimization methods [14]. They can also be constructed in an inverse configuration to directly design metasurfaces from their optical responses [18–20], though this approach can be limited by the many-to-one problem for many metasurface classes. However, ANNs tend to require much larger datasets than many gradient-based and evolutionary optimization methods, potentially mitigating any benefits gained from their implementation [9,11,17,19,21].
Other optimization methods used to design metasurfaces and other optical devices include genetic algorithms (GA) [22], covariance matrix adaptation evolutionary strategy (CMA-ES) [9], Bayesian optimization [23], gradient-based adjoint methods [24], coupled mode theory [25], and topology optimization (TO) [26]. Disadvantages of some of these techniques stem from their complexity, with numerous hyperparameters which may require additional optimization. Moreover, the user must be diligent when selecting hyperparameters, lest the method fail to converge to an optimal result. While there is no single ideal algorithm for any arbitrary problem, particular methods are better suited depending on context. We applied the heuristic found in Fig. 1 of Ref. [9] to identify particle swarm optimization (PSO) as a strong candidate for metasurface design. PSO is a gradient-free method that is well-suited for design problems with a limited number of design parameters and where a good starting point is not known. Its implementation is relatively simple and intuitive, and it can be generalized to problems involving discretized parameter spaces and non-differentiable cost functions.
Particle swarm optimization is a type of swarm-intelligence algorithm that mimics the collective behavior seen in flocks of birds and schools of fish [27]. The decentralized system can self-organize to discover the globally optimal solution to a problem, with the proper interplay of exploration and exploitation [28]. While some swarm-intelligence algorithms, such as predator-prey optimization, accomplish this through a competitive nature [29,30], PSO uses a cooperative nature, similar to those behaviors seen in ant colony optimization [31,32] and artificial bee colony optimization [33,34]. PSO contains a collection of particles in a swarm, where the position of each particle is an $N$-dimensional vector that represents all $N$ parameters of the design space to be optimized. As the particles collectively move around the design space, their positions will converge to the optimal design. We can guarantee convergence by selecting hyperparameters that satisfy the inequality in Ref. [27]. This, as well as the algorithm’s simplicity and effectiveness, helped us decide to use particle swarm optimization (PSO) to discover our durable and scalable metasurface reflectors [14,35,36].
In this manuscript, we use PSO to find a broadband metasurface reflector centered at $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$ with a reflectance of $R \geq 99\%$ over the band $\Delta \lambda / \lambda _0 = 20\%$. Beyond applications in displays, such broadband capability has applications for lasers, photodiodes, and spectral filters [3,37–42]. We compare our results to an optimized silicon (Si) open-cylinder array on silica (SiO$_2$) metasurface discovered through use of ANNs in Ref. [19]. In Ref. [19] the Si open-cylinders are arranged in a square lattice with spacing $\Lambda = {0.758}\mathrm{\mu} \rm {m }$, outer and inner radii of $r_\textrm {o} = {0.279}\mathrm{\mu} \rm {m }$ and $r_\textrm {i} = {0.077}\mathrm{\mu} \rm {m }$, respectively, and open-cylinder height of $h = {0.593}\mathrm{\mu} \rm {m }$. However, the structure of our metasurface designs differ from Ref. [19] for two significant reasons: (1) to improve the durability of the device and (2) to increase the scale of manufacturing to potential production levels. We have observed considerable damage to freestanding Si-based metasurfaces after cleaning with solvent and lens tissue. Embedding nanostructures should enable conventional cleaning processes, improve robustness to external physical damage, and prevent moisture within and between meta-atoms, while also allowing the stacking of multiple metasurface layers. In an effort to allow scalable manufacturing of these designs via roll-to-roll fabrication [43–46], we also replace the silica (SiO$_2$) substrate with polydimethylsiloxane (PDMS). PDMS can be stamped or embossed with the inverse of a metasurface design and then backfilled with Si, rather than using the conventional, scale-limited methods of Si deposition, lithography, and etching [47–50]. We investigate three potential Si metasurface designs embedded in PDMS: (1) a 1D grating, (2) a 2D open-cylinder array, and (3) two stacked 1D gratings arranged in a cross-hatched pattern. We also introduce a decoupled method approximation (DMA) that can be used to decrease simulation times for multilayer stacks.
2. Optimization and modeling of metasurfaces
2.1 Metasurface designs and simulation method
Three proposed metasurface designs are shown in Fig. 1, including (a) a 1D Si grating embedded in PDMS (GE1), (b) a 2D Si open-cylinder array embedded in PDMS (OCE), and (c) two Si gratings embedded orthogonal to each other at the two surfaces of a thin PDMS film (GEX). OCE allows us to directly compare embedded open-cylinder spectral performance to that of non-embedded open-cylinders (OC) given by Ref. [19], while GE1 decreases computation time at the cost of polarization dependence, and GEX takes advantage of attributes from both OCE and GE1. To compare the effectiveness of these embedded designs to a non-embedded metasurface design, we use the Si-in-air open-cylinder array on SiO$_2$ mentioned in Ref. [19] as a reference.
Fig. 1. Metasurface designs with (a) a 1D Si grating embedded in PDMS, abbreviated as GE1, (b) a 2D Si open-cylinder array embedded in PDMS, abbreviated as OCE, and (c) two orthogonal (cross-hatched) 1D Si gratings embedded in the top and bottom of a PDMS film, abbreviated GEX.
For our simulations, light is incident normal to the patterned surface (z-direction), with the electric field aligned along either the x- or y-direction for p- or s-polarized light, respectively (Fig. 1). The material dispersion curves of Si, SiO$_2$, and PDMS were taken from Refs. [51], [52], and [53], respectively. The metasurface design parameters shown in Fig. 1 are the lattice constant, $\Lambda$; the grating width, $w$; the open-cylinder outer and inner radii, $r_\textrm {o}$ and $r_\textrm {i}$, respectively; the metasurface height, $h$; and the grating spacer thickness, $d$. We simulated the reflectance and transmittance using the Stanford Stratified Structure Solver (S$^4$), which uses rigorous coupled wave analysis (RCWA) [54]. This simulation method is the natural choice for periodic metasurfaces since it works in Fourier space and expands the plane wave in terms of spatial harmonics. We use subpixel smoothing as our Fourier modal method (FMM) formulation, with the number of harmonics set to $n_\textrm {G} = 121$ and 625 for the 1D (GE1) and 2D (OC, OCE, GEX) cases, respectively. The total number of spatial harmonics used in this study was based on Ref. [19] and determined from additional 1D and 2D convergence tests, where we increased $n_\textrm {G}$ until successive RCWA simulated spectra no longer differed.
2.2 Optimization procedure for GE1 and OCE
In order to find the optimal metasurface designs with broadband reflectance around $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$, we used a modified form of PySwarms, which is a free and open-source Python package with a general, high-level implementation of PSO and an API that enables extensive customization [55]. A flowchart of the optimization and design modeling procedures is shown in Fig. 2. The blue box on the left illustrates the PSO procedure, which was used to optimize each of the GE1 and OCE designs independently. The GEX design was not directly optimized via PSO, but instead used the optimized GE1 solution for both of its grating layers. The green box on the right shows additional exploration of this multilayer structure, including variation of the spacer thickness and application of the decoupled method approximation (DMA) for obtaining the properties of a multilayer device from single-layer calculations. For the rest of Sec. 2.2 we will describe the methodology of PSO shown in the blue box, while the procedures shown for GEX simulations (green box) and DMA (red box) will be discussed further in Sec. 2.3.
Fig. 2. Design optimization and simulation flowchart. The blue box shows the PSO algorithm used for both GE1 and OCE designs. The optimal GE1 design is used for several GEX designs with various spacer thickness (procedures shown in the green box). GEX designs with $d \in \{ {2}\mathrm{\mu} \rm {m }, {10}\mathrm{\mu} \rm {m }, {100}\mathrm{\mu} \rm {m }\}$ spectra are simulated via the DMA (red box) and RCWA methods. Additionally, GEX designs with $d$ normally distributed are also simulated.
In PSO, each particle’s position, given by the vector $\mathbf {X}$, represents a unique metasurface design, where each vector component is a parameter of the design space to be optimized. The algorithm starts by first randomly initializing a collection of particle positions bounded within our desired design space. For GE1 $\mathbf {X} = (\Lambda, w/\Lambda, h)$ is bounded by the intervals $\Lambda \in [{0.3}\mathrm{\mu} \rm {m }, {2}\mathrm{\mu} \rm {m }]$, $w/\Lambda \in [0, 1]$, and $h \in [{0.1}\mathrm{\mu} \rm {m }, {2}\mathrm{\mu} \rm {m }]$. For OCE $\mathbf {X} = (\Lambda, r_\textrm {o}/\Lambda, r_\textrm {i}/r_\textrm {o}, h)$ is bounded by the intervals $r_\textrm {o}/\Lambda \in [0, 0.5]$, $r_\textrm {i}/r_\textrm {o} \in [0, 1]$, and uses the same boundaries for $\Lambda$ and $h$ used for GE1. For each time step in PSO, the particle’s position is updated by a velocity vector ($\mathbf {V}$), with $\mathbf {X}$ and $\mathbf {V}$ given by
The inertia weight $\omega = 0.707$ controls the exploration aspect of the PSO algorithm. Parameters $u_1, u_2 \sim \mathcal {U}(0, 2.04)$ are randomly chosen from a uniform distribution, where $u_1$ and $u_2$ are the cognitive and social parameters, respectively, and control the exploitation aspect of the PSO algorithm. The value of $\omega$ and the limits of $u_1$ and $u_2$ were taken from Ref. [27]. The vector $\mathbf {X}_\textrm {pbest}$ is the particle’s personal best position, and $\mathbf {X}_\textrm {gbest}$ is the swarm’s global best position, both of which are determined by an objective cost function that depends on the particle’s position (i.e., metasurface geometry). The subscript $i$ represents the i$^\textrm {th}$ particle, while the subscript $j$ represents the j$^\textrm {th}$ design parameter, and $t$ is the current time step index in the iteration. The velocity is clamped such that $\mathbf {V} \leq 0.3(\mathbf {X}_\textrm {upper} - \mathbf {X}_\textrm {lower})$. Updated particle positions that fall outside of the design space boundaries are moved back into the search space via the reflection method [56,57] and the velocity is set to zero. We use swarm sizes of 100 particles and run each PSO through 50 iterations, simulating 5000 total spectra for each design (significantly fewer than the tens of thousands required to provide training data for an ANN).
In order to determine the goodness of a design that most closely matches our desired performance, which is $R \geq 99\%$ over a band of $\Delta \lambda / \lambda _0 = 20\%$ centered at $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$, we scalarize our multi-objective cost function by summing an inverted 2D Gaussian function given by
This is the objective cost ($C$) that we seek to minimize for our design problem. Each spectrum consists of wavelengths ${1.395}\mathrm{\mu} \rm {m } \leq \lambda \leq {1.705}\mathrm{\mu} \rm {m }$ with a step size of ${0.01}\mathrm{\mu} \rm {m }$. From Eq. (3), wavelengths closer to $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$ contribute more to the cost and thus play a bigger role in the minimization process. Likewise, reflectances closer to the maximum value of unity also contribute more to the minimization. The standard deviation of our 2D Gaussian objective function is determined by our desired performance metrics such that $\sigma _{\lambda _0} = 0.5$, $\Delta \lambda = 0.1 \lambda _0$ and $\sigma _{R} = 1\%$, with the latter softly enforcing $R \gtrapprox 99\%$. After each iteration in the PSO procedure (Fig. 2), each particle’s personal best position and cost, $\mathbf {X}_i^\textrm {pbest}$ and $C_i^\textrm {pbest}$, are replaced by the particle’s current $\mathbf {X}_i$ and $C_i$ if the latter outperforms the former. The swarm’s global best $\mathbf {X}^\textrm {gbest}$ and $C^\textrm {gbest}$ are also updated accordingly.
2.3 Direct and approximate modeling of GEX
The advantage of simulating GE1 (1D grating) over OCE (2D array of open cylinders) is that GE1 is invariant in the y-direction, allowing us to compute the electromagnetic interaction across the xz-plane only. This is in contrast to OCE, which requires simulations of the entire xyz-volume. Consequently, GE1 requires fewer spatial harmonics in the RCWA simulation, and since the computation time scales with the cube of the harmonics used, GE1 can be simulated roughly $140 \times$ faster than OCE with the same level of accuracy. However, the reflectance spectrum of the 1D grating GE1 exhibits a strong polarization dependence due to its asymmetrical shape in the xy-plane (Fig. 1(a)) [58]. In contrast, while the 2D open-cylinder array OCE takes significantly longer to reach a converged solution, its reflectance can be polarization independent.
Thus, to take advantage of the fast computation speeds for 1D designs while also obtaining a polarization-independent response, we simulate the effects of adding a second 1D grating embedded in the bottom interface of the PDMS film, with the periodicity perpendicular to the top grating (Fig. 1(c)). This is similar to Refs. [59,60], but has embedded features that can be fabricated in roll-to-roll processing [43–45]. As mentioned in Sec. 2.2, GEX is not directly optimized using PSO, but instead uses the PSO-optimized GE1 design for both the top and bottom gratings. Since the optimized GE1 design reflects p-polarized light at $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$ across $\Delta \lambda / \lambda _0 = 20\%$, a second identical grating orthogonal to the first will exhibit a similar reflection profile, but instead for s-polarized light. No additional optimization is necessary.
With GEX we simulate the effects of various spacer thicknesses, $d$, on the reflectance profiles. These steps are shown in the GEX portion of Fig. 2. For these simulations, the three spacer thicknesses between the two embedded gratings are $d \in \{ {2}\mathrm{\mu} \rm {m }, {10}\mathrm{\mu} \rm {m }, {100}\mathrm{\mu} \rm {m } \}$. These values were chosen so that we could observe the effects of $d \gtrapprox \lambda$ and $d \gg \lambda$. We expect GEX to behave similarly to GE1, except with polarization independence. In practice, there will be deviations from the expected reflectance spectrum due to thin film Fabry-Perot interference caused by the spacer layer. Nevertheless, we will show in Sec. 3.3 that the total reflectance of GEX can be treated as a product of individual layer stacks via the decoupled method approximation (DMA).
The decoupled method approximation ignores the phase that light acquires as it propagates through each layer. DMA calculates the total reflectance from (or transmittance through) two interfaces by summing all internal reflections between (or transmissions through) the individual interfaces, while also accounting for the absorption of the material. The sum creates a geometric series that converges to
DMA is typically only valid when $d \gg \lambda$, since this will lead to high-frequency Fabry-Perot oscillations in the spectrum that cannot be resolved experimentally due to limitations in detector resolution. Nevertheless, this method can be valid for even smaller values of $d$ due to light incoherence, which can result from small imperfections in the metasurface and/or layer interface roughness. To mimic incoherence caused by layer roughness, we randomly sample 1000 spacer thicknesses from a normal distribution described by $d \sim \mathcal {N}(d_\textrm {avg} \! = \! {100}\mathrm{\mu} \rm {m }, \, \sigma _d \! = \! {0.175}\mathrm{\mu} \rm {m })$, and average over the response reflectances. We show in Sec. 3.3 that this destroys the high-frequency Fabry-Perot oscillations seen in the reflectance and converges to the decoupled method approximation.
3. Results and discussion
3.1 Optimization convergence
Convergence of the PSO algorithm for the GE1 and OCE optimizations is illustrated in Fig. 3, which shows the decrease in cost (i.e., the value obtained by evaluating Eq. (3) at a particle’s position in the design space) during the optimization. The cost distribution for all particles in each swarm is shown as quartiles, with the solid colored line showing the median ($Q_2$) and the dark and light shaded regions showing the interquartile (IQR: $Q_3 - Q_1)$ and $\textrm {max} - \textrm {min}$ ranges, respectively. The plots in Fig. 3(a) show the current (top row) and personal best (bottom row) costs of the particles at each iteration for GE1, while the plots in Fig. 3(b) show the same quantities for OCE.
Fig. 3. The objective costs of the swarm for (a) GE1 and (b) OCE. The top and bottom rows show the current and personal best costs of the particles at each time step, $C_i$ and $C_i^\textrm {pbest}$, respectively. The solid colored line shows the median cost (3) of all particles in the swarm, with the dark shaded region showing the interquartile range and the light shaded region showing the entire cost range. The dashed black line shows the swarm’s current global best cost, $C^\textrm {gbest}$.
Trends observed for the GE1 and OCE costs are similar. For both swarms most particle costs start at $C \approx 0$, the maximum cost allowed by Eq. (3), giving little variance in the distribution of costs. As time progresses, the variance grows as particles begin to explore and find designs with costs other than zero. The noise seen in $C_i$ is due to the stochastic nature of the particles’ velocities, specifically from the $u_1$ and $u_2$ terms in Eq. (2). Nevertheless, it is clear that the current costs $C_i$ on average decrease over time, as particles tend to move toward their personal and the swarm’s global optimal positions. While the median and IQR of the personal best costs, $C_i^\textrm {pbest}$ (bottom panels), show a smoother trend, the minimum values show step features, indicating that the swarm’s global best cost remains unchanged for several iterations until a particle’s $C_i < C^\textrm {gbest}$ and the global best is updated; similar step features can be seen for the maximum values representing the worst of the personal bests.
As the optimization continues, both the distributions of $C_i$ and $C_i^\textrm {pbest}$ show a shift in skewness, from larger to lower costs. For $t \leq 5$ there is a large spread of cost values between the minimum and the median costs ($> 5$), but a very small spread between the median and maximum costs ($< 0.1$). However, for $t \gtrapprox 40$ the distribution around the median becomes more symmetric for $C_i$, and even skews to the lower side for $C_i^\textrm {pbest}$, with a smaller spread of cost values between the minimum and median costs than that between the median and maximum costs. In addition, convergence is seen in multiple ways. As the IQR grows, but the $\textrm {max} - \textrm {min}$ range does not, the tailedness of the distribution decreases, suggesting fewer outliers and a larger concentration of similar cost values. Also, the $C_i^\textrm {pbest}$ curve flattens as it approaches the global best, $C^\textrm {gbest}$, implying that each particle’s personal best is approaching the global best. Finally, there is also a flattening of the $C^\textrm {gbest}$ curve as time increases, further indicating convergence to a global optimum.
3.2 Reflectance of optimized designs
The PSO-optimized GE1 and OCE designs are shown in Table 1, and their associated reflectance spectra are shown in Figs. 4(a) and (b), respectively. Figure 4(c) shows the reflectance for the GEX design comprising two optimized GE1 gratings stacked in an orthogonal configuration with a spacer separation of $d = {2}\mathrm{\mu} \rm {m }$. The left column shows the reflectance spectra on a linear scale, illustrating the strong broadband reflectance in an intuitive fashion. The right column shows the same reflectances, but instead $R$ is plotted on an inverted log scale. This logarithmic scale enables a more quantitative evaluation of the device performance, to assess whether the optimized designs meet the objective criteria for high broadband reflectance (i.e., reflectance values above the horizontal dotted line, $R \geq 99\%$, over the entire gray shaded region of the spectrum, $\Delta \lambda /\lambda _0 = 20\%$). These plots also clearly show the individual resonances responsible for creating the broadband features.
Fig. 4. Simulated reflectance profiles of the (a) GE1, (b) OCE, and (c) GEX PSO-optimized metasurface designs embedded in PDMS. The left columns show the reflectance spectra $R$ on a linear scale, whereas the right columns show the same $R$ on an inverted log scale. For all embedded designs, the yellow solid line represents s-polarization, the blue solid line represents p-polarization, and the yellow-blue dashed line represents s- and p-polarization overlap (i.e., polarization-independent reflectance). For reference, the polarization-independent OC reflectance is shown in all plots as a solid gray line. The horizontal and vertical black dotted lines show the desired minimum reflectance ($R \geq 99\%$) and the center wavelength of interest ($\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$), respectively. The gray shaded region shows the bandwidth $\Delta \lambda / \lambda _0 = 20\%$.
Table 1. Particle swarm optimized designs.
Figure 4(a) shows the reflectance spectra of the optimized GE1 design as solid colored lines; as a reference, the spectrum of the Si-in-air OC design from Ref. [19] is shown as a solid gray line. For p-polarized light, the broadband reflectance profile of GE1 (solid blue line) outperforms OC in terms of the high reflectance bandwidth ($R \geq 99\%$) and in-band average reflectance. Specifically, GE1 shows an average reflectance of $R = 99.789\%$ over the range $\Delta \lambda / \lambda _0 = 28.1\%$, while OC shows an average reflectance of $R = 99.775\%$ over the range $\Delta \lambda / \lambda _0 = 24.3\%$. However, the broadband reflectance is completely destroyed in the GE1 case for s-polarized light (solid yellow line), while OC is unaffected due to its azimuthal symmetry. From the right column of plots, three main resonances contribute to the broadband reflectance at $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$. It can be seen that while GE1 has a wider high reflectance band than the OC design for p-polarization, its peak reflectance values are slightly smaller. This is due to the smaller refractive index contrast between the Si grating and PDMS background in the GE1 device compared to OC with Si open-cylinders in air. This may lead to weaker field confinement within the metastructures and a more uniform field distribution throughout the entire metasurface film. Consequently, since PDMS is more absorptive than either Si or air, we would expect a drop in peak reflectance values for the Si-in-PDMS compared to the Si-in-air designs.
While the reflectance of the 1D Si grating in PDMS (GE1) shows a strong polarization dependence, polarization independence can be enabled through azimuthal symmetry, as in the Si open-cylinder array on SiO$_2$. Figure 4(b) compares the s- and p-polarization reflectances of the optimized OCE design (dashed yellow-blue lines) to that of the OC design (solid gray line). The OCE optimized design shows an average reflectance of $R = 99.757\%$ over the range $\Delta \lambda / \lambda _0 = 19.9\%$. As in the GE1 case, the lower index contrast between Si and PDMS (compared to Si and air) negatively affects the reflectance. Here, this difference results in both decreased resonant magnitude and smaller bandwidth. Compared to the GE1 device, the OCE device displays greater absorptance because it contains a larger ratio of PDMS to Si. The PDMS fill factor is $73.7\%$ for OCE and $37.1\%$ for GE1. Nevertheless, OCE still shows a bandwidth of $\Delta \lambda / \lambda _0 \approx 20\%$ with $R \geq 99\%$, and unlike GE1, its response is polarization independent. Compared to the OC design, OCE should exhibit enhanced durability and offers potential for facile, production-scale manufacturing.
Nonetheless, it is also possible to achieve the the ultra-wide broadband response displayed by GE1, while also incorporating the polarization independence seen in OCE. This can be achieved by embedding two 1D Si gratings in the top and bottom surfaces of a PDMS thin film, with both gratings orthogonal to each other (e.g., the top is periodic in x and the bottom is periodic in y). This result is demonstrated in Fig. 4(c), where the distance between the top and bottom gratings was set to $d = {2}\mathrm{\mu} \rm {m }$. The cross-hatched grating (GEX) exhibits a nearly polarization-independent broadband reflection (left plot of 4(c) on a linear scale), where the distinction between s- and p-polarization is only apparent on the inverted log scale. GEX shows a broader reflectance peak for both s- and p-polarizations (yellow and blue lines, respectively) compared to the dielectric-in-air OC design (solid gray line).
The spectrum of the GEX device is directly related to the response of the individual grating, GE1. In the inverted log scale plot in Fig. 4(a), when the light polarization is orthogonal to the GE1 grating’s lattice, the reflectance shows a broadband feature centered around $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$ with three resonances (GE1 p-polarization blue line). In contrast, when the light polarization is parallel to the GE1 grating’s lattice, the reflectance shows a single resonance around $\lambda = {1.25}\mathrm{\mu} \rm {m }$ (GE1 p-polarization yellow line), outside of the target band. Thus, for p-polarized light incident on GEX, the top grating creates three resonances that are responsible for the broadband reflectance and the bottom grating is responsible for the single-resonance reflectance, shown in Fig. 4(c) (GEX p-polarization blue line). This relationship is inverted for s-polarized light incident on GEX, with the top grating responsible for the single-resonance reflectance, and the bottom grating responsible for the three-resonance broadband reflectance, shown in Fig. 4(c) (GEX s-polarized yellow line). However, in the latter case, most of the s-polarized light must pass through both the first grating and the PDMS spacer layer before it is reflected from the bottom grating. This results in additional absorption and Fabry-Perot effects, explaining the differences between s- and p-polarization seen on the inverted log scale of Fig. 4(c).
Based on the simulated spectra shown in the inverted log scale plot of Fig. 4(c), it should be possible to re-optimize GEX so the resonance at $\lambda = {1.25}\mathrm{\mu} \rm {m }$ (polarization parallel to grating orientation) is closer to the three in-band resonances (polarization orthogonal to grating orientation), broadening the high reflectance band even further. With the current design, obtained from the GE1 PSO, GEX shows an average reflectance of $R = 99.716\%$ and $R = 99.769\%$ over the range $\Delta \lambda / \lambda _0 = 31.2\%$ for s- and p-polarizations, respectively, excluding the $\lambda = {1.25}\mathrm{\mu} \rm {m }$ s-like resonance. Including the additional resonance, GEX has an average reflectance of $R = 99.711\%$ and $R = 99.761\%$ over the range $\Delta \lambda / \lambda _0 = 38.8\%$ for s- and p-polarizations, respectively, though this value drops to $R \approx 78\%$ around $\lambda = {1.3}\mathrm{\mu} \rm {m }$.
3.3 Approximation of the multilayer metasurface
Figure 5 compares the spectra of the GEX multilayer design obtained from full-device RCWA simulations and from the DMA approach (using RCWA simulations of the individual layers). These results show the validity of using DMA, even when the spacing between metasurfaces is comparable to $\lambda$. Figure 5(a) shows the spectra of the same GEX design as in Fig. 4(c), comparing the results obtained from RCWA and DMA. Although $d = {2}\mathrm{\mu} \rm {m }$, a value only five times greater than $\lambda _0 / n_\textrm {PDMS}$, the predicted decoupled reflectance still closely approximates the true reflectance profile of this device. As the spacing between metasurfaces increases, the Fabry-Perot oscillations increase in frequency and become more visible in the RCWA spectra, though these are not included in the DMA spectra.
Fig. 5. The simulated reflectance spectra of GEX with interface grating spacer thicknesses of $d =$ (a) ${2}\mathrm{\mu} \rm {m }$, (b) ${10}\mathrm{\mu} \rm {m }$, and (c) ${100}\mathrm{\mu} \rm {m }$ are shown on an inverted log scale. The left column is for s-polarized incident light (warm colors), while the right column is for p-polarized (cool colors). The red and blue lines show the RCWA simulation in both cases, which includes Fabry-Perot interference from coherent light. The yellow and light blue lines show the DMA in both polarization cases, which assumes incoherent light. The black dashed and dot-dashed lines in the (c) ${100}\mathrm{\mu} \rm {m }$ spacer case illustrate the average $R$ of 1000 RCWA simulated spectra, where the spacer thicknesses $d$ were taken from the normal distribution shown in the inset. The horizontal and vertical black dotted lines show the desired minimum reflectance ($R \geq 99\%$) and the wavelength of interest ($\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$), respectively. The gray shaded region shows the bandwidth $\Delta \lambda / \lambda _0 = 20\%$.
Changing the spacer thickness also has some effect on the device performance, seen in both the RCWA and DMA reflectance. For $d = {10}\mathrm{\mu} \rm {m }$ the GEX design still exhibits $R \geq 99 \%$ over the desired minimum band $\Delta \lambda / \lambda _0 = 20\%$. As the spacing is increased even further to $d = {100}\mathrm{\mu} \rm {m }$, the cross-hatched design still meets the objective criteria for p-polarization, but no longer for s-polarization. This is due to the top grating primarily reflecting p-polarized light, but transmitting s-polarized light, which is reflected from the orthogonal bottom grating. Consequently, s-polarized light has a round-trip optical path length of $2 d n_\textrm {PDMS} \approx {278}\mathrm{\mu} \rm {m }$ through the PDMS spacer, which increases absorptance.
The Fabry-Perot oscillations seen in the simulated RCWA spectra for large spacer thicknesses are not typically observed experimentally due to layer roughness, device imperfections, and detector resolution limits. To mimic incoherence caused by layer roughness, we incorporated the effects of small variations in the spacer thickness into an additional set of RCWA simulations. We randomly sampled 1000 spacer thicknesses from a normal distribution described by $d \sim \mathcal {N}(d_\textrm {avg} \! = \! {100}\mathrm{\mu} \rm {m }, \, \sigma _d \! = \! {0.175}\mathrm{\mu} \rm {m })$, and averaged the RCWA-simulated reflectance spectra for all designs. In the Fig. 5(c) inset, the normal distribution, $\mathcal {N}$, is shown as a black line and the random sampling taken from $\mathcal {N}$ is shown as a green histogram. The mean spacer thickness is $d_\textrm {avg} \! = \! {100}\mathrm{\mu} \rm {m }$ and the standard deviation is $\sigma _d \! = \! {0.175}\mathrm{\mu} \rm {m }$. The black dashed and dot-dashed lines in Fig. 5(c) show the average of all 1000 reflectance spectra. We can see that the high-frequency Fabry-Perot oscillations are destroyed and the average reflectance is in very good agreement with the DMA predictions. The distribution of spacer thicknesses results in spectral shifts of the high-frequency constructive and destructive interference fringes for each $d$, so that while the Fabry-Perot interference effects are no longer seen, the overall spectral shape of $d_\textrm {avg} = {100}\mathrm{\mu} \rm {m }$ is preserved. By ignoring these Fabry-Perot interference effects, DMA provides a relatively accurate approximation even for small spacer thicknesses ($d \gtrapprox \lambda$), and especially for large spacer thicknesses ($d \to \infty$) and incoherent light caused by device imperfections.
4. Conclusions
Using particle swarm optimization, three potential polymer-embedded metasurface designs were identified, exhibiting broadband high reflectance of $R \geq 99\%$ over the range $\Delta \lambda / \lambda _0 \geq 20\%$ centered at the wavelength $\lambda _0 = {1.55}\mathrm{\mu} \rm {m }$. In addition, all designs exhibit at least three individual resonances with $R > 99.9\%$ and an average reflectance of at least $99.7\%$ across the target spectral band. A single 1D Si grating embedded in PDMS showed a broader reflectance peak for p-polarized light than a 2D Si open-cylinder array embedded in PDMS. As expected, however, the 1D grating performance suffered greatly for s-polarization, whereas the 2D open-cylinder case exhibited polarization independence. Thus, we also demonstrated that two Si gratings can be embedded in a PDMS film—one at each of the two surfaces—orthogonal to each other to maintain a broader reflectance profile and provide polarization independence. Additionally, one can take the spectral response of a single metasurface layer, stack them, and then use a decoupled method approximation to predict the total spectral response from the multilayer stack. This approach could be used with a library of single-layer metasurfaces to efficiently generate the optical responses for all possible combinations of layers with a much lower computational cost than full RCWA simulations of all possible multilayer devices. All of these proposed polymer-embedded metasurface designs have the benefits of improved durability, elimination of condensation problems within or between meta-atoms, and the ability to stack numerous designs. Most notably, this work reveals that, despite reduction in index contrast, embedded metasurfaces can meet high performance metrics conventionally set by high index metasurfaces in air, and the embedded form factor enables scalable (e.g., roll-to-roll) manufacturing of mechanically robust and environmentally stable metasurfaces.
Funding
Funding for this work was provided by the Air Force Research Laboratory (AFRL) under MaRSS contract task order #13 FA8650-16-D-5404-0013.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. Y. Ding and R. Magnusson, “Resonant leaky-mode spectral-band engineering and device applications,” Opt. Express 12(23), 5661 (2004). [CrossRef]
2. R. Magnusson and M. Shokooh-Saremi, “Physical basis for wideband resonant reflectors,” Opt. Express 16(5), 3456 (2008). [CrossRef]
3. C. J. Chang-Hasnain and W. Yang, “High-contrast gratings for integrated optoelectronics,” Adv. Opt. Photonics 4(3), 379 (2012). [CrossRef]
4. Y. Ding and R. Magnusson, “Band gaps and leaky-wave effects in resonant photonic-crystal waveguides,” Opt. Express 15(2), 680 (2007). [CrossRef]
5. O. Quevedo-Teruel, H. Chen, A. Díaz-Rubio, G. Gok, A. Grbic, G. Minatti, E. Martini, S. Maci, G. V. Eleftheriades, M. Chen, N. I. Zheludev, N. Papasimakis, S. Choudhury, Z. A. Kudyshev, S. Saha, H. Reddy, A. Boltasseva, V. M. Shalaev, A. V. Kildishev, D. Sievenpiper, C. Caloz, A. Al, Q. He, L. Zhou, G. Valerio, E. Rajo-Iglesias, Z. Sipus, F. Mesa, R. Rodríguez-Berral, F. Medina, V. Asadchy, S. Tretyakov, and C. Craeye, “Roadmap on metasurfaces,” J. Opt. 21(7), 073002 (2019). [CrossRef]
6. X. Han, Z. Fan, Z. Liu, C. Li, and L. J. Guo, “Inverse design of metasurface optical filters using deep neural network with high degrees of freedom,” InfoMat 3(4), 432–442 (2021). [CrossRef]
7. M. Y. Shalaginov, S. D. Campbell, S. An, Y. Zhang, C. Ríos, E. B. Whiting, Y. Wu, L. Kang, B. Zheng, C. Fowler, H. Zhang, D. H. Werner, J. Hu, and T. Gu, “Design for quality: Reconfigurable flat optics based on active metasurfaces,” Nanophotonics 9(11), 3505–3534 (2020). [CrossRef]
8. S. An, B. Zheng, M. Y. Shalaginov, H. Tang, H. Li, L. Zhou, J. Ding, A. M. Agarwal, C. Rivero-Baleine, M. Kang, K. A. Richardson, T. Gu, J. Hu, C. Fowler, and H. Zhang, “Deep learning modeling approach for metasurfaces with high degrees of freedom,” Opt. Express 28(21), 31932 (2020). [CrossRef]
9. S. D. Campbell, D. Sell, R. P. Jenkins, E. B. Whiting, J. A. Fan, and D. H. Werner, “Review of numerical optimization techniques for meta-device design,” Opt. Mater. Express 9(4), 1842 (2019). [CrossRef]
10. S. So, T. Badloe, J. Noh, J. Rho, and J. Bravo-Abad, “Deep learning enabled inverse design in nanophotonics,” Nanophotonics 9(5), 1041–1057 (2020). [CrossRef]
11. M. M. Elsawy, S. Lanteri, R. Duvigneau, J. A. Fan, and P. Genevet, “Numerical optimization methods for metasurfaces,” Laser and Photonics Reviews 14(10), 1900445 (2020). [CrossRef]
12. R. S. Hegde, “Deep learning: a new tool for photonic nanostructure design,” Nanoscale Adv. 2(3), 1007–1023 (2020). [CrossRef]
13. W. Ma, Z. Liu, Z. A. Kudyshev, A. Boltasseva, W. Cai, and Y. Liu, “Deep learning for the design of photonic structures,” Nat. Photonics 15(2), 77–90 (2021). [CrossRef]
14. Z. Ma and Y. Li, “Parameter extraction and inverse design of semiconductor lasers based on the deep learning and particle swarm optimization method,” Opt. Express 28(15), 21971 (2020). [CrossRef]
15. I. U. Idehenre and M. S. Mills, “Multi-directional beam steering using diffractive neural networks,” Opt. Express 28(18), 25915 (2020). [CrossRef]
16. C. C. Nadell, B. Huang, J. M. Malof, and W. J. Padilla, “Deep learning for accelerated all-dielectric metasurface design,” Opt. Express 27(20), 27523 (2019). [CrossRef]
17. J. R. Thompson, J. A. Burrow, P. J. Shah, J. Slagle, E. S. Harper, A. Van Rynbach, I. Agha, and M. S. Mills, “Artificial neural network discovery of a switchable metasurface reflector,” Opt. Express 28(17), 24629 (2020). [CrossRef]
18. Z. Liu, D. Zhu, S. P. Rodrigues, K.-T. Lee, and W. Cai, “Generative model for the inverse design of metasurfaces,” Nano Lett. 18(10), 6570–6576 (2018). [CrossRef]
19. E. S. Harper, E. J. Coyle, J. P. Vernon, and M. S. Mills, “Inverse design of broadband highly reflective metasurfaces using neural networks,” Phys. Rev. B 101(19), 195104 (2020). [CrossRef]
20. E. S. Harper and M. S. Mills, “Bayesian optimization of neural networks for the inverse design of all-dielectric metasurfaces,” in Emerging Topics in Artificial Intelligence 2020, G. Volpe, J. B. Pereira, D. Brunner, and A. Ozcan, eds., (SPIE, 2020), p. 34.
21. Y. Kiarashinejad, M. Zandehshahvar, S. Abdollahramezani, O. Hemmatyar, R. Pourabolghasem, and A. Adibi, “Knowledge Discovery in Nanophotonics Using Geometric Deep Learning,” Advanced Intelligent Systems 2(2), 1900132 (2020). [CrossRef]
22. S. Jafar-Zanjani, S. Inampudi, and H. Mosallaei, “Adaptive genetic algorithm for optical metasurfaces design,” Sci. Rep. 8(1), 11040 (2018). [CrossRef]
23. M. M. Elsawy, S. Lanteri, R. Duvigneau, G. Brière, M. S. Mohamed, and P. Genevet, “Global optimization of metasurface designs using statistical learning methods,” Sci. Rep. 9(1), 17918 (2019). [CrossRef]
24. M. Mansouree, A. McClung, S. Samudrala, and A. Arbabi, “Large-scale parametrized metasurface design using adjoint optimization,” ACS Photonics 8(2), 455–463 (2021). [CrossRef]
25. M. Zhou, D. Liu, S. W. Belling, H. Cheng, M. A. Kats, S. Fan, M. L. Povinelli, and Z. Yu, “Inverse design of metasurfaces based on coupled-mode theory and adjoint optimization,” ACS Photonics 8(8), 2265–2273 (2021). [CrossRef]
26. E. B. Whiting, S. D. Campbell, L. Kang, and D. H. Werner, “Meta-atom library generation via an efficient multi-objective shape optimization method,” Opt. Express 28(16), 24229 (2020). [CrossRef]
27. S. Sengupta, S. Basak, and R. Peters, “Particle swarm optimization: A survey of historical and recent developments with hybridization perspectives,” Machine Learning and Knowledge Extraction 1(1), 157–191 (2018). [CrossRef]
28. A. Silva, A. Neves, and T. Gonçalves, “An heterogeneous particle swarm optimizer with predator and scout particles,” in Lecture Notes in Computer Science, vol. 7326M. Kamel, F. Karray, and H. Hagras, eds. (Springer Berlin Heidelberg, 2012), pp. 200–208.
29. A. Silva, A. Neves, and E. Costa, “An empirical comparison of particle swarm and predator prey optimisation,” in AICS, M. O’Neill, R. F. E. Sutcliffe, C. Ryan, M. Eaton, and N. Griffith, eds. (Springer-Verlag, Berlin, Heidelberg, 2002), pp. 103–110.
30. A. Silva, A. Neves, and T. Gonçalves, “Using scout particles to improve a predator-prey optimizer,” in Adaptive and Natural Computing Algorithms, vol. 7824M. Tomassini, A. Antonioni, F. Daolio, and P. Buesser, eds. (SpringerBerlin Heidelberg, 2013), pp. 130–139.
31. D. Z. Zhu, P. L. Werner, and D. H. Werner, “Design and optimization of 3-D frequency-selective surfaces based on a multiobjective lazy ant colony optimization algorithm,” IEEE Trans. Antennas Propag. 65(12), 7137–7149 (2017). [CrossRef]
32. D. Z. Zhu, E. B. Whiting, S. D. Campbell, P. L. Werner, and D. H. Werner, “Inverse design of three-dimensional nanoantennas for metasurface applications,” in 2019 International Applied Computational Electromagnetics Society Symposium, (IEEE, 2019), pp. 1–2.
33. B. Yuce, M. S. Packianather, E. Mastrocinque, D. T. Pham, and A. Lambiase, “Honey bees inspired optimization method: The bees algorithm,” Insects 4(4), 646–662 (2013). [CrossRef]
34. K. D. Donda and R. S. Hedge, “Rapid design of wide-area heterogeneous electromagnetic metasurfaces beyond the unit-cell approximation,” Prog. Electromagn. Res. M 60, 1–10 (2017). [CrossRef]
35. S. Kessentini and D. Barchiesi, “Particle swarm optimization with adaptive inertia weight,” International Journal of Machine Learning and Computing 5(5), 368–373 (2015). [CrossRef]
36. F. Paiva, J. Costa, and C. Silva, “A serendipity-based pso approach to delay premature convergence using scout particles,” International Journal of Innovative Computing, Information and Control 12, 1141–1163 (2016).
37. S. Jahani and Z. Jacob, “All-dielectric metamaterials,” Nat. Nanotech. 11, 23–36 (2016). [CrossRef]
38. V. Raghunathan, J. Deka, S. Menon, R. Biswas, and A. S. Lal Krishna, “Nonlinear optics in dielectric guided-mode resonant structures and resonant metasurfaces,” Micromachines 11(4), 449 (2020). [CrossRef]
39. S. Zhang, Y. H. Ko, and R. Magnusson, “Broadband guided-mode resonant reflectors with quasi-equilateral triangle grating profiles,” Opt. Express 25(23), 28451 (2017). [CrossRef]
40. Y. H. Ko and R. Magnusson, “Wideband dielectric metamaterial reflectors: Mie scattering or leaky Bloch mode resonance?” Optica 5(3), 289 (2018). [CrossRef]
41. V. Karagodsky, F. G. Sedgwick, and C. J. Chang-Hasnain, “Theoretical analysis of subwavelength high contrast grating reflectors,” Opt. Express 18(16), 16973–16988 (2010). [CrossRef]
42. M. Shokooh-Saremi and R. Magnusson, “Wideband leaky-mode resonance reflectors: Influence of grating profile and sublayers,” Opt. Express 16(22), 18249 (2008). [CrossRef]
43. X. Yin, J.-S. Park, K. K. Stensvad, R. L. Brott, M. B. Wolk, F. Capasso, and N. Rubin, “Roll-to-roll dielectric metasurfaces,” in Metamaterials, Metadevices, and Metasystems 2020, N. Engheta, M. A. Noginov, and N. I. Zheludev, eds. (SPIE, 2020), p. 34.
44. H. B. Van Lengerich, K. K. Stensvad, E. L. Kusilek, M. S. Stay, C. T. Nelson, C. S. Lyons, M. M. David, J. L. Solomon, M. B. Wolk, N. C. Erikson, J. Zhu, and J. M. Nelson, “Materials and methods for forming nano-structures on substrate,” WIPO patent 2020/095258 A1 (14 May 2020).
45. M. B. Wolk, J. M. Nelson, K. K. Stensvad, H. B. Van Lengerich, C. S. Lyons, M. M. David, J. L. Solomon, N. C. Erikson, and C. T. Nelson, “Nanostructured optical films and intermediates,” WIPO patent 2020/097319 A1 (14 May 2020).
46. J. Chang, S. Lee, K. B. Lee, S. Lee, Y. T. Cho, J. Seo, S. Lee, G. Jo, K.-y. Lee, H.-S. Kong, and S. Kwon, “Overlay accuracy on a flexible web with a roll printing process based on a roll-to-roll system,” Rev. Sci. Instrum. 86(5), 055108 (2015). [CrossRef]
47. S. C. Ambhire, S. Palkhivala, A. Agrawal, A. Gupta, G. Rana, R. Mehta, D. Ghindani, A. Bhattacharya, V. G. Achanta, and S. S. Prabhu, “"Pattern and Peel" method for fabricating mechanically tunable terahertz metasurface on an elastomeric substrate,” Opt. Mater. Express 8(11), 3382–3391 (2018). [CrossRef]
48. X. Y. Han, Z. L. Wu, S. C. Yang, F. F. Shen, Y. X. Liang, L. H. Wang, J. Y. Wang, J. Ren, L. Y. Jia, H. Zhang, S. H. Bo, G. Morthier, and M. S. Zhao, “Recent progress of imprinted polymer photonic waveguide devices and applications,” Polymers 10(6), 603 (2018). [CrossRef]
49. D. K. Oh, T. Lee, B. Ko, T. Badloe, J. G. Ok, and J. Rho, “Nanoimprint lithography for high-throughput fabrication of metasurfaces,” Frontiers of Optoelectronics 14(2), 229–251 (2021). [CrossRef]
50. C. Zhang, H. Subbaraman, Q. Li, Z. Pan, J. G. Ok, T. Ling, C. J. Chung, X. Zhang, X. Lin, R. T. Chen, and L. J. Guo, “Printed photonic elements: Nanoimprinting and beyond,” J. Mater. Chem. C 4(23), 5133–5153 (2016). [CrossRef]
51. H. H. Li, “Refractive index of silicon and germanium and its wavelength and temperature derivatives,” J. Phys. Chem. Ref. Data 9(3), 561–658 (1980). [CrossRef]
52. I. H. Malitson, “Interspecimen Comparison of the Refractive Index of Fused Silica,” J. Opt. Soc. Am. 55(10), 1205 (1965). [CrossRef]
53. X. Zhang, J. Qiu, X. Li, J. Zhao, and L. Liu, “Complex refractive indices measurements of polymers in visible and near-infrared bands,” Appl. Opt. 59(8), 2337 (2020). [CrossRef]
54. V. Liu and S. Fan, “S4: A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. 183(10), 2233–2244 (2012). [CrossRef]
55. L. J. V. Miranda, “PySwarms: A research-toolkit for particle swarm optimization in python,” Journal of Open Source Software 3(21), 433 (2018). [CrossRef]
56. N. Padhye, K. Deb, and P. Mittal, “Boundary handling approaches in particle swarm optimization,” in Advances in Intelligent Systems and Computing, vol. 201J. Bansal, P. Singh, K. Deep, M. Pant, and A. Nagar, eds. (Springer, India, 2013), pp. 287–298.
57. S. Helwig, J. Branke, and S. Mostaghim, “Experimental analysis of bound handling techniques in particle swarm optimization,” IEEE Trans. Evol. Computat. 17(2), 259–271 (2013). [CrossRef]
58. Y. H. Ko, M. Niraula, K. J. Lee, and R. Magnusson, “Properties of wideband resonant reflectors under fully conical light incidence,” Opt. Express 24(5), 4542 (2016). [CrossRef]
59. Y. H. Ko, K. J. Lee, and R. Magnusson, “Experimental demonstration of wideband multimodule serial reflectors,” Opt. Express 25(8), 8680 (2017). [CrossRef]
60. D. Zhao, H. Yang, Z. Ma, and W. Zhou, “Polarization independent broadband reflectors based on cross-stacked gratings,” Opt. Express 19(10), 9050–9055 (2011). [CrossRef]
61. M. C. Troparevsky, A. S. Sabau, A. R. Lupini, and Z. Zhang, “Transfer-matrix formalism for the calculation of optical response in multilayer systems: from coherent to incoherent interference,” Opt. Express 18(24), 24715 (2010). [CrossRef]
62. C. C. Katsidis and D. I. Siapkas, “General transfer-matrix method for optical multilayer systems with coherent, partially coherent, and incoherent interference,” Appl. Opt. 41(19), 3978 (2002). [CrossRef]
63. J. Zhang, S. Shi, H. Jiao, X. Ji, Z. Wang, and X. Cheng, “Ultra-broadband reflector using double-layer subwavelength gratings,” Photonics Res. 8(3), 426 (2020). [CrossRef]
