Fabrication robustness in BIC metasurfaces

1 Introduction

Nanophotonics has been a vibrant topic of research in the last two decades, focusing on the interactions of subwavelength objects with light and enabling the nanoscale concentration of electromagnetic fields beyond the optical diffraction limit [1, 2]. Controlling light at the nanoscale is traditionally associated with surface plasmon polaritons supported on metallic nanoresonators [3, 4]. However, since plasmonic resonators suffer from strong dissipative losses due to the intrinsic damping of the constituent metals, the field of resonant all-dielectric photonics provides novel possibilities for tuning the electric and magnetic responses with low intrinsic material losses [5], [6], [7], [8], [9], [10], [11]. In particular, photonic bound states in the continuum (BICs) have recently emerged as an exciting toolkit for the design and experimental realization of all-dielectric platforms with sharp resonances and strong near-field enhancements [12], [13], [14], [15], [16].

In this context, BIC-based metasurfaces composed of unit cells with broken in-plane inversion symmetry have unlocked precise control over the linewidth and hence the Q-factor of their resonances [17]. Combined with tunability over a wide frequency range [18], this approach has launched a broad spectrum of applications ranging from biochemical sensing [19], [20], [21], [22], [23], [24] to higher harmonic generation [25], [26], [27], and lasing [28]. While extremely high Q-factors are straightforward to achieve in numerical simulations [29, 30], the experimental realization of such metasurfaces through nanofabrication invariably introduces geometrical perturbations that can significantly broaden the resonance linewidths and therefore greatly reduce the experimental Q-factors, limiting practical applications.

Here, we propose a combined numerical/experimental methodology to assess and compare the reduction of the radiative Q-factor of different BIC unit cell designs for a range of fabrication tolerances. Specifically, we investigate and compare three established all-dielectric BIC-based unit cell designs [17]: tilted ellipses, asymmetric rods, and split rings. Our method provides invaluable insights on the robustness of BIC-based resonant metasurfaces prior to fabrication, enabling more rapid prototyping and structure discovery as well as improved optical performance and durability for practical applications such as field-deployed sensing of pathogens [31], environmental monitoring [32], or ultrathin optical elements [33].

2 Results and discussion

To establish our fabrication robustness model, three different BIC-based metasurface unit cells are designed numerically to each exhibit a sharp quasi-BIC resonance at 750 nm with a Q-factor of 200. The resonators are composed of silicon, which is chosen for its proven nanophotonic properties and low losses above 650 nm. We select the visible wavelength range for the validation of our methodology, since the effect of errors in the fabrication process play a critical role there.

For each unit cell design, two geometric parameters with strong influence on the BIC resonance are identified by estimating how fabrication-induced variations will likely affect the resonator shape. As an example, in the asymmetric double rod geometry, the width w of the rods and the length l of the longer rod are selected (Figure 1A). The set of these two parameters can be used to approximate the resonance changes, since we expect that dose variations during the structure patterning process will lead to similar length modifications for both rods. For the tilted ellipse design, we select the major and minor axes of the ellipses while keeping the tilting angle constant. In the case of the split ring geometry, the width of the ring and the size of the two gaps are chosen.

Figure 1:

Numerical modelling of resonance broadening. (A) Sketch of a BIC unit cell design and identification of two geometric parameters with strong influence on the resonance. (B) A 5 by 5 grid of simulation points with variation of the parameters by ±7.5% and ±15% is interpolated in 1000 steps, resulting in a finely gridded parameter space of geometric variation. (C) Selection of points out of Gaussian random distribution for two different standard deviations σ ₁ and σ ₂ > σ ₁. (D) Spectral averaging yields a broadened resonance spectrum corresponding to a fabricated metasurface with random geometric variation.

One straightforward way of modelling the fabrication-induced effects on the optical response of the metasurface is to introduce random geometrical variations for each unit cell independently. However, this approach would lead to a prohibitively large simulation domain for common metasurface arrangements, which often require more than 100 unit cells to support their respective resonant modes [10, 21]. Instead, our methodology simulates the spectra of full unit cell arrays with varying geometrical parameters and uses statistical averaging to predict the final optical response of the fabricated metasurfaces. In order to implement this method numerically, the two chosen geometric parameters are varied by ±7.5% and ±15%, respectively, to construct a 5 by 5 grid of equidistant points, where the initial parameter values of the design are located at the center (Figure 1B), drastically reducing the simulation overhead. Reflectance spectra for all 25 parameter combinations are then numerically simulated and linearly interpolated using n = 1000 steps. This process creates a parameter space with a fine grid of simulated spectra, allowing us to capture even small variations of the geometric parameters. To model the geometric variations in a fabricated sample, N = 5000 random points are selected from a Gaussian random distribution over the parameter space, mimicking the experimental distribution of slightly different geometrical parameters (Figure 1C). The spectra corresponding to the randomly drawn points of the parameter space are the naveraged to obtain an approximation of the optical response of a fabricated metasurface (Figure 1D).

Notably, even low standard deviations as low as σ = 1 nm or σ = 3 nm can have a strong influence on the resonance amplitude and the Q-factor. For example, the resonance amplitude of the ellipse design for σ = 1 nm is at around 0.8, whereas it is at 0.7 for the rods, illustrating the markedly different impact of the same amount of fabrication tolerance on BIC unit cell designs with the same resonance position and Q-factor (Figure 2A–C). This behavior clearly shows that the magnitude of resonance broadening is correlated with each individual unit cell design. Figure 2D displays the Q-factor of the different BIC designs for standard deviations of the geometrical variation ranging from σ = 0 nm to σ = 5 nm. Notably, the ellipse design retains the highest Q-factor over the entire variation range, while the other designs fall off more quickly (Figure 2D). Fitting an inverse quadratic function to the Q-factors shows good agreement and confirms the proportionality of Q(σ) ∼ 1/σ ² previously found for photonic crystal cavities [34]. Interestingly, the ordering of the different designs by their Q-factors is maintained for all standard deviations as well, highlighting that the fabrication robustness is an underlying property of each BIC design.

Figure 2:

Broadening model predictions for common BIC geometries.

(A–C) Simulated quasi-BIC resonance spectra and modelled spectra for the different resonator designs and two standard deviation values for the geometric variation. For σ = 3 nm both the resonance intensity and the Q-factor are significantly reduced compared to the simulation. (D) Calculated Q-factor against standard deviation values of the geometrical parameters up to σ = 5 nm (points). Inverse quadratic fits (lines) showing good agreement with the Q-factors of the broadened resonances. The tilted ellipse design exhibits higher resonance sharpness over the entire variation range compared to the other designs.

The three BIC metasurface geometries are then realized experimentally in amorphous silicon through electron-beam lithography and directional reactive ion etching (for details see Section 4). To quantify the fabrication tolerance experimentally, we capture large sets of scanning electron microscopy (SEM) images for each BIC unit cell design and perform image processing to assess the distribution of the geometric parameters. Specifically, in order to achieve a high level of uniformity and reproducibility for the parameter extraction, we implement an image fitting algorithm to overlap the designs with the SEM images. Images are first up-scaled eightfold for higher precision and undergo a multistep denoising process through applying a series of filters (Figure 3, top). First, a bilateral filter removes fine noise from the image while preserving the edges of larger objects. The next step is a set of erosion and dilation filters, resulting in overall brighter and smoother features in the images, bringing out the resonator geometry in more detail. Subsequently, a threshold is applied to transform the image into filled black and white image areas. Finally, the design can be fitted iteratively by creating a mask with varying geometrical parameters and maximizing the number of overlapping pixels between mask and image to obtain the optimal parameter values. For better statistics, this fitting algorithm is applied to 10 images containing 12 unit cells for each design, providing higher accuracy as well as consistency and reproducibility compared to a manual measurement of the parameters from the SEM images. The resulting distributions of extracted geometric parameters clearly show Gaussian profiles for all unit cell designs, confirming the assumption of the theoretical model (Figure 4).

Figure 3:

Methodology for comparison of theory and experimental realization. Evaluation of fabricated metasurface geometric parameters by fitting of resonator designs, which are used as input for the broadening model to calculate a resonance that mimics random fluctuations in a fabricated sample. The modelled resonance is compared to the measured spectra of the fabricated metasurface.

Figure 4:

Evaluation of fabrication-induced variation of resonator geometry.

Structural analysis output of resonator dimensions for the tilted ellipse (A), split ring (B) and double rod (C) unit cell designs, all of which show a Gaussian distribution of their geometric parameters.

Comparing the mean dimensions from the fabricated structures to the design values (Figure S1, Supplementary material), only the asymmetric rod design shows a notable discrepancy from the as-designed geometry with a 12 nm difference in rod width (Figure 4C), with all other designs closer to the desired parameters. Since the metasurface arrays for all three BIC unit cell designs are fabricated on one substrate, we average the standard deviations of all designs to obtain a unified measure of the geometric tolerances in fabrication, which yields σ _avg = 2.4 nm. However, it is important to note that if metasurfaces from different fabrication runs are to be compared, the standard deviation of each sample should be determined separately. Interestingly, in the split ring design, the fitted parameters show a correlation between gap size and width, where overall up- and downscaled structures are more common. This effect is absent in the other BIC designs due to a simpler unit cell structure of two resonator elements with larger spatial separation.

The average geometric variation σ _avg determined from the SEM images is used as an input for our numerical broadening model (Figure 3, bottom), allowing us to predict realistic BIC resonance properties from first-principle electromagnetic simulations. In order to improve the accuracy of the prediction, we take into account some intrinsic sources of error related to the SEM imaging. For example, the imaging quality is affected by a small amount of lens aberration, which leads to a magnification of objects at the image corners and thus to an increased observed size, while objects in the center remain almost unchanged. In turn, this increases the standard deviation of geometrical parameters independent of the actual variation of the fabricated metasurface, since slight blurring at the edge of objects can already influence the extracted parameter on the order of a few nanometers. Furthermore, because the SEM images are taken at a resolution of 2048 by 1536 pixels and one pixel corresponds to around 1 nm at 60,000× magnification, there is a natural limit to the accuracy of this methodology, which we estimate at σ _im = 1 nm. Electrostatic charging of the metasurface during imaging can similarly lead to slight variations of the exposure between different images taken for the same BIC design. To quantify this effect, we calculate the standard deviation of the mean geometrical parameters between different images and obtain a correction factor of σ _c  = 0.5 nm. These considerations are important so that the true fabrication tolerance is not overestimated, which would lead to a very broad resonance (see Figure 2A–C for σ = 3 nm).

Combining these effects, we obtain a total standard deviation of σ _tot = σ _avg − σ _im − σ _c  = 0.9 nm. The modelled resonances with this standard deviation value show a reduction in Q-factor of around 40–110 compared to the as-designed value of Q = 200 for all investigated unit cell geometries (Figure 5A), highlighting the crucial importance of taking these broadening effects into account. Significantly, the tilted ellipse design stands out as the BIC unit cell that retains the highest Q-factor, while the other designs have broader resonances and thus lower Q-factors. In experimental transmission measurements (Figure 5B, for details see Section 4), both the resonance shapes and amplitudes are in excellent agreement with the predictions of our broadening model apart from a slightly lower resonance wavelength. This offset is most likely related to a small discrepancy between the dielectric function of silicon as used in the simulations and the actual optical properties of the thin film as deposited in experiments. Nevertheless, the resonance Q-factors of the fabricated metasurfaces and our model show good agreement and demonstrate the possibility to assess the influence of random geometric variation during fabrication on the resonance, enabling the accelerated discovery, and development of robust BIC-based nanophotonic systems.

Figure 5:

Comparison between broadening model output and measured optical response of the BIC metasurfaces.

(A) Modelled BIC resonance spectra for the standard deviation value of geometric variation derived from SEM image analysis. (B) Measured spectra of the fabricated metasurfaces showing excellent agreement with the model.

3 Conclusion

In conclusion, we have developed a method for modelling the effects of unavoidable fabrication tolerances on the optical performance of BIC-based metasurfaces. For each given unit cell design, our methodology first maps the parameter space of fabrication-induced geometrical variations using numerical simulations of the metasurface reflectance. Subsequently, a statistical approach is used to predict the modified reflectance response as a function of the standard deviation of the geometrical parameters. Comparing three prominent designs (tilted ellipses, asymmetric double rods, and split rings), we find a strong disparity in fabrication robustness behavior, with the tilted ellipse design showing more than 50% higher Q-factors compared to the split rings, demonstrating the crucial impact of such geometrical broadening effects. The numerical predictions provide important guidelines for the design of BIC-based unit cells and are in excellent agreement with optical measurements on fabricated metasurface samples, where the experimental variation of the geometrical parameters is extracted using sophisticated analysis of sets of SEM images. Since our methodology only requires simulated optical spectra as input, it is not limited to BIC-driven geometries, but can be applied to all nanophotonic systems where high Q-factors are needed for critical applications, enabling the a-priori design and discovery of new classes of fabrication-robust resonant systems in fields ranging from biomolecular sensing and environmental monitoring to tailored light/matter interaction.

4 Methods

Numerical simulations: Simulated reflectance spectra of the different metasurface designs were obtained using the finite-element frequency-domain Maxwell solver included in CST Studio Suite 2020 (Dassault Systèmes), with periodic boundary conditions applied at both sides of the unit cell.

Metasurface fabrication: The all-dielectric resonator geometries were fabricated from a 100 nm thick layer of amorphous silicon deposited onto a fused silica substrate by plasma-enhanced chemical vapor deposition (PECVD). For the lithography step, the sample was first spin-coated with a layer of PMMA 950K A2 photoresist followed by a highly conducting polymer layer (ESpacer 300Z) to avoid charging effects during exposure. Subsequently, the three BIC unit cell designs were patterned using electron-beam lithography (Raith eLINE Plus) with an acceleration voltage of 20 kV and a 7.5 μm aperture. After development for 50 s in 3:1 MIBK:IPA solution, a hard-mask consisting of 14 nm SiO₂ and 20 nm Cr was deposited by evaporation, followed by lift-off in Microposit Remover 1165 for 2 h at 80° C. The pattern was then transferred into the silicon layer by reactive ion etching (RIE). Lastly, the chrome layer was removed by wet etching and the SiO₂ layer with an additional RIE step. The obtained Si resonator metasurfaces on a silica substrate were then imaged for geometrical analysis using the same Raith eLINE Plus system.

Optical measurement: The spectral characterization of the fabricated metasurface samples was carried out with a WiTec optical microscope in transmission mode using a collimated Thorlabs OSL2 white light source for illumination. The light is collected with an Olympus 10× objective (NA = 0.25), normalized to a background spectrum on the sample and converted to reflectance by 1 − transmittance.

SEM image fitting and geometrical analysis: For evaluation of the SEM images and the dimensions of the unit cells, a Python script was written based on the OpenCV library for computer vision. In order to detect individual resonators, the input image was filtered, starting with a bilateral filter for noise reduction followed by a set of dilation and erosion filters. This series of filters allows for thresholding of the image to isolate the features from noise and facilitate corner recognition, which is needed to get the resonators position in the image as well as its size to function as starting parameters for the geometry fitting. For this process, a mask of the target resonator was generated, placed on top of the processed resonator image, and the overlapped area calculated by counting the overlapping pixels, the number of which was then returned to the optimizer. Finally, the fit parameters were converted from pixels to nanometers by reading out the length of the scale bar in the SEM image.