1. Introduction

Metasurfaces (MSs) operate by superposing the emitted waves of individual nanoscale elements placed in a large-scale array, enabling complex shaping of optical beams in compact devices [1]. Recent works have demonstrated the use of MSs as lenses [2] to achromatically focus light on either planar substrates [3] or optical fiber [4], to create, manipulate, and detect new optical states, such as bound states in the continuum [5] or non-classical states of light [6]. Moreover, MSs have been used for highly efficient nonlinear frequency conversion [7], and for the simultaneous control of the intensity, polarization and phase of optical beams [8,9]. A common embodiment of a MS is the dielectric metasurface array, consisting of a periodic array of high refractive index (RI) contrast cylindrical nanopillars placed on a supporting substrate [10–14]. Here, the individual nanopillar acts as a waveguide-based antenna with tailored radiation properties that, when combined in a large-scale lattice, provide the desired MS functionality. In this geometry, the transmission properties can be controlled by changing the physical dimensions (diameter, height) of the cylinders, as well as by changing the properties of the array (i.e., lattice pitch and geometry) in which they are arranged; we note that while changing the spectral response by tuning the cylinders is a well-established approach, tuning the response by changing the lattice properties has been much less commonly exploited, even though it represents a important degree of freedom in designing MSs.

To model how such a MS will behave in a particular configuration, numerical tools such as the Finite Element (FE) method [15] or Finite-Difference Time Domain (FDTD) method [16] are usually used. However, it is often useful to supplement these full numerical simulations with analytical or semi-analytical models that can quickly and effectively predict important quantities such as accumulated phase or transmission. These models can be used directly as a tool in guiding fabrication, and can also lead to physical insight which further informs the design. The most commonly-used of these analytical models is to approximate the transmission through the array by the transmission of the excited mode through a single isolated cylinder. This approach is directly related to the well-understood physics of modes in step-index optical fibers, in which the modal dispersion can be found by solving a transcendental equation [17,18]. This model is explicitly limited to sparse arrays in which the effect of the array-type arrangement is negligible. However, as shown in the context of photonic crystal fibers [19,20], arranging cylindrical waveguides in a lattice is known to have a substantial effect on the modes in terms of phase and transmission, even for dilute systems of elements [21,22]. The fact that the cylinders are finite in height further complicates the situation, because the 3D layer of cylinders gives rise to discrete resonances that can dramatically change the transmission [23]. Semi-analytic approaches that can handle this behaviour include the transfer matrix method [24], modal layer formulations [25] and layer impedance modelling [26]; these methods take into account both the geometry of the cylinder array as well as the finite cylinder height, but become very similar in scale and complexity to full numerical simulations of the entire structure.

Here, we present a semi-analytical model, which we denote the lattice dipole model, that can be used as an extension of the single isolated cylinder model to include the effects of the array as a whole. This model takes into account not only the coupling between adjacent cylinders but also includes coupling between all the cylinders in the array. This is achieved via the inclusion of lattice sums into the dispersion equation for a single cylinder; these sums depend solely on the array parameters, and can be computed in a straightforward way [27]. Without limiting the generality of the model, we focus in this work on square arrays of different pitches, and show how the lattice parameters influence the phase transmission through the MS. In comparing with full 3D FE simulations, we find that the phase predicted by the lattice dipole model is much more accurate over a broad range of geometric parameters compared to the isolated cylinder. Furthermore this model can be used in conjunction with a Fabry-Perot model to design MS geometries that have close-to-unity transmission and a full $2\pi$ phase shift. We also explore the physics of the transmission through the MS using a combination of a full 3D simulation and a 2D mode-based model, and use this to gauge the limits of our analytic model. We find that the lattice dipole model is accurate up to the appearance of sharp resonant features associated with the lattice resonances of the 3D MS layer; we show that the regime where these resonances appear can be predicted using a just one transcendental equation. The physical insight gained from the combination of the analytical model and the FE simulations gives firm guidelines for the design of nanopillar-based MSs.

2. Approximate models

We consider an array of dielectric nanocylinders with refractive index $n_1$, height $h$ and radius $a$, aligned parallel to the $z$-axis and arranged in a two-dimensional array with lattice spacing $\Lambda$ (see Fig. 1(a)). Generally, the lattice can have any rectangular or triangular arrangement, with square and hexagonal lattices being the most common; throughout this work, we consider a square array of silicon nanopillars surrounded by air and located on a glass substrate. We denote the refractive index of the material in the MS layer surrounding the cylinders by $n_2$, and that of the substrate $n_s$. The cover material above the layer of cylinders we consider to be infinite in extent and of refractive index $n_c$. Light of frequency $f = \omega /2\pi$ or wavelength $\lambda =c/f$ is incident normally from the substrate side of the array; unless otherwise specified, we assume an incident free-space wavelength of $\lambda =1.0\mu m$, which defines the indices for the different materials at the operating wavelength $\lambda$ ($n_{\rm s}=1.5$, $n_{\rm c}=1.0$, $n_{\rm 1}=3.6$, $n_{\rm 2}=1.0$).

 

Fig. 1. (a) Illustration of the metasurface structure, consisting of a 2D array of cylinders of refractive index $n_1$ placed on a substrate with refractive index $n_s$. (b) Illustration of the lattice dipole model, in which the accumulated phase and transmission through the layer is determined by the guided mode of the central cylinder (blue) of the array.

Download Full Size | PDF

2.1 Isolated cylinder model

In the isolated cylinder model, the transmission through the periodic array is computed under the approximation that the light couples completely into one waveguide mode while passing through the array layer. The usual assumption, derived from the symmetry of the fields, is that the mode couples to the fundamental cylinder mode (hybrid $\rm {HE}_{11}$-mode), which is quasi linearly polarized and dipolar in both axial field components $E_z$ and $H_z$ (one example of such a dipolar distribution is shown in Fig. 2(b)). The effective mode index can be found for a given frequency by finding the propagation constant $\beta$ that solves the transcendental equation

Here $J_1(z)$ is the first-order Bessel function of the first kind, and $K_1(z)$ is the first-order modified Bessel function of the second kind. The parameters $U$ and $W$ are the dimensionless transverse wavenumbers defined as [17]

 

Fig. 2. Details of the model discussed here. (a) Working principle of the individual nanopillar element (light blue). The purple lines correspond to incident and scattered light and the excited mode. The plot shows the phase accumulated while the mode propagates through the cylinder. The yellow part on the right illustrates the Fabry-Perot oscillation. The two right-handed contour plots show selected spatial field distributions of the fundamental mode (top: longitudinal magnetic field component $H_z$. bottom: Poynting vector $S_z$). Note the dipolar-type distribution of the longitudinal component. (b) Geometry of the periodic array used in the lattice dipole model. Here a square lattice is shown, which forms the underlying geometry used in the discussion of this work; more general arrays can be generated by adjusting the lattice vectors $\mathbf {R}_p$.

Download Full Size | PDF

2.2 Lattice dipole model

The isolated cylinder model neglects contributions to the effective index from the cylinders in the surrounding lattice. One can take these into account by computing the effective index of the entire 2D array of cylinders infinitely extended along the $z$-direction. For such a computation higher-order multipoles need to be included in the field expansion, which accounts for the light scattered by neighboring cylinders. The general expansion for the longitudinal fields (along the $z$-direction) in the medium external to the cylinders is [28]

Here the functions $I_m(z)$ and $K_m(z)$ are the modified Bessel functions of the first and second kind; the expansion in terms of modified Bessel functions is needed for fields that are evanescent in the medium outside the cylinders. The $I_m(z)$ functions, unlike the $K_m(z)$ functions, are non-singular at the origin; in general, field sources are given by singularities at the source point, and so the singular part of the field expansion in (3) given by the $K_m$ terms represents sources from the central cylinder, whereas the non-singular terms given by the $I_m$ represent fields from all the other cylinders in the 2D array. The coefficients $a_m^{\rm E,H}$ can therefore be thought of as representing, in polar coordinates, a transversely-incident field on the central cylinder, with the $b_m^{\rm E,H}$ representing the resulting scattered field.

The coefficients $a_m^{\rm E,H}$ and $b_m^{\rm E,H}$ are related in two ways. First, they are connected via the boundary conditions at $r = a$ on the central cylinder. Second, because all the cylinders outside the central unit cell are sources of scattering, the non-singular field in the vicinity of the central cylinder can be written as a sum of singular terms arising from the 2D lattice. This relation takes the form (see Supplement 1):

Such an identity is known as a Rayleigh identity [29,30]. The terms $\tilde {S}_m$ are lattice sums, and can be written in dimensionless form as

To obtain a first approximation to the fields that incorporates the effect of the lattice, we assume, as in the isolated cylinder case, that the dominant form of transmission through the MS layer occurs via a mode that has the same symmetry as the fundamental ${\rm HE}_{11}$ mode of the isolated cylinder, and can be regarded as a perturbation from it resulting from the presence of the lattice. The validity of this assumption will be numerically verified in Section 4; this approximation is expected to hold exactly in the limit that the filling fraction of the cylinders becomes small. For such a mode we neglect all multipole orders $m\neq 1$, and write Eq. (4) as

Equation (7) contains a dipolar correction to the isolated cylinder equation (1) in the form of the term $\tilde {S}_0 I_1$, which contains the influence of the lattice and forms a central part of the discussion presented in this work. In principle it would be possible to incorporate higher-order multipole corrections which would give better accuracy for more closely-packed cylinders, however the inclusion of these additional terms (containing lattice sums $\tilde {S}_2, \tilde {S}_4,$ etc.) requires solving increasingly complicated transcendental equations, and begins to approach a fully numerical treatment in terms of complexity. We therefore confine ourselves to the dipole approximation, containing only $\tilde {S}_0$, from this point forward. The equation (7) contains real quantities and can be solved numerically to find the effective index $n_{\rm eff}$ of the layer. The evaluation of the lattice sum $\tilde {S}_0$ requires a sum over the lattice for each value of $W$ via equation (5); provided $W$ remains real, this sum is rapidly convergent. Because the equation is of the same form as (1) the correction can be readily incorporated into existing algorithms. It can be seen that $\tilde {S}_0\rightarrow 0$ as $\Lambda /a$ becomes large, reducing the equation to that of the isolated cylinder case.

2.3 Phase and transmission

Once $n_{\rm eff}$ is found, the accumulated phase shift of the mode after propagating through the periodic layer is determined by the equation $\Phi = n_{\rm eff}k_0h$ (Fig. 2(a)) This equation shows that, in addition to the element’s height, the geometric dimensions of the nanocylinder, as well as the properties of the lattice, can be used to adjust the phase of the emitted light through modulating $n_{\rm eff}$. This, for instance, was used in the work of H. Ren et al. [4] to generate an MS-based lens that is achromatic across the entire telecommunications range.

The effective index $n_{\rm eff}$ also gives way of estimating the transmission through the MS: by approximating the finite-height nanocylinder with a Fabry-Perot thin film with a central medium having a RI equal to the effective mode index ($n_{\rm central}=n_{\rm eff}$) (Fig. 2(b)), we find that the complex amplitude transmission through the grating is given by

3. Comparison of models

We now compare the accumulated phase $\Phi$ through a silicon MS using the two approximate models, and compare these with a full 3D finite element (FEM) simulation (details in Supplement 1). We consider square lattices of different pitch at an operating wavelength $\lambda =1000$ nm, with fixed cylinder heights ($h=500$ nm). While direct comparison of the speed and computational requirements between numerical and analytical models should be used with caution, for this configuration the FEM model requires approximately 67000 degrees of freedom and 11GB of RAM and solves a single frequency point in about 5s on a workstation (2 Xeon 8 core E5-2665 processors, 3GHz clock speed). The analytic models are of course practically instantaneous, taking less than 50ms to find each point on a laptop-based MATLAB installation. Figures 3(a)-(c) show the accumulated phase $\Phi$ transmitted through the MS by applying the isolated cylinder model (black dashes, from solving 1) compared to the lattice dipole model (blue line, from solving 7), for three different lattice pitches $\Lambda$. Note that the phase from the isolated cylinder model does not change with the lattice pitch. It can be seen that the lattice dipole model gives a far more accurate value for the phase over a greater range of cylinder radii than the isolated cylinder model. In particular, for cylinders with small radius, the lattice dipole model and 3D simulation overlap completely. The agreement continues for larger cylinder fill fractions, until the appearance of sharply resonant features in the 3D results. These features we associate with lattice resonances of the slab as a whole, incorporating scattering from the array in combination with reflections from the cylinder end-facets along the $z$-direction; these resonances cannot be accounted for in either of the approximate models. More details about the occurrence of these resonances will be shown in the following section. It is important to note here that for the correct operation of a MS, only phase shifts within the domain $0\,<\,\Phi \,<\,2\pi$ are required, which is completely covered by the lattice dipole model, whereas the isolated cylinder model fails to accurately predict the phase.

 

Fig. 3. Comparison of models for phase difference and transmission across a square metasurface array for different lattice pitches (red dots: 3D simulations, solid blue lines: lattice dipole model, dashed-dotted black lines: isolated cylinder). (a)-(c) show the phase difference $\Phi$ as a function of cylinder radius. (d)-(f) show the amplitude and phase of the transmission $T$, with the radius plotted on the vertical axis. The projection indicates the phase interval that can be covered by changing the radius.

Download Full Size | PDF

We also show the transmission (i.e., the real and imaginary parts of the complex amplitude transmission) as functions of cylinder radius $a$ through the MS in Fig. 3(d)-(f), comparing the full 3D simulations with the lattice dipole model combined with the Fabry-Perot equations stated above (Eq. 8). The transmission is predicted by the lattice dipole model to within 10% accuracy up until the appearance of the lattice resonances, and remains close to unity for the entire $0\,<\,\Phi \,<\,2\pi$ range.

Within the small radius regime ($a/\Lambda \leq {\rm 0}.3$) the accumulated phase depends linearly on the height of the silicon cylinder $h$, as confirmed by full numerical simulations (Fig. 4(a)-(c)). We see that the slope, i.e., the constant of proportionality is given by the effective index $n_{\rm eff}$, which is accurately predicted by the lattice dipole approach in this regime. By contrast, the accuracy of the isolated waveguide method varies with the radius of the cylinders and does not match the predictions of full 3D simulations. It can also be seen in Fig. 4(d)-(f) that the magnitude of the transmission is well predicted by both approximate models, with the lattice dipole model giving a better estimate of the magnitude of the Fabry-Perot fringes. Note that the presence of oscillations visible in Fig. 4(a)-(c) (particularly for $a/\Lambda$=0.4) overall justifies the use of the Fabry-Perot model, i.e., the model confirms that both boundary surfaces of the nanocylinder must be considered to fully describe the transmission behavior. The approximation of the entire MS structure by a uniform slab neglects effects such as fringing fields that may become important at longer wavelengths and as the slab thickness becomes small. However, the close agreement between the lattice dipole model and the FEM simulations (e.g., positions of minima and maxima and the slopes in Fig. 4(c)) indicates that this effect is practically negligible. This argument is supported by the error analysis shown in Supplement 1 (Sec. S.5, Fig. S.6), also showing a good match between model and FEM results.

 

Fig. 4. Phase difference for different nanopillar heights $h$, for a square array of pitch $\Lambda = 300$nm, for three different cylinder radii in the low radius regime (red dots: 3D simulations, solid blue lines: lattice dipole model, dashed-dotted black lines: isolated cylinder). The top graph in each case shows the computed power transmission, while the bottom graph shows the linear relationship with the layer thickness. The lattice dipole model gives an accurate prediction of the phase in this region. The reference phase $\Delta \Phi$ is the phase at a height of 200nm, which is the minimum used for the Finite Element simulation.

Download Full Size | PDF

4. Modal investigation and lattice resonances

The two assumptions underlying the lattice dipole model are that the fields within each unit cell of the metasurface layer remain dipolar, and that the fields of the finite MS layer can be approximated by a single mode of the infinite structure. It can be seen in Fig. 3 that one of these assumptions breaks down for large cylinder radii, as resonant feature appear. We therefore employ 3D Finite Element modeling to investigate the modal behavior of the fields within the MS layer, with the aim of understanding the transmission spectrum as a whole, as well as obtaining a quantitative description of where the approximate models of the previous section can be applied. We approach this as follows: first, 3D simulations are performed to obtain the fields within the MS layer. Next, we obtain the modes of the layer by Maxwell’s equations for an array of infinitely-extended cylinders. We then compute the numerical overlap between the 3D field and the 2D modes to give us the magnitude of each mode field in the 3D simulation (for details see Supplement 1).

The results of this numerical modal decomposition are shown in Fig. 5 for a square array of lattice pitch $\Lambda = 400$nm and layer thickness of $h=0.5\mu$m. A sharp resonance occurs at cylinder radius $a = 0.154\mu$m, with a smaller resonant feature at $a = 0.12\mu$m. In Fig. 5(b) we show the modal decomposition for six configurations, i.e., small radius (A), moderate radius (B), as well as radii approaching and in close proximity to a resonance (C-F). We see that for smaller radii (A-B) a single mode (labeled as $\pm 1$) appears, which entirely propagates in the forward direction (A, labeled by $+$) or includes a small fraction of the same mode in the backward direction (B, labeled by $-$). This mode is the dipolar-like mode of the lattice, and it has an effective index close to that of the fundamental ($\textrm{HE}_{11}$) mode of the isolated waveguide. For larger radii (C) a second mode (labeled as $\pm 2$) appears, which grows in relative strength towards the main resonance at $a=0.154\mu m$ (D). Close to that resonant peak, both modes are strongly excited in both forward and backward propagating direction (E-F). The existence of the second mode indicates the presence of a resonance of the lattice itself, in which the fields are strongly enhanced in the MS layer. Such lattice resonances have been studied previously [23] and are a known feature of MS arrays.

 

Fig. 5. (a) Power transmission computed using the Finite Element method for different cylinder radii, for a square array of pitch $\Lambda = 0.4\mu$m and height $h = 0.5 \mu$m. The lattice resonance is visible at $a = 0.154 \mu$m. (b) Numerically-extracted mode amplitudes at the points labelled in (a). The modes are ordered on the horizontal axis by increasing effective index; with positive modes being forward propagating. The vertical axis gives the power in each mode relative to the power of the incident plane wave. A second mode appears at point C, and becomes equally dominant with the first mode at the lattice resonance (E).

Download Full Size | PDF

The second mode that gives rise to the lattice resonances has an effective index below that of the light line of the matrix material (RI: $n_1$) - unlike a waveguide mode, this new mode has fields that are propagating, rather than being evanescent, in the $x$-$y$ plane (defined in Fig. 1). This mode arises from the second branch of the photonic band structure of the infinite 2D array, and, like the fundamental waveguide mode, is also strongly dipolar [31]. These transversely propagating lattice modes have been noted previously as being key to maximizing absorption in nanostructured solar cell arrays [32]. We show the dependence of the effective index of these modes in nanopillar radius $a$ for different pitches in Fig. 6(a). Note that the dipolar mode (red) does not exhibit a cut-off and approaches the index of the surrounding medium $n_2$ for smaller radii, while the second mode (blue) has a clear cut-off at which $n_{\rm eff}=0$. This is a clear indication of a mode that exists only because of the lattice and does not occur for the isolated waveguide. This effect is very prominent in hollow-core waveguide optics (e.g., [33]) and is responsible for the resonances in the transmission spectra of these waveguides. Note that in contrast to isolated systems, propagating modes in a periodic environment take the form of Floquet-Bloch modes which, apart from a phase contribution, have identical field patterns in each unit cell. These modes exist in regions below the RI of the cladding material. Thus, there is a clear connection between the behavior of nanopillar-based MS and the modes in waveguides with transversely periodic structures. It is evident from Fig. 6(a) that there is a nanopillar radius at which this mode cuts off that depends inversely on the lattice pitch–at this cutoff the effective index (i.e., the propagation constant) is equal to zero. We can obtain a quantitative prediction of where the cutoff of the propagating dipole mode occurs by adapting the lattice dipole method, with the wavenumber $W$ replaced by

Care must be taken while evaluating this sum since it is only conditionally converges; however, it can be converted to a rapidly convergent sum in reciprocal space [27] (see Eq. S.10 in Supplement 1). Figure 6(b) shows the behavior of the cut-off radius $a_c$ of the second mode as a function of the pitch $\Lambda$ for different wavelength $\lambda$, computed by solving (11). Note that the curves are bound on the small pitch side (grey area in Fig. 6(b)) by the restriction that the cylinders should not overlap. We observe that, for each operating wavelength, $a_c$ decreases as the pitch increases: this implies that lattices with smaller pitch will possess a larger range of cylinder radii $a\,<\,a_c$ where only one mode exists — it is within this region that sharp resonances can be avoided. We also observe that the cut-off radii calculated by (11) agree well with the predictions of the lattice dipole model (7) and with the full numerical simulations. The accuracy of the lattice dipole model and the isolated cylinder model can be determined by comparing the phase predicted by these models with the numerical results of the FEM (see Supplement 1 Sec. S.5, Fig. S.6), assuming that the latter has converged to the "true" solution.

 

Fig. 6. (a) Computed refractive indices of for pitches $\Lambda =0.35\mu$m, $\Lambda = 0.4 \mu$m and $\Lambda = 0.5\mu$m. The radius at which the second mode appears (the cutoff radius $a_c$) is labelled in each case. (b) Cutoff radius $a_c$ as a function of pitch $\Lambda$, computed using equation (11) for different incident wavelengths.

Download Full Size | PDF

5. Discussion of the model and future directions

The model presented in this work allows explanation of the transmission properties of MSs in terms of modes of periodically-arranged waveguides. This represents an extension of the isolated waveguide approach used so far (see e.g., [34]) and allows detailed understanding of mode formation in periodically arranged cylinders, resulting in modes that are different from those of isolated single cylinders. This sheds light on the operating principle of MSs from a mostly overlooked perspective, while also providing fundamental principles for the design of MSs, which involves making various trade-offs.

From the application perspective, it is often crucial for a MS to cover the entire phase space from 0 to 2$\pi$ through structural modification that are straightforward to implement, while maintaining transmission close to unity. In this context, the discussed model offers a direct and straightforward approach to explore parameter spaces over a wide range.

Another crucial consideration is ensuring that the operating range of the MS remains far from resonances, thereby guaranteeing correct device operation even in the presence of manufacturing inaccuracies. The present study reveals that lattice resonances occur when the cylinder radius falls within the two-mode region. To avoid these resonances, which are absent in isolated systems, the radius should be selected well below the cutoff radius specified by the cut-off condition (11). In the case of the square lattice examined here (with an operating wavelength of 1 µm), this implies that the radius must be smaller than 280 nm. Choosing a small cylinder radius, with designed effective index given by (7), might appear to be the best option. However, one should also balance the need to achieve the required phase shift, as a low index will require unfeasibly tall cylinders.

The engineering of modal dispersion is a vital aspect in numerous photonic applications. In this regard, the model presented here enables the direct calculation of group delay and group velocity dispersion of the dipolar mode within the lattice geometry, without the need to consider the height of the nanopillars. Such capability becomes particularly significant for applications involving spectral dependence, including broadband achromatic metalenses, dispersion compensation devices, and manipulation of broadband ultrashort optical pulses.

It is important to note that the waveguide approach used in this work is not necessarily limited to on-axis, cylindrical systems and can potentially be extended to off-axis incidence and more complex cross-sections. A non-normally-incident beam can, in principle, be specified by including a phase term (in the form of a non-zero Bloch vector) in the lattice sum expression of (5). This may lead to enhanced coupling to the resonances of the slab. In addition, square and rectangular cross-sections, in addition to coating the nanopillars with dielectric films, can be effectively treated using waveguide models of the type presented here [35]. In cases where multiple waveguides are present within a single unit cell, the well-established coupled mode theory, widely used in the field of fiber optics, can be used to compute the waveguide modes [17].

Furthermore, our model is not restricted to a specific material system, thus enabling investigations into MSs with lower refractive index contrast [4], or even, depending on the geometry, those composed of entirely different materials such as metals [36]. For metallic inclusions, we expect that the model can be applied for long cylinders that are capable of supporting long-range Surface Plasmon Polariton (SPP) modes, since these are very similar to the dipolar waveguide modes discussed here [37]. Our model would not, however, be applicable to geometries (such as short cylinders or plates) supporting localised resonances such as localized surface plasmon resonances (LSPR), since these structures consist of subwavelength resonators for which there is no propagating mode.

An interesting question arises when considering the inversion of refractive indices: In the field of fiber optics, such a structure can form waveguides which exhibit only one optical mode for all wavelengths: the endlessly single mode fibers [38]. This effect primarily relies on the control of the fundamental space-filling mode, which, in the context of MSs, has not yet been thoroughly explored. The presented model serves as a valuable tool to enhance our understanding of the formation of this type of mode.

6. Conclusion

We have introduced a mathematical model that elucidates the operational principles of metasurfaces in terms of periodically-arranged dielectric waveguides modes. Unlike the isolated cylinder model, our model encompasses the effects of the lattice as a whole, thereby providing a more accurate representation for various practical scenarios. The model relies on lattice sums, which are computable straightforwardly and solely depend on the array’s geometry.

Our investigation focuses on the transmission properties of a silicon nanopillar-based metasurface. By understanding mode formation in conjunction with a Fabry-Perot approach, we are able to accurately compute both the radiated phase and transmitted power. A key aspect of our findings is a semi-analytic expression for the dominant dipolar mode existing within the nanopillars, which greatly influences the behavior of the metasurface within the relevant parameter range. We demonstrate that the dipolar approach remains valid below a certain threshold radius, allowing to cover a full 2$\pi$ phase space. Beyond this threshold, an additional mode, with a substantial portion of power confined within the inter-cylinder region, emerges. This mode induces distinct resonances in the transmission spectrum and corresponds to a higher-order photonic crystal mode of the lattice.

Overall, this model provides a waveguide perspective for understanding metasurface properties, establishing a direct connection to the field of fiber optics. Furthermore, it represents a novel and practical tool that can be readily employed for the design of future metasurface devices with complex functionalities.