Finite-Size and Illumination Conditions Effects in All-Dielectric Metasurfaces

Abstract

Dielectric metasurfaces have emerged as a promising alternative to their plasmonic counterparts due to lower ohmic losses, which hinder sensing applications and nonlinear frequency conversion, and their larger flexibility to shape the emission pattern in the visible regime. To date, the computational cost of full-wave numerical simulations has forced the exploitation of the Floquet theorem, which implies infinitely periodic structures, in designing such devices. In this work, we show the potential pitfalls of this approach when considering finite-size metasurfaces and beam-like illumination conditions, in contrast to the typical infinite plane-wave illumination compatible with the Floquet theorem.

1. Introduction

All-dielectric metasurfaces have recently emerged as a valid and promising alternative to their plasmonic counterparts. Their lower dissipative losses and thermal heating at optical frequencies allow working at higher pump power for nonlinear frequency generation and the performance of non-invasive sensing of biological samples [1,2,3,4,5,6,7,8,9,10,11,12]. Moreover, in all-dielectric resonators, electric and magnetic resonances can be used to engineer near and far fields, thus providing more degrees of freedom in the design of optical devices [1,13,14,15,16]. In this framework, III-V compounds have been used for second-harmonic generation, due to their non-centrosymmetric structure and their high second-order nonlinear coefficients [2,3,4,5,17,18,19,20,21,22,23]. Strong third-harmonic generation has been demonstrated instead with Si-based platforms [24,25]. Lately, great efforts have been made to improve the design of optical nanoresonators, either in the limit of single isolated nanoantennas or in the opposite scenario where perfect periodicity is assumed and thus the finite-size real-world structure is approximated by an infinite structure [26,27,28]. The two opposite situations are required in order to limit the numerical burden in the modelling when resorting to the full-wave solution of Maxwell’s equations. Despite the fact that this has given very helpful guidelines, these approaches fail in the accurate design of a device, because they unavoidably neglect the effects that are due to the size of a real device, which is necessarily finite and is often pumped by input beams exciting a limited portion of the device (see Figure 1). Such limitations become particularly relevant when considering metasurfaces supporting bound states in the continuum (BIC), which are non-radiating resonant modes in an open system that are not coupled to the radiating channels and exist only in infinite structures [29]. A common approach to obtain extremely high Q-factor resonances is to introduce a defect into a BIC-supporting metasurface (e.g., by breaking its symmetry), thus realizing a so-called quasi-BIC. This approach is mainly dictated by the fact that currently available simulation tools allow only infinite structures to be inspected, and thus it is necessary to break the symmetry in order to couple to such modes. Quasi-BICs, however, are expected to arise naturally by simply considering the finite extent of a real metasurface.

The necessity of using a finite-size model arises from the fact that the currently adopted models simulate an elementary unit cell and exploit the periodicity of the structure through Floquet’s theorem to obtain the optical parameters of an infinite extended metasurface. The results of these models are very good, but they do not take into account the finiteness of the structure and the border effects. Furthermore, when Floquet’s theorem is employed, the only possible excitation is a plane wave, unlike in experiments, where a Gaussian beam is usually employed. The combination of the two types of models can help to better foresee the behaviour of an actual metasurface in the experimental phase.

The differences between an infinite model and a finite one are highlighted in [30,31], where the situation is analysed with a dipole coupling model for a 2D array of nanoparticles. However, dielectric resonators often support higher-order resonances which cannot be described in the dipole approximation. In order to overcome this limitation, it is necessary to account for higher-order multipolar components, which are of particular interest for dielectric metasurfaces and the study of metamaterials in general [26,32,33,34,35]. Babicheva and Evlyukhin [36] provided an analytical model to include a quadrupolar contribution, but the authors employed an infinite structure approximation. In [37], the authors employed the superposition T-matrix scheme with the Ewald sum formalism to describe periodic structures with large unit cells. In [38], the authors studied finite-size effects in plasmonic metasurfaces, highlighting the full-width half-maximum dependence of resonances upon the number of resonators.

Our study aims to highlight the importance of the spatial finiteness of a metasurface supporting high-order multipole components. For this purpose, we employed SMUTHI [39], which implements a linear space-limited model, and tested its performance considering a metasurface composed of silicon nanopillars on a silicon oxide substrate. We analysed the total reflectivity behaviour as a function of the incident wavelength and compared it with the results obtained with the COMSOL infinite periodic model. The reflectivity of the metasurface strongly depends on the number of elements and converges to the infinite array approximation for a large array size, thus highlighting the importance of considering the finiteness of the metasurface. Furthermore, this approach allows the metasurface to be illuminated with a Gaussian beam of a different size under the paraxial approximation, which is not possible when employing Floquet’s theorem.

2. Methods

We simulated an all-dielectric metasurface composed of Si cylinders on a ${\mathrm{SiO}}_2$ substrate, both with a periodic model using COMSOL (exploiting the periodicity via Floquet’s theorem) and a space-limited model using SMUTHI, which allowed us to simulate very large finite structures due to the lower computational power required. To test the finite model, we calculated the reflectivity (R) for both models, varying the wavelength of the incident wave $\lambda$ from 1250 nm to 1800 nm, in steps of 10 nm.

The finite structure was composed of a maximum of 169 cylinders lying on the substrate in a 13 × 13 configuration, as seen in Figure 2. The refractive indices were 1 for the air, 3.45 for the cylinders, and 1.45 for the substrate, the last two values being the silicon and silicon oxide indices at the working frequencies. The cylinder heights were 428 nm and their radii were 214 nm. The distance between two adjacent cylinders was 856 nm. These dimensions were chosen to create a structure that supports resonant modes at the working frequencies on two levels: inside single cylinders and globally due to the periodicity. The layers consisted of a substrate with thickness $2\lambda /n_{\mathrm{substrate}}$ and an air layer $2\lambda$ thick, and the domain on the plane of the cylinders was a square with a side length of 15,500 nm. Although SMUTHI allows the simulation of lossy media, in this study we employed a fixed real value for the refractive indices, because in the spectral region of interest the silicon and silicon dioxide absorption coefficients and the dispersion were very small and thus could be neglected.

The incident wave used for the COMSOL periodic model is always a plane wave polarized along the y-axis, while with SMUTHI it is also possible to simulate an incident Gaussian beam. We set an incident wave propagating from the air layer towards the substrate at normal incidence, at a wavelength which varied as previously described. The polarization of the wave was linear, with the electric field oscillating along the y-axis and an initial intensity of $E_0=1 \mathrm{V}{\mathrm{m}}^{-1}$. The maximum value of the electric field was at the air–substrate interface. To calculate the reflectivity, we wrote a MATLAB script, which takes as an input the SMUTHI results related to the scattered field from the cylinders and the reflected one from the substrate. Then, we calculated the flux of the Poynting vector through three surfaces $xy$, $yz$, and $xz$ at the border of the domain described above, inside the air layer.

The formulas used to evaluate the Poynting vector components were the following [40]:

$$\begin{array}{c}\hfill p_x=\frac{|E_y{|}^2+{\left|E_z\right|}^2}{2Z}{p}_y=\frac{|E_x{|}^2+{\left|E_z\right|}^2}{2Z}\\ \\ \hfill p_z=\frac{|E_x{|}^2+{\left|E_y\right|}^2}{2Z}{p}_0=\frac{|E_{0y}{|}^2}{2Z}=\frac{|E_0{|}^2}{2Z}\end{array}$$

where $E_x$, $E_y$, and $E_z$ are the field components given by the combinations of the scattered and the reflected waves, $E_0$ is the initial intensity, and Z is the characteristic impedance of vacuum. Integrating these values along the respective surfaces, it is possible to calculate the reflectivity:

$$\begin{array}{c}\hfill R=\frac{P_z+2P_x+2P_y}{P_0}\end{array}$$

where capital $P_x$, $P_y$, and $P_z$ are the powers relative to the power densities with the same subscript and are obtained through surface integration over the simulation domain walls. This procedure was repeated for each wavelength to obtain the full spectrum.

3. Validation

To develop the space-limited model we used SMUTHI, a recently developed open-source Python package based on numerical methods, to speed up the computation of light scattering problems [41,42] (an analytical treatment is reported in Appendix A).

To determine how the software performs in our situation, we developed two simple models of the metasurface, one using COMSOL and the other using SMUTHI, to compare the outputs for a structure composed of nine cylinders arranged in a 3 × 3 configuration around the origin. The parameters were the same as those set in the previous section, except for the domain width in the plane of the cylinders, which was changed to 5000 nm, and the incident plane wave that was fixed at 1200 nm. In COMSOL, we simulated a larger domain: a cube with a side length of 6500 nm with a perfectly matched layer all around, in order to prevent unwanted reflections from the side. We performed the integration on the same surface as that employed in SMUTHI. We employed the iterative GMRES solver in COMSOL and set a mesh size of ${\lambda }_e/7$, where ${\lambda }_e$ is the effective refractive index of the material.

Figure 3 shows the electric fields at a distance of 1430 nm from the last cylinder, obtained in COMSOL (left) and SMUTHI (right), respectively.

A good agreement was found using a multipolar expansion with $l=m=3$ for the cylinders, which resulted in very similar field distributions and integrated reflectance values of:

$$\begin{array}{c}\hfill R_{\mathrm{COMSOL}}=0.0649\\ \hfill R_{\mathrm{SMUTHI}}=0.0542\end{array}$$

Then, for the transmission coefficient:

$$\begin{array}{c}\hfill T_{\mathrm{SMUTHI}}=0.9482\end{array}$$

and thus the condition $R+T=1$ is almost respected (this must hold because the refractive indices are real) [41,42].

We underline the times required for the two simulations: for the structure reported in this chapter, COMSOL required 16 min with a reduced geometry, obtained by exploiting symmetry considerations, and 74 min with the full structure, to solve the problem; on the other hand, SMUTHI required 3 min (PC specifications: CPU Intel^© Xeon^© E5-2650 0, 2 GHz, 8 core. RAM DIMM-DD3 1600 MHz, 64 GB. GPU NVIDIA Quadro K2000). Although COMSOL provides a more complete simulation by computing the field in the full volume of the structure, we only need the field on three specific surfaces, and thus SMUTHI is much more time-efficient.

4. Results

Figure 4a shows the plots of R, calculated on the same surface, as a function of the wavelength, for different metasurface sizes. When the number of cylinders increases, the spectral features related to the multipolar resonances inside the cylinders and their interplay become more evident. We note that for smaller arrays, the covered area inside the integration surface is much smaller, thus resulting in smaller reflectivity values.

In Figure 4b–d we compare the reflectivity as a function of the wavelength obtained with the two models.

For the finite-size model, the values are relative to the 13 × 13 cylinder structure and are normalized to have a maximum equal to one. The normalization is only for comparison purposes. Indeed, in the finite structure case, when the integration domain is enlarged, the unpatterned surface increases, thus reducing the reflectivity. However, the relative size of the peaks does not change as this is only due to the cylinders. The normalization of the finite structure to the infinite structure allows the relative intensities of the peaks to be compared, and thus the finite-size effect to be evaluated. We compared the spectra with an incident plane wave for both structures and then changed the wave for the finite model into a Gaussian beam 15 μm and 6 μm wide. Our simulations showed that as long as the paraxial simulation holds and the beam covers the metasurface, the plane wave approximation deviates slightly from the real situation.

The normalized behaviours are very faithful to the COMSOL behaviours, either with an incident plane wave or with a Gaussian beam. In particular, we see two peaks at 1610 nm and 1330 nm, which are at exactly the same wavelengths for the two models. The resonance at 1610 nm is related to the magnetic dipole inside the cylinder and is also present in the scattering cross section in the isolated structure. The different values around the shorter wavelength peak are due to the finite extension of the structure and the excitation beam; therefore, this feature can only be highlighted through a finite-size model. Indeed, a Gaussian beam can be described by considering a larger amount of k-vectors, which is not possible with the plane-wave excitation employed in COMSOL. This is further evident from Figure 5, where we report also the results obtained considering the flux through the walls around the central cylinder in the finite-size model and the reflectivity of the infinite periodic structure. It is interesting to note that the characteristic of the bigger structure is more similar to the periodic model. As expected, the relative weight of the resonant peak due to the interaction between the cylinders (around 1330 nm) and the single-pillar magnetic dipole (1610 nm) increases with the array size. We emphasize that the resonance at 1610 nm is related to the magnetic dipole in the isolated pillar, thus it is well reproduced even for small arrays and beam sizes. Thus, these results show that it is essential to consider the size of the system and the excitation, to correctly describe a real finite-size system. For this purpose, it should be noted that the size of the structure compared to the wavelength is more important than the number of unit cells considered in the simulation. However, the size of the unit cells plays a key role in determining the interaction strength between the resonators. Indeed, if the cylinders are well separated from each other, they act as isolated sources, and the scattering properties can be evaluated by computing the electromagnetic field of an isolated pillar and considering the system as an array of antennas.

To further prove the importance of finite-size effects in conjunction with resonant behaviour, we consider a structure supporting bound states in the continuum. Figure 6 shows the reflectivity as a function of the wavelength for a suspended array of cylinders, which was previously reported for sustaining a BIC [43]. The Gaussian beam waist was kept fixed at 10 μm, while the number of resonators was increased. It is clear that increasing the number of resonators reduces the line width of the Fano resonance associated with the quasi-BIC. We highlight that in [43], the authors introduced a gap in the middle of the cylinders in order to couple to the resonant mode, while in our case this was not necessary since both the illumination and the structure were finite and the modes could couple to the radiating channels.

5. Conclusions

In conclusion, we built a space-limited model with SMUTHI to simulate a real metasurface, to determine the border effects. First, we tested the software performance in the presence of the metasurface and verified the conservation of energy in the simulated domain ($R+T=1$), as well as the agreement with the COMSOL simulation. The resulting reflection coefficients R were similar for the two structures, and the produced scattered fields were very similar.

We developed the finite model by introducing an approximation for the scattered field and the refractive indices. Moreover, with the finite-domain model, it was possible to highlight a difference in the higher-frequency peak related to collective excitation. We could also simulate an incident Gaussian beam, which is not possible in the periodic case. SMUTHI allowed us to simulate the wavelength sweep for a large finite structure with 13 × 13 cylinders, which is impractical to replicate in COMSOL, and also allowed us to simulate the configuration with 3 × 3 cylinders 5.3 times faster. This approach also allowed us to employ a non-periodic distribution of resonators, which is not possible when exploiting the Floquet theorem. Furthermore, it allowed us to study real situations involving BIC, which are intimately related to the finiteness of the system. The proposed approach may be further extended by analogy with [44], where the authors employed multipolar decomposition to study free-standing single nanoparticles, to simulate nonlinear frequency generation in finite-size metasurfaces.