Over the past few years, metasurfaces have gained increasing attention due to their low form factor and incredible light-control capabilities. As a consequence, we have seen a plethora of promising applications emerge, such as light refraction,1 polarization control,2 and holography.3 Naturally, this has further sparked an increasing interest towards ever more advanced applications, such as optical analog processing using nonlocal interactions4^,5 or sensing, detection, and polarization multiplexing using chirality.6^–10

These advances have prompted the need for appropriate modeling approaches able to predict the electromagnetic behavior of a metasurface so as to leverage the adequate available degrees of freedom that they offer for an optimal implementation of the desired specifications or simply to achieve higher performance. Conveniently, several metasurface modeling techniques have been proposed throughout the years,11^–16 which essentially all revolve around the same concept, i.e., replacing the metasurface by a homogeneous sheet of effective material parameters such as impedances, polarizabilities, or susceptibilities. These models are all based on a dipolar description of the metasurface electromagnetic response, and they typically include bianisotropic responses that are necessary to model chiral effects.17^–19

However, despite being generally powerful, these techniques have been shown to be limited when it comes to modeling the angular scattering properties of a metasurface.20^,21 This is particularly a problem for applications pertaining to optical analog signal processing, where accurate angular scattering control over a large range of incidence angles is crucial.4^,5 One of the main reasons for such a limitation is the fact that optical metasurfaces have unit cells that are relatively large compared to the wavelength, especially those based on dielectric resonators, implying that multipolar contributions beyond the dipolar regime start contributing significantly to their scattering response.21 Since multipolar components, such as dipoles and quadrupoles, have different angular scattering behaviors,22^–27 it is clear that metasurface modeling techniques solely based on dipolar components cannot properly model the response of a metasurface under an arbitrary illumination, such as an oblique plane wave. Thankfully, extensions to the conventional dipolar metasurface modeling techniques have been recently proposed.21^,28 They include higher-order multipolar contributions as well as nonlocal spatially dispersive responses that are necessary to adequately connect the multipolar components to the fields exciting the metasurface.29^–32

While multipolar metasurface modeling techniques exhibit promising features in terms of improved modeling accuracy and better physical insight into the scattering process, they also pose a challenging problem. Compared to dipolar modeling, the inclusion of the required nonlocal third and fourth rank effective material tensors introduces a significant number of additional degrees of freedom and makes the practical application of the modeling significantly more difficult and cumbersome. Nevertheless, this problem may be mitigated by considering that many of these degrees of freedom do not play a role in most conventional metasurfaces and can thus be neglected. For instance, parameters that are related to polarization rotation should not be taken into account for metasurfaces that do not induce polarization rotation. Knowing which parameters should or should not be included in the model may be decided, in the most simple cases, just by intuitive reasoning, which has been a common practice in the case of dipolar metasurfaces.16 For slightly more complicated cases, for instance, those that require bianisotropic responses, it is possible to determine which are the dominant parameters that should be considered using a series of numerical simulations with a complex scheme of illumination conditions.33^–35 While such approaches are powerful tools to reduce the complexity of the modeling problem, they require a complicated simulation setup and an important computational cost. Finally, it is also possible to individually test the parameters within the homogenized model to determine their angular scattering behavior and check whether this matches that which corresponds to the system and therefore may be present.21 However, this is cumbersome, and furthermore, requires intuition of the symmetries of the scattered fields. This therefore begs the need for a simple, fast, rigorous, and effective method to determine which multipolar parameters should be considered and which ones should be excluded from the model.

On the other hand, the use of chiral responses in metasurfaces for sensing, polarization control, or asymmetric transmission applications has led to fascinating works investigating the origins of such responses. It was, for instance, shown that 3D geometrically chiral structures are not necessary to achieve chiral responses.36^–50 This is possible because a chiral response may be obtained either using 2D chiral scattering particles placed on a substrate,36^–41 which breaks the symmetry of the system in the third dimension, or by illuminating asymmetric particles along specific oblique directions: a phenomenon commonly referred to as extrinsic chirality (also called pseudo-chirality).42^–46 Due to their dependence on the spatial symmetries of the metasurface or the direction of wave propagation of the illumination, these types of exotic effects turn out to be particularly difficult to grasp from intuition alone, which is deemed to be even more arduous in the case of a multipolar metasurface. This therefore suggests the need for a method to effortlessly predict the existence of a chiral response in a given metasurface structure and for a specific illumination condition.

In this work, our goal is to address these two requirements, namely, to devise a method to determine the existence of multipolar components and that of chiral responses (or other peculiar scattering effects) in the case of an arbitrary metasurface. To this end, we propose to establish a connection between the spatial symmetries of a metasurface and its effective material parameters, as well as its scattering response.

Clearly, the idea of drawing such a connection is not new, as many works have already proposed using spatial symmetries to model metamaterials. For instance, concepts pertaining to group theory have been leveraged to associate the point group of a metasurface to its dipolar material parameters51^,52 or to predict the existence of certain responses based on the symmetries of the metamaterial structure.53^–56 Spatial symmetries have also been used in several other contexts, such as the design of photonic crystals,57^,58 their relationships with reciprocity and chirality and bianisotropy,59^–66 the existence of bound states in the continuum,67^,68 or even the implementation of asymmetric nonlinear responses.69 Finally, several works have also shown a connection between the scattering response of a metasurface and its spatial symmetries.70^–77 These works typically consider the symmetries associated with the superposition of the metasurface structure with the fields that interact with it, leading to the formulation of corresponding Jones, Mueller, or scattering matrices.

Based on the existing literature, we aim at providing a simple, straightforward, and coherent framework that connects the symmetries, material tensors, and scattering response of a metasurface in a way that does not require an extensive knowledge of group theory. Thus, the proposed approach only requires listing the spatial symmetries of a metasurface, which makes it accessible to the largest possible audience. It also does not require performing any complicated numerical simulation, which makes it fast and simple to use. In addition, this framework extends the existing methods by including the presence of quadrupolar and nonlocal responses, which is crucial for properly modeling certain types of chiral responses, as we shall see. Finally, we will not restrict ourselves to the study of chiral responses but instead investigate the complete scattering response of a metasurface by computing its full scattering matrix. To facilitate the use of the proposed framework, we also provide a Python script, which is accessible on GitHub,78 and that simply requires listing the spatial symmetries of the metasurface. Finally, it should be noted that while we concentrate our attention on the topic of metasurfaces, the formalism developed thereafter also applies to 3D metamaterials, providing that they can be homogenized so that their electromagnetic response may be model via effective material parameters.

This paper is organized as follows. In Sec. 2.1, we review the general approach for transforming the dipolar material parameters of a metasurface according to an arbitrary spatial symmetry, which is a crucial step for ultimately obtaining the material parameters corresponding to the metasurface. In Sec. 2.2, we extend that approach to multipolar responses and, in Sec. 2.3, finally show how to connect spatial symmetries and multipolar material parameters. In Sec. 3, we connect the spatial symmetries of a metasurface and the fields that interact with it to its corresponding scattering matrix. Finally, in Sec. 4, we provide three examples illustrating the application of the proposed framework.

The general approach for connecting the effective material tensors of a metasurface to its spatial symmetries is based on Neumann’s principle.79^–81 This principle states that if a system, such as a crystal or an electromagnetic structure, is invariant under certain symmetry operations, then its physical properties should also be invariant under the same symmetry operations. It follows that the relationships between spatial symmetries and material tensors may be established by deriving symmetry invariance conditions that apply to the multipolar tensors describing the effective electromagnetic response of a metasurface.

We shall next review the conventional method used to derive such invariance conditions in the case of a medium described by dipolar responses. Then, we will extend these conditions to the case of multipolar responses.

It should be noted that throughout this work, we shall consider that a metasurface is an electrically thin array consisting of a subwavelength periodic arrangement of reciprocal scattering particles. The period of the array is considered small enough compared to the wavelength so that no diffraction orders exist (besides the zeroth orders in reflection and transmission), irrespective of the direction of wave propagation and the refractive index of the background media. Under these assumptions, it follows that the electromagnetic response of such a metasurface can be modeled as that of a homogeneous and uniform sheet of effective material parameters.16^,82

Let us consider a uniform and homogeneous bianisotropic metasurface whose constitutive relations are given as83[DB]=[ϵ¯¯ξ¯¯ζ¯¯μ¯¯]·[EH],(1)where ϵ¯¯ is the permittivity matrix, μ¯¯ is the permeability matrix, and ξ¯¯ and ζ¯¯ are magnetoelectro-coupling matrices. In what follows, we will consider that this metasurface is reciprocal, implying that ϵ¯¯=ϵ¯¯T,μ¯¯=μ¯¯T,ζ¯¯=−ξ¯¯T,(2)where T is the transpose operation.83

To obtain the invariance conditions that apply to parameters ϵ¯¯, μ¯¯, ξ¯¯, and ζ¯¯ in Eq. (1), we must first understand how the electric and magnetic fields, E and H, transform under a given symmetry operation. For this purpose, consider the transformation matrix Λ¯¯, that corresponds to an arbitrary symmetry operation such as those described in Appendix. It follows that E and H, respectively, transform into E′ and H′ as51^,52E′=Λ¯¯·E,(3a)H′=1jωϵ0(Λ¯¯·∇)×(Λ¯¯·E)=det(Λ¯¯)Λ¯¯·H,(3b)where we have used the fact that the magnetic field H, being a pseudo-vector, transforms as the curl of electric field E. Note that all polar vectors would transform in the same way as E in Eq. (3a), whereas all pseudo-vector vectors transform as H in Eq. (3b). For the system in Eq. (1), we thus have D′=Λ¯¯·D,(4a)B′=det(Λ¯¯)Λ¯¯·B.(4b)

We next use Eqs. (3) and (4) to obtain the transformation relations that apply to the material parameters in Eq. (1). To do so, we reverse Eqs. (3) and (4) to express E, H, D, and B in terms of E′, H′, D′, and B′, respectively and then substitute the resulting relations into Eq. (1). By association, this readily yields51ϵ¯¯′=Λ¯¯·ϵ¯¯·Λ¯¯T,(5a)μ¯¯′=Λ¯¯·μ¯¯·Λ¯¯T,(5b)ξ¯¯′=det(Λ¯¯)Λ¯¯·ξ¯¯·Λ¯¯T,(5c)ζ¯¯′=det(Λ¯¯)Λ¯¯·ζ¯¯·Λ¯¯T,(5d)where we have used the fact that Λ¯¯ is an orthogonal matrix, implying that Λ¯¯−1=Λ¯¯T and that det(Λ¯¯)=det(Λ¯¯)−1, since det(Λ¯¯)=+1 for rotation symmetries and det(Λ¯¯)=−1 for reflection symmetries (refer to Appendix).

The relations in Eq. (5) represent how dipolar material parameters change under a particular spatial transformation defined by Λ¯¯. We shall see in Sec. 2.3 how such relations may be transformed into symmetry invariance conditions but, first, we shall investigate how these relations may be extended to an arbitrary multipolar order, which is the topic of the next section.

As we shall see in Sec. 4, some electromagnetic effects cannot be described solely using a purely dipolar model. This motivates the need to extend the dipolar framework discussed in the previous section to include higher-order multipolar components and their associated spatially dispersive responses. For this purpose, we combine concepts from the multipolar theory typically used for modeling metamaterials22^–24^,26^,27 along with concepts associated to spatial dispersion.29^–31 For the sake of conciseness, we shall next restrict ourselves to dipolar and quadrupolar responses, while higher-order multipolar responses may be considered in future works by following a procedure identical to the one discussed thereafter.

It follows that the quadrupolar constitutive relations read21^,31^,64D=ϵ0E+P+∇·Q¯¯,(6a)B=μ0(H+M+∇·S¯¯),(6b)where P and M are electric and magnetic polarization densities, and Q‾‾ and S‾‾ are irreducible (symmetric and traceless) electric and magnetic quadrupolar density tensors, respectively. The reason for considering irreducible tensors is that they provide a physically consistent description of the metasurface response.25 These quantities may be related to the fields via the spatially dispersive relations21^,64[PiMiQilSil]∝[χeeijχemijχee′ijkχem′ijkχmeijχmmijχme′ijkχmm′ijkQeeiljQemiljQee′iljkQem′iljkSmeiljSmmiljSme′iljkSmm′iljk]·[EjHj∇kEj∇kHj],(7)where χeeij, χemij, χmeij, and χmmij are related to the parameters in Eq. (1), whereas all the other terms are there to fully model the quadrupolar response of the metasurface. Note that reciprocity imposes relationships between parameters in Eq. (7), with the relevant reciprocity relations presented in Ref. 64.

As was done in Eq. (5), we may express the transformation relations that apply to the parameters in Eq. (7). To do so, we make use of tensor notation (in what follows, we omit the summations over repeated indices for convenience), since it applies more conveniently to the third- and fourth-order tensorial parameters in Eq. (7). It follows that an arbitrary tensor T‾‾, corresponding to any of the tensorial parameters of orders 1 to 4 in Eq. (7), transforms under the transformation Λ¯¯ as84Ti′=aΛijTj,(8a)Tij′=aΛimΛjkTmk,(8b)Tijk′=aΛilΛjmΛknTlmn,(8c)Tijkl′=aΛimΛjnΛkoΛlpTmnop,(8d)where a=det(Λ¯¯)n with n=0 for the “ee” or “mm” tensors in Eq. (7) and n=1 for the “em” or “me” tensors (in the case of polar vectors, n=0, whereas n=1 for pseudo-vectors), respectively. It is clear that Eqs. (8a) and (8b) are the tensor notation counterparts of Eqs. (3) and (5), respectively.

Extending this formalism to higher-order multipolar contributions is rather trivial. To do so, one first needs to include these additional contributions in Eq. (7) along with their corresponding spatial derivatives, which increases the total number of spatially dispersive material parameter tensors. For instance, taking into account octupolar contributions would require adding 20 new material tensors of orders 4 to 6 in Eq. (7), since second-order spatial derivatives must be included. Then, the relations in Eq. (8) would also need to be extended to apply to the additional material tensors of orders 5 and 6, which can be easily achieved by logical extrapolation.

We have just established that under a given symmetry operation, the material parameters in Eq. (7) transform according to the relations in Eq. (8). We are now interested in connecting the material parameters in Eq. (7) to the spatial symmetries of a metasurface. This may be achieved by considering that, according to Neumann’s principle, if a given structure is invariant under a symmetry operation, then its material parameters should also be invariant under the same operation.79^–81 Such an invariance condition may be mathematically expressed from Eq. (8) as Ti=aΛijTj,(9a)Tij=aΛimΛjkTmk,(9b)Tijk=aΛilΛjmΛknTlmn,(9c)Tijkl=aΛimΛjnΛkoΛlpTmnop,(9d)which implies that the tensor T‾‾ remains equal to itself after being transformed by Λ¯¯.

To obtain the material parameters that correspond to a given metasurface structure, we shall now describe a technique directly based on Eq. (9). This technique differs from those described in the literature, such as those that consist in expressing the material parameters for various symmetry groups51^,52 or those using the orthogonality theorem to find the irreducible representation of the structure.53^–56 Instead, we propose an approach that consists in recursively solving Eq. (9) for each symmetry operation for which the metasurface structure is invariant. By starting with a full material parameter tensor T‾‾ (one that contains all parameters), each iteration of this procedure leads to a system of equations formed by Eq. (9) that, when solved, reduces the complexity of T‾‾ by connecting some of its components or by setting others to zero. At the end of this process, one is left with a tensor that is precisely invariant under all symmetry operations that define the metasurface structure and thus properly models its effective response.

To illustrate this process, let us consider two different metasurfaces formed by periodically arranging in the x−y plane either the unit cell shown in Fig. 1(a) or the one shown in Fig. 1(b). These two unit cells are composed of the same scattering particle and only differ in the geometry of their lattice.

Inspecting the symmetries of these two metasurfaces, it is obvious that both are invariant under reflections through the x, y, and z axes, which corresponds to the symmetry operations Λ¯¯=σx,σy,σz (see Appendix). However, although they are composed of the same scattering particle, the metasurface unit cell in Fig. 1(a) possesses a square lattice, whereas the one in Fig. 1(b) possesses a rectangular lattice. It follows that these two metasurfaces do not exhibit the same rotation symmetry along the z axis. Indeed, a metasurface composed of a periodic arrangement of the unit cell in Fig. 1(a) would have a Λ¯¯=C4,z rotation symmetry [corresponding to 90 deg rotation symmetry along z; see Eq. (33)], whereas the one composed of the unit cell in Fig. 1(b) would have a Λ¯¯=C2,z rotation symmetry (corresponding to 180 deg rotation symmetry along z).

Now that we have found the symmetries corresponding to the two metasurfaces in Fig. 1, we can apply the process described above to find the permittivity matrix of these structures. Here, we restrict ourselves to finding ϵ¯¯ for simplicity but without loss of generality, since the process to obtain the other material parameters in Eq. (7) would be identical. We provide more complete examples that include all material parameters in Sec. 4.

The process of finding the permittivity matrix corresponding to the metasurfaces in Fig. 1 starts by considering the full permittivity matrix given as ϵ¯¯=[ϵxxϵxyϵxzϵxyϵyyϵyzϵxzϵyzϵzz],(10)where the reciprocity conditions in Eq. (2) have been considered. Next, Eq. (9b) is solved using Eq. (10) and Λ¯¯=σx with σx given in Eq. (31a), which leads to the simplified permittivity matrix, ϵ¯¯=[ϵxx000ϵyyϵyz0ϵyzϵzz].(11)

This process is then repeated by solving Eq. (9b) again, but this time with Eq. (11) and Λ¯¯=σy, and then once again with Λ¯¯=σz to obtain ϵ¯¯=[ϵxx000ϵyy000ϵzz].(12)

We are now left with the two rotation symmetries, C4,z and C2,z. However, it turns out that the permittivity matrix in Eq. (12) already corresponds to the structure in Fig. 1(b). Indeed, we do not even need to apply a further iteration of the process with Λ¯¯=C2,z because a structure possessing reflection symmetries along the x and y axes necessarily and automatically possesses a C2,z rotation symmetry. (However, the opposite is not necessarily true. Indeed, a structure possessing a C2,z symmetry, such as an S-shaped structure lying in the x−z plane, would not possess a reflection symmetry along the x axis.) This would therefore make a further application of the process with Λ¯¯=C2,z redundant. On the other hand, the rotation symmetry C4,z is not redundant and applying Eq. (19) with Λ¯¯=C4,z and Eq. (9b) leads to ϵ¯¯=[ϵxx000ϵxx000ϵzz],(13)which is the permittivity matrix that corresponds to the metasurface in Fig. 1(a).

Obviously, for such simple structures as those in Fig. 1, the results obtained in Eqs. (12) and (13) could have been guessed simply based on intuition. However, for more complex structures, intuition alone is usually not sufficient to guess the correct shape of the material parameter tensors, especially those related to quadrupolar responses in Eq. (7). The usefulness of the approach in such cases is emphasized using examples found later in Sec. 4.

Note that the background medium is assumed to be identical on both sides of the metasurfaces in Fig. 1. However, in practice, the scattering particles are usually deposited on top of a substrate instead of being embedded into a uniform background medium. In this case, the presence of the substrate would actually break the symmetry of the system in the z direction, meaning that the metasurface response would not be invariant under Λ¯¯=σz. While this would not change the shape of the permittivity matrices in Eqs. (12) and (13), it would have significant influence on the other material parameters in Eq. (7). For instance, it has been shown that the presence of a substrate is sufficient to lead to 3D chiral responses for metasurfaces composed of only 2D chiral scattering particles such as Gammadion structures.41

In the previous sections, we have seen how the effective material tensors of a metasurface may be predicted by considering the spatial symmetries of the metasurface unit cell. We shall now investigate the relationships between the scattering properties of a metasurface and its spatial symmetries.

It turns out that the spatial symmetries of the metasurface unit cell are not the only parameters to consider when investigating the metasurface scattering response. Indeed, one needs to also take into account the influence of the waves interacting with the metasurface. This is because the combined system composed of the superposition of the metasurface and the incident and scattered waves exhibits spatial symmetries that are not the same as those of the metasurface alone.70^–76 Such considerations have already been well discussed in the literature, specifically in the pioneering works by Dmitriev for the case of normally incident waves72 and oblique propagation.73

In what follows, we shall review the general concepts discussed in Refs. 72 and 73 for the case of reciprocal metasurfaces and develop a procedure to obtain the scattering matrix of a metasurface for a given direction of propagation that is similar to the procedure described in Sec. 2.3.

Let us first consider the case of a reciprocal metasurface illuminated by obliquely propagating plane waves. Since the metasurface is considered to be uniform and its lattice period is deeply subwavelength, it scatters the incident waves according to Snell’s law. Such a situation is depicted in Fig. 2, where two input and output “ports” transmit and receive, respectively, TE- and TM-polarized waves. Here, the plane of incidence is arbitrarily chosen to be the x−z plane. The scattering matrix associated with the interactions in Fig. 2 is given as72^,73S¯¯=[R11R12T13T14R21R22T23T24T31T32R33R34T41T42R43R44],(14)where Rab and Tab are reflection and transmission coefficients, respectively, with a,b=1,2,3,4 corresponding to the polarization state of the different waves where even and odd numbers are associated to TE and TM polarizations, respectively.

We are now looking for an invariance condition that would apply to S‾‾ in a similar way that the invariance conditions in Eq. (9) applied to the material tensor T‾‾. For this purpose, we need to relate S‾‾ to the symmetries of the combined system (metasurface and direction of wave propagation) defined by Λ¯¯. From that prospect, we start by defining a rotated version of Λ¯¯ that is aligned with the orientation of the plane of incidence, since the latter is generally not restricted to the x−z plane. This rotated version of Λ¯¯ is given as Λ¯¯′=Rz(ϕ)·Λ¯¯·Rz(ϕ)T,(15)where the rotation matrix Rz is defined in Eq. (32c) and ϕ represents the angle between the plane of incidence and the x−z plane. Consequently, if the plane of incidence coincides with the x−z plane, then ϕ=0deg or ϕ=180deg, whereas if it is aligned with the y−z plane, then ϕ=90deg or ϕ=270deg.

Now, since we are considering a given plane of incidence, it follows that only the spatial symmetries that map an input/output port to another one are to be considered for obtaining a symmetry-invariant scattering matrix.73 The other symmetry operations, such as the C4,z rotation symmetry, should not be considered, since they would be inconsistent with a fixed plane of incidence. Therefore, the only spatial symmetries that make sense to be taken into account are Λ¯¯=σi and Λ¯¯=C2,i, where i=x,y,z as well as the inversion symmetry operation defined by Λ¯¯=−I¯¯.

Among this set of relevant symmetry operations, we may further classify them according to their effects on the polarization and the direction of wave propagation of the incident and scattered waves. For instance, applying the symmetry operation Λ¯¯=σx on the system in Fig. 2 would flip the sign of the TM polarizations while leaving the TE polarizations unaffected, and it would also flip the direction of the propagation vectors k. On the other hand, the operation Λ¯¯=σy would flip the sign of the TE polarizations while leaving both the TM polarizations and the k-vectors unchanged.

Finally, we need to define a transformation operation that would apply to S‾‾ in the same way that the tensor T‾‾ was transformed by the operation Λ¯¯ in Eq. (8). This may be achieved by defining the transformed scattering matrix, S¯¯′, as S¯¯′=M¯¯·S¯¯, where M‾‾ is a transformation matrix that is related to Λ¯¯.73 Note that in Ref. 73 a specific transformation matrix M‾‾ was defined for each of the relevant symmetry operations specified previously, based on the intuitive effects that these operations have on the polarizations and propagation directions of the waves. In what follows, we shall instead provide a technique that directly connects M‾‾ with Λ¯¯.

Noting that the invariance condition may be expressed as S‾‾=M‾‾·S‾‾·M‾‾−1 for the symmetry operations that do not affect the direction of the k vectors, whereas it is given by S‾‾=M‾‾·S‾‾T·M‾‾−1 for the symmetry operations that flip the k vectors,73 we define the general invariance condition as S¯¯=M¯¯·[bS¯¯+(1−b)S¯¯T]·M¯¯−1,(16)where b=1−c2(1−Λxx′), with Λxx′=x·Λ¯¯′·x and c=c1c2 and c1 and c2 given as c1={1,ifΛxz′=Λyz′=0,0,otherwise,(17a)c2={1,ifΛ¯¯′·y=±R¯¯z(ϕ)T·R¯¯z(ϕ)·y,0,otherwise.(17b)

The transformation matrix M‾‾ in Eq. (16) is now obtained by generalizing the results provided in Ref. 73, which leads to M¯¯=[a¯¯+a¯¯−a¯¯−a¯¯+],(18)where the submatrices a‾‾+ and a¯¯− are defined as a¯¯+=c2(1+Λzz′)[Λyy′00Λxx′]+(1−c)I¯¯t,(19a)a¯¯−=c2(1−Λzz′)[Λyy′00Λxx′],(19b)where I‾‾t is a two-dimensional identity matrix.

The procedure to obtain the scattering matrix of a given metasurface under oblique incidence is quite similar to the approach described in Sec. 2.3 for finding the metasurface material parameter tensors. It consists in first defining the orientation of the plane of incidence with ϕ, and then recursively solving Eq. (16) for each of the relevant symmetry operations expressed using Eq. (15).

For the special case of normally incident plane waves impinging on the metasurface, it is required to slightly modify this procedure.72 First, at normal incidence, the orientation of the plane of incidence loses its meaning and Eq. (15) should be bypassed by setting ϕ=0deg. Then, the TE and TM polarizations should now be associated with x and y polarizations, respectively. Finally, it now actually makes sense to consider the C4,z rotation symmetry as it correctly maps TE to TM polarizations together. For this specific symmetry operation, the invariance condition Eq. (16) reduces to S‾‾=M‾‾·S‾‾·M‾‾−1 with M‾‾ given as72M¯¯C4,z=[0−1001000000−10010].(20)

For all the other symmetry operations, relations Eqs. (16)–(19) should be used.

Before looking at examples illustrating the application of the method described above, which will be presented in the next section, we would like to comment on how to deduce the polarization effects induced by a metasurface directly from its scattering matrix. Considering the general scattering matrix [Eq. (14)], we see that it may be reduced to four internal submatrices composed of two sets of reflection and transmission matrices. Each one of these submatrices may be associated with a Jones matrix, since they describe how different polarizations are scattered by the metasurface.85 From a given Jones matrix, it is then possible to deduce the general effect that the metasurface has on different polarization states.75^–77^,86 For completeness, we provide in Table 1 the relationships between some common Jones matrices and their corresponding polarization effects.

We shall now look at three examples that illustrate the application of the techniques described in Sec. 2 and in Sec. 3. In these examples, we will represent the material parameter tensors in way that is reminiscent of how they are depicted in Ref. 34, meaning that we will make use of the fact that the quadrupolar tensors in Eq. (6) are assumed to be irreducible tensors (symmetric and traceless), which allows us to greatly reduce the number of independent components in Eq. (7). This also means that the derivative operators in Eq. (7) can be simplified such that, for instance, the electric field gradient becomes34$$\begin{gathered} ⋄E={(∂xEy+∂yEx)/2,(∂xEz+∂zEx)/2, \\ (∂yEz+∂zEy)/2,∂xEx,∂yEy,∂zEz}, \end{gathered}$$(21)and similarly for the magnetic field gradient.

Extrinsic chirality corresponds to an electromagnetic effect in which a medium (in our case a metasurface) exhibits a chiral response even though it is not composed of geometrically 3D chiral scattering particles.43^–46 This chiral response emerges only when the metasurface is illuminated along a specific direction, hence the reason why it is referred to as “extrinsic”: it depends on the direction of wave propagation.

Consider the metasurface unit cell depicted in Fig. 3 composed of a T-shaped metallic particle embedded in a uniform background medium. The metasurface is illuminated by an obliquely incident plane wave propagating either in the x−z plane, as in Fig. 3(a), or in the y−z -plane, as in Fig. 3(b).

The spatial symmetries of the metasurface, assuming that it has a square lattice, are directly found to be σx and σz. Note that it also exhibits a C2,y rotation symmetry; however, it is redundant, since we have already considered its σx reflection symmetry, as explained in Sec. 2.3. Applying the procedure outlined in Sec. 2.3, we find the multipolar components of this metasurface and present them in Fig. 4, using Eq. (21), in way that is similar to the representations in Ref. 34.

In Fig. 4, the vertical and the horizontal axes, respectively, correspond to the multipolar densities and the excitation vector, on the left- and right-hand side of Eq. (7), where the term ⋄xye corresponds to the first element in Eq. (21), and so on. The numbers and the colors in the figure are purely arbitrary and only help visualize which components are allowed to exist due to the symmetries of the structure and which components are related to each other (components that have the same number and color).

The information provided in Fig. 4 reveals the complexity of the metasurface multipolar response, notably its bianisotropic dipolar nature due to the existence of the components numbered 4 and 5 in the matrix. Comparing the results in Fig. 4 with those of the dolmen structure provided in Ref. 34 reveals a perfect match between the two methods. It should be noted that results in Fig. 4 only correspond to the components that are allowed to exist, and it does not provide any information about the magnitude of these components. On the other hand, the results presented in Ref. 34 provide the numerical value of the multipolar components, since the method used to obtain them consists in numerically simulating the structure with a complex set of different illumination conditions, which results in an important computational cost.

Now, we use the approach described in Sec. 2 to compute the scattering matrix for the obliquely propagating waves of Fig. 3. Using Eqs. (15)–(19) with Λ¯¯=σx, Λ¯¯=σz, and ϕ=0deg yields the scattering matrix, for the case in Fig. 3(a), given as S¯¯=[R11R12T13T14−R12R22−T14T24T13T14R11R12−T14T24−R12R22].(22)

Similarly, using ϕ=90deg, we obtain the scattering matrix for the case in Fig. 3(b) as S¯¯=[R110T1300R220T24T130R1100T240R22].(23)

Comparing Eqs. (22) and (23), it is clear that illuminating the metasurface along different directions of wave propagation leads to quite different polarization effects. When the metasurface is illuminated in the x−z plane, as in Fig. 3(a), we see that the internal Jones matrices in Eq. (22) exhibit circular birefringent properties typical of chiral responses (see Table 1). On the other hand, illuminating the metasurface in the y−z plane, as in Fig. 3(b), leads to a linear birefringent response, as evidenced by the scattering matrix in Eq. (23). This analysis confirms that the metasurface is extrinsically chiral, since its chiral response is dependent on the direction of wave propagation.

One way to understand the emergence of a chiral response in this structure is to consider its bianisotropic dipolar response. For an illumination in the x−z plane, both TE and TM polarizations excite the components 4 and 5 in Fig. 4, which correspond to the components χemxz and χemzx in Eq. (7), respectively, or to their reciprocal counterparts χmezx and χmexz. Since these two components have different values that thus cannot cancel each other, they lead to a net chiral response.10 However, for an illumination in the y−z plane, both TE and TM polarizations excite either the set of parameters χemxz and χmezx or the set χemzx and χmexz, respectively. Since by reciprocity [see Eq. (2)], χemxz=−χmezx and χemzx=−χmexz, these two contributions cancel each other, leading to an absence of chiral response.

We are now interested in designing a metasurface that exhibits an asymmetric angular transmittance when illuminated in the x−z plane. While it has been shown that asymmetric angular transmittance may be achieved using spatially varying metasurfaces with phase-gradient modulations,87 we are instead interested in designing a uniform metasurface whose asymmetric scattering response stems from the asymmetry of its structure.

Let us consider the two metasurfaces depicted in Fig. 5, composed of a periodic array of L-shaped scattering particles. In Fig. 5(a), the structures lie flat on the plane of the surface, whereas they stand vertically in Fig. 5(b).

Since both structures are asymmetric in the x direction, one may a priori expect that they would both lead to an asymmetric scattering response for waves propagating in the x−z plane in opposite directions, like the waves Ψ1 and Ψ2 in Fig. 5. However, that is not the case, as we shall next demonstrate.

The spatial symmetries of the structure in Fig. 5(a) are σz and σxy, whereas those of the structure in Fig. 5(b) are σy and σxz, where σxy and σxz refer to symmetries with respect to the 45 deg diagonal between the x and y axes or between the x and z axes, respectively. These symmetries may be defined by combining those provided in the Appendix so that σxy=Rz(−45deg)·σx·Rz(45deg) and σxz=Ry(−45deg)·σx·Ry(45deg). Applying the approach discussed in Sec. 2.3 to these two sets of symmetries yields their corresponding material parameter tensors, which are provided in Fig. 6. As can be seen, both structures are bianisotropic and, if one only considers their permittivity matrix, we see that ϵxx=ϵyy and ϵxy≠0 for the flat L-shaped structures of Fig. 5(a), whereas ϵxx=ϵzz and ϵxz≠0 for the vertical L-shaped structures of Fig. 5(b), as may be expected intuitively.

We next compute their respective scattering matrices. For the case in Fig. 5(a), and the waves Ψ1 and Ψ2 that propagate in the x−z plane (ϕ=0deg⁡), we have S¯¯=[R11R12T13T14−R12R22−T14T24T13T14R11R12−T14T24−R12R22].(24)

Interestingly, we see that this structure exhibits a form of extrinsic chirality similar to the one in Fig. 3(a). However, this metasurface does not possess an asymmetric angular transmittance, since the two transmission submatrices in Eq. (24) are equal to each other, implying that Ψ1 and Ψ2 transmit through the metasurface identically.

Note that if the plane of incidence coincides with the in-plane axis of symmetry along which σxy is satisfied, as it is the case for the wave Ψ3 in Fig. 5(a), the extrinsic chiral response disappears, as evidenced by the scattering matrix (ϕ=45deg), S¯¯=[R110T1300R220T24T130R1100T240R22].(25)

For the metasurface in Fig. 5(b), and the waves Ψ1 and Ψ2 that propagate in the x−z plane (ϕ=0deg), we have S¯¯=[R110T1300R220T24T310R3300T420R44].(26)

This time, not only is there no cross-polarized scattering, but the two transmission submatrices are different from each other (T13≠T31 and T24≠T42), implying that the waves Ψ1 and Ψ2 interact differently with the metasurface. This indicates that asymmetric angular transmittance is possible in this situation. For the incident wave Ψ3 that propagates in the y−z plane, there is, however, no co-polarized angular transmittance asymmetry, as evidenced by its scattering matrix (ϕ=90deg), S¯¯=[R11R12T13T14−R12R22T23T24T13−T23R33R34−T14T24−R34R44].(27)

Nevertheless, illuminating the metasurface under this direction of propagation should induce a form of asymmetric transmittance when illuminated in the +z or the −z directions. This therefore does not correspond to the sought-after asymmetric angular transmittance effect.

The last example investigates the scattering response of the two metasurfaces depicted in Fig. 7 under normal illumination. The metasurface in Fig. 7(a) consists of a simple periodic array of square-shaped scattering particles, whereas the one in Fig. 7(b) is made of Gammadion structures. In both cases, the metasurface is embedded within a uniform medium.

While it is intuitive to guess that a normally incident plane wave impinging on the metasurface in Fig. 7(a) would not undergo polarization conversion, it is less trivial to intuitively predict the response of the metasurface in Fig. 7(b). Referring to the literature, it was for instance suggested in Ref. 41 that such an ideal Gammadion array should not exhibit a chiral response unless the scattering particles are asymmetric in the z direction (or the metasurface lies on a substrate). In what follows, we shall see that the problem is in fact more complex than it appears.

Let us first compute the material parameter tensors of these two metasurfaces. The symmetries associated with the metasurface in Fig. 7(a) are σx, σy, σz, and C4,z, whereas those for the metasurface in Fig. 7(b) are only σz and C4,z. The corresponding material parameter tensors are given in Fig. 7. Inspecting the differences between Figs. 8(a) and 8(b) reveals a striking feature: they both exhibit identical purely dipolar responses, i.e., the four 9×9 submatrices relating p and m to E and H are the same for both metasurfaces. This means that a modeling approach solely based on Eq. (1) would fail to predict a difference in the scattering response of these metasurfaces. It is only by considering quadrupolar responses and higher-order spatially dispersive effects that the differences between these metasurfaces become evident.

Following the approach described in Sec. 3 for normally incident waves, the scattering matrix of the metasurface in Fig. 7(a) is given as S¯¯=[R110T1300R110T13T130R1100T130R11],(28)whereas the one for the metasurface in Fig. 7(b) reads S¯¯=[R11R12T13T14−R12R11−T14T13T13T14R11R12−T14T13−R12R11].(29)

As may be expected, the metasurface in Fig. 7(a) is isotropic and does not induce any polarization rotation or conversion at normal incidence. On the other hand, the shape of the scattering matrix of Eq. (29) reveals that the metasurface in Fig. 7(b) should exhibit a chiral response. This may come as a surprise, considering that the metasurface is not made of geometrically 3D chiral structures and that it is illuminated at normal incidence. However, a detailed inspection of the material parameter tensors in Fig. 8(b) helps us understand why such a chiral response exists.

Consider an x-polarized normally incident plane wave impinging on the metasurface in Fig. 7(b) with an electric field defined by E=x^ei(ωt−kz) and a magnetic field given by H=y^ei(ωt−kz)/η, with k and η being the wavenumber and impedance of the background medium. This wave excites the components 5 and 19 in Fig. 8(b) via the field derivatives ⋄xze and ⋄yzh, which, in this case, correspond to ∂zEx and ∂zHy, respectively. Since these spatial derivatives are not zero for the considered excitation, it follows that these two material components, which do not exist for the metasurface in Fig. 7(a), respectively induce an x-polarized magnetization of the metasurface, mx, and a y-polarized polarization, py, suggesting rotation of polarization. By reciprocity,64 the components 5 and 19 also induce the quadrupolar responses Qyz and Sxz that are excited via the field components Hy and Ex, respectively, and that also contribute to rotation of polarization.

In order to verify that the multipolar components py, mx, Qyz, and Sxz are indeed excited when illuminating such a structure with an x-polarized normally propagating plane wave, we have numerically simulated an isolated Gammadion particle. From its scattered fields, we have extracted its spherical multipolar components following the approach provided in Ref. 88. The resulting simulations are plotted in Fig. 9, where the solid and dashed lines correspond to the co- and cross-polarized multipolar components, respectively. It is thus clear that a cross-polarized response is achieved, although it is small compared to the co-polarized one. This confirms that the metasurface in Fig. 7(b) exhibits a chiral response, as suggested by its scattering matrix, Eq. (29).

This type of chiral response slightly differs from the one discussed in Sec. 4.1, where the chirality was due to an obliquely propagating wave interacting with the bianisotropic effective material parameters of the structure. In the case of the metasurface in Fig. 7(b), the chiral response is due to multipolar and spatially dispersive components excited by a normally propagating plane, thus indicating the presence of multipolar extrinsic chirality. The extrinsic nature of this chiral response is here not directly related to the direction of wave propagation, as it was the case in Sec. 4.1, but rather to the gradient of the fields.

We have established a relationship between the spatial symmetries of a metasurface and its corresponding material parameter tensors and scattering matrix. This relationship has been obtained based on a simple approach that consists in the recursive application of invariance conditions, instead of relying on complicated concepts pertaining to group theory. This makes this approach versatile and accessible to a large audience. Moreover, to facilitate the application of the proposed approach, we have implemented a Python script that conveniently computes the form of the material parameter tensors and the scattering matrix directly from a list of specified symmetries. This Python script is freely accessible on GitHub.78

Based on this approach, we have shown how easy it is to investigate responses such as extrinsic chirality or asymmetric angular transmittance. We have also demonstrated the possibility of multipolar extrinsic chirality, where chiral responses may be achieved in achiral structures even for normally incident waves due to the excitation of multipolar components; this was shown to result from the sensitivity of certain structures to field gradients. These examples demonstrate that intuition alone is insufficient to predict the effective material tensors and rich scattering behavior that can be achieved, and thus highlight the usefulness of our method.

Reflection symmetries may be expressed using the Householder equation89Λ¯¯=I¯¯−2nn,(30)where n corresponds to the reflection axis. It follows that the symmetry operations σx, σy, and σz are defined as σx=[−100010001],(31a)σy=[1000−10001],(31b)σz=[10001000−1].(31c)

On the other hand, the 3D rotation matrices are defined as Rx(θ)=[1000cosθ−sinθ0sinθcosθ],(32a)Ry(θ)=[cosθ0sinθ000−sinθ0cosθ],(32b)Rz(θ)=[cosθ−sinθ0sinθcosθ0001].(32c)

These matrices may be used to define the symmetry operation CN,i as CN,i=Ri(2πN).(33)

Karim Achouri received his BSc and his MSc degrees in photonics from Swiss Federal Institute of Technology in Lausanne (EPFL), in 2010 and 2012, respectively. In 2017, he received his PhD in electromagnetics from Ecole Polytechnique de Montréal. In 2017, he joined the Nanophotonics and Metrology Laboratory at EPFL as a postdoctoral fellow. His research interests are in the areas of electromagnetics, metamaterials, applied optics, photonics, and nanotechnologies.

Ville Tiukuvaara received his BEng. degree in engineering physics and his MASc. from Carleton University, Ottawa, Canada, in 2019 and 2021, respectively. Currently, he is pursuing a PhD in photonics at Swiss Federal Institute of Technology in Lausanne (EPFL), Switzerland, with a research interest in metamaterials, along with analytical and computational electromagnetics.

Olivier J. F. Martin conducts a comprehensive research that combines the development of numerical techniques, advanced nanofabrication, and experiments on plasmonic systems. Applications of his research include optical antennas, metasurfaces, nonlinear optics, optical nano-manipulations, heterogeneous catalysis, security features, and optical forces at the nanoscale. He has authored over 300 journal articles and holds a handful of patents and invention disclosures.