1. Introduction

In conventional optics, light manipulation is achieved by relying on propagation effects in refractive components. According to their trajectories, the light rays propagate along different optical paths δ, thus accumulating different optical phase shifts, ϕ = kδ, where k is the free space wavevector to mold the wavefront in a desired way. This refractory optics approach uses bulky materials, which, for a vast majority, have fixed optical properties. In order to fully address the direction of light, the phase shifts accumulated along the different ray trajectories, traversing across one or several materials of various shapes, should cover mod(2π), meaning that the longitudinal dimension of the optical elements have to be several orders thicker than the wavelength, l > λ. Leveraging on nanostructured materials and arranging tiny features in subwavelength arrays, it is possible to break away from our dependence on refractive materials and propose innovative optical components capable of controlling the wavefront of light in an almost arbitrary manner. These artificial components, which do not exist naturally, are called metasurfaces [1–3] and are closely related to the high-contrast subwavelength dielectric grating (HCG) [4,5], studied in details in the nineties. The interested reader will find more details for example in these papers, reviews and reading references therein [6–11]. HCGs have unconventional properties obtained by varying the material composition, i.e. by variation of the duty cycle of the grating. Beside the similarities with HCG, considerable progress in nanofabrication has recently boosted this field of research to the point that various sorts of new metasurfaces, achieving multiple optical functionalities going beyond simple phase addressing, are revolutionizing optical designs. Today, metasurfaces find numerous important applications in imaging [12,13], in nonlinear and quantum optics, in digital coding [14–17], polarimetry [18–22] and life-science, i.e. medical applications [23]. For these applications to bear fruits at the industrial levels, the device performance and their reliability should be improved and maximized. Being aware of the requirement for real-world applications of metasurfaces, we are considering suitable material systems for their realization which could be cost-effective, enabling large-area fabrication techniques and which also maintain the required material properties after nanoprocessing. The material system for realizing the devices should also meet the stringent requirements of industrial manufacturing processes such as C-MOS compatibility. Gallium Nitride (GaN), is a promising industrially relevant semiconductor. Here we have selected this material for the realization of our optimized metasurface subwavelength gratings because it offers multiple advantages: GaN is highly transparent in the entire visible range and components for various wavelengths can be designed; its refractive index (n > 2) is relatively high and ensures that GaN nanostructures can possess strong scattering resonances, while the dispersion is sufficiently small as to be considered fixed within the wavelength of interest; its considerably high thermal and chemical stability make this material suitable even for extreme applications; from the industrial point of view, it is a mature material, which is also highly compatible with the semiconductor fabrication techniques, as shown by its numerous implications in electronics and optoelectronics. To date, the most common way of designing metaoptics components usually considers the optical response of a unique scatterer (or a subwavelength array of identical elements). The proximity effects between adjacent nanostructures and the impact of various sorts of nanofabrication uncertainties are generally neglected. Near-field coupling plays a detrimental role in the scattering response of closely packed nanostructures. This effect is sufficiently strong such that it can be used as a design parameter to realize multiwavelength and broadband metasurfaces [24–29]. However, considering the proximity effects of closely packed antennas is not a simple task and it often requires tedious and time consuming numerical iteration procedures to optimize the overall metasurface response. Several promising works have been recently realized in order to tackle this optimization problem using different optimization methods including for example particle swarm optimization algorithms [30], gradient-based inverse design method [31], adaptive generic [32] and genetic algorithms [33], evolutionary [34] and deep learning models [35]. Although all of these optimization methods could lead to highly efficient designs, their experimental realizations hardly reach the expected performances. The fabrication errors and imperfection with respect to the designed structures make these optimization tools excessively complicated, going way beyond the experimental feasibility. In parallel of the development of metasurface technology, significant progresses have recently been achieved in the area of uncertainty quantification (UQ), where the aim is to estimate the influence of fabrication and measurement errors on physical quantities of interest. Most notably, higher-order stochastic spectral methods [36–40] have been developed, analyzed and applied with great success to many different application areas. Contributions in electromagnetics are reported in [41–43] among others. Relying on surrogate or meta-modeling, with the aid of modern computational architectures, addressing complex computational models with a large number of uncertain inputs comes into reach. Yet, a key prerequisite is a smooth input-to-output behavior of the underlying model, which is not always present in applications. In this case, Monte Carlo techniques and their contemporary multilevel and multifidelity variants [44,45] are still the methods of choice. In the area of nanoplasmonics, recent work on uncertainty quantification has been reported in [46]. In this work we have studied and quantified the robustness of optimized metasurfaces, considering the impact of structural errors within given experimental uncertainties. The performance of fabricated optimized structures show relatively good transmission and deflection performances in agreement with our numerical simulations. As we will discuss in Section 3, polynomial surrogate models are not readily applicable for the class of metasurfaces considered in this work and hence, a Monte Carlo strategy is adopted. Starting from several optimized designs obtained by computing the electromagnetic response of arrays of subwavelength ridges using rigorous coupled wave analysis software (RCWA) [47], we performed an UQ analysis and explored the impact of technologically relevant geometrical parameters such as ridge widths δx_i and relative positions x_i. A Monte Carlo ensemble of a million of random realizations has been numerically simulated to evaluate the sensitivity of the optimal design with respect to manufacturing imperfections through statistical indicators like mean deflection efficiency, standard deviation to the mean value and confidence intervals. The results show that the number of elements per phase gradient period plays a considerable role in the reliability of the structures with respect to small ±5 nm uncertainties in the widths and positions. We conclude our discussion with a study on non-periodic perturbations, i.e. we include more and more gratings with independent uncertain parameters into a periodic unit cell. These last results are of capital importance for end-users and industrial users interested in understanding the reproducibility and the overall mean performances of a large number of metaoptics replica.

2. Simulation, optimization and experimental realization

The overall aim of this work is to design an ensemble of optimized, high contrast, GaN subwavelength gratings with close to unity blazed response, while seeking robust structures, insensitive to small manufacturing imperfections. Dedicated tools for robust optimization are available in the literature [48], where uncertainty typically enters the objective function and constraints via a penalization. Here, we decouple optimization and uncertainty analysis, for simplicity. In particular, we first optimize several structures at nominal parameter values before performing Monte Carlo analysis on the optimized design to carry out UQ. The main reason for this simplification is that our sampling approach is numerically too expensive to be considered within the optimization routine. Constructing more efficient UQ methods and performing robust optimization for this application is subject of ongoing research efforts. Nevertheless, the uncertainty analysis at the optimum gives important insights into the efficiency of the manufactured structure. The numerical methods used for simulation and optimization are described below, followed by a report on experimental results.

2.1. Methodology

GaN metasurfaces and a sketch of the optimization principle are depicted in Fig. 1(a). In the following, all metasurfaces are optimized for an incident plane wave with a wavelength λ, impinging onto the structure at normal incidence and deflecting the light in the transmission mode at desired angles θ, given by the grating law mλ = Γ sin(θ), where Γ > λ represents the width of the unit cell. According to this grating law, only one fixed wavelength with a fixed period can provide the first order at a desired deflection angle. Due to the periodicity of the structure (Fig. 1(b)), the deflection angle follows the grating law, but because of the subwavelength nanostructuration, we expect a strong blazing effect as depicted in Fig. 1(c). The nanoridges in a unit cell have different widths to introduce a gradual phase retardation ranging from 0 to 2π in a period, thus concentrating the light flux into one selected diffraction order only. The purpose of our optimization procedure and the uncertainty analysis is to reveal the best arrangement of subwavelength features to adjust the transmission and phase gradient within a given period, similarly as it is done in echelette gratings. The latter refractive devices which work by accumulating propagation phase shift across the tiny prisms, from 0 to 2π going from the apex to the base of the prism as shown in see Fig. 1(b), are difficult to manufacture and have high deflection efficiency only around a narrow angular range. An accurate electromagnetic model of the GaN metasurface is thus crucial for obtaining satisfactory optimization results. For this purpose and in the scope of this work, we have used two complementary simulation tools, i.e. two numerical solvers for Maxwell’s equations. First, RCWA is chosen for its extremely efficient monochromatic simulations, which is very well suited for this type of 1D optimization calculations. For broadband calculations, the Discontinuous Galerkin Time-Domain (DGTD) method [49] gives accurate broadband results and serves as a reference. Here, we have used the RCWA implementation reticolo within the optimizer. The use of RCWA can be justified by a numerical cross-comparison with the DGTD solver from the DIOGENeS software suite [50], as depicted in Fig. 1(d). The comparison is carried out over a wide spectrum and shows excellent agreement. Generally, if not indicated otherwise, all simulations have been performed in the transverse magnetic (TM) polarization, i.e. for an electric field oscillating perpendicular to the ridges axis.

 

Fig. 1 (a) Illustrative schematic of the angular deflection property of a phase gradient metasurface. Metasurfaces based on arrays of unit cells distributed in a periodic manner with Γ > λ, allowing discrete diffracted modes only. In the most general case, and assuming an incident light at normal incidence, in transmission, the light will be deflected into the different orders (permitted by the periodicity of the microscopic structure). A single diffracted beam is achievable by controlling the light diffracted by each unit cell. In (b), the device works essentially as conventional echelette blazed grating. Replacing the periodic echelette with a subwavelength array of nanoridges (c), we create a metasurface. A strong influence of the nanoridge response as a function of their width is obtained for TM polarized incident light. The width of a periodic unit cell is denoted by Γ, the centers of the ridges by x = [x₁, . . ., x_i, . . ., x_N], where i = 1, . . ., N is the index of the nanoridge and N the total number of nanoridges. In the same manner, we define the widths of all ridges by δx = [δx₁, .., δx_i, .., δx_N]. (d) Typical broadband response of the transmission efficiency for an optimized metasurface obtained using two different electromagnetic simulation solvers, the Discontinuous Galerkin Time-Domain solver (DGTD) solver and the Rigorous Coupled Wave Analysis (RCWA) solver in orange dashed.

Download Full Size | PDF

2.2. Optimization results

For the optimization we apply a gradient-free pattern search algorithm, implemented in Matlab’s optimization toolbox. Pattern search is well suited for this kind of setups due to its simplicity and robustness. The optimizer seeks to maximize the transmission efficiency of the first diffracted order. In a first step, we optimize the one-dimensional grating structure at a fixed design angle and investigate the broadband and geometric robustness of the transmission efficiency η₋₁ of the first order diffracted mode. The height of each grating element, the grating period and the minimum feature size are given by h = 1000 nm, Γ = 1322 nm and 𝔡 = 90 nm, respectively. The grating period Γ implies a deflection angle of θ = 27°, i.e. we seek for an optimal geometry resulting in maximum light intensities at θ = 27°. All grating elements are made of GaN, epitaxially grown on a double side polished Al₂O₃ (Sapphire) substrate using molecular beam epitaxy.

We define the geometry, i.e. the positions and widths of the grating elements, through abstract geometry parameters α ∈ [0, 1]^N and a ∈ [0, 1]^N, which ensure that we obtain a grating with N ∈ ℕ non-overlapping nanoridges inside of the grating period Γ of size bigger or equal to the feature size. The concrete connection between α_i, a_i and the grating element geometry x_i, δx_i is given in the appendix. We define the objective function as f_N,m(α, a) = 1 − η_N,m(α, a), where η_N,m : ℝ^2N → ℝ refers to the transmission efficiency of the mode m ∈ ℤ at a fixed wavelength λ_design ∈ ℝ. In our case, all materials are passive and hence η_N,m ∈ [0, 1].

We formulate the nonlinear constrained optimization problem as

Figure 2 depicts the broadband transmission behavior of the first diffraction order η₋₁ of all five designs N ∈ {2, 3, 4, 5, 6}, after optimization. Details of the geometric data is given in Table 2 in the appendix. For the TM mode, the transmission efficiencies are increasing with the number of nanoridges N. Also, for N = 6, the efficiency is rather independent in a small vicinity of λ = 600 nm. Figure 3 outlines the impact of a varying height on the grating. Again, for N = 6, η₋₁ is very robust with respect to changes of h around the nominal value h = 1000 nm. This robustness is less pronounced for decreasing N.

 

Fig. 2 Diffraction efficiency spectra for θ = 27°, h = 1000 nm, d = 1322 nm, 𝔡 = 90 nm, and varying N. The design wavelength is λ_design = 600 nm. The upper and lower plots correspond to a TM and TE polarization, respectively.

Download Full Size | PDF

 

Fig. 3 Diffraction efficiency spectra for θ = 27°, 𝔡 = 90 nm, different N and varying grating height h. The design wavelength and height are λ_design = 600 nm and h = 1000 nm, respectively.

Download Full Size | PDF

Table 1. UQ results for η₋₁ in the single frequency case and for uniform input distributions together with measured η₋₁ of the fabricated structure.

View Table | View all tables in this article

Table 2. Grating geometry. Numerical values of the grating’s geometric parameters for a deflection angle θ = 27°, grating period Γ = 1322 nm and different numbers of nanoridges per period N. The optimizations were carried out using the larger minimal feature size 𝔡 = 90 nm.

View Table | View all tables in this article

The same optimization strategy as outlined above has now been employed for different design angles θ, and hence a changing grating period Γ(θ). Figure 4 depicts the transmission efficiency η₋₁ as a function of the design angle θ on the bottom and their broadband behavior on the top. Again, the concrete geometric data is provided in the Tables 3 and 4.

Table 3. Grating geometry. Numerical values of the grating’s geometric parameters for different deflection angles θ and hence grating periods Γ. The optimizations were carried out using the smaller minimal feature size 𝔡 = 50 nm. Here, we only show the setting with the optimal number of nanoridges per period.

View Table | View all tables in this article

Table 4. Grating geometry. Numerical values of the grating’s geometric parameters for different deflection angles θ and hence grating periods Γ. The optimizations were carried out using the larger minimal feature size 𝔡 = 90 nm. Here, we only show the setting with the optimal number of nanoridges per period.

View Table | View all tables in this article

 

Fig. 4 Top: Diffraction efficiency spectra for h = 1000 nm, 𝔡 = 50nm (left) and 𝔡 = 90 nm (right), and varying design angles θ and Γ(θ), respectively. The design wavelength is λ_design = 600 nm. Bottom: Angle dependence of η₋₁. The number of nanoridges N varies depending on the angle. The design wavelength is again λ_design = 600 nm.

Download Full Size | PDF

2.3. Experimental realization and characterization of optimized metagratings

In this section, we discuss the experimental investigations of optimized metagratings fabricated using standard nanofabrication processing and their optical characterizations. The latter are performed using a homemade k-space microscopy setup to measure the angular deflection efficiency as a function of the incident wavelength. A schematic of the experiment is presented in Fig. 5(a).

 

Fig. 5 Experimental setup and k-space microscopy images of metagrating deflection properties. (a) the broadband light transmitting through the metasurface is collected with a high microscope objective (N A = 0.9) of an inverted microscope. A modified 4f lensing system relays the k-space information to the entrance slit of an imaging spectrometer. (b) shows a typical k-space image of the light transmitted through the metagrating when the spectrometer is set in imaging mode, i.e. by looking at the back focal plane image reflected toward the CCD camera by the diffraction grating oriented to the specular reflection angle. The deflected fields at m = −1 diffraction order appears as an elongated spot around 45°. The red circle indicates the NA of the microscope objective. (c) The spectral response of the components is characterized by closing the spectrometer entrance slit and spectrally dispersing the vertical k-space information. The image displays both spectral information (along the x dimension) and transverse k_y momentum information. The m = −1 has higher intensity indicating strong blazing behaviour as expected from the design. Note also the typical k-space angular dispersion properties of the metagrating as function of the wavelength.

Download Full Size | PDF

The excitation of the metasurface from the substrate side is realized using a supercontinuum broadband laser (NKT photonics) coupled to a spectral filter to select 100 nm bandwidth around 600 nm wavelengths. The beam is then filtered in polarization and focused on the metasurface with a spot size of about 200 μm. After interacting with the metasurface, a high numerical aperture (NA = 0.9) microscope objective coupled to a tube lens collects the transmitted signal and images the metasurface plane at the exit port of an inverted microscope Nikon eclipse TE. A pinhole is inserted at this image plane to select part of the transmitted signal. This way, we are characterizing the transmitted fields passing exclusively across the metasurfaces, thus avoiding unwanted additional straight light at normal incidence which would change the deflection efficiency. A modified 4f setup with f₁ = 2 f₂ is added to magnify the objective back focal plane on the entrance slit of an Andor Shamrock 500i spectrograph (500 mm focal distance), which is in turn coupled to an intensified Istar CDD camera. When the spectrometer is operating in imaging mode at the zero order, i.e. in the case of specular reflection on the diffraction mirror, the 2D intensity map gives full k-space information as presented in Fig. 5(b). By closing the slit and spectrally dispersing the vertical k-space information, we obtain images that contain both spectral information (along the x dimension) and transverse k_y momentum information as depicted in Fig. 5(c). It is important to underline that the metasurfaces are oriented such that the deflection, related to direction along which the phase gradient acts, is chosen in agreement with k_y imaging capability of our microscopy system. Using a non-optimized calibration sample, designed to deflect normally incident 600 nm light at an angle of 45 deg, we observe an intense spot at the m = −1 order and weaker signals at m = 0, 1 indicating that the metasurface is efficiently refracting light at the designed refracted angle of 45° for 600 nm wavelength. Moreover, we observe an increasing deflection angle as a function of the incident wavelength in agreement with the abnormal dispersion given by the generalized Snell law $\text{sin}({\theta }_t)=\frac{\lambda }{2\pi }\frac{d\varphi }{dy}$. The k-space calibration of our experimental Fourier plane images is performed by measuring the k_ymin and k_ymax of the objective lens, bringing the spectrometer diffraction grating to the zero order and calculating the transverse k momentum as a function of the number of pixels filling the NA of the objective. The latter is delimited by the red circle in Fig. 5(b).

For the fabrication of the structures, we start the process by growing a GaN layer on a double side polished (111) Sapphire substrate using Molecular Beam Epitaxy (MBE) reactor. The state of the art MBE technique can provide accurate thickness control as well as large scale uniformity. A Hydrogen silsesquioxane (HSQ) resist is then spin-coated onto the GaN layer for electron beam lithography process. The clearly perforated resist pattern is used as an etching mask for the RIE process. After etching the GaN layer, the rest of the resist is removed using chemical native oxide removal by dipping the patterned films in a BOE etch. The results of the fabrication are shown in Fig. 6(a). A scanning electron micrograph of the resulting structures is shown in Fig. 6(b) showing clearly defined GaN ridges on sapphire substrate.

 

Fig. 6 Example of a fabricated metasurface. (a, top) The schematics indicate the fabrication work flow, starting from GaN growth on Sapphire substrate, followed by eBeam lithography on HSQ resist which is cured and used as a hard mask for RIE etching. The excess of resist is then washed away using BOE etch. (b) A typical scanning electron micrograph of a fabricated structure. (c) Zoomed SEM images of the different nanoridge gratings with increasing number of elements per period from 2 to 6. In the large scale figure, the assemblies of nanoridges highlighted in yellow represent unit cells of the metasurfaces. In the inset, a zoomed of unit cells revealing the sub-unit nanoridges chosen to discretize the 2π phase ramp. In this example, the unit cells are composed of 2 to 6 nanoridges indicated by the colors ranging from red to violet. (d) Fourier plane optical transmission curves. The red curves represent the metasurface transmission as a function of the deflection angles for various nanoridges per periods. The data have been measured at 600 nm. (e) Summary of the deflection efficiencies as function of the number of elements per period.

Download Full Size | PDF

The experimentally measured transmission spectra as a function of the deflection angles for an interface deflecting at 27° are presented in Fig. 6(d) for different numbers of nanoridges per period. A summary of the deflection efficiency is reported in (e) in which we compare the η₋₁ and η₀ by normalizing the signal directly from the bare substrate transmission (blue curves). We observe that the structure is highly efficient (> 80%) at the designed wavelength of 600 nm for 5 nanoridges per period. The performances decrease by decreasing the number of phased elements, as expected from numerical simulations. We also indicated the ratio between −1 and 0 diffraction orders to clearly show the blazing efficiency of the components. The nanoridges have important structural birefringence and are therefore sensitive to the state of polarization. We have calculated and experimentally confirmed that good blazing effect occurs only for TM incident polarization as presented in Fig. 7.

 

Fig. 7 Measured diffraction efficiency spectra for TM and TE incident polarizations. The structure is composed of 5 nanoridges per unit cell as presented in Fig. 6 for θ = 27°, h = 1000 nm, 𝔡 = 90 nm, and varying polarization. The design wavelength is λ_design = 600 nm. (a) and (b) are raw experimental k space images obtained with a microscope objective (0.9 NA). (c) are normalized efficiencies for S and P polarizations.

Download Full Size | PDF

3. Uncertainty quantification

We conduct a number of uncertainty quantification (UQ) studies with the aim of quantifying the impact of random geometrical deviations in the considered structures upon the transmission efficiency η₋₁. As explained in the introduction, this analysis is of critical importance to understand the limitation and to evaluate the mean performance of a large number of components. We believe that such a study will comfort practical engineers who are interested in having meaningful numbers to evaluate the general metasurface performances. We consider structures with N ∈ {2, 3, 4, 5, 6} nanoridges per period and model the geometrical uncertainties as random shifts in the x-coordinates of the bounds of the ridges, thus obtaining a total of 2N uncertain parameters, equivalently, random variables (RVs). This stochastic modeling results in randomly varied grating element widths and positions.

A realization of a random x-coordinate is given by $X_i=X_i^{\text{nom}}+{\xi }_i(\theta )$, i = 1, . . ., 2N, where $X_i^{\text{nom}}$ denotes the x-coordinate in the nominal geometry of the structure, ξ_i the corresponding RV, and θ a random outcome. We use capital letters to distinguish the nanoridge boundaries X_i from its centers x_i. We model the RVs ξ_i to be independent and identically distributed (iid), following either uniform distributions with support within the range [−5 nm, 5 nm], or beta [−7.5 nm, 7.5 nm] and with shape parameters α = β = 4. With the selected shape parameters, the beta distribution is a good approximation to the normal distribution with a mean value μ = 0 nm and a standard deviation σ = 2.5 nm, while having a bounded support. An illustration is provided in Fig. 8(a). The effect of the distribution on the results will be analyzed below, yet, as we did not observe large differences for several computed quantities, the results for the beta distribution are often omitted.

 

Fig. 8 Input and output probability density functions.

Download Full Size | PDF

Denoting the probability density function (PDF) of a given distribution with ϱ_i(ξ_i), under the independent and identically distributed random variables assumption, the joint probability function is given by $\rho (\xi )={\prod }_{i=1}^{2N}{\rho }_i({\xi }_i)$, where ξ = (ξ₁, . . ., ξ_2N) is a random vector.

Our numerical experiments indicate that advanced UQ techniques based on (adaptive) spectral distributions with support in the range approaches cannot be used for the given problem, due to the fact that the input-to-output map, i.e. the behaviour of the underlying model, is not sufficiently smooth. However, since a single call of the RCWA solver requires less than 1 s of execution time, the classic Monte Carlo (MC) sampling method remains an attractive option. Exploiting parallel computing resources, even several millions of solver calls is a feasible task.

Using the MC sampling method, we generate a sample of M_MC realizations of the input random vector, ${\left\{{\xi }^{(m)}\right\}}_{m=1}^{M_{\text{MC}}}$. The realizations are drawn randomly from the joint input PDF ϱ(ξ). The solver is called for each one of the realizations and we obtain the corresponding solutions ${\eta }_{-1}^{(m)}={\eta }_{-1}\left({\xi }^{(m)}\right)$, m = 1, . . ., M_MC. Then, the expected (mean) value of the transmission efficiency, 𝔼[η₋₁], can be estimated as

A further attractive feature of the MC method is that it is unbiased with a root-mean-squared-error (RMSE) that can be estimated as

Finally, the sample of evaluations ${\left\{{\eta }_{-1}^{(m)}\right\}}_{m=1}^{M_{\text{MC}}}$ can be used for the estimation of other statistical measures, such as confidence intervals, quantiles, or the PDF of the now random η₋₁. For the latter, we use here a kernel density estimation approach, based on an Epanechnikov kernel [51].

3.1. Single frequency and periodic deformations

The simulations performed to obtain the results presented in this section have been executed for a fixed wavelength of 600 nm and a fixed deflection angle of θ = 27°. We consider the cases of N ∈ {2, 3, 4, 5, 6} nanoridges. For each configuration we conduct a MC experiment with M_MC = 10⁶ random realizations of the geometry, yielding a negligible sampling error. We calculate the MC estimates for the mean value and the standard deviation of η₋₁, as described above. We further compute the maximum and minimum efficiency values that occur in the generated output sample. The results for uniform input distributions are presented in Table 1. Similar results have been obtained for beta input distributions, but are here omitted. In this table, we also report the measured efficiencies. Simulation and experiment are in good agreement since, except for N = 3, the measured values do not deviate significantly more than one standard deviation of the numerical prediction. Note that the uncertainty analysis is based on estimations of the manufacturing uncertainty solely and further improvements seem possible by taking into account geometry measurements. The PDFs of η₋₁ for uniformly and beta-distributed RVs are presented in Figs. 8(b) and 8(c).

3.2. Broadband and periodic deformations

We consider wavelengths in the range 500 nm to 700 nm. The frequency range is discretized in an equidistant way with 41 discrete frequency points, i.e. with a step size of 5 nm. We generate 10⁵ random input realizations for each design, i.e. for N = 2, 3, 4, 5, 6 nanoridges per period. For each input realization, the computational model is evaluated at all frequency points, thus resulting in 41 transmission efficiency values. We also consider both uniform and beta input distributions regarding the 2N geometrical uncertainties.

The combinations between input distributions and nanoridges per period, result in 10 distinct settings. For each one of those settings, we have calculated expected (mean) values, standard deviations, 10%, 20%, . . ., 90% percentiles, and best/worst case efficiency values, over the given frequency range. As an example, we present in Fig. 9 the results for the case of uniform inputs and N = 5 nanoridges per period. In the same plot, the 10% percentile per frequency point is presented, indicating the percentage of realizations resulting in transmission efficiency values below the corresponding (black) curve.

 

Fig. 9 Broadband UQ results for uniform input distributions and N = 5 nanoridges per period.

Download Full Size | PDF

3.3. Non-periodic deformations

In the previous subsection, the simulation was confined to a single grating with periodic boundary conditions, implying periodic geometric deformations as well. Insight on the structure of the geometric uncertainties can be gained from measurement data by applying Bayesian parameter estimation techniques, for instance. Such a study is beyond the scope of the present work. Instead we analyze a second case for comparison, where the influence of the periodic boundary condition is systematically minimized. To this end, we proceed as illustrated in Fig. 10: instead of one, we consider two gratings per unit cell with identical geometry but independent geometric uncertainties and repeat the uncertainty analysis. This procedure is repeated with more and more gratings per unit cells, yielding an ever-increasing number of random variables, until changes in the mean value become negligible. Such a procedure is tractable in a Monte Carlo setting, as the computational accuracy only mildly depends on the number of parameters. The accuracy is controlled for each simulation by choosing the sample size according to (4). A detailed description of this algorithm is given in the appendix. For simplicity, in this study, we focus on the case N = 5 and uniform random inputs.

 

Fig. 10 Sketches of considered numerical models with N = 5 nanoridges per period. Black solid lines illustrate periodic boundary conditions. Black dots indicate uncertain parameters. Top: unit cell consisting of n = 1 grating as considered in Section 3.1 and 3.2. Middle, bottom: increased number of gratings per unit cell (top: n = 2; bottom: n = 4) in order to systematically reduce the influence of the periodic boundary conditions (Section 3.3).

Download Full Size | PDF

Figure 11 depicts the results for the mean value and variance over an increasing number of gratings per unit cells (random variables). The mean efficiency decreases slightly, before saturating after n = 16. The standard deviation in turn decreases continuously with a slope of 0.5. Hence, in absence of systematic fabrication offsets the problem features a more and more deterministic character, which can be physically attributed to interference effects between the individual variations. More precisely, the individual grating contributions add up to an overall efficiency and the central limit theorem implies that the distribution of this overall efficiency is normally distributed with a standard deviation proportional to $1/\sqrt{n}$, which is clearly visible in Fig. 11. In summary, the impact of manufacturing imperfections is reduced for robust designs due to averaging interference effects between neighboring structures. Small deviations from the ideal gradient phase within neighboring unit cells create a diffuse background signal, which decreases slightly the overall diffraction efficiency. As shown in this last figure, small variations do not drastically affect the optical performances of robust metasurfaces designs, indicating that optimized designs would have reproducible performances. These results are of significant practical importance, showing that multiple interference effects help to stabilize the overall performance of a large number of components within a few percent below the design mean value.

 

Fig. 11 Effect of periodic boundary conditions on UQ results. Mean value (left) and standard deviation (right). The periodic part consists of an increasing number of gratings n. The mean value decreases for an increasing number of gratings before saturating after roughly n = 15. The standard deviation and hence, the uncertainty due to the manufacturing are decreasing continuously.

Download Full Size | PDF

4. Conclusion and outlook

The development of metasurfaces and their deployment in real world applications require new numerical tools to design reliable and high efficient devices. In view of processing a large number of metasurfaces at once, these performances have to be sufficiently robust with respect to minute fabrication errors generally occurring in high-throughput wafer-scale manufacturing processes. In this work, we have analyzed both theoretically and experimentally the deflection efficiency of 1D GaN phase gradient metasurfaces. We could observe that the optimized structures are rather robust against manufacturing errors, in particular for a larger number of nanoridges, even though robustness was not yet taken into account during optimization. By quantifying uncertainties for the optimized design we obtained confidence intervals for the transmission efficiencies, which are in good agreement with the experimental results. Evaluating the potential and the physical limits of GaN metaoptics is also of importance to push further toward large market applications such as imaging, lighting and displays and hybrid optical components [52,53]. The results discussed in this manuscript comfort the initial assumption about the potential impact of GaN metamaterial and its implication in designing novel and efficient passive and active meta-optical components for practical applications.

Appendix A: Geometrical data

We present the numerical values of the optimized geometric parameters, i.e. centers x and widths δx of the nanoridges, for different gratings in Table 2, 3 and 4. In all cases the design wavelength is λ_design = 600 nm and the grating height is h = 1000 nm. Note that the respective grating period Γ corresponding to a deflection angle θ can be obtained using λ_design = Γ sin(θ).

Appendix 2: Algorithmic UQ aspects

Results presented in Section 3.3 are obtained using Fig. 12 with M_min = 100 and ∊ = 10⁻⁴. Table 5 shows that, due to the ever decreasing standard deviation, the needed number of MC samples decreases with respect to n, too. This partially mitigates the rapidly increasing computational cost for an increasing number of gratings per unit cell n.

Table 5. UQ results for an increasing number of nanoridges n per unit cell. N_par = 2nN denotes the number of uncertain parameters and M^MC the number of Monte Carlo samples used to estimate the mean value E_n and the standard deviation σ_n.

View Table | View all tables in this article

 

Fig. 12 Uncertainty quantification in quasi-periodic structures.

Download Full Size | PDF

Appendix 3: Details on geometric parametrization

The positions of the nanoridges can be derived from the abstract geometry parameters α_i, a_i ∈ [0, 1], i = 1, . . ., N as follows: the width of the i-th nanoridge is given as

$$\delta x_i=𝔡+{\alpha }_i(\mathrm{\Gamma }02N𝔡).$$

The corresponding center position $x_i\in \left(-\frac{\mathrm{\Gamma }}{2},\frac{\mathrm{\Gamma }}{2}\right)$ can then be obtained as

$$x_i=\frac{-\mathrm{\Gamma }+\delta x_i}{2}+\{\begin{array}{ll}0\hfill & i=1,\hfill \\ {\sum }_{j=1}^{i-1}\delta x_j+𝔡+a_j\left(\mathrm{\Gamma }-N𝔡-{\sum }_{k=1}^N\delta x_k\right)\hfill & i>1.\hfill \end{array}$$

Funding

European Research Council (ERC) (639109); German Research Foundation (DFG) (RO4937/1-1)

References

1. P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, “Recent advances in planar optics: from plasmonic to dielectric metasurfaces,” Optica 4, 139–152 (2017). [CrossRef]  

2. A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science 339, 1232009 (2013). [CrossRef]   [PubMed]  

3. H.-T. Chen, A. J. Taylor, and N. Yu, “A review of metasurfaces: physics and applications,” Reports on Prog. Phys. 79, 076401 (2016). [CrossRef]  

4. C. J. Chang-Hasnain and W. Yang, “High-contrast gratings for integrated optoelectronics,” Adv. Opt. Photonics 4, 379–440 (2012). [CrossRef]  

5. A. Arbabi, Y. Horie, A. J. Ball, M. Bagheri, and A. Faraon, “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nat. Commun. 6, 7069 (2015). [CrossRef]   [PubMed]  

6. M. Moharam and T. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” JOSA 71, 811–818 (1981). [CrossRef]  

7. P. Lalanne, J. P. Hugonin, and P. Chavel, “Optical properties of deep lamellar gratings: a coupled bloch-mode insight,” J. Light. Technol. 24, 2442–2449 (2006). [CrossRef]  

8. C. F. Mateus, M. C. Huang, L. Chen, C. J. Chang-Hasnain, and Y. Suzuki, “Broad-band mirror (1.12–1.62/spl mu/m) using a subwavelength grating,” IEEE Photonics Technol. Lett. 16, 1676–1678 (2004). [CrossRef]  

9. S. Collin, “Nanostructure arrays in free-space: optical properties and applications,” Reports on Prog. Phys. 77, 126402 (2014). [CrossRef]  

10. D. Fattal, J. Li, Z. Peng, M. Fiorentino, and R. G. Beausoleil, “Flat dielectric grating reflectors with focusing abilities,” Nat. Photonics 4, 466 (2010). [CrossRef]  

11. U. Levy, H.-C. Kim, C.-H. Tsai, and Y. Fainman, “Near-infrared demonstration of computer-generated holograms implemented by using subwavelength gratings with space-variant orientation,” Opt. Lett. 30, 2089–2091 (2005). [CrossRef]   [PubMed]  

12. S. Colburn, A. Zhan, and A. Majumdar, “Metasurface optics for full-color computational imaging,” Sci. Adv. 4, eaar2114 (2018). [CrossRef]   [PubMed]  

13. M. Khorasaninejad, A. Y. Zhu, C. Roques-Carmes, W. T. Chen, J. Oh, I. Mishra, R. C. Devlin, and F. Capasso, “Polarization-insensitive metalenses at visible wavelengths,” Nano Lett. 16, 7229–7234 (2016). [CrossRef]   [PubMed]  

14. S. M. Kamali, E. Arbabi, A. Arbabi, Y. Horie, M. Faraji-Dana, and A. Faraon, “Angle-multiplexed metasurfaces: Encoding independent wavefronts in a single metasurface under different illumination angles,” Phys. Rev. X 7, 041056 (2017).

15. H. Zhao, B. Quan, X. Wang, C. Gu, J. Li, and Y. Zhang, “Demonstration of orbital angular momentum multiplexing and demultiplexing based on a metasurface in the terahertz band,” ACS Photonics 5, 1726–1732 (2017). [CrossRef]  

16. D. Wen, F. Yue, G. Li, G. Zheng, K. Chan, S. Chen, M. Chen, K. F. Li, P. W. H. Wong, K. W. Cheah, E. Y. B. Pun, S. Zhang, and X. Chen, “Helicity multiplexed broadband metasurface holograms,” Nat. Commun. 6, 8241 (2015). [CrossRef]   [PubMed]  

17. M. Mehmood, S. Mei, S. Hussain, K. Huang, S. Siew, L. Zhang, T. Zhang, X. Ling, H. Liu, J. Teng, A. Danner, S. Zhang, and C. Qiu, “Visible-frequency metasurface for structuring and spatially multiplexing optical vortices,” Adv. Mater. 28, 2533–2539 (2016). [CrossRef]   [PubMed]  

18. F. Ding, Y. Chen, and S. I. Bozhevolnyi, “Metasurface-based polarimeters,” Appl. Sci. 8, 594 (2018). [CrossRef]  

19. E. Arbabi, S. M. Kamali, A. Arbabi, and A. Faraon, “Full stokes imaging polarimetry using dielectric metasurfaces,” arXiv preprint arXiv:1803.03384 (2018).

20. A. Pors, M. G. Nielsen, and S. I. Bozhevolnyi, “Plasmonic metagratings for simultaneous determination of stokes parameters,” Optica 2, 716–723 (2015). [CrossRef]  

21. J. B. Mueller, K. Leosson, and F. Capasso, “Ultracompact metasurface in-line polarimeter,” Optica 3, 42–47 (2016). [CrossRef]  

22. W. T. Chen, P. Török, M. R. Foreman, C. Y. Liao, W.-Y. Tsai, P. R. Wu, and D. P. Tsai, “Integrated plasmonic metasurfaces for spectropolarimetry,” Nanotechnology 27, 224002 (2016). [CrossRef]   [PubMed]  

23. H. Pahlevaninezhad, M. Khorasaninejad, Y.-W. Huang, Z. Shi, L. P. Hariri, D. C. Adams, V. Ding, A. Zhu, C.-W. Qiu, F. Capasso, and M. J. Suter, “Nano-optic endoscope for high-resolution optical coherence tomography in vivo,” Nat. Photonics 12, 540 (2018). [CrossRef]  

24. F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science 347, 1342–1345 (2015). [CrossRef]   [PubMed]  

25. E. Arbabi, A. Arbabi, S. M. Kamali, Y. Horie, and A. Faraon, “Controlling the sign of chromatic dispersion in diffractive optics with dielectric metasurfaces,” Optica 4, 625–632 (2017). [CrossRef]  

26. E. Arbabi, A. Arbabi, S. M. Kamali, Y. Horie, and A. Faraon, “Multiwavelength polarization-insensitive lenses based on dielectric metasurfaces with meta-molecules,” Optica 3, 628–633 (2016). [CrossRef]  

27. M. Khorasaninejad, Z. Shi, A. Y. Zhu, W.-T. Chen, V. Sanjeev, A. Zaidi, and F. Capasso, “Achromatic metalens over 60 nm bandwidth in the visible and metalens with reverse chromatic dispersion,” Nano Lett. 17, 1819–1824 (2017). [CrossRef]   [PubMed]  

28. W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nat. Nanotechnol. 13, 220 (2018). [CrossRef]   [PubMed]  

29. S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, M.-K. Chen, H. Y. Kuo, B. H. Chen, Y. H. Chen, T.-T. Huang, J.-H. Wang, R.-M. Lin, C.-H. Kuan, T. Li, Z. Wang, S. Zhu, and D. P. Tsai, “A broadband achromatic metalens in the visible,” Nat. Nanotechnol. 13, 227 (2018). [CrossRef]   [PubMed]  

30. J. R. Ong, H. S. Chu, V. H. Chen, A. Y. Zhu, and P. Genevet, “Freestanding dielectric nanohole array metasurface for mid-infrared wavelength applications,” Opt. Lett. 42, 2639–2642 (2017). [CrossRef]   [PubMed]  

31. A. Zhan, T. K. Fryett, S. Colburn, and A. Majumdar, “Inverse design of optical elements based on arrays of dielectric spheres,” Appl. Opt. 57, 1437–1446 (2018). [CrossRef]   [PubMed]  

32. S. I. Jafar-Zanjani, Samad, and H. Mosallaei, “Adaptive genetic algorithm for optical metasurfaces design,” Sci. Reports 8, 11040 (2018). [CrossRef]  

33. V. Egorov, M. Eitan, and J. Scheuer, “Genetically optimized all-dielectric metasurfaces,” Opt. Express 25, 2583–2593 (2017). [CrossRef]  

34. K. D. Donda and R. S. Hegde, “Rapid design of wide-area heterogeneous electromagnetic metasurfaces beyond the unit-cell approximation,” Prog. In Electromagn. Res. M 60, 1–10 (2017). [CrossRef]  

35. W. Ma, F. Cheng, and Y. Liu, “Deep-learning-enabled on-demand design of chiral metamaterials,” ACS Nano 12, 6326–6334 (2018). [CrossRef]  

36. R. G. Ghanem and P. D. Spanos, “Stochastic finite element method: Response statistics,” in Stochastic Finite Elements: A Spectral Approach, (Springer, New York, NY, 1991).

37. D. Xiu, Numerical methods for stochastic computations: a spectral method approach (Princeton University Press, 2010). [CrossRef]  

38. I. Babuska, R. Tempone, and G. E. Zouraris, “Galerkin finite element approximations of stochastic elliptic partial differential equations,” SIAM J. on Numer. Analysis 42, 800–825 (2004). [CrossRef]  

39. H. G. Matthies and A. Keese, “Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations,” Comput. Methods Appl. Mech. Eng. 194, 1295–1331 (2005). [CrossRef]  

40. F. Nobile, R. Tempone, and C. G. Webster, “A sparse grid stochastic collocation method for partial differential equations with random input data,” SIAM J. on Numer. Analysis 46, 2309–2345 (2008). [CrossRef]  

41. C. Chauviere, J. S. Hesthaven, and L. Lurati, “Computational modeling of uncertainty in time-domain electromagnetics,” SIAM J. on Sci. Comput. 28, 751–775 (2006). [CrossRef]  

42. R. Gaignaire, S. Clénet, B. Sudret, and O. Moreau, “3-d spectral stochastic finite element method in electromagnetism,” IEEE Transactions on Magn. 43, 1209–1212 (2007). [CrossRef]  

43. A. C. Austin and C. D. Sarris, “Efficient analysis of geometrical uncertainty in the fdtd method using polynomial chaos with application to microwave circuits,” IEEE Transactions on Microw. Theory Tech. 61, 4293–4301 (2013). [CrossRef]  

44. B. Peherstorfer, K. Willcox, and M. Gunzburger, “Survey of multifidelity methods in uncertainty propagation, inference, and optimization,” SIAM Rev. 60, 550–591 (2018). [CrossRef]  

45. M. B. Giles, “Multilevel monte carlo methods,” Acta Numer. 24, 259–328 (2015). [CrossRef]  

46. N. Georg, D. Loukrezis, U. Römer, and S. Schöps, “Uncertainty quantification for an optical grating coupler with an adjoint-based leja adaptive collocation method,” arXiv preprint arXiv:1807.07485 (2018).

47. J. P. Hugonin and P. Lalanne, “Reticolo software for grating analysis,” https://www.lp2n.institutoptique.fr. Online; accessed 29 January 2014.

48. A. Ben-Tal, L. El Ghaoui, and A. Nemirovski, Robust optimization, vol. 28 (Princeton University Press, 2009).

49. S. Lanteri, C. Scheid, and J. Viquerat, “Analysis of a generalized dispersive model coupled to a DGTD method with application to nanophotonics,” SIAM J. Sci. Comp. 39, A831–A859 (2017). [CrossRef]  

50. “DIOGENeS — a DG-based software suite for nano-optics,” (2018). [Online; accessed 28-December-2018].

51. V. A. Epanechnikov, “Non-parametric estimation of a multivariate probability density,” Theory Probab. & Its Appl. 14, 153–158 (1969). [CrossRef]  

52. I. Kim, G. Yoon, J. Jang, P. Genevet, K. T. Nam, and J. Rho, “Outfitting next generation displays with optical metasurfaces,” ACS Photonics 5, 3876–3895 (2018). [CrossRef]  

53. R. Sawant, P. Bhumkar, A. Zhu, P. Ni, F. Capasso, and P. Genevet, “Mitigating chromatic dispersion with hybrid optical metasurfaces,” Adv. Mater. 31, 1805555 (2018).