Synthesis and Applications
of Birefringent Metasurfaces
Birefringent metasurfaces are two-dimensional structures capable of independently controlling the amplitude, phase and polarization of orthogonally polarized incident waves. In this work, we propose a in-depth discussion on the mathematical synthesis of such metasurfaces. We compare two methods, one that is rigorous and based on the exact electromagnetic fields involved in the transformation and one that is based on approximate reflection and transmission coefficients. We next validate the synthesis technique in metasurfaces performing the operations of half- and quarter-wave plates, polarization beam splitting and orbital angular momentum multiplexing.
Birefringence, also called double refraction, is the physical property of an anisotropic medium to exhibit an angle dependent refractive index . This phenomenon, first observed in crystals more than 300 years ago , has already lead to the realization of several major optical components and applications. More recently, metasurfaces [3, 4, 5, 6, 7, 8], the two dimensional counterparts of metamaterials, have seen an important rise of interest due to their rich potential in the transformation of electromagnetic fields. Combined with the orthogonality property of and polarized waves, birefringent structures, and especially metasurfaces, have the capability to independently control the amplitude, phase and polarization of two orthogonal electromagnetic waves, leading to a wealth of possible applications at optical and microwave frequencies.
In this work, we propose an in-depth discussion on the synthesis of birefringent metasurfaces. This discussion is based on the general bianisotropic metasurface synthesis technique developed in . We compare two different synthesis methods, one that is rigorous and one that is based on paraxial approximate transmission and reflection coefficients. We also report the implementation details of several metasurfaces that were only briefly presented in [9, 10]. These metasurfaces perform the operations of half-wave plate [11, 12, 4], quarter-wave plate [13, 14], polarization beam splitting [15, 16, 17] and orbital angular momentum [18, 19, 20, 21, 22] multiplexing.
The paper is organized as follows. In the second section, the mathematical synthesis as well as the physical realization of metasurfaces are addressed. Two different approaches for the mathematical synthesis are discussed and compared. In the third section, we present the implementation of the four metasurfaces introduced above.
Ii Metasurface Design
Ii-a Mathematical Synthesis
The metasurface synthesis technique in  stems from the continuity relations, initially derived by Idemen , which express the discontinuities of the electromagnetic fields in the presence of a spatial discontinuity, such as a metasurface. In the simple case of a monoanisotropic (zero magnetoelectric coupling coefficients) birefringent metasurface given in terms of its transverse diagonal susceptibility tensors, and , the continuity relations read
where it is assumed that the metasurface lies in the plane at . In (1), the differences of the electric and magnetic fields between both sides of the metasurface (expressed by the operator ) are related to the metasurface susceptibilities excited by the average electric and magnetic fields (denoted by the subscript ‘av’). The system of (1) can be easily solved to express the susceptibilities in terms of the specified fields. Due to the orthogonality between and polarizations, the solutions split into two independent sets which are respectively given by
At this stage, the metasurface is completely defined by the susceptibilities in (2) and performs the required transformation between the incident, reflected and transmitted waves [8, 10]. The birefringent operation leverages the property of orthogonality between (2a) and (2b).
In order to realize the metasurface, one has to find the appropriate shapes of the scattering particles. Therefore, it is convenient to find the transmission and reflection coefficients of the metasurface which could then be related to the transmission and reflection coefficients of each scattering particle obtained via full-wave simulations. To obtain these coefficients, it is assumed that the metasurface is illuminated by a normally incident plane wave and that it reflects and transmits normally propagating plane waves (either or polarized). The corresponding electric and magnetic fields are thus, respectively, given by and , and and . These fields are then substituted into relations (2), which are then solved for the coefficients and , which read 
where is the free space wave number. Using the susceptibilities from (2a) or (2b) will yield the transmission and reflection coefficients for or polarizations, respectively. Since a monoanisotropic metasurface is necessarily symmetric with respect to the direction , the reflection coefficients of our metasurface are the same from both sides of the structure.
Relations (3) can be reversed to express the susceptibilities in terms of the transmission and reflection coefficients as
where and are found assuming that and are the coefficients of -polarized waves. Alternatively, and are found when and are the coefficients of -polarized waves. Although Eqs. (2) can be rigorously used to synthesize metasurfaces, relations (4) suggest an alternative synthesis technique which would consists in specifying the required transformation in terms of transmission and reflection coefficients instead of the tangential electromagnetic fields, as would be done in (2).
In what follows, we will compute the responses of metasurfaces synthesized using the methods based on relations (2) and (4) and compare them. Let us consider a reflection-less metasurface that refracts at a normally incident -polarized plane wave. The electromagnetic fields of the incident wave are, at , given by , while the fields of the transmitted wave are . The first metasurface synthesis method consists in substituting these fields into (2a), which results in the following susceptibilities
The second synthesis method, based on relations (4), seems a priori unsuitable for such a kind of transformation (i.e. refraction) since relations (3) and (4) were obtained assuming that all waves impinging on or scattered by the metasurface are propagating normally to the structure, which is obviously contradictory with the concept of refraction. Indeed, this second synthesis technique rigorously applies only to normally propagating waves, but it may also be used as an approximation to synthesize refractive metasurfaces in the case of small refraction angles, i.e. paraxial approximation, as will be shown next. In fact, this metasurface synthesis technique allows one to obtain the material properties, through relations (4), from the transmission and reflection coefficients that would be initially defined using the complex amplitude transmittance method  (in the case of zero reflection) or, more generally, the momentum transformation technique introduced in .
respectively. Because this metasurface is assumed to be reflection-less, the complex amplitude transmittance method can be used to define the parameter in (7), that is to say
where and are the phase profiles, projected on the metasurface plane, of the incident and the transmitted waves, respectively. Since the incident wave is normally impinging on the metasurface, we have that . The transmission coefficient is then simply defined by , which transforms (7) to
Now, let us compare the susceptibilities in (5), obtained from the first synthesis method, which is rigorous, and the susceptibilities in (9), derived from the second method, which is approximate. The real and imaginary parts of the electric and magnetic susceptibilities are plotted at the top and bottom of Fig. 1, respectively. The left-hand side of the figure corresponds to the first method while the right-hand side corresponds to the second method. They are two main differences between the two methods. Firstly, the electric and magnetic susceptibilities are different in (5) whereas they are equal in (9). Secondly, the susceptibilities in (5) are complex with negative imaginary parts corresponding to absorption whereas the susceptibilities in (9) are purely real and thus correspond to a lossless and passive structure.
The reason for the imaginary parts in (5), and the resulting absorption, is the unequal normal power flow between the incident wave and the transmitted wave [26, 24]. Indeed, the transmitted wave, which propagates under , has a lower normal transmitted power than the normally incident wave. This translates into a reduced transmission efficiency where the excess energy of the incident wave is absorbed by the metasurface.
Another interesting comparison to establish between the two methods is the differences of the transmission and reflection coefficients in (3) and (8). Substituting the susceptibilities found in (5) into (3) yields the transmission and reflection coefficients, for the first synthesis method, that are plotted on the left-hand side of Fig. 2 in solid blue lines and in dashed red lines, respectively. The transmission and reflection coefficients of the second synthesis method are simply and , and are plotted on the right-hand side of Fig. 2.
It appears that the metasurface synthesized with the first method may be seen as an equivalent amplitude and phase grating in transmission and reflection, while the other metasurface is a simple transmission phase gradient structure. The non-zero reflection coefficient that is plotted in Fig. 2 (top-left) seems, a priori, contradictory with the prescription of zero reflection specified to obtain relations (5). In fact, no propagating reflected wave is produced by the metasurface because the -vector of the reflection phase (bottom-left in Fig. 2), defined as , where is the period of the reflection phase, is larger than the free space wave number, . This means that the reflected wave is an evanescent wave and thus does not propagate. Moreover, the non ideal transmission efficiency that was discussed above and that is responsible for the imaginary parts of the susceptibilities in (5) also contributes to the non-zero reflection coefficient in Fig. 2 and the loss evidenced by the fact that .
Finally, let us see how these two synthesis techniques compare by performing full-wave simulations. We have made two slightly different kinds of simulations. One using an home made Finite-Difference Frequency Domain (FDFD) code  that simulates an exactly zero-thickness metasurface. The other simulations were made using COMSOL by assuming a thin metasurface of thickness . Note that, in COMSOL, the susceptibilities were introduced by assigning relative permittivity and permeability to a slab which are respectively defined by and where the division by dilutes the effect of the susceptibilities over the thickness of the slab . Simulation results showing the real part of are plotted in Fig. 3, where the top and bottom rows correspond to the FDFD simulations and COMSOL simulations, respectively. The left- and the right-hand sides correspond to the first and second synthesis techniques, respectively.
As can be seen, the best result (Fig. 3 top left) was obtained with the FDFD code and using the rigorous relations in (2) (first synthesis technique). This is not surprising because, in this simulation, the metasurface has exactly zero-thickness and, therefore, it can be rigorously described by the continuity equations (1). The diffraction efficiency, defined as the transmitted power density at divided by the incident power density, is about for the first method and for the second one which, surprisingly, works fairly well except for undesired reflection. The simulations with COMSOL give worse results due to the thickness of the metasurface. For both synthesis techniques, several diffraction orders appear, either in reflection or in transmission. Moreover, due to the thickness, some modes are trapped (guided modes) inside the metasurface, which contributes to further lower the diffraction efficiency into the desired direction to and for the first and second methods, respectively. The rigorous method gives, in that case, worse results than the approximate method which might be explained by the fact that the susceptibilities in (2) are lossy and thus the wave is more attenuated by propagating through the thickness of the metasurface whereas the metasurface obtained with the second method is not lossy at all.
To conclude this section, it must be noted that, while the first synthesis technique is the most rigorous one and gives the best results in FDFD simulation, it remains much more complicated to implement than the second method. This is because the physical realization of these metasurfaces necessarily requires a mapping between the susceptibilities, given in Fig. 1, and the scattering parameters, given in Fig. 2. And, as can be seen from the scattering parameters, the realization of the metasurface synthesized with the first method would require implementing non-uniform reflection and transmission coefficients that present different phase gradients, moreover this metasurface would also be lossy. Compare this now to the metasurface synthesized with the second method and that presents a uniform unity transmission coefficient and a phase gradient, it is clear that this second method is much easier to realize and considering the excellent results in FDFD simulation, it is the one that is usually preferred for the realization of most metasurfaces. The metasurfaces presented in the following section of this paper are synthesized based on the second synthesis technique.
Ii-B Physical Realization
In order to fabricate the metasurfaces, the required scattering parameters are discretized with subwavelength resolution. At each lattice site, a scattering particle (or unit cell) is realized such that it exhibits the required scattering behavior. The unit cells are simulated in a commercial software and assuming periodic boundary conditions. The shape of the unit cells are optimized such that the scattering matrices obtained by simulation correspond to the transmission and reflection coefficients given in (3).
To implement each unit cell, we have used a cascade of three metallic layers (with identical outer layers) held together by two dielectric spacers. This type of unit cells has been shown to present full transmission (assuming lossless material) and a complete -phase coverage .
Each metallic layer of the unit cell consists in a Jerusalem cross, as shown in Fig. 4 with all its variable dimensions. All the metasurfaces discussed thereafter have been realized with unit cells of size mm which corresponds to at 10 GHz. The dielectric substrates used are Rogers RO3003 () with a thickness of 1.52 mm for each spacer leading to a total metasurface thickness of 3.04 mm ().
In this section, several types of birefringent metasurfaces are discussed and demonstrated experimentally. The realized metasurface are, in order of appearance, a half-wave plate, a quarter-wave plate, a polarization beam splitter and an orbital angular momentum generator. These four operations are illustrated in Fig. 5.
Iii-a Electromagnetic Wave Plates
Electromagnetic wave plates are birefringent structures that exhibit specific transmission phase difference between - and -polarizations defined as . Here, we present the two most common wave plates: a half-wave plate and a quarter-wave plate which, respectively, correspond to and . The half-wave plate performs a rotation of polarization for linear polarization or changes the handedness of circular polarization. The quarter-wave plate transforms linear polarization into circular polarization and vice-versa.
For such electromagnetic transformations, the metasurface is uniform since there is no variation in the direction of propagation of the waves and, therefore, no phase gradient is required. This make these metasurfaces very easy to design since only one unit cell has to be implemented and repeated periodically to form the metasurface.
Iii-A1 Half-Wave Plate
The fabricated metasurface is shown in Fig. 6 while the dimensions of its unique unit cell are given in Table I. The metasurface is made of unit cells corresponding to an aperture of about 5. Note that the two holes on both sides of the metasurface are used to attach it to the measurement setup.
The metasurface has then been measured. Two horn antennas, placed on both sides of the metasurface, have been used to measure the normal transmission from a normally incident wave. The measured transmissions for and polarizations are, respectively, plotted in Figs. (a)a and (b)b, where the red solid lines correspond to measurements with the metasurface and the blue dashed lines are the reference line of sight measurements of the horn antennas. The phase difference, , is plotted in Fig (c)c.
As can be seen, the metasurface transmission is almost unity around the specified frequency of operation of 10 GHz. The ideal phase shift difference of is obtained at GHz. In order to verify the rotation of polarization capability of this structure, cross-polarization measurements with and without the metasurface have been performed and the result is shown in Fig. 8 only for the metasurface transmission case for convenience.
The result in Fig. 8 confirms that the metasurface behaves almost as a perfect half-wave plate with a power transmission efficiency of at 10.2 GHz and a -3-dB bandwidth of about .
Iii-A2 Quarter-Wave Plate
The quarter-wave plate metasurface was designed and realized following exactly the same procedure as that of the half-wave plate metasurface. The fabricated metasurface is shown in Fig. 9 and the dimensions of its unit cell are given in Table II.
As was done for the half-wave plate metasurface, the measurements of the quarter-wave plate metasurface, corresponding to polarization transmission, polarization transmission and phase difference, are plotted in Figs. (a)a, (b)b and (c)c, respectively.
The metasurface exhibits very good transmission (near unity) for both and polarization around the frequency of operation. The phase difference reaches the required value of at the specified frequency of 10 GHz. Finally, the linear-to-circular conversion efficiency has been estimated from the and polarization amplitude and phase measurements and has been plotted in Fig. 11.
As can be seen, the linear-to-circular conversion efficiency is very good reaching about at 10 GHz with a -3-dB bandwidth of about .
Iii-B Generalized Birefringent Reflection and Refraction
The concept of generalized birefringent reflection and refraction consists in controlling independently the reflection and transmission angles and amplitudes of orthogonally polarized waves. To illustrate this concept, we have realized a polarization beam splitting (PBS) reflection-less metasurface that refracts in opposite directions normally incident and polarized waves. The synthesis of this metasurface essentially follows the procedure elaborated in the introduction of this paper and which corresponds to the second synthesis technique that was discussed. Accordingly, the transmission coefficients for and polarization are, respectively, given by and , where is the specified transmission angle. Note that separation of both polarizations was initially specified to be along the direction. The phase gradient, corresponding to and , has a period that is given by . For the realization of that metasurface, we decided to sample the period with 8 unit cells of lateral size . Consequently, the transmission angle is determined by the unit cell size and the number of unit cells and is, thus, given by .
Each unit cell has then been implemented to realize a specific phase shift for and polarizations, respectively, and . The transmission phases, for each unit cell, are given in Table III. Note that the absolute phase shift of a single unit cell is irrelevant, only the phase shift difference between adjacent unit cells (here ) matters.
As can be seen in Table III, the supercell (formed by the 8 unit cells) has an asymmetric phase progression meaning that the unit cells number 5, 6, 7 and 8 have opposite and phase shifts as the unit cells number 4, 3, 2 and 1, respectively. This means that the 4 lasts unit cells are simply the rotated version of the first 4 unit cells. Consequently, the realization of this metasurface is greatly simplified since only 4 unit cells need to be implemented.
After designing the supercell, we performed full-wave simulations to verify the beam splitting behavior of the metasurface. The polarization refraction yielded good result but, unfortunately, the polarization transformation was not good, which can be explained by the presence of spurious coupling between adjacent unit cells. While the coupling affected both polarizations, it turns out that it was more damaging to the polarization than the polarization. It has to be noted that the metasurface is non-uniform only in the direction while being perfectly uniform in the direction, this asymmetry in the structure was hypothesized to be the cause of the different behavior of the two polarizations. To overcome this difficulty, we modified the metasurface such that the same non-uniformity was present in both and directions. Consequently, the supercell is now made of unit cells instead of . The realized metasurface is shown in Fig. 12, note the sinusoidally varying pattern in the diagonal direction indicating the direction of the phase gradients. The metasurface is made was a repetition of 9 supercells.
The dimensions of each of the first 4 unit cells are given in Table IV. As said above, the 4 remaining unit cell are just the rotated version of the first 4 ones.
Because the metasurface has now a period in the diagonal direction, the dimension of the phase gradient period is reduced to . This changes the transmission angle to at 10 GHz.
The metasurface has then been measured. A horn antenna was used to generate the normally incident waves while a probe was scanning the near-field over a plane parallel to the metasurface in the transmission region. Near-field to far-field transformation  was then used to evaluate the transmission response of the metasurface. The measured and polarization transmissions, in the diagonal plane of the metasurface, are plotted in Fig. 13 as a dashed blue line and a solid red line, respectively. Note that the curves have been normalized with respect to the -polarized transmission.
As can be seen, the metasurface effectively separates the two polarizations which are refracted, with almost identical amplitude, at about and from broadside, respectively. The frequency corresponding to the results in Fig. 13 is about 10.4 GHz and the transmission efficiency, defined has the ratio between the transmitted electric field and the incident electric field, is about . The efficiency of the metasurface versus frequency is plotted in Fig. 14. The reasons for which the metasurface efficiency does not exceed can be explained partly by the presence of loss in the dielectric layers but mostly from undesired refraction orders (either in reflection or in transmission) that are due to the spurious coupling of the unit cells. For instance, zeroth diffraction orders are clearly visible in the measurements shown in Fig. 13.
Iii-C Orbital Angular Momentum Multiplexing
The last metasurface that was realized is a structure that generates waves possessing orbital angular momentum (OAM) of different topological charges depending on the polarization of the incident wave. The OAM wave that we have chosen as the transmitted wave is a Hypergeometric Gaussian (HyG) wave that corresponds to the solution of the paraxial Maxwell equations in cylindrical coordinates. The reason why the HyG wave has been chosen over the more common Bessel wave  is that the HyG has the advantage of being linearly polarized whereas the Bessel wave is either radially of longitudinally polarized making the multiplexing of two OAM waves with a single metasurface more difficult for Bessel waves.
The metasurface is thus required to transform an -polarized normally incident wave into an HyG wave of topological charge and -polarized normally incident wave into an HyG wave of topological charge . The electric field of an HyG wave reads 
where is the confluent hypergeometric function, is the gamma function, is the OAM topological charge, is a real parameter, and where , , with being the beam waist and the Rayleigh range given by .
The amplitude and phase of the HyG wave, for , are plotted in Figs. (a)a and (b)b, respectively. It is clear that this kind of wave has a non-periodic phase pattern, as evidenced in Fig. (b)b, contrary to the oblique transmitted plane waves that were specified for the polarization beam splitting metasurface in the previous section. This means that a larger number of unit cells has to be implemented because of the aperiodicity of the transformation. Moreover, the amplitude of the HyG wave is non-uniform which further complicates the realization of the unit cells. However, these difficulties may by overcome by assuming that the amplitude of the and transmission coefficients are instead of following the profile in Fig. (a)a. Despite the fact that this approximation might a priori seem extreme, it turns out that the main properties of the HyG wave may be obtained by only implementing its phase evolution. For instance, the null amplitude at the center of the wave is achieved by destructive interferences due to the phase rotation around the center. Moreover, the orbital angular momentum information is contained not in the amplitude but rather in the phase of the wave. These considerations justify the assumption that only the phase profile of the transmitted waves should be implemented while their respective magnitude can be assumed to be uniform and equal to 1. Additionally, the phase profiles of the two OAM waves were discretized by four phase samples each. Consequently, the metasurface is made of unit cells with phase shifts for and polarizations as given in Figs. (c)c and (d)d, respectively.
When combining together the phase shifts in Figs. (c)c and (d)d, it turns out that the total number of different unit cells composing the metasurface is 16. The Fig. 16 represents the unit cells of the metasurface. In that figure, each color corresponds to a specific unit cell having unique phase shift for and polarizations. Interestingly, the unit cells in the highlighted regions 1, 2 and 3 are quarter-wave plates, half-wave plates and isotropic wave plates (where ), respectively. Consequently, 6 out of the 10 correspond to the same unit cells but rotated in the plane of the metasurface with respect to each other.
The metasurface has then been simulated and measured and the results are reported in Figs. 18 and 19 which respectively correspond to and polarizations. In these two figures, plots (a) and (b) are the simulated amplitude and phase transmissions. These simulations were obtained by first measuring the radiated reference field of the exciting horn antenna at the position of the metasurface. Then, the expected scattered field of the metasurface was calculated using the antenna reference field and assuming ideal transmission of flat and unity amplitude and phase profiles as in Figs. (c)c and (d)d. This simulation technique allows fair comparison between the expected scattered fields and the measured scattered fields from the metasurface that are shown in plots (c) and (d).
The measured results are in good agreements with the expected simulated results. The topological charges of and are achieved with a transmission efficiency near at 10 GHz. Finally, the transmission efficiency of the metasurface was evaluated for a frequency band between 8 and 12 GHz. The result is reported in Fig. 20.
In the first part of this work, we have discussed and elaborated two different approaches for the synthesis of birefringent metasurfaces. The two methods yield the metasurface electric and magnetic susceptibilities either when the exact electromagnetic fields on both sides on the metasurface are specified or when the transmission and reflection coefficients are specified. The first synthesis technique is rigorous while the second one is an approximation. However, we have seen that both methods lead to similar results in terms of the metasurface scattering response. One of the main differences between the two techniques lies in the difficulty of the physical realization of the metasurfaces. While the first technique requires the implementation of complicated non-uniform amplitude and phase transmission and reflection coefficients, the second usually only requires the implementation of non-uniform phase profiles. Consequently, although less rigorous, the second synthesis method is generally preferred over the first one.
In the second part, we have presented the synthesis and realization of four different birefringent metasurfaces performing the operation of half-wave plate, quarter-wave plate, polarization beam splitting and orbital angular momentum multiplexing, respectively. These metasurfaces were synthesized based on the second synthesis technique and their measurements were found in good agreement with the expected scattering responses.
-  B. Saleh and M. Teich, Fundamentals of Photonics, ser. Wiley Series in Pure and Applied Optics. Wiley, 2007.
-  E. Bartholin, Experimenta crystalli islandici disdiaclastici quibus mira & infolita refractio detegitur. Phil. Trans. of the Royal Society of London, 1669.
-  N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nature Mater., vol. 13, no. 2, pp. 139–150, 2014.
-  D. Lin, P. Fan, E. Hasman, and M. L. Brongersma, “Dielectric gradient metasurface optical elements,” science, vol. 345, no. 6194, pp. 298–302, 2014.
-  C. Pfeiffer and A. Grbic, “Metamaterial Huygens’ surfaces: tailoring wave fronts with reflectionless sheets,” Phys. Rev. Lett., vol. 110, p. 197401, May 2013.
-  C. Pfeiffer, C. Zhang, V. Ray, L. J. Guo, and A. Grbic, “High performance bianisotropic metasurfaces: asymmetric transmission of light,” Phys. Rev. Lett., vol. 113, p. 023902, Jul 2014.
-  M. Kim, A. M. Wong, and G. V. Eleftheriades, “Optical huygensâ metasurfaces with independent control of the magnitude and phase of the local reflection coefficients,” Physical Review X, vol. 4, no. 4, p. 041042, 2014.
-  K. Achouri, M. A. Salem, and C. Caloz, “General metasurface synthesis based on susceptibility tensors,” IEEE Trans. Antennas Propag., vol. 63, no. 7, pp. 2977–2991, July 2015.
-  ——, “Birefringent generalized refractive metasurface,” in 2015 IEEE International Symposium on Antennas and Propagation USNC/URSI National Radio Science Meeting, July 2015, pp. 878–879.
-  K. Achouri, B. A. Khan, S. Gupta, G. Lavigne, M. A. Salem, and C. Caloz, “Synthesis of electromagnetic metasurfaces: principles and illustrations,” EPJ Applied Metamaterials, vol. 2, p. 12, 2015.
-  F. Ding, Z. Wang, S. He, V. M. Shalaev, and A. V. Kildishev, “Broadband high-efficiency half-wave plate: a supercell-based plasmonic metasurface approach,” ACS nano, vol. 9, no. 4, pp. 4111–4119, 2015.
-  A. Pors, M. G. Nielsen, and S. I. Bozhevolnyi, “Broadband plasmonic half-wave plates in reflection,” Optics letters, vol. 38, no. 4, pp. 513–515, 2013.
-  N. Yu, F. Aieta, P. Genevet, M. A. Kats, Z. Gaburro, and F. Capasso, “A broadband, background-free quarter-wave plate based on plasmonic metasurfaces,” Nano letters, vol. 12, no. 12, pp. 6328–6333, 2012.
-  Y. Zhao and A. Alù, “Manipulating light polarization with ultrathin plasmonic metasurfaces,” Physical Review B, vol. 84, no. 20, p. 205428, 2011.
-  A. Pors, O. Albrektsen, I. P. Radko, and S. I. Bozhevolnyi, “Gap plasmon-based metasurfaces for total control of reflected light,” Scientific reports, vol. 3, 2013.
-  M. Farmahini-Farahani and H. Mosallaei, “Birefringent reflectarray metasurface for beam engineering in infrared,” Optics letters, vol. 38, no. 4, pp. 462–464, 2013.
-  J. H. Lee, J. W. Yoon, M. J. Jung, J. K. Hong, S. H. Song, and R. Magnusson, “A semiconductor metasurface with multiple functionalities: A polarizing beam splitter with simultaneous focusing ability,” Applied Physics Letters, vol. 104, no. 23, p. 233505, 2014.
-  C. Pfeiffer and A. Grbic, “Controlling vector bessel beams with metasurfaces,” Phys. Rev. Applied, vol. 2, p. 044012, Oct 2014.
-  E. Karimi, S. A. Schulz, I. De Leon, H. Qassim, J. Upham, and R. W. Boyd, “Generating optical orbital angular momentum at visible wavelengths using a plasmonic metasurface,” Light: Science & Applications, vol. 3, no. 5, p. e167, 2014.
-  N. Yu, P. Genevet, M. A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: generalized laws of reflection and refraction,” Science, vol. 334, no. 6054, pp. 333–337, 2011.
-  C.-F. Chen, C.-T. Ku, Y.-H. Tai, P.-K. Wei, H.-N. Lin, and C.-B. Huang, “Creating optical near-field orbital angular momentum in a gold metasurface,” Nano letters, vol. 15, no. 4, pp. 2746–2750, 2015.
-  W. Wang, Y. Li, Z. Guo, R. Li, J. Zhang, A. Zhang, and S. Qu, “Ultra-thin optical vortex phase plate based on the metasurface and the angular momentum transformation,” Journal of Optics, vol. 17, no. 4, p. 045102, 2015.
-  M. M. Idemen, Discontinuities in the Electromagnetic Field. John Wiley & Sons, 2011.
-  V. Asadchy, M. Albooyeh, S. Tcvetkova, Y. Ra’di, and S. Tretyakov, “Metasurfaces for perfect and full control of refraction and reflection,” arXiv preprint arXiv:1603.07186, 2016.
-  M. A. Salem and C. Caloz, “Manipulating light at distance by a metasurface using momentum transformation,” Opt. Express, vol. 22, no. 12, pp. 14 530–14 543, Jun 2014.
-  B. O. Zhu and Y. Feng, “Passive metasurface for reflectionless and arbitary control of electromagnetic wave transmission,” IEEE Transactions on Antennas and Propagation, vol. 63, no. 12, pp. 5500–5511, Dec 2015.
-  Y. Vahabzadeh, K. Achouri, and C. Caloz, “Simulation of metasurfaces in finite difference techniques,” arXiv preprint arXiv:1602.04086, 2016.
-  C. Pfeiffer and A. Grbic, “Bianisotropic metasurfaces for optimal polarization control: Analysis and synthesis,” Phys. Rev. Applied, vol. 2, p. 044011, Oct 2014.
-  C. A. Balanis, Antenna theory: analysis and design. John Wiley & Sons, 2016.
-  E. Recami and M. Zamboni-Rached, “Localized waves: a review,” Advances in Imaging and Electron Physics, vol. 156, pp. 235–353, 2009.
-  E. Karimi, G. Zito, B. Piccirillo, L. Marrucci, and E. Santamato, “Hypergeometric-gaussian beams,” Opt. Lett., vol. 32, no. 21, pp. 3053–3055, 2007.