Optical metasurfaces have drawn significant attention in recent years due to their ability to control various properties of light [1,2]. They have especially proved to be very efficient in tailoring the polarization states of light [3–7]. Changes in the polarization states are governed by the polarimetric properties of the metasurface that can be modeled through different formalisms. In particular, the Jones matrix, coherency matrix, Poincaré sphere, and Stokes-Mueller matrix formalism have been developed to define these polarization properties. The Stokes-Mueller matrix offers an unambiguous and complete representation of the polarimetric properties of optical systems, especially for rough or heterogeneous metasurfaces [8–10]. This powerful technique can deal with polarized, non-polarized, or partially polarized light and allows for a complete characterization of the polarization properties of optical systems [11,12]. While Mueller matrix spectroscopy has played an important role in the experimental measurement of polarimetric properties of metasurfaces [13], the numerical computations of the Stokes-Mueller matrix elements of metasurfaces remain relatively limited [14,15].

In recent years, several numerical and experimental studies have explored the azimuth-resolved polarimetric properties of achiral metasurfaces when illuminated at specific polar angles of incidence [14,16]. The polarimetric properties of pseudo-chiral metasurfaces under specific polar and azimuthal angles of illumination have also been investigated [17,18]. On the other hand, Arteaga et al. have conducted several studies on the Mueller matrices of anisotropic and bianisotropic materials and their symmetries [19,20]. However, a comprehensive study of optical metasurfaces with specific macroscopic chiral properties illuminated across a full range of polar and azimuthal incidence angles, and a comparative analysis of corresponding symmetries of their Mueller matrices remains unexplored. In 2010, Demésy et al. presented a scattered field formulation using the Finite Element Method to calculate the diffracted vector amplitudes for arbitrarily shaped two-dimensional (or crossed or bi-periodic) gratings [21,22]. The advantage of this numerical method lies in its complete independence from the geometry, material, number of resonators, and illumination conditions. In this work, we benefit from this versatile approach to model the polarimetric properties of achiral, pseudo-chiral, and chiral metasurfaces through the calculation of their Mueller matrices. As the chiral metasurface, we introduce a novel design using achiral building blocks.

The resulting model, based on open-source finite element codes [23,24] is available online as a template model [25,26]. In this formulation, the metasurface is considered as a two-dimensional periodic structure composed of nano-resonators bearing on a substrate. The diffraction problem is treated as a radiation problem whose sources are localized within the diffracting resonators. The work presented here extends the model to rigorously calculate the complex Rayleigh coefficients in reflection and transmission, for both linear and circularly polarized light, followed by the rigorous calculation of every Mueller matrix element from the specularly reflected order. We take advantage of a parallelization technique to compute the polarimetric properties for illumination conditions across a wide range of polar and azimuthal angles of incidence. We also assess and compare the polarization-control capabilities of these three metasurfaces using the Poincaré sphere representation.

1. Computation of the Mueller matrix with the finite element method

A monochromatic light beam incident on a sample surface is conventionally represented by its electric field vector given by

 

Fig. 1. (a) 3D configuration of the unit cell of a metasurface. (b) 3D configuration of the semi-infinite metasurface with directions of orientation of the incident wave vector and the electric field. The green plane corresponds to the plane of incidence containing the metasurface normal along the z-axis, and the incident wavevector $\mathbf {k^+}$, and the wavevector is orthogonal to the red plane. (c) Schematic description of the numerical formulation.

Download Full Size | PDF

The objective of this work is to calculate the complete Jones and Mueller matrices of meta-surfaces with various shapes of scatterers using full-wave simulations (see Numerical code). This formulation has been able to accurately calculate the electromagnetic fields diffracted by bi-periodically arranged scatterers of arbitrary shape, embedded in multilayered structures, illuminated by a plane wave with arbitrary polar and azimuthal angles of incidence and arbitrary polarization [21]. We extend this model to post-process the reflection amplitudes constituting the Jones matrix, followed by the calculation of the Mueller matrix.

Figure 1(a) depicts a schematic configuration of the addressed unit cell in the Cartesian coordinate system. It consists of two homogeneous, isotropic, lossless layers called the substrate ($z<0$), and the superstrate ($z>0$). The unit cell is extended into the semi-finite metasurface shown in Fig. 1(b) by applying Bloch periodic conditions on the lateral surfaces $x_{1}$, $x_{2}$, $y_{1}$ and $y_{2}$. The size of the unit cell in the $x$ and $y$ directions, $P_{x}$ and $P_{y}$, define the periodicity of the metasurface. To truncate the infinite superstrate and substrate regions, they are extended by a set of Perfectly Matched Layers (PML) along the $z$-axis. Therefore, the diffracted field decays exponentially within the PMLs along the $z$-axis and propagates in the infinite $Oxy$ plane. The superstrate, substrate, and scatterers are characterized by their complex-valued relative permittivities and relative permeabilities $\epsilon _+$, $\mu _+$, $\epsilon _-$, $\mu _-$, $\epsilon _{scat}$ and $\mu _{scat}$, respectively. The metasurface is illuminated from the superstrate, with an incident field given by

In this expression, $k_{0}=\frac {\omega }{c}$ is the incident wave vector in free space while $\epsilon _r$ and $\mu _r$ represent the relative permittivity and permeability functions (of space and the driving frequency $\omega$) that entirely define the opto-geometric parameters of the metasurface. The wavenumbers in the superstrate and substrate are given by $k^{+}=k_{0}\sqrt {\epsilon _+\mu _+}$ and $k^{-}=k_{0}\sqrt {\epsilon _-\mu _-}$, respectively. A scattered field formulation [21] is used to convert the vector diffraction problem with faraway sources into a radiation problem with localized sources, the latter being more compliant with the bounded nature of the finite element computational box. The numerical approach is described in Fig. 1(c) with a simple schematic block diagram and the working principle is recalled hereafter. The problem of diffraction in the metasurface is carried out by addressing separately the configurations with and without resonators. The total field can indeed be expressed as :

Equation 6 can be readily solved analytically since it relies on the Fresnel coefficients of the interface. More precisely, the annex field is the superposition of two plane waves in the superstrate (the incident field $\mathbf {E}_{inc}$, the reflected field $\mathbf {E}_1^r$) and the transmitted field denoted $\mathbf {E}_1^t$ in the substrate.

Now, making use of the linearity of Eqs. (4,6) together with Eq. (5) allows to establish the propagation equation for the unknown field $\mathbf {U}$ :

The right-hand side of the above equation can be interpreted as sources localized in the scatterers, since outside the scatterers $\epsilon _r=\epsilon _\pm$ and $\mu _r=\mu _\pm$. The problem of diffraction has now essentially changed to a problem of radiation. Equation (8) can be solved through its weak formulation as elaborated in Ref. [28] using periodic conditions in its lateral sides and perfectly matched layers along with boundary conditions on its top and bottom sides.

From the numerical point of view, all geometries and conformal meshes are performed with the Gmsh software [23], and the finite element computations are implemented thanks to the flexibility of the finite element software GetDP [24]. Open-source models used in this study are provided in [25,26] (see Numerical code). The size of each edge element is determined by $\lambda _0/(N\sqrt {Re(\epsilon _r)})$, where $N$ is the number of tetrahedral elements per wavelength of the electromagnetic field inside a given material ($\lambda _0/\sqrt {Re(\epsilon _r)}$), characterized as the mesh refinement factor, which will allow performing convergence tests reported in the next sections. The unknown field $\mathbf {U}$ is discretized on the tetrahedral mesh using high-order Webb hierarchical edge elements [29–31] with 26 Degrees Of Freedom (DOFs) per tetrahedron (3 DOFs per edge, 2 DOFs per face). The FEM sparse matrix inversion is carried out by the direct solver MUMPs [32] which is natively interfaced in the finite element software GetDP [24].

Finally, in order to define the Rayleigh coefficients let us call the total diffracted field $\mathbf {E}_d$ such that

In other words, the field diffracted by the metasurface is a summation of the outgoing reflected and transmitted fields from the annex problem, and the outgoing unknown field $\mathbf {U}$, i.e., $\mathbf {E}_d=\mathbf {E}-\mathbf {E}_{inc}=\mathbf {U}+\mathbf {E}_1^r$ in the superstrate and $\mathbf {E}_d=\mathbf {E}=\mathbf {U}+\mathbf {E}_1^t$ in the substrate. Once the diffracted field is computed, the vector amplitudes $\mathbf {e}^\pm _{i,j}$, known as Rayleigh coefficients, of the reflected (superscript $^+$) and transmitted (superscript $^-$) diffraction orders can be computed by integrating cuts of the diffracted field in the superstrate and substrate [21] :

The calculation of these quantities allows for the determination of the reflected and transmitted diffraction efficiencies defined as:

For the specular orders, whose wavevectors remain in the plane of incidence, the unit vector $\mathbf {s}$ normal to the plane of incidence and defined as $[\sin (\phi _0),-\cos (\phi _0),0]=\hat {\mathbf {z}}\times \mathbf {k}^+/|\hat {\mathbf {z}}\times \mathbf {k}^+|$ remains unchanged. However, the unit vector $\mathbf {p}$ takes different directions for the transmitted and reflected specular orders, owing to the different wavevectors for these two cases, and are defined as $\mathbf {p}_{0,0}^\pm =\mathbf {k}^\pm /|\mathbf {k}^\pm |\times \mathbf {s}$. Note that the non-specular orders lie outside the plane of incidence and the unit vector $\mathbf {s}$ loses its meaning. A workaround could be to define extra “outgoing planes” with the wavevector of the considered order and the normal to the surface: $\mathbf {s}_{i,j}^\pm =\hat {\mathbf {z}}\times \mathbf {k}_{i,j}^\pm /|\hat {\mathbf {z}}\times \mathbf {k}_{i,j}^\pm |$, then $\mathbf {p}_{i,j}^\pm =\mathbf {k}_{i,j}^\pm /|\mathbf {k}_{i,j}^\pm |\times \mathbf {s}_{i,j}^\pm$.

In this paper, we focus on the reflected specular order, and the simulations are performed for both $s$ and $p$-polarized incident plane waves. For a specific illumination condition, the specular reflection Rayleigh coefficient $\textbf {e}^{+}_{0,0}$ is calculated at the top surface of the superstrate (i.e. the interface with the PML) and projected on the $s$ and $p$ basis (scalar product with $\mathbf {s}$ or $\mathbf {p}_{0,0}^+$), constituting the elements of the Jones matrix labeled by $r_{mn}$. For the sake of clarity, if the incident beam is for instance $p$-polarized ($n=p$), the reflection coefficient along the direction of $s$-polarization ($m=s$) will be given by $r_{sp}$. Finally, all the elements of the Mueller matrix are deduced from the Jones matrix, using the theoretical framework provided in Ref. [14], given by:

2. Design of metasurfaces and analysis of their polarimetric properties

We now illustrate the generality and relevance of the presented numerical model by designing three kinds of metasurfaces with specific polarimetric properties. The accuracy of this approach is assessed for both Mie-resonant and plasmonic scatterers. The numerical method presented in the previous section is implemented to characterize the polarimetric properties of achiral, pseudo-chiral, and chiral metasurfaces. More precisely, we consider (i) an achiral metasurface composed of Si-based nano-disks, (ii) a plasmonic metasurface composed of gold U-shaped resonators, and (iii) the case of a chiral structure composed of four Mie-resonant scatterers per cell with heights varying on left or right circular directions. In the latter case, we benefit from the versatility of the numerical approach to design left or right-handed Si-based quadrumers exhibiting strong circular dichroism.

2.1 Achiral metasurface with Mie-resonant Si nano-disks

The first metasurface consists of an array of dielectric nano-disks composed of silicon (see Fig. 2(a)). The building blocks of this metasurface have proven to exhibit chiral density in their nearfield but no circular dichroism of their own [33–36]. The radius, height, and periodicity of the pillars are given by $r = {80}{\textrm {nm}}, h = {130}{\textrm {nm}}$ and $P_x = P_y = {420}{\textrm {nm}}$. The pillars are bearing on a glass substrate ($\epsilon _- = 2.25$) while the superstrate is made of air ($\epsilon _+ = 1$). The metasurface is illuminated by $s$ and $p$ polarized incident plane waves with an incidence defined by the angles $\theta$ and $\phi$.

 

Fig. 2. (a) Top view of the metasurface with Si nano-disks (top) with the top view of the 3D mesh (right) and 3D view of the unit cell (bottom) comprising one nano-disk bearing on the substrate. (b) Computation accuracy with respect to the mesh refinement by calculating and plotting the relative differential error $m_{D}$ for non-zero Mueller matrix elements $M_{11}$, $M_{12}$, $M_{33}$ and $M_{34}$ in logarithmic scale (top) and the global energy balance for $s$ and $p$ polarization (bottom) as a function of $N$. The frequency-dependent dielectric parameters of Si in the mentioned wavelength range have been taken from the data by Schinke et al. [37]. (c) Polar plots of the Mueller matrix elements for $\phi =90^{\circ}$ and $\theta$ varying from $0^{\circ}-45^{\circ}$ for the wavelength range $400-{800}{\textrm {nm}}$. (d) Azimuthal plots of the Mueller matrix elements for $\theta =45^{\circ}$ and $\phi = 0^{\circ}-360^{\circ}$ in the same spectral range.

Download Full Size | PDF

2.1.1 Convergence test

The Finite Element Method is of high interest for modeling the polarization properties of metasurfaces since the mesh can be very easily adapted to scatterers of arbitrary shape. The crucial point is to assess its convergence when refining the mesh size. In order to identify an appropriate value of the mesh refinement parameter $N$, a convergence test is performed for oblique incidence ($\theta =45^{\circ}$) at $\lambda = {600}{\textrm {nm}}$ for a range $N = \{2,\ldots,12\}$. A parameter was defined to quantify the relative differential error of the Mueller matrix elements, as $m_{D}=|\frac {M_{N}-M_{N-1}}{M_{N_{max}}}|$, where $M_{N}$ is an element of the Mueller matrix $M_{ij}$ for a particular value of the mesh refinement parameter $N$. Figure 2(b) represents the convergence of $m_{D}$ for the non-zero elements $M_{11}$, $M_{12}$, $M_{33}$ and $M_{34}$ in logarithmic scale (top) and the global energy balance for $s$ and $p$ polarization (bottom) as a function of $N$. The total energy efficiency for both $s$ and $p$ polarizations tends toward 1 as the mesh refinement parameter $N$ increases. As, for example, the $|M_{12}|$ element obtained for $N=7$ $(0.0370)$ and $N=12$ $(0.0368)$ leads to a relative difference of $0.5{\% }$ only, a trade-off between accuracy and computation time is found by taking $N=7$ for all calculations. However, in order to maintain a uniform meshing throughout the spectrum, the mesh parameter was fixed at $\lambda _{min}/(N\sqrt {Re(\epsilon _r)})$ instead of $\lambda _{0}/(N\sqrt {Re(\epsilon _r)})$ for the dielectric parts of all three systems.

2.1.2 Dependence on θ, φ and λ

To compute the Mueller matrix elements, numerical simulations are performed by varying two of the three illumination parameters ($\theta$, $\phi$, and $\lambda$) while keeping the third one ($\theta$ or $\phi$) fixed. All the elements are normalized with respect to $M_{11}$. Figure 2(c) displays the polar plot spectra for the 16 elements ($m_{ij} = M_{ij}/M_{11}$) of the normalized Mueller matrix, with $m_{11}$ replaced by a diagram showing the varying parameters. The spectra are obtained by illuminating the metasurface in the wavelength range of $\lambda =400\!-\!{800}{\textrm {nm}}$, at a fixed azimuthal angle of incidence $\phi =90^{\circ}$ and varying polar angles from $\theta =0^{\circ}\!\!-\!45^{\circ}$. The plots have a clockwise angle increment, where the wavelength increases radially from the inner circle of ${400}{\textrm {nm}}$ to the outer circle of ${800}{\textrm {nm}}$. It can be observed that at normal incidence, the matrix is diagonal, and consists of two distinct non-zero elements, $m_{22}$ and $m_{33}$. On the other hand, as the incidence becomes oblique, two other non-zero unique elements, $m_{12}$ and $m_{34}$, are introduced, and the matrix transforms into a block-diagonal form for each $\theta$,

This characteristic can be traced back to a diagonal Jones matrix ($r_{ps}=r_{sp}=0$), that results from the isotropic, achiral nature of the resonators, and the mirror symmetry of the metasurface with respect to the plane of incidence [19,38]. Among these elements, $m_{12}$, $m_{33}$, and $m_{34}$ tend to display a dispersive signature with the angle of incidence deviating from normal, while the others remain unchanged at zero. The azimuthal plots of the metasurface, obtained by varying $\phi$ from $0^{\circ}\!\!-\!360^{\circ}$ at $\theta = 45^{\circ}$ are shown in Fig. 2(d). It is interesting to note that the matrix symmetry for each $\phi$ has changed and now appears in the form:

The off-block-diagonal elements of the matrix exhibit significant non-zero values for all azimuthal angles except at multiples of $45^{\circ}$. These non-zero values arise due to the oblique incidence on the metasurface which leads to cross-polarization conversion ($p$-$s$ or $s$-$p$). As a result, the Jones matrix becomes non-diagonal but antisymmetric ($r_{ps}=-r_{sp}$) [38]. Under oblique incidence, the metasurface satisfies only one symmetry condition which is the $180^{\circ}$ rotational symmetry around its normal, for any $\phi$ except multiples of $45^{\circ}$. For those specific angles, the metasurface satisfies an additional mirror symmetry condition with respect to the incidence plane, due to the $45^{\circ}$ symmetry of the square lattice. Therefore, the matrix simplifies at these points and retains the same form as in Eq. (14). Furthermore, it can be observed that under oblique incidence, the off-block-diagonal elements display a $180^{\circ}$ change of phase every $45^{\circ}$ interval of rotation, which comes from the $\pi$ phase transition of $r_{ps}$ and $r_{sp}$ at these angles. So, with respect to an axis cutting through them every $45^{\circ}$ interval, the off-block-diagonal elements are antisymmetric, whereas the block-diagonal elements are symmetric. This has been observed previously in plasmonic metasurfaces, where the overlapping between LSPR and Rayleigh anomalies played an important role [14,16,39], but not in dielectric metasurfaces.

Element $m_{14}$ is often attributed to the presence of circular dichroism (CD) generated by the metasurface. However, there should be no CD observed in the case of Si nanodisks. It has been established that the element $m_{14}$ could also contain information about linear dichroism and birefringence which could be mistaken for CD [40,41]. In the case of a Mueller matrix determined in reflection at oblique incidence, the so-called polar decomposition [42] allows for the determination of the different polarimetric contributions to the Mueller matrix. In our calculation, there is no depolarization, so the analysis was performed by assuming that the Mueller matrix of the Si nanodisks could be represented as the product of a pure retardance and a pure diattenuation matrix. The retardance matrix was then decomposed as the product of a pure retardator followed by a pure optical rotator [43]. The optical rotator would describe the contribution of CD to the Mueller matrix of the nanodisks. The parameters of the retardator (retardance and fast axis orientation) and of the optical rotator (rotation angle) were fitted to the Mueller matrix of Fig. 2(d). Figure 3(c) presents the rotation angle associated with the optical rotator extracted from the fitting. It can be seen that except for a few numerical noises, the optical rotation was zero indicating that the non-zero values of the elements $m_{14}$ and $m_{23}$ presented in Figs. 3(a) and (b) were arising from linear dichroism and birefringence.

 

Fig. 3. Azimuthal variation (first quadrant) of the elements (a) $m_{14}$ associated with circular dichroism and (b) $m_{23}$ associated with circular birefringence, and (c) the corresponding optical rotation extracted via Lu-Chipman decomposition.

Download Full Size | PDF

2.2 Pseudo-chiral metasurface with U-shaped resonators

Let us now investigate a more complex shape associated with a dispersive and metallic material. The metasurface under study is a pseudo-chiral metasurface composed of a two-dimensional square array of U-shaped gold resonators on a glass substrate. In recent years, they have emerged as a unique class of metasurfaces in exhibiting optical activity. The term ’Pseudo-chirality’ refers to the optical activity observed in metasurfaces that leverage simple achiral building blocks to achieve polarization-sensitive responses that closely resemble those of their chiral counterparts. This form of chiral behavior is referred to as extrinsic chirality, first reported by Verbiest et al. [44]. Over the past decades, several studies have numerically and experimentally confirmed the presence of circular dichroism in individual and arrays of U-shaped nanostructures arising from magneto-electric coupling [17,45–48].

The dimensions of the structure, shown in Fig. 4(a), are set to $P_x=P_y={400}{\textrm {nm}}$, $L_x=L_y={200}{\textrm {nm}}$, $d_x=d_y={70}{\textrm {nm}}$, and $d_z={40}{\textrm {nm}}$, respectively. The sharp corners were rounded with a radius of ${20}{\textrm {nm}}$ in order to make the design more realistic. The frequency-dependent dielectric parameters of gold were taken from Rakic et. al. [49]. The meshing parameter was set to $N=7$ for the substrate and superstrate. However, due to the metallic nature of gold, the size of the tetrahedral element for the scatterer was fixed at ${15}{\textrm {nm}}$ instead of $\lambda _0/(N\sqrt {Re(\epsilon _{scat})})$ to effectively capture the decaying electric field within the skin depth.

 

Fig. 4. (a) 3D configuration of the metasurface with U-shaped scatterers (top) and an individual scatterer from a top angular view with the directions of orientation of the incident light beam. The inset shows the geometry of the U-shape scatterer in a top view. (b) Polar plots of the Mueller matrix elements as functions of wavelength varying from $400$ to ${1200}{\textrm {nm}}$ and polar angle $\theta$ from $0^{\circ} - 45^{\circ}$. The azimuthal angle of illumination is set to $\phi =90^{\circ}$. (c) Azimuthal plots of the Mueller matrix elements as functions of the azimuthal angle $\phi = 0^{\circ} - 360^{\circ}$ in the same spectral range. The polar angle of illumination is fixed at $\theta =45^{\circ}$. (d) The azimuthal plots of the real parts of the Jones matrix elements in the same configuration. The blue-filled spectra show the spectral variation of $r_{sp}$ and $r_{ps}$ at $90^{\circ}$ and $270^{\circ}$, respectively.

Download Full Size | PDF

The polar plots of the 16 elements $m_{ij}$ are illustrated in Fig. 4(b) and Fig. 4(c) for polar angles ranging from $\theta =0^{\circ}-45^{\circ}$ at increments of $3^{\circ}$, where the azimuthal angle is fixed at $\phi =90^{\circ}$ and for an azimuthal angle ranging from $\phi =0^{\circ}-360^{\circ}$ at $5^{\circ}$ intervals, where the polar angle is fixed at $\theta =45^{\circ}$, respectively. For both cases, each angle was swept in the wavelength range $400-{1200}{\textrm {nm}}$ by ${5}{\textrm {nm}}$ intervals. It can be observed from Fig. 4(b) that at normal incidence, the metasurface behaves like an isotropic and achiral photonic component, with its Mueller matrix satisfying the same symmetry conditions seen in Eq. (14). However, for oblique incidence at $\phi =90^{\circ}$, i.e. along the ($Oyz$) plane, unlike the isotropic case, the matrix exhibits a different type of symmetry of the following form:

This property results from the symmetry of the Jones matrix at this angle $\phi =90^{\circ}$ [17] and for $\phi =270^{\circ}$ ($r_{ps}=r_{sp}$, see Fig. 4(d)). It can be attributed to the mirror symmetry of the metasurface with respect to the plane ($Oxz$) perpendicular to the plane of incidence ($Oyz$) [38]. Contrary to the isotropic disks, magneto-electric coupling plays an important role in this resonator geometry, giving rise to a cross-polarization conversion [17,45,47]. On the other hand, even in oblique incidence, all the off-block-diagonal elements are null for $\phi =0^{\circ}$ as shown in Fig 4(c), and the Mueller matrix again follows the symmetry of Eq. (14). It hints towards the absence of cross-polarization conversion at $\phi =0$, which can be confirmed from the azimuthal polar plots of the Jones matrix, shown in Fig 4(d). For all other azimuthal angles, any mirror symmetry is broken and the Mueller matrix takes the general form:

However, in general, all the elements demonstrate a symmetry condition of $m_{ij}(\phi )=\pm m_{ji}(\phi +\pi )$, regardless of the azimuthal angles. Figure 4(d) also reveals that the diagonal elements of the Jones matrix ($r_{pp}$ and $r_{ss}$) are symmetric with respect to the vertical axis connecting $\phi =0^{\circ}- 180^{\circ}$ ($Oxz$ plane), whereas the off-diagonal elements are antisymmetric. This feature translates to the phenomenon that the signs of circular and $\pm 45^{\circ}$ linear dichroism (birefringence) are reversed as the incidence plane is rotated by $180^{\circ}$, whereas linear dichroism (birefringence) remains the same. Since all the elements are non-zero, the specific linear and circular polarimetric properties can be extracted by employing one of the several decomposition methods, such as polar decomposition [42], reverse decomposition [50], symmetric polar decomposition [51] etc.

2.3 Height-modulated chiral quadrumer metasurface

Achiral dielectric resonators exhibiting strong chiral nearfields have already proven to be important building blocks to create chiral metasurfaces [52–55]. We now benefit from the versatility of the FEM formulation to design an all-dielectric chiral metasurface composed of height-varying nano-discs exhibiting strong circular dichroism. The configuration is illustrated in Fig. 5(a). The metasurface features four nano-disks per unit cell, with gradually decreasing heights in clockwise and counterclockwise orientations. The disks have a fixed radius of ${60}{\textrm {nm}}$ and are arranged in a two-dimensional array along the $x$ and $y$ directions with a fixed periodicity of $P_x = P_y = {420}{\textrm {nm}}$. The height of disk-1 is constant at ${130}{\textrm {nm}}$, while the heights of the neighboring disks decrease by $h_d$ in the right-circular (RC) or left-circular (LC) direction compared to their respective previous disks. The distances between neighboring pillars in the $x$ and $y$ directions are fixed at $d_x = d_y = {40}{\textrm {nm}}$. The objective is to maximize the circular dichroism of the metasurface by optimizing the $m_{14}$ elements as a function of the height difference $h_d$.

 

Fig. 5. (a) Configuration of the unit cell of a metasurface with four nano-pillars of decreasing height in a right-circular (top-left) and left-circular orientation (bottom-left), top view of the metasurface with the unit cell highlighted with a black square (top-right), and side-view of the right-circular (RC) metasurface showing the height difference between three disks (bottom-right). (b) Normalized spectra of the Mueller matrix elements for the LC (orange) and RC (green) metasurfaces, and with equally tall disks (red) at normal incidence. (c) Polar plots of the Mueller matrix elements for LC metasurface as functions of the polar angle $\theta = 0^{\circ} - 45^{\circ}$ in the spectral range $400-{1000}{\textrm {nm}}$, with the azimuthal angle fixed at $\phi =0^{\circ}$. (d) Full azimuthal polar plots of the off-block-diagonal Matrix elements for LC metasurface at normal incidence in the spectral range $400-{800}{\textrm {nm}}$. The black dotted circles correspond to Rayleigh anomalies.

Download Full Size | PDF

Figure 5(a) illustrates the schematic design of the metasurfaces of opposite handedness with top and side views of the unit cell. Figure 5(b) showcases the Mueller matrix spectra of these two metasurfaces (with $h_d = {15}{\textrm {nm}}$) along with the case of self-similar disks ($h_d = 0$). The metasurface with uniformly tall multi-disks mimics similar behavior to the achiral metasurface with a single disk in normal incidence discussed in subsection 2.1, whereas, significant changes can be observed when modifying the heights of the disks. It can also be noted that $m_{12}$ ($m_{34}$), $m_{14}$ ($m_{23}$) elements, representing linear and circular dichroism (birefringence) respectively, exhibit exact and opposite responses upon reversing the handedness of the metasurface. On the other hand, the $m_{13}$ ($m_{24}$) element, associated with $45^{\circ}$ dichroism (birefringence), is the same for both right- and left-handed circular (RC and LC) metasurfaces. This phenomenon can be explained by the fact that when the handedness is reversed, disks 2 and 4 exchange positions, while disks 1 and 3 remain in the same place. As a result, the effective ratio of matter-volume encountered by horizontal and vertical polarization components on both sides of the unit cell center changes from RC to LC. However, the $\pm 45^{\circ}$ components encounter the same volume distribution, leading to similar dichroism for both handednesses. The opposite circular dichroism results from the inherent chirality present in the twist that interacts with the chirality of circular polarization components. The polar angle-resolved Mueller matrix of the LC metasurface is shown in Fig. 5(c). It can be observed that while at normal incidence the matrix takes the form of Eq. (15), it tends to lose symmetry as the incidence becomes oblique. The multiple lines visible in the block-diagonal elements as compared to the case of a single disk in oblique incidence result from the coupling effect between the disks. Figure 5(d) presents the $\phi$-resolved polar plots of the off-block-diagonal elements. All the elements obey $m_{ij}=\pm m_{ij}(\phi +\pi )$, which is true for all the block-diagonal elements as well. It can be observed that the circular dichroism and birefringence responses are very close to the Rayleigh anomalies shown in black dotted lines. Note that circular dichroism is constant for any value of $\phi$, whereas circular birefringence changes phase after every $45^{\circ}$ interval. This phenomenon can be explained by the changing phase difference between $r_{sp}$ and $r_{ps}$ at different $\phi$ locations. It is interesting to note that the $m_{13}$, $m_{23}$, $m_{24}$ elements and their symmetric counterparts exhibit antisymmetry with respect to the axes which are about $\sim 10^{\circ}$ rotated from the vertical and horizontal axes. A possible explanation could come from the symmetry breaking introduced by the height modulation, which results in the rotation of this axis.

The circular dichroism can be maximized through parametric optimization, in particular in terms of the inter-pillar distances $d_x$, $d_y$, and the height modulation $h_d$. The dependence of $m_{14}$, associated with the circular dichroism, is shown in Fig. 6(a)-(c). Figure 6(a) depicts the influence of the inter-pillar distance on the $m_{14}$ element, with $d_x = d_y = d$ varying between $10-{170}{\textrm {nm}}$. The plots reveal that significant resonances occur within a narrow range of inter-pillar distances, specifically between $d_x = d_y = d = 40-{60}{\textrm {nm}}$, with a pronounced peak resonance observed at $d={50}{\textrm {nm}}$. This indicates that the coupling effect between the disks within a single unit cell reaches its maximum at this distance, subsequently decreasing as the inter-pillar distance is further increased. However, interestingly, when $d$ is further increased, $m_{14}$ experiences another enhancement, reaching its highest value at $d={130}{\textrm {nm}}$. This can be attributed to the coupling between disks in neighboring cells, which form a unit cell (marked with a black box) of the same handedness. As a result, the resonances in the $m_{14}$ spectra exhibit a symmetric pattern about $d={90}{\textrm {nm}}$, representing the distance at which coupling between the disks is minimized. Figure 6(b) and 6(c) illustrate the dependence of $m_{14}$ on $d_x$ and $d_y$, respectively, while keeping the other one fixed at $d = {50}{\textrm {nm}}$. Figure 6(b) reveals that as $d_x$ increases, the peaks suffer redshifts and eventually go through a sign change near the Rayleigh line. For further increase of $d_x$ beyond ${90}{\textrm {nm}}$, the peaks get reversed. This arises from the coupling between disks of the neighboring cells in the $x$-direction which now effectively form a cell of opposite-handedness. The highest magnitude of $m_{14}$ is achieved for $d_x={50}{\textrm {nm}}$ and ${130}{\textrm {nm}}$. A similar sign-reversal phenomenon can be observed for $d_y$ variation as well, except for a smaller redshift of peaks. The highest $m_{14}$ is achieved for $d_y={60}{\textrm {nm}}$ and ${120}{\textrm {nm}}$.

 

Fig. 6. (a) Variation of the $m_{14}$ element as a function of $d_x=d_y$ and wavelength, with $h_d$ fixed at ${15}{\textrm {nm}}$. The white dotted line corresponds to Rayleigh anomaly. (b) Spectral dependence of $m_{14}$ with $d_x$ and (c) $d_y$ while the inter-pillar distances in the respective orthogonal directions are kept fixed at $d={50}{\textrm {nm}}$. (d) $m_{14}$ variation with the height modulation $h_d$. (e) The optimized $m_{14}$ spectra associated with circular dichroism. The Rayleigh line is shown with the black dotted line.

Download Full Size | PDF

Figure 6(d) illustrates the dependence of $m_{14}$ on the height modulation $h_d$. It is evident that the intrinsic chirality resulting from the arrangement of the unit cell intensifies as $h_d$ increases, reaching its highest point at $h_d={15}{\textrm {nm}}$, and then diminishes as $h_d$ continues to increase. The Rayleigh anomaly close to the resonance is shown with the white dotted line. Figure 6(e) shows the $m_{14}$ spectrum corresponding to the optimized $d_x$, $d_y$ and $h_d$ parameters. The peak of the spectrum is observed at ${619.5}{\textrm {nm}}$, while the Rayleigh line is at ${630}{\textrm {nm}}$.

2.3.1 Analysis of polarization-modulation using the Poincaré sphere

It is worth integrating the Stokes-Mueller matrix representation with the Poincaré sphere to extract the polarization states of the output beam for the aforementioned metasurfaces. Poincaré sphere is a well-established approach to portray polarization states as points on its surface, with their coordinates corresponding to the Stokes parameters S$_1$, S$_2$, and S$_3$ in the Cartesian system [56].

The wavelength-resolved trajectory of the reflected polarization for the three metasurfaces under study illuminated in normal incidence with a left-circular polarization is shown in Fig. 7. The associated wavelength for the trajectory is indicated by the color bar. It can be observed that the achiral Si-based metasurface maintains the incident polarization across all wavelengths, and induces a $\pi$ phase shift, resulting in a right-circularly polarized reflected beam. In contrast, the gold U-shaped metasurface transforms the incident circular polarization into an elliptical polarization that depends on the wavelength. The resulting trajectory is confined to the upper section of the sphere, suggesting a small degree of polarization modulation. The chiral multi-disk metasurface, on the other hand, has a much broader polarization trajectory. The reflected polarization trajectory is not only confined to the upper half but rather crosses the equatorial line twice to reach the lower half, meaning that the metasurface behaves as a quarter waveplate at those two wavelengths.

 

Fig. 7. Poincaré sphere representing the (a) incident left-circular polarization, (b) the wavelength-resolved polarization trajectory of the light reflected from the achiral metasurface with single nanodisk (Subsection 2.1), (c) pseudo chiral metasurface with U-shaped resonators (Subsection 2.2), and (d) 3D chiral metasurface with height modulated multi-nanodisks at normal incidence (Subsection 2.3). The color bar represents the wavelengths of the trajectory.

Download Full Size | PDF

3. Conclusion

To conclude, we presented a numerical method based on open-source modules developed to compute the polarimetric properties of periodic metasurfaces with the Finite Element Method. The versatility and capability to model intricate shapes and rounded edges serve as the driving factors behind the adoption of the finite element method. We benefit from this numerical method to model and analyze the polarimetric properties of achiral, pseudo-chiral, and chiral metasurfaces. The computation of the complete Mueller matrices reveals the polarimetric properties of the associated metasurfaces under normal and oblique incidences. Additionally, the symmetries of the Mueller matrix elements are explained in relation to the metasurfaces’ symmetries. We designed and optimized the element representing the circular dichroism of a chiral metasurface composed of 4 silicon nanodiscs of varying height per cell. The circular dichroism is optimized in normal incidence with respect to the interparticle distance and twist of the 4 Si pillars. The spectral distribution around the Poincaré sphere shows that this Si-based metasurface exhibits the properties of a quarter waveplate at two different frequencies. We conclude our study with a comparison between the trajectories of the reflected polarizations retrieved from the metasurfaces.

Numerical code

We provide a bash script [26] allowing to reproduce the polar plots (cf. Figures 4(b) and 4(c)) given for the U-shaped resonator array. This script successively calls Gmsh [23] for the geometry definition and meshing, GetDP [24] for solving the electromagnetic problem for a specific $(\lambda,\theta,\phi,\psi )$ incident plane wave, and Python for plotting the spectrally and angularly resolved Mueller matrix elements. For computational efficiency, looping over the angles and wavelengths is performed in an embarrassingly parallel way using the “GNU parallel” library [57]. This script has been successfully tested on various Linux and macOS machines.