Do bilayer metasurfaces behave as a stack of decoupled single-layer metasurfaces?

Flat optics or metasurfaces have opened new frontiers in wavefront shaping and its applications. Polarization optics is one prominent area which has greatly benefited from the shape-birefringence of metasurfaces. However, flat optics comprising a single layer of meta-atoms can only perform a subset of polarization transformations, constrained by a symmetric Jones matrix. This limitation can be tackled using metasurfaces composed of bilayer meta-atoms but exhausting all possible combinations of geometries to build a bilayer metasurface library is a very daunting task. Consequently, bilayer metasurfaces have been widely treated as a cascade (product) of two decoupled single-layer metasurfaces. Here, we test the validity of this assumption by considering a metasurface made of TiO2 on fused silica substrate at a design wavelength of 532 nm. We explore regions in the design space where the coupling between the top and bottom layers can be neglected, i.e., producing a far-field response which approximates that of two decoupled single-layer metasurfaces. We complement this picture with the near-field analysis to explore the underlying physics in regions where both layers are strongly coupled. Our analysis is general and it allows the designer to efficiently build a multi-layer metasurface, either in transmission or reflection, by only running one full-wave simulation for a single-layer metasurface.


Introduction
Metasurfaces or flat optics refer to subwavelength-spaced arrays of scatterers with spatially varying geometries (shape, size, and orientation) and have been widely utilized as a compact wavefront shaping tool [1][2][3][4][5][6][7][8][9].Metasurfaces made of shape-birefringent nanofins have unlocked many possibilities in polarization optics ranging from vectorial structured light and holography to imaging and polarimetry [10][11][12][13][14][15][16][17].While many wavefront shaping capabilities have been demonstrated using plasmonic metasurfaces [18][19][20][21][22][23][24], dielectric metasurfaces have gained much attention recently, especially for phase and polarization control, due to their relatively low loss.Dielectric form birefringent meta-atoms can be locally represented by a spatial arrangement of linearly birefringent wave-plate-like elements, each mathematically described in Cartesian coordinates by the Jones matrix formalism [15].A general Jones matrix has four components, each with two degrees of freedom, amplitude and phase.Hence, there is a maximum of eight degrees-of-freedom, or free parameters, that can be tuned in any linear material.The more free parameters that can be varied in a Jones matrix, the more diverse the capabilities of the structure to manipulate light's polarization [10,11,14,15].
The Jones matrix that can be physically realized with a single-layer dielectric metasurface is subject to being unitary and symmetric [15,25].This sets a fundamental constraint on the retardance and di-attenuation values (i.e., possible polarization transformations) which can be imparted on incoming light.Retardance and di-attenuation respectively refer to the relative phase and amplitude imparted on two input orthogonal polarizations.Notably, matrix symmetry is a fundamental constraint that cannot be surmounted by design.It originates from the linear shapebirefringence of dielectric metasurfaces which fails to realize circular or elliptical birefringence.Hence, a single-layer metasurface cannot be used to build a circular polarizer/retarder -its eigen polarizations must be linear.This limitation exists in any single-layer meta-atom with vertical sidewalls (regardless of its geometry) as long as it is reciprocal [26].Surmounting this constraint requires breaking the in-plane symmetry; either using slanted or a bilayer stack of meta-atoms.Unitarity, on the other hand, is a less fundamental constraint that stems from the lossless nature of dielectric meta-atoms and the typical choice to operate off-resonance to realize higher efficiencies.Nevertheless, a unitary metasurface can still modulate both the amplitude and phase of incoming light by dumping light onto the diffraction orders (which behave as loss channels) as is standard in holography [27,28].
By combining the propagation phase (which arises from varying the dimensions of the nanofins/meta-atoms) and geometrical phase (related to the rotation angle of the nanofin about its axis), it is possible to implement a symmetric Jones matrix with 3 independent DOF-namely, two different phase terms on the diagonal and two identical off-diagonal phase terms.With this functionality, a single metasurface can impart two independent phase [11,29] or amplitude [30] profiles on any two input orthogonal polarizations (linear or elliptical), with the caveat that input elliptical polarizations will flip their handedness at the output.Using clever arrangements of metasurface unit cells, a Jones matrix with more DOFs can be realized.For instance, super-cell based single-layer metasurfaces can offer six DOFs for the Jones matrix, enabling complex amplitude modulation on two orthogonal polarization bases in the far field [27,28].Another variation of super-cell metasurfaces have recently been used to generate multiple polarizationsensitive holograms (exceeding 10 chanels) by exploiting higher diffraction orders as energy loss channels [31].In all these works, however, one cannot freely decouple the input and output polarization states, hindered by matrix symmetry.To tackle this limitation, and access all 8 DOFs of the Jones matrix, the in-plane symmetry must again be broken by constructing a multi-layer system [32].In principle, a bilayer metasurface can impart arbitrary and independent amplitude and phase control on any set of two orthogonal polarizations, while completely decoupling the input and output polarization states [33].Multi-layer metasurfaces also suggest new ways for achieving lossless polarization transformations [34] as well as non intuitive ways for wavefront shaping by utilizing more general configurations of the Pancharatnam-Berry phase [35].
With recent advances in nanofabrication, multi-layer metasurfaces have now become feasible.For instance, by cascading two single-layer metasurfaces made of silicon at design wavelength of 808 nm, each with six DOFs in their Jones matrix, a spatially varying Jones matrix with full parameters of eight DOFs has been realized [36].The latter has been utilized in polarization-selective holography with 16 different channels.In addition, compound meta-optics made of -silicon in the near infrared (in combination with inverse design) have been utilized in spatial mode multiplexing, optical mode conversion, and vectorial holography with very high diffraction efficiencies [37,38].On the other hand, the design of compound metasurfaces (made of either bilayer or cascaded meta-atoms) is computationally intensive compared to single-layer devices.For example, building an exhaustive bilayer library is not a straight-forward process due to the huge size of the parameter space which requires varying the nanofins dimensions (  and   , assuming a rectangular geometry) for both the top and bottom layers while performing the simulation at different angular rotations between the two.Furthermore, assuming a library is built, mapping a target 2-by-2 Jones matrix profile, pixel by pixel, to this massive library is a challenging and time consuming task.Accordingly, bilayer dielectric metasurfaces have been widely treated as a stack of two decoupled metasurfaces.A rigorous validation of this general assumption, however, has not been presented in the literature.To address this gap, here we use full-wave simulations to build a bilayer metasurface library and we study the coupling between the top and bottom nanofins by varying the full design space and analyzing the far-field and near-field responses.Our aim is to provide a systematic recipe that allows the designer to narrow down regions in the parameter space where the top and bottom nanofins are effectively decoupled while avoiding the geometries that suffer from strong bilayer coupling.This viewpoint simplifies the construction of a bilayer metasurface as it only requires one full-wave simulation for a single-layer nanofin in transmission.Using simple Jones calculus, the output response of the former can then be cascaded (for e.g., using matrix multiplication) to build a multi-layer structure.We perform our analysis for bilayer metasurfaces operating in transmission and reflection.Although we studied TiO 2 at the visible wavelength range because of the limited literature on bilayer metasurfaces in that regime, our analysis provides a holistic guideline to the metasurface community which can also be applied to various wavelengths and material platforms as will be shown.

Operating in transmission
We start by creating a library for a single-layer metasurfaces using the finite-difference timedomain (FDTD) and represent its response in terms of a Jones Matrix [11].A model of the simulated structure is shown in Figure 1(a).It is made of a titanium dioxide (TiO 2 ) rectangular nanofin on top of a fused silica substrate and can impart two different phase delays along its major and minor axes; hence the shape-birefringence.The boundary condition applied at the edges of the simulation box is the Periodic Boundary Condition (PBC) which emulates an infinitely periodic array of the same rectangular nanofins.By sweeping the nanofin dimensions and recording the phase and transmission response in the far-field, a single-layer metasurface library can be built (as depicted in Supplementary Figure 1).The response of each nanofin to -and -polarized light can then be mathematically cast in a 2-by-2 Jones matrix.
As a next step, we extend our analysis to the case of a bilayer metasurface.We explore the possibility of computing the Jones matrix of the bilayer as the product of two decoupled (single-layer) Jones matrices.This requires that the coupling between the two nanofins is negligible regardless of the dimensions of the two nanofins and their relative rotation.To test this assumption, we compare the results of the FDTD simulations of the bilyer with the "analytical" results obtained from the product of two single-layer metasurfaces.The model of the simulated structure is shown in Fig. 1(b).Initially, we fixed the dimension of the bottom nanofin at 134 nm x 202 nm, an arbitrarily chosen geometry which emulates a quarter-wave plate.The Jones matrix of this geometry was then extracted from the single-layer library reported in Supplementary Figure 1.Afterwards, we performed a parameter sweep for the dimensions of the top nanofin without introducing a rotation angle between the nanofins; i.e., keeping the angular orientation of both nanofins at 0 • .This simulation helps us verify if the two layers can be considered decoupled for all the geometries.If that is true, then the Jones matrix of the bilayer could be expressed as: The transmission and phase responses for this bilayer structure under -polarized illumination are depicted in Figs.1(c-d), respectively.From these plots, one can observe transmission dips which can be attributed to resonances.Here, we refer to all geometries whose transmission is lower than 80 % to be in the regime of resonance.In Fig. 1(e), we plot the error between the phase response (from the full-wave simulation above) and the analytical phase response calculated using Eq.(1).Since the rotation angle between the two nanofins is zero ( = 0) the Jones matrix that describes the structure is diagonal.Hence, in the plots we only show the results related to the element  11 .The error plots for element  22 (under incident y-polarization) are similar due to the rectangular geometry and are thus not included.The average absolute phase error here is less than 3 • .For most of the geometries, the simulation results coincide with the analytical ones.The few geometries that exhibits a larger error (only 2% of the geometries have an error higher than 10 • ) are the same ones that correspond to the resonance lines in the top right of Fig. 1(d).Therefore, for these geometries and others with large phase error, a full-wave simulation of the bilayer structure is needed to accurately capture its Jones matrix.However, given that full 0-2 phase coverage can be achieved with enough geometries away from resonance, these resonant elements can simply be filtered out from the library.In the next section, we will provide a more in depth investigation of these resonances (and their types) and demonstrate cases in which the coupling effects in a bilayer meta-atom due to resonance can be neglected.
To study the case of a non-diagonal Jones matrix, we introduce a relative rotation between the top and bottom nanofins.In Supplementary Note 2, we considered two cases of bilayer structures made of identical top and bottom nanofins and we rotated the top nanofin at increments of 15 • to study the effect of angular orientation on coupling.The two bilayer structures are made of identical nanofin with dimensions 134 nm × 202 nm and 114 nm × 154 nm, respectively.For each of the two structures, we tabulated the amplitude and phase error for each element in their corresponding Jones matrices.Our analysis suggests an average phase error -between FDTD simulations and analytical calculation of Eq. (1) -on the order of 5 • .The phase error plots are shown in Supplementary Figure 2 for one of the two geometries confirming that 0 to 2 phase coverage can be achieved with acceptable errors.
As a next step, we extend our analysis by performing a full sweep for the dimensions of both the top and bottom nanofins.We then examine the effect of bilayer coupling by comparing these FDTD simulation results with the analytical calculation of Eq. ( 1), as before.Figures 2(b-c) depict the phase shift and power transmission obtained from FDTD.Each figure exhibits a grid that is made of 121 subplots (cells).Each cell represents one parameter sweep where the top nanofin has the dimensions reported on the horizontal/vertical axes while the dimensions of the bottom nanofin are swept from 50 nm to 250 nm.Interestingly, the response of each subplot mimics the behaviour of the entire plot, as if it is of fractal nature.This behavior can be reconciled with sampling theorem since the figure is compiled from a Jones matrix product, between the top and bottom nanofin, that is reminiscent of a convolution between their respective phase plots.To quantify the bilayer coupling effects, we plot the phase and transmission errors (with respect to the analytical prediction of Eq. ( 1)) as shown in Figs.2(d-e), respectively.The axes of the two plots can be interpreted the same as described above.This analysis confirms the possibility of expressing the Jones matrix of a bilayer dielectric metasurface as the product of the single-layer Jones matrices over most of the geometries included in the considered parameter space.The average phase error in this case is 3.19 • .Therefore, for a large number of geometries, the hypothesis that the two layers can be treated as decoupled is a valid one.Notably, the phase redundancy afforded by the metasurface library offers a sufficiently large number of decoupled geometries with full 0 to 2 phase coverage.As a rule of thumb, the geometries that exhibit larger errors fall within two main categories: a) structures in which one of the two nanofins is operating near resonance, and b) structures with significant reflections where the top nanofin is much larger than the bottom one.The latter can specially be inferred by examining the top right part of Fig. 2(c).This can be attributed to the multiple (Fabry-Perot like) reflections that occur between the substrate and the lower base of the top nanofin.To gain more insights into the bilayercoupling, we consider different geometries and perform a near-field analysis to examine the cases with large errors.The underlying physics of this problem is discussed next.

Near-field analysis
In this section, we complement the far-field Jones matrix calculations presented before with a near-field analysis.Our goal is: a) to evaluate the accuracy of our Jones matrix analysis in predicting the response of the bilayer structure, and b) to gain more insights into cases where the far-field response obtained from FDTD deviates from the Jones matrix analysis.Towards these aims, we adopted the scattering matrix (S-matrix) and transmission matrix (T-matrix) approaches while calculating the fields using the eigenmode expansion (EME) method [39][40][41][42].Since we are interested in understanding the bilayer coupling in a transmissive metasurface, our ports are defined at the input of the bottom nanofin and at the output of the top nanofin, respectively.
To assess the effects of guided and evanescent modes as well as back reflections, separately, we adopt four different calculations: a) The full S-martrix which captures the contributions of guided, evanescent, and back reflected fields (hence, it is the most accurate calculation or ground truth).b) The full T-matrix is similar to (a) but without recording the back reflections; only forward propagating guided and evanescent fields.(c) The (2 × 2) S-matrix which only Comparing between these four calculations will help quantify the effects of bilayer coupling (which can be inferred from the strength of evanescent fields) and back reflections for various meta-atom geometries.To further examine the effects of evanescent coupling and impedance mismatch, we introduce a small air gap between the top and bottom nanofins.By varying the gap size between the two nanofins and recording the amplitude variation one can identify the regimes where bilayer meta-atoms no longer behave as a stack of two decoupled single-layer metasurfaces.
We consider four main categories of bilayer dielectric meta-atom geometries under x-and y-polarizations, respectively.Figure 3(a) depicts the first case: a bilayer meta-atom composed of two off-resonance nanofins * with a larger nanofin at the bottom.By looking at the amplitude response under the input polarizations,   (11) and   (00), one can observe that the S-matrices transmission oscillate only slightly around a mean value that matches the T-matrices result.Here, the deviation in amplitude between the S-and T-matrices calculations is on the order of 0.01%.This implies that the (2 × 2) T-matrix (and so the Jones matrix) correctly describe the nanofin dynamics accurately; the evanescent coupling and back reflections in this type of geometry are small enough to be ignored.Hence, under x-polarized illumination, light will be reflected due to the size mismatch between the two nanofins.In this case, a Fabry-Perot like cavity will be created between the substrate and the top nanofin.As these back reflections are not captured by the T-matrices, a large deviation in the S-and T-matrix amplitude response is observed.On the other hand, when the same nanofin is illuminated by y-polarized light (blue curves), the amplitude response evaluated by the T-matrices and S-matrices are in good agreement; the back reflections and evanescent coupling in this geometry is minimal.These observations can be reconciled with waveguide theory.Both geometries involve impedance mismatch between the top and bottom nanofins (i.e., the modes are characterized by two different effective indices).However, since the small nanofin will have a smaller effective index (close to the cladding material-air), its placement above the large nanofin is already captured by the Jones matrix of the large nanofin.This is not the case when the small nanofin is terminated by the large one on top.The impedance mismatch in the latter is more significance and is thus not fully captured by the Jones matrix of the small nanofin.
Figure 3(c) shows the case of two identical nanofins operating at resonance.In this case it is not expected that back reflection between the top nanofin and the substrate can occur (given the agreement in cross sectional area) and indeed the discrepancy between the full S-matrix and full T-matrix is insignificant.However, the large discrepancy between the full S-and full T-matrices versus the (2 × 2) S-and (2 × 2)T-matrices suggest that the higher order evanescent modes play a fundamental role in the amplitude response.Here, the evanescent field coupling between the two nanofins is very significant (due to the operation at resonance) and, hence, cannot be neglected.Additional categories of nanofins that feature back reflections and that exhibit coupling through evanescent fields are shown in Supplementary Figure 3.The former is typically observed when the top nanofin is larger in size whereas the latter occurs when at least one of the two nanofins is at resonance.
Figure 3(d) highlights another case of resonant geometries.Here, only the bottom nanofin is at resonance whereas the top nanofin is of much smaller dimensions.In this case, due to the small spatial overlap between the two nanofins, the evanescent coupling is not significant.The smaller dimensions of the top nanofin also suggest that the back reflections are minimal.These two behaviors-i.e., weak evanescent coupling and reflections-are indeed suggested by the response of the (2 × 2) T-matrix which captures the amplitude transmission fairly accurately.Hence, the full S-matrix and full T-matrix are in close agreement.
In short, our near-field analysis confirms that (under some circumstances) the S-and T-matrices responses coincide.Cases with disagreement correspond to geometries that involve large back reflections (small spatial overlap between top and bottom nanofins) or geometries that lie close to resonance.By avoiding these regions in the parameter space, one can populate a set of bilayer meta-atoms which effectively behave as a stack of two decoupled single-layer metasurfaces.
Thanks to the redundancy afforded by the metasurface library such a task is possible to achieve. Figure 4(a) depicts a schematic which visualizes the nanofin categories.The set of allowable (decoupled) nanofin geometries are represented in the first row (green zone).After filtering out all the geometries with bilayer coupling, the remaining meta-atoms (around 50% of the parameter space) still densely span from 0 to 2 phase shift as shown in Fig. 4(b).
Note that thus far we have chosen a specific material, namely titanium dioxide, as a platform for conducting our analysis.Nevertheless, we expect the physical dynamics associated with evanescent coupling and back reflection, to be universal across other dielectric metasurface libraries using different material platforms or design wavelengths.For instance, silicon nitride (SiN) is another widely used material in the visible, albeit with less index contrast [43].Therefore, due to the less mode confinement, the coupling strength can be more significant causing all resonant geometries to be considered coupled (including for e.g., the category in Fig. 4(a)-iii.
Other effects, such as back reflections, should preserve their behavior.In addition, other material platforms also exist for the near infrared and telecom regimes including, for e.g., silicon.To validate the generality of our analysis, we study bilayer dielectric metasurfaces based on silicon, at a design wavelength of 1550 nm, and show that the coupling effects are governed by the same physical dynamics.We summarize the results of this analysis in Supplementary Section 7.

Operating in reflection
Thus far we have shown that we can evaluate the Jones matrix of a transmissive bilayer metasurface starting from the Jones matrix of a single-layer metasurface under some constraints.This allows the designer to build a bilayer metasurface by only utilizing the single-layer metasurface library, thereby simplifying the design process.In this section, we investigate the validity of this assumption for reflective dielectric metasurfaces.We examine if a birefringent metasurface operating in reflection can be expressed as the product of four Jones matrix (each describing a single-layer nanofin).In this case, the product of one unit cell (pixel) on the bilayer metasurface, assuming no rotation, is given by: Equation ( 2) describes the path that light makes when impinging on the metasurface as it passes through the two nanofins, gets reflected by the mirror, before traversing the two nanofins again, in the reverse order.To test the validity of Eq. ( 2), we start by simulating a unit cell consisting of a single-layer in reflection and then we build on it by considering the full bilayer metasurface in reflection.Following the same approach used for a transmissive metasurface, we consider a single-layer reflective metasurface first to verify if, in presence of a mirror, one can express the Jones matrix of the system as the following product: Supplementary Figure 4 shows a comparison between the phase responses of the the reflective (single-layer) metasurface obtained from FDTD and and the one calculated analytically using Eq.
(2), suggesting significant discrepancy between the two.Note that placing the nanofin on top of a mirror perturbs the phase response of the structure, possibly introducing standing wave patterns (similar to terminating a waveguide with a complex load).The mirror-dielectric interface effects cannot be accurately captured by Eq. ( 2).Therefore, in contrast to the case of transmissive bilayer metasurfaces, where the analytical approach was fairly accurate in predicting the amplitude and phase responses, the requirements are more stringent when operating in reflection.To bypass this challenge, we introduce a simple modification to the structure by inserting a layer of silica (named "spacer") between the mirror and the nanofin, essentially matching the impedance between the two, as depicted in Fig. 5(a).To compensate for the effective index perturbation of the TiO 2 fin which arose from the mirror-dielectric interface, the spacer dimensions need to be optimized depending on nanofin size.The goal is obtain an overall response for the reflective structure that matches the behaviour of the transmissive nanofin as if it was in direct contact with the substrate.This allows us to build a reflective bilayer metasurface using the Jones matrix product by starting from the same single-layer metasurface library in transmission.
To optimize the spacer, we sweep its thickness between 20 nm and 165 nm, and for each thickness, we perform full sweep on the dimensions of the nanofin.For each nanofin geometry, we find the optimized spacer dimension that minimizes the phase error compared to the analytical one of Eq.  (4). Figure 5(b) shows the phase shift obtained for each geometry when its optimized spacer thickness is selected.The error between the simulated results and the analytical ones is 0.5 • , on average, provided that the spacer thickness is optimized, as depicted in Fig. 5(c).In addition, Fig. 5(d) shows the optimum spacer thickness as a function of the nanofin dimensions.Here, the vertical axis depicts the phase error between FDTD and the analytical product of Eq. ( 4).Therefore, it is now possible to express a single-layer reflective metasurface as the Jones matrix product of single-layer decoupled nanofins.This intermediate step is essential as it enables us to design a bilayer metasurface by only relying on a single-layer library.However, the insertion of the matching spacer obviously imposes a constraint on the design of the metasurface; as it is not possible (at least with our current fabrication techniques) to build a device with spatially varying spacer thickness.Instead, we envision the final device to be made of a bilayer meta-atom where the bottom nanofin is of fixed dimensions while only sweeping the top nanofin.A metasurface unit cell of this kind can still realize an asymmetric (yet unitary) Jones matrix, point-by-point, across the structure.By making use of super-cell metasurfaces or one can break the unitarity constrainst as well.
To design a bilayer metasurface in reflection we repeat the previous analysis performed in transmission.We fix the dimension of the bottom nanofin at 134 nm × 202 nm.This is justified because a stack of two nanofins provides 6 degrees-of-freedom, whereas an arbitrary Jones matrix requires only 4 [36].Hence, by fixing the bottom nanofin and varying the top one, all 4 degrees-of-freedom can still be accessed.We set the spacer thickness 100 nm which is the optimized value for that selected geometry as suggested by Fig. 5(d).We performed a parameter sweep of the dimensions of the top nanofin without introducing a rotation angle between the fins.The previous analysis on the spacer optimization was needed to select the spacer dimension that optimizes the response of the fixed bottom geometry.This allow us to verify the assumption of decoupling in reflection for all the geometries.We also confirm that the spacer dimension from the single-layer is valid for the design of a bilayer metasurface in reflection; in which case the Jones matrix of the bilayer can be written as shown in Eq. ( 2).The results of the simulation are shown in Fig. 5(f) and Fig. 5(g).For most of the geometries, the spacer selection rule defined above allows to accurately build a reflective bilayer starting from the library of a transmissive single-layer metasurface given the low phase error reported in the color map of Fig. 5(i).The latter represents the phase difference between the FDTD simulations and the analytical product of Eq. ( 2).The average absolute phase error in this case is 5.3 • .The cases that show larger errors are the ones composed of a top nanofin that is much larger than the the bottom one.This can be due to the reflections that occur between the mirror and the base of the top nanofin.This is reminiscent of the observation we made for the case of transmissive bilayer metasurfaces.

Conclusion
We showed that a bilayer dielectric metasurface operating in transmission can be expressed as the product of two decoupled single-layer metasurfaces under some constraints.In this process, we distinguished regions in which the bilayer coupling is governed by resonance versus back reflections and we provided systematic recipes to avoid operation in both regimes.Furthermore, we demonstrated that it is also possible to express a reflective bilayer metasurface as the product of five matrices which describe the nanofins composing the structure, the reflective mirror, and a matching spacer in between.By combining our near and far-field analysis, we narrow down the design space to a smaller subset of geometries that are essentially decoupled.Notably by excluding the meta-atoms with strong coupling from the design library-either by avoiding resonant structures or bilayer geometries with very large top nanofins-one can efficiently build a multi-layer metasurface as a cascade of single-layer meta-atoms.In this case, fitting a target profile to a library will entail decomposing it into a product of two matrices and fitting each one following the same selection criteria of a single-layer nanofin.We validated the applicability of our approach to a wide range of libraries by considering titanium dioxide platform at a design wavelength of 532 nm in addition to silicon at 1550 nm, showing the generality of our approach.

Generating a single-layer metasurface library
In this section, we present a systematic strategy for designing a single-layer metasurface.We start by building a metasurface library based on the results and ideas discussed in [1].To the aim of obtaining a direct relation between the phase shift and the dimension of the nanofins, a The dimension of the meta-atom, namely   and   , range from 50 to 250 nm spanning a total of 2601 geometries which cover the phase range from 0 to 2.The center-to-center separation , which defines the unit cell size, is 420 nm.The structure is illuminated by an x-polarized plane wave source at 532 nm that has been placed at a distance of 600 nm from the bottom of the substrate.At this wavelength the refractive index of the fused silica and of TiO 2 are set to 1.46 and 2.48, respectively.For the reflective metasurfaces, we set the refractive index of the Aluminum mirror to 0.811 + 0.366.A monitor is placed a few wavelengths above the nanofin in order to evaluate the phase response in the far-field and the percentage of transmitted power.The mesh size of this simulation is set to 2.5 nm × 2.5 nm × 2.5 nm.
The obtained results are shown in Fig. 1 showing an excellent agreement with those reported in Ref [1].Note that for the case of a -polarized source, the plots of the phase (  ) and power transmission (  ) would be identical to those shown in Fig. 1 but with the  and  axes exchanged

Bilayer simulations: the effect of relative rotation
In this section, we consider a bilayer metasurface and we analyze the effect of relative rotation between its top and bottom nanofins.Two cases are tested: in the first case we set the dimensions of both the bottom and top nanofins to be 134 nm × 202 nm while in the second case we set the dimensions of both nanofins to be 114 nm × 154 nm.The rotation angle of the bottom nanofin around its geometrical center is set to zero ( = 0) and we rotate the top nanofin ( ′ ) between 0 • and 90 • with increments of 15 • .Other rotation angles for the top nanofin beyond  ′ > 90 are naturally covered by this sweep (following a simple symmetry argument).In Tables 1 -4, we report element-by-element comparison between the Jones matrix of the simulated bilayer (J s ) and the "analytical Jones matrix (J a )" obtained as follows: where ( ′ ) is the 2 × 2 rotation matrix.
Table 1.Element-by-element comparison between the magnitude of the Jones matrix of the simulated bilayer J s and the elements of the analytical Jones matrix J a for the bilayer composed of two nanofins of dimensions 134 nm × 202 nm.
Table 3. Element-by-element comparison of the magnitude of the Jones matrix elements of the simulated bilayer J s and the elements of the analytical Jones matrix J a for the bilayer composed of two nanofins of dimensions 114 nm × 154 nm.
The discrepancy that arises between the Jones matrix of the simulated bilayer and the "analytical Jones matrix" is on average not so large and can be considered acceptable for both the case of phase and amplitude.This implies that, at least for the tested geometries, the output response can be calculated using Eq. ( 5).To reconcile the discrepancies, recall that we apply the Periodic Boundary Conditions (PBC) when simulating both the single and bilayer metasurfaces.The PBC emulates a periodic array composed of identical nanofins with zero rotation angle everywhere producing an electric field at the borders of the unit cell that is different from the case of a rotated nanofin.Although the nanofins are rotated, the square unit cells comprising the simulated structure are always fixed.Hence, the lattice symmetry is different from the case reported in Fig. 1.Accordingly, applying the rotation matrix on the simulation data of the aligned single-layer nanofins-calculating Eq. ( 5)-does not entirely capture their response under rotation and is thus expected to slightly differ from the full-wave simulations.
The data reported in Tables 1-4 suggest that the discrepancy between the analytical and simulated cases is higher for the bilayer geometry of dimensions 134 nm × 202 nm.This is likely due to the larger mismatch between the major and minor axes (compared to the 114 nm × 154 nm geometry) which in turn reduces the overlap region between the top and bottom nanofins, thereby introducing more reflections and making this geometry more sensitive to the rotation-dependent variation of the electric field at the boundary.

Phase coverage of the bilayer metasurface operating in transmission
Here, we report the results obtained from the simulation of the structure shown in Fig. 1(b).The dimension of the bottom nanofin are fixed (134 nm × 202 nm) while the dimensions of the top nanofin are swept between 50 nm and 250 nm. Figure 2 depicts the phase shift of the diagonal elements of the Jones matrix of the simulated structure,  1,1 and  2,2 .These results suggest that both elements (in the analyzed parameter space) can cover a phase range between 0 and 2 as reported on the vertical axis of the plots.The color bars depict the error between the simulated results and the analytical ones for both Jones matrix elements.Similar results were obtained for the case of a bilayer metasurfce operating in reflection with fixed bottom nanofin dimensions.

The effect of resonance and of back reflections
There are two main sources of errors when describing the bilayer as a stack of two decoupled single layers: a) Fabry-Perot-like back reflections that typically occur between the top nanofin and the substrate if the top nanofin is larger than the bottom one, and b) coupling through the evanescent fields when either nanofin operates near resonance.In this section, we study these sources of error in more detail.Figure 3 shows the field profile in the xz-plane of all possible geometries that cover these two sources of error.The results were obtained using full wave (FDTD) simulations.By varying the gap between the top and bottom nanofins and recording the amplitude response, we gain some insights into the bilayer coupling strength.
In Fig. 3(a) we consider a bilayer meta-atom composed of two nanofins at resonance.In this case, as it has been pointed out in the near-field section, the two nanofins are coupled through the evanescent modes (which are more significant when the nanofins are in proximity).Consequently, by varying the gap between the two nanofins, the field profiles and their confinement is altered.
Such behavior is consistent with the S-and T-matrix analysis (of Fig. 3(c)) where a discrepancy between full S-and T-matrices versus the (2 × 2) matrices was observed.
Figure 3(b) depicts another case in which the top nanofin is at resonance and is larger than the bottom one.Interference fringes are observed between the top nanofin and the substrate due to the back reflections.This is consistent with the oscillations in transmission amplitude observed in Fig. 3(b).The mode shape in the top nanofin is unperturbed.Hence, the effect of back reflections is more significant than resonance, causing the two nanofins to be strongly coupled.
To complement this picture, we consider another geometry in which the top nanofin is still at resonance but is smaller than the nanofin bottom; hence, the back reflections are mitigated.In this case, the field profile still does not change as a function of gap size and no interference fringes are observed.Therefore, unlike the previous case, the two nanofins can be considered decoupled even in the presence of resonance.
Resonance becomes much more significant when it occurs in the bottom nanofin; the bilayer coupling becomes strong due to the interaction of the evanescent fields.If the top nanofin is too large (considering our parameter space, if its dimensions roughly exceed 130 nm × 130 nm), these evanescent fields are strongly coupled.This is confirmed by the FDTD simulations shown in Fig. 3(d On the other hand, in the limit when the top nanofin is much smaller (while the bottom is still at resonance), the coupling is weak and the field profile does not change by varying the gap size as shown in Fig. 3(e).
Figure 3(f) emphasizes the effect of back reflections in the absence of resonance (both nanofins are off-resonance).One can observe the Fabry-Perot-like interference fringes that arise between the top nanofin and the substrate.This effect was not captured by the (2 × 2) T-matrix but rather by the S-matrices.This confirms that this geometry cannot be treated as a stack of two decoupled single layers.
Lastly, in Fig. 3(g) we consider a more straight forward case with neither resonances nor back reflections.Here, both nanofins operate off-resonance.The top nanofin is smaller than the bottom one.As it can be noticed, the field profile remains identical regardless of the gap size.No fringes arise and the two nanofins are effectively decoupled.

Single-layer metasurface in reflection
Here, we show the simulation results of the single-layer metasurface operating in reflection without introducing a spacer between the mirror and the nanofin.As mentioned in the main text, looking at the results of the analytical calculations (Fig. 4(a)) and those of the FDTD simulations (Fig. 4(b)), one can immediately notice the significant differences over the full parameter space.
The average phase error between the two cases is on the order of 10 • .

The effect of a relative rotation of 45 • in reflection: sweeping the dimensions of the top nanofin while fixing the dimensions of the bottom one
To verify if the assumption of decoupling is also valid when a relative rotation between the metasurface nanofins is introduced, we simulated a bilayer structure while fixing the bottom nanofin and sweeping the dimension of the top nanofin.Here, the top nanofin is rotated with respect to the bottom one by an angle of 45 • .In this case it is necessary to introduce a rotation matrix () that sandwiches the Jones matrix of the top nanofin to take into account the rotation of the top nanofin.The Jones matrix of this structure can be analytically written as follows: where () is the 2 × 2 rotation matrix that is defined as: Hence, by introducing a relative rotation, it is possible to obtain a metasurface where the four Jones matrix elements are all different.In Supplementary Figures 5 and 6 we show the magnitude and phases of the four elements of the Jones matrix describing the bilayer obtained via FDTD simulations.Interestingly, the geometries that now exhibit larger phase error -for elements  1,1 and  2,2 -are not the ones in the top right of the plot (as was the case with no rotation) but rather the ones in the top left and bottom right.
The average phase error is 1.14 • for the element  1,1 and 1.02 • for the element  2,2 , exhibiting less error compared to the case where no rotation angle was introduced.This is because the rotation increases the overlap region between the top and bottom nanofins, thereby reducing multiple reflections between the top nanofin and the substrate.For the off diagonal Jones matrix elements, the largest phase error is in both cases on the diagonal of the plot.On the other hand, the phase of these element is not relevant since, as it can be noticed from Fig. 5, the magnitude of these elements is zero.The average absolute phase error for the elements  1,2 and  2,1 , excluding the points on the diagonal is around 3 • .

Near-field analysis for a Silicon bilayer metasurface
In this section, we expand our analysis by considering silicon as another material platform at a design wavelength,  = 1550 nm.Our aim is to show that the physical dynamics associated with back reflections and coupling in a bilayer metasurface are, to some extent, independent of the material platform, its parameter space, and wavelength.The phase and transmission library for a single layer silicon nanofin is depicted in Fig. 8. Using this library, we then performed a near-field analysis for a silicon bilayer metasurface at 1550 nm, considering a unit-cell size of 500 × 500 nm and 1-m tall nanofins.
We adopted the same scattering and transmission matrix approaches used for the TiO 2 platform in the main text.We tested the same four categories considered before, assuming x-polarized incident light.Figure 9(a) shows the first case: a bilayer metasurface composed of two off-resonance nanofins with a smaller nanofin at the top.It is observed that the Full S-matrix transmission oscillates slightly (0.01%) around a mean value which matches the (2x2) T-matrix.In this case, the Jones matrix is able to correctly describe the nanofin dynamics as if they were decoupled.This is consistent with the results obtained for TiO 2 .Furthermore, in the same figure, we show the field profile in the xz-plane obtained using FDTD simulations.In this case, since neither the evanescent coupling nor back reflections play a significant role, the field profile remains identical regardless of the gap size.Here, the coupling trough evanescent modes causes a large discrepancy between the Full S-(and T-) matrices, compared to the (2x2) T (and S) matrices.By looking at the field profiles, we observe that the field confinement due to resonance is highly dependent on the gap size.On the other hand, if the top nanofin is much smaller, while the bottom is still at resonance, the evanescent coupling becomes weak and field confinement exhibit less dependence on the gap size.This is depicted in 9(d) where the full S-matrix and the (2x2) T-matrix responses are in very good agreement.In summary, this analysis shows that the rules governing the evanescent coupling and back reflections in bilayer metasurfaces are somewhat universal, regardless of the material platform and design wavelength.The underlying physics dictating if a bilayer meta-atom can be described as two decoupled nanofins thus apply to different metasurface libraries.

Figure 1 .
Figure 1.Bilayer dielectric metasurface with fixed bottom nanofin, operating in transmission.(a) Model of the unit cell of a transmissive single-layer metasurface.(b) Model of the unit cell of the simulated bilayer.The dimensions of the bottom nanofin are fixed and set to be 134 nm × 202 nm.The dimensions of the top nanofin are swept from 50 to 250 nm.The Power transmission and phase response of the bilayer metasurface are shown in (c) and (d), respectively.(e) Absolute error in phase shift between the simulation results shown in (d) and the results analytically obtained from the matrix product of two cascaded single-layer metasurfaces.

Figure 2 .
Figure 2. Bilayer dielectric metasurface with top and bottom nanofins of variable dimensions, operating in transmission.(a) Model of simulated bilayer metasurface: the arrows refer to sweeping the dimensions of top and bottom nanofins.The phase shift and power transmission of the structure in (a) is recorded using FDTD simulation and is shown in (b) and (c), respectively.Each cell in these grids refers to a separate simulation in which the bottom nanofin's dimensions are swept while fixing the dimensions of the top nanofin.The top nanofin's dimensions are depicted on both axes.(d) Error in phase shift between the simulation results shown in (b) and the results analytically obtained from the single-layer library.(e) Error in power transmission between the simulation results shown in (c) and the results analytically obtained from the single-layer library.

Figure 3 .
Figure 3. Scattering and transmission matrix analysis for transmissive bilayer metasurfaces with aligned top and bottom nanofins and variable air gap in between.Four cases are considered: (a) Two nanofins operating off-resonance.The top nanofin is smaller than the bottom one.The transmission amplitude response from the S-and T-matrices under two orthogonal polarizations,   (11, blue curve) and   (00, red curve), is plotted.(b) A birefringent bilayer meta-atom with different dimensions along x and y.The two nanofins are off-resonance and the corresponding amplitude response under x-and y-polarizations is shown.(c) Bilayer meta-atom composed of two identical nanofins operating at resonance.(d) Only the bottom nanofin of the bilayer meta-atom is at resonance and is larger than the top nanofin.

Figure 3 (
Figure 3(b) represents the second case which is anisotropic: the two nanofins have different dimensions along the x and y directions but neither of the two is at resonance.The top nanofin is larger along the x direction (D ,top > D ,bottom ) and is smaller along y (D ,bottom > D ,top ).Hence, under x-polarized illumination, light will be reflected due to the size mismatch between the two nanofins.In this case, a Fabry-Perot like cavity will be created between the substrate and the top nanofin.As these back reflections are not captured by the T-matrices, a large deviation in the S-and T-matrix amplitude response is observed.On the other hand, when the same nanofin is illuminated by y-polarized light (blue curves), the amplitude response evaluated by the T-matrices and S-matrices are in good agreement; the back reflections and evanescent coupling in this geometry is minimal.These observations can be reconciled with waveguide theory.Both geometries involve impedance mismatch between the top and bottom nanofins (i.e., the modes are characterized by two different effective indices).However, since the small nanofin will have a smaller effective index (close to the cladding material-air), its placement above the large nanofin is already captured by the Jones matrix of the large nanofin.This is not the case when the small nanofin is terminated by the large one on top.The impedance mismatch in the latter is more significance and is thus not fully captured by the Jones matrix of the small nanofin.

Figure 4 .
Figure 4. (a) Schematic of different categories of bilayer dielectric meta-atoms.The first row corresponds to the geometries for which the coupling between the nanofins is negligible: (i) neither the bottom nor the top nanofins are at resonance and bottom nanofin is larger than the top, (ii) only the bottom nanofin is at resonance and is larger than the top, (iii) only the top nanofin is at resonance but the bottom nanofin is much larger than the top.The second row depicts the cases for which the two nanofins in the bilayer are strongly coupled: (iv) neither the bottom nor the top nanofins are at resonance and the top nanofin is larger than the bottom, (v) only the top nanofin is at resonance and larger than the bottom, (vi) only the top nanofin is at resonance and slightly smaller than the bottom, (vii) both nanofins are at resonance.These cases are detailed more fully in Supplementary Section 4. (b) Complex transmission of the dashed region in (a).The electric field amplitude       is plotted on the complex plane demonstrating the 0-2 phase coverage afforded by the considered geometries while maintaining low loss.The red circle is the unit circle.The black circle is the average of   over all considered geometries.

Figure 5 .
Figure 5. Operation of a reflective metasurface.(a) Model of the unit cell of a single-layer metasurface operating in reflection.A layer of silica has been inserted between the mirror and the nanofin.(b) Phase shift obtained from the FDTD simulation of the structure shown in (a).For each geometry the spacer thickness has been optimized so that the phase difference between the simulation results and the analytical product is minimized.(c) Difference between the simulated and analytical phase shifts for the structure shown in (a) while optimizing the spacer thickness for each geometry.(d) Optimum spacer dimension as a function of the nanofin geometry.The vertical axis depicts the phase error for each optimum spacer thickness.(e) Model of the unit cell of bilayer metasurface working in reflection.(f) Phase shift obtained from the FDTD simulation of the structure shown in (e).(g) Phase shift error for the structure shown in (e) with optimized spacer thickness for each geometry.(h) Power transmission obtained from FDTD simulation for the structure shown in (e).i Error in power transmission for the structure in (e) when the spacer thickness is optimized for each geometry.
library has been generated by performing a parameter sweep with finite-difference time-domain (FDTD) simulations.In a simplified picture, each subwavelength structure can be considered as a truncated waveguide or a low-quality-factor Fabry-Perot resonator.Nanofins with different dimensions (length and width) will induce different confinement of the field impinging on the structure.This confinement provides an effective refractive index that differs along the two polarization components.The FDTD software used is Ansys Lumerical.It solves Maxwell's equations on a discrete spatial and temporal grid in complex geometries such as the analyzed case.A model of the simulated structure is the one shown in Fig.1(a); it is composed of a TiO 2 rectangular nanofin on a fused silica substrate.The boundary conditions applied at the edges of the simulation box are the Periodic Boundary Conditions (PBC) which emulate the existence of an infinitely periodic array of these rectangular fins so that the simulated structure is an infinite array of the same metasurface unit cell.

(Figure 1 .
Figure 1.Simulation data for two-dimensional parameter sweeps of TiO 2 rectangular fins (ℎ = 600 nm).(a) Phase shift   on x-polarized light.The phase shift has been computed as the ratio between the phase collected in the center of the far field projection of a monitor above the structure when the nanofin is present on top of the substrate and the phase at the same monitor when only the silica substrate Units in radians.(b) Power Transmission   for x-polarized light.The total power passing through a monitor above the structure relative to the source.(c) Complex Transmission.The blue dots represent the electric field plotted on the complex plane for all the simulated geometries.The black circle corresponds to the averaged transmission.The red circle is the unit circle.

Figure 2 .
Figure 2. Phase error between the Jones matrix elements of the full-wave simulated bilayer and the corresponding analytically computed ones.Full 0 to 2 phase coverage can be achieved with low phase error.Here the bottom nanofin dimensions are fixed.(a) Difference between the phase of element  1,1 extracted from the simulation and the phase of  1,1 analytically computed from the data of the single layer.The vertical axis depicts the phase shift of each geometry, suggesting 0 to 2 phase coverage.(b) Difference between the magnitude of  2,2 extracted from the simulation and analytical calculations.Full phase coverage can again be achieved, as before.

Figure 3 .
Figure 3. Simulated field profiles of different examples bilayer meta-atoms (in the xz-plane).The input polarization is   (in plane).(a)Both top and bottom nanofin are at resonance; exhibiting strong coupling.(b) Top nanofin is at resonance and is larger than the bottom one; causing back reflections.(c) Top nanofin is at resonance but is smaller in size than the bottom one; exhibiting very weak coupling.(d) Bottom nanofin is at resonance.The top nanofin dimensions are fairly large (170 nm × 170 nm); showing strong coupling.(e) Bottom nanofin at resonance.The top nanofin dimensions are fairly small (90 nm × 90 nm); effectively decoupled.(f) Neither bottom nor top nanofin at resonance.Top nanofin is larger than the bottom one; introducing back reflections.(g) Neither the bottom nor the top nanofins are at resonance.Top nanofin is smaller than the bottom one; a perfectly decoupled geometry.

Figure 4 .
Figure 4. single-layer metasurface working in reflection.(a) Phase shift obtained from the analytical product.(b) Phase shift obtained from FDTD simulation of a single-layer metasurface working in reflection without the introduction of the spacer.

Figure 5 .
Figure 5. Magnitude of the four elements of the Jones matrix of a bilayer metasurface whose bottom nanofins dimensions are set to be 134 nm × 202 nm while the dimensions of the top nanofin are sweeping between 50 and 250 nm.

Figure 6 .
Figure 6.Phase of the four elements of the Jones matrix of a bilayer metasurface whose bottom nanofins dimensions are set to be 134 nm × 202 nm while the dimensions of the top nanofin are sweeping between 50 and 250 nm.

Figure 7 .
Figure 7. Error in phase of the four elements of the Jones matrix of a bilayer metasurface whose bottom nanofins dimensions are set to be 134 nm × 202 nm while the dimensions of the top nanofin are sweeping between 50 and 250 nm.

Figure 8 .
Figure 8. Simulation data for two-dimensional parameter sweeps of Si rectangular fins (ℎ = 1000 nm).(a) Power Transmission   for x-polarized light.The total power passing through a monitor above the structure relative to the source.(b) Phase shift   on x-polarized light.The phase shift has been computed as the ratio between the phase collected in the center of the far field projection of a monitor above the structure when the nanofin is present on top of the substrate and the phase at the same monitor when only the silica substrate is present.Units in radians.

Figure 9 (
b) shows the case when the top nanofin is larger than the bottom one.Due to the size mismatch between the bottom and top nanofins, the effect of back reflections result in a large deviation between the responses of the S-matrices and the T-matrices.More specifically, the Fabry-Perot cavity effect created between the top nanofin and the substrate is emphasized here by looking at the field profile and the associated interference fringes.Figure9(c) represents the third case where two identical nanofins operate at resonance.