## Abstract

Modeling meta-surfaces as thin metamaterial layers with continuously varying bulk parameters, we employed a rigorous mode-expansion theory to study the scattering properties of such systems. We found that a meta-surface with a linear reflection-phase profile could redirect an impinging light to a non-specular channel with nearly 100% efficiency, and a meta-surface with a parabolic reflection-phase profile could focus incident plane wave to a point image. Under certain approximations, our theory reduces to the local response model (LRM) established for such problems previously, but our full theory has overcome the energy non-conservation problems suffered by the LRM. Microwave experiments were performed on realistic samples to verify the key theoretical predictions, which match well with full-wave simulations.

© 2013 Optical Society of America

## 1. Introduction

Last decade has witnessed tremendous progresses in using metamaterials (MTMs) to manipulate light on a subwavelength scale. MTMs are artificial materials composed by manmade functional electromagnetic (EM) microstructures, typically in deep-subwavelength sizes and exhibiting tailored electric and/or magnetic responses. Early studies were largely conducted on *homogenous* MTMs, with discovered light-manipulation phenomena including negative refraction [1], super and hyper lensing [2–5], polarization control [6], and so on. Recently, with the help of transformation optics (TO) theory [7,8], *inhomogeneous* MTMs with *slowly* changing material properties were widely studied, based on which more fascinating effects were discovered, such as invisibility cloaking [9,10], illusion optics [11–13], lensing [14–17] and beam bending [18]. These works revealed that a carefully designed “*inhomogeneity*” could manipulate the local propagation phase inside the medium so as to *adiabatically* guide light travelling along a desired trajectory dictated by the TO theory, leading to new physics and phenomena.

Very recently, much attention were paid to *inhomogeneous* systems with *abruptly* changing materials properties, in particular, ultra-thin MTMs (also called meta-surfaces) constructed by carefully designed planar subwavelength components with tailored EM responses [19–29]. Rather than modulating the *propagation phase* inside a bulk medium, these planar systems explored another degree of freedom to modulate the lateral distribution of the *abrupt phase change* of reflected/transmitted light across the systems. It was shown that the transmission/reflection of light follow a generalized Snell’s law at the interface between air and a carefully designed gradient meta-surface, in which the parallel momentum of light is not conserved at the interface but rather gain an additional term contributed by the lateral gradient of the transmission/reflection phase change [19,20]. In particular, the meta-surface can perfectly convert an impinging propagating wave to a bounded surface wave under certain conditions [20]. Other fascinating wave-manipulation effects include focusing by planar lens [27–29], generalization of optical vortex [19,23], and so forth.

In sharp contrast to the exciting achievements on experimental side, the theoretical developments for studying such complex systems are far behind. Apparently, standard approaches for homogeneous MTMs are not applicable here. Furthermore, the TO theory [7,8], being a powerful tool for inhomogeneous MTMs with *slowly varying* properties, is also not suitable for present inhomogeneous meta-surfaces exhibiting *abruptly* changing material properties [19,20]. In previous works, theoretical understandings were either based on heuristic Fermat-Huygens wave interference arguments or based on full-wave simulations [19,21–29]. The former is very intuitive but is far from rigorous and thus cannot fully explain the rich physical phenomena observed (say, the multi-channel reflections/refractions discovered in [19]). The latter is basically a computational experiment, which is rigorous enough but suffers the limitation of physically less illuminating. We notice that a local response model (LRM), which assumed each local part of the inhomogeneous meta-surface to response *locally* to the incident wave, was established previously for such systems [30, 31]. However, we will show in the following sections that the LRM suffers severe energy *non-conservation* problems in many cases. Therefore, a general theoretical approach to study light scatterings by such inhomogeneous meta-surfaces, which can yield energy-conserved results, is still lacking and is highly desired.

Here, we develop such a theoretical framework. In contrast to previous theories assuming zero thicknesses for the studied systems, we model the meta-surfaces as thin MTM layers with finite thickness and with laterally varying bulk EM parameters. Although some limited results of the theory (obtained under restricted conditions) have been reported in [20], here we present all details of the theoretical developments under general conditions (Sec. 2) and then apply the theory to study two specific examples, one exhibiting a linear reflection-phase distribution and another a parabolic one (Sec. 3). Section 4 is devoted to highlighting the key advantage of our theoretical approach. Specifically, we show that our theory can reduce to the LRM under several approximations, but the full theory (without taking these approximations) has overcome the energy non-conservation problems faced by the LRM. Finally, we design and fabricate a series of realistic samples and perform microwave experiments to verify the theoretical predictions (Sec. 5). We conclude our paper in Section 6.

## 2. Developments of the mode-expansion theory

We chose a particular system as a concrete example to illustrate the developments of our theory. The model system is schematically shown in Fig. 1(a), which is an inhomogeneous MTM layer of thickness *d* (much smaller than wavelength *λ*) put on top of a perfect electric conductor (PEC). For simplicity, we assume that the MTM layer is inhomogeneous only along *x* direction, but is invariant along both *y* and *z* directions. Adding a PEC substrate makes the entire system totally reflecting for EM waves so that we do not need to worry about the transmitted signals, which significantly simplifies the theoretical developments.

Although the model looks ideal, we emphasize that actually it can be realized by various means in practice. For example, one can put powders of high-index materials [32] onto a PEC with laterally different densities to realize this model. Also, it has been proven in [6, 33, 34] that the usual high-impedance structures (HIS) can be well represented by such a double-layer model (in a homogeneous version), and thus it is straightforward to design an (laterally) inhomogeneous HIS to practically realize the model presented in Fig. 1.

As shown in Fig. 1, the entire space is divided into three regions, where region I denotes the vacuum, region II denotes the inhomogeneous MTM with permittivity and permeability matrices given by

*x*. We are mostly concerned on the scattering properties of such an inhomogeneous system under arbitrary light illuminations.

Consider first the case that the excitation is a transverse-electric (TE) polarized incident plane wave, with field components explicitly given by (in region I)

*ω*being the working frequency and

*c*the speed of light, ${k}_{x}^{in}$ and ${k}_{z}^{in}={({k}_{0}^{2}-{({k}_{x}^{in})}^{2})}^{1/2}$ are the parallel and vertical components of the

**k**vector for the incident wave, and

*Z*

_{0}= (

*μ*

_{0}/

*ε*

_{0})

^{1/2}is the impedance of vacuum. Here, we use the physics convention and omit the common time oscillation term $\mathrm{exp}[-i\omega t]$ throughout this paper. Due to lacking of translational invariance on the

*xy*-plane, the reflected beam does not necessarily exhibit a single

**k**vector, but must in principle be a linear combination of plane waves with all allowed

**k**vectors (each defined as a reflection channel). In general, there is no restriction on choosing the value of

*k*for each channel. However, it is more convenient in practical computations to introduce periodic boundary conditions at the two ends of considered system (with length denoted by

_{x}*L*), i.e., (

**E**,

**H**)

_{x}_{= 0}= (

**E**,

**H**)

_{x}_{=}

*. Introducing a super periodicity makes computations more tractable and will not affect the final results when*

_{L}*L*is large enough. In addition, we note that some of the designed/fabricated meta-surfaces already exhibited super periodicities to which our method is naturally applicable [19–22].

For each reflection channel, the corresponding EM fields can be written as

*N*for

*n*so that the final number of reflection channels is 2

*N*+ 1. To obtain reasonable results, convergence tests against both

*N*and

*L*(in case of a finite system without super periodicity) should be carefully performed. With these notations, EM fields in region I can be formally expressed as a sum of the incident plane wave and the reflected waves,

Let us turn to consider the EM fields in region II. We need to first calculate the eigen wave functions of EM waves travelling inside this region, which are governed by the following equation

*E*field component. Since the MTM is homogenous along

_{y}*z*direction, we can perform variable separation to assume thatwhere superscript $\pm $ stands for forward ($+$) and backward ($-$) propagating waves with wave-vector

*q*being a positive number. The parameter

_{z}*q*is used to label the eigenmodes inside region II and will be determined later. Inserting Eq. (6) to (5), we find that the $G({q}_{z},x)$ function satisfies

_{z}Once *E _{y}* is known by solving Eq. (7),

**H**fields can be derived from Maxwell’s equations. Explicitly, we have

It is hard to solve Eq. (7) analytically, so that we now develop a numerical approach. As shown in Fig. 1(b), discretizing a super cell of the inhomogeneous MTM slab into 2*N* + 1 sub-cells [36], each with length *h* = *L*/(2*N* + 1), we can rewrite the differential Eq. (7) as the following 2*N* + 1 linear equations

**sub cell with position located at$${H}_{mm\text{'}}=({k}_{0}^{2}{\mu}_{M,m}^{x}{\epsilon}_{M,m}^{y}-2{\mu}_{M,m}^{x}\gamma ){\delta}_{mm\text{'}}+{\mu}_{M,m}^{x}\gamma {\delta}_{m,m\text{'}-1}+{\mu}_{M,m}^{x}\gamma {\delta}_{m,m\text{'}+1}.$$Diagonalizing the **

*x*= (2_{m}*m*- 1)*h*/2, and ${\epsilon}_{M,m}^{y}$, ${\mu}_{M,m}^{y}$ and $G({q}_{z},m)$ are the values of functions ${\epsilon}_{M}^{y}(x)$, ${\mu}_{M}^{x}(x)$ and $G({q}_{z},x)$ taken at the position*x*=*x*. It is worth mentioning that we take the periodic boundary condition so that${\epsilon}_{M,m+2N+1}^{y}={\epsilon}_{M,m}^{y}$, ${\mu}_{M,m+2N+1}^{x}={\mu}_{M,m}^{x}$ and $G({q}_{z},m+2N+1)=G({q}_{z},m).$ We can further rewrite Eq. (10) as the following matrix formwhere ${H}_{mm\text{'}}$ is a $(2N+1)\times (2N+1)$ matrix with elements defined by [37]_{m}**H**matrix, we can obtain 2*N*+ 1 eigenvalues labeled as ${q}_{z,j}^{2}$. The eigen vector corresponding to the*j-*th eigenvalue is just [$G({q}_{z,j},1)$,…, $G({q}_{z,j},i)$,…, $G({q}_{z,j},2N+1)$]*, which gives the wave function of $G({q}_{z},x)$ in a discretized manner. The discretized versions of*^{T}*E*,_{y}*H*,_{x}*H*can be easily obtained from the $G({q}_{z},x)$ function based on Eqs. (6)–(9)._{z}Knowing all non-vanishing field components for every eigenmode, we can then formally write the EM fields in region II as linear combinations of these eigenmodes. Thus, in general we have

*N*+ 1) eigenvalues of

*q*and [${C}^{+}({q}_{z,j})$, ${C}^{-}({q}_{z,j})$] is another set of coefficients to be determined.

_{z}We now match the boundary conditions at two interfaces. On the MTM/PEC interface located at *z* = −*d*, EM fields should follow the PEC boundary condition ($\overrightarrow{n}\times \overrightarrow{E}\equiv 0$) so that

*z*= 0. The tangential EM fields (

*E*and

_{y}*H*) should be continuous crossing the interface, i.e.,Put the explicit forms of fields (Eqs. (4) and (13)) into Eq. (16), we get that

_{x}*N*+ 1) linear equations with (2

*N*+ 1) unknowns $\left\{\rho ({k}_{x}^{r,n})\right\}$ and (2

*N*+ 1) unknowns$\left\{{C}^{+}({q}_{z,j})\right\}$. Solve these equations by standard linear algebra, we get the final form of reflection coefficients as

*n*th-order plane wave (in region I) and the eigenmodes in region II, and ${A}^{-1}({k}_{x}^{in},{q}_{z,j})$ is the inverse matrix of $A({q}_{z,j},{k}_{x}^{in})$. Therefore, all reflection coefficients of the EM waves scattered to different channels can be calculated from the above equations.

The same technique can be easily extended to the case of a transverse-magnetic (TM) plane wave excitation (i.e., ${\overrightarrow{H}}^{inc}\left|\right|\widehat{y},{\overrightarrow{k}}^{in}={k}_{x}^{in}\widehat{x}-{k}_{z}^{in}\widehat{z}$). After tedious calculations, we found that the reflection coefficients $\rho ({k}_{x}^{r,n})$ can still be calculated by Eq. (18), except that the two overlapping integrals are now defined as

**H**matrix originally defined in Eq. (12) should now be defined as

We note that the coefficients in front of exp[2*iq _{z}*

_{,}

*] exhibit different signs in Eqs. (20) and (22), due to different boundary condition requirements for TE and TM cases. Besides this, in fact we can interchange ${\overleftrightarrow{\epsilon}}_{M}$ and ${\overleftrightarrow{\mu}}_{M}$ to derive the TM formulas from the TE case, thanks to the excellent symmetry properties of EM fields.*

_{j}dWe mention three points before concluding this section. First, after knowing the scattering properties of the system for both TE and TM excitations, we can in principle obtain all information of the scattered field under an *arbitrary* excitation (not necessarily a plane wave). Second, the developed technique is so general that there are no difficulties to extend it to more complicated cases, say, the inhomogeneous MTMs with materials properties depending on both *x* and *y*. Third, so far the developed formulas and the model adopted are directly applicable only to those meta-surfaces *with* ground planes [6, 20, 22, 33, 34], but extensions of the theory to the cases *without* ground planes (e.g., single-layer meta-surfaces [19,21,23–28]) are straightforward.

## 3. Applications of the theory

The developed theoretical approach can be applied to many inhomogeneous meta-surfaces. Below we present two explicit examples.

#### 1. Meta-surfaces with linear reflection-phase profiles

In this subsection, we design a gradient system (based on the model depicted in Fig. (1)) working for the TE polarization, and then employ the newly developed mode-expansion theory to study its scattering properties. To determine the model parameters of the system, we adopt a local-phase argument similar to that taken in [20]. Although obviously such a designing scheme neglected the diffraction effects, our mode-expansion theory will take all such effects into account, and therefore can serve as a serious justification on such a designing scheme. Specifically, we fix the model parameters (i.e., ${\overleftrightarrow{\epsilon}}_{M}(x),{\overleftrightarrow{\mu}}_{M}(x)$) by letting the whole structure exhibit a linearly varying reflection-phase profile

for normally incident EM wave with polarization $\overrightarrow{E}\left|\right|\widehat{y}$ (instead of $\overrightarrow{E}\left|\right|\widehat{x}$ assumed in [20]). We note that there are multiple solutions for ${\overleftrightarrow{\epsilon}}_{M}(x),\text{\hspace{0.17em}}{\overleftrightarrow{\mu}}_{M}(x)$ making Eq. (24) satisfied, and here we take two specific solutions to illustrate the applications of our theory. The first model is an ideal impedance-matched model, where we assume that ${\overleftrightarrow{\epsilon}}_{M}(x)={\overleftrightarrow{\mu}}_{M}(x)$. A simple calculation shows thatwith*κ*≡

*ξ*/

*2k*

_{0}

*d*. Realizing the difficulties in matching the impedance at every local point, in the second model we set ${\epsilon}_{M}^{y}=1$ and let ${\mu}_{M}^{x}$ vary as a function of

*x*. The ${\mu}_{M}^{x}(x)$ distribution can be easily obtained by solving Eq. (24) with local reflection phase determined by the following equation [20]

As explicit illustrations, we show in Figs. 2(a) and 2(b) the distributions of material properties for two models with different values of *ξ*. We employed the mode-expansion theory to study the scattering coefficients of meta-surfaces constructed by two different models with different *ξ* under normal excitations with TE polarizations. The spectra depicted in Figs. 2(c) and 2(d) show that $|\rho ({k}_{x}^{r}){|}^{2}$ [38] take maximum values at ${k}_{x}^{r}=\xi $ for all the cases studied, indicating that the incident wave is indeed redirected to the desired anomalous channel after reflections.

We also employed the mode-expansion theory to study the cases of oblique incident excitations. Figure 3(a) shows the $|\rho ({k}_{x}^{r}){|}^{2}$ spectra for a *ξ* = 0.8*k*_{0} meta-surface (based on the impedance-matched model) under illuminations with different oblique angles specified by the values of ${k}_{x}^{in}$. Different spectra are maximized at different values of ${k}_{x}^{r}$, but it is interesting to note that the relation

*ξ*is always provided by the meta-surface. We performed a systematic study on three different meta-surfaces with

*ξ*= 0, 0.4

*k*

_{0}, 0.8

*k*

_{0}, respectively, under TE excitations with different incident angles. The results depicted in Fig. 3(b) show that Eq. (27) holds perfectly for all the cases studied.

We now identify the conversion efficiency for such anomalous reflection. Since the anomalous reflection beam (with parallel wave vector ${k}_{x}^{r}$) travels along an *off-normal* direction, its beam width is reduced by a factor of $\mathrm{cos}{\theta}_{r}={(1-{({k}_{x}^{r}/{k}_{0})}^{2})}^{1/2}$ as compared to the incident beam along the *normal* direction [20, 39]. Therefore, the final expression for the conversion efficiency should be

*anomalous*reflections. One may easily verify that the conversion efficiency

*R*→ 1 for most cases studied, indicating that almost all incoming energies are converted to these non-specular channels after reflections by the meta-surfaces. However, in the case of the [${\epsilon}_{M}^{y}=1,{\mu}_{M}^{x}(x)$] model with

_{c}*ξ*= 0.8

*k*

_{0}, we found that $|\rho (0.8{k}_{0}){|}^{2}=1.58$ [see Fig. 2(d)] so that the corresponding conversion efficiency can be easily calculated as

*R*~0.95 based on Eq. (28). We note that a small peak appears at

_{c}*k*= −0.8

_{x}*k*

_{0}in the spectrum [see blue line in Fig. 2(d)], which means that some of the incoming energy is converted to other channels due to the Bragg scatterings, so that the conversion efficiency to the desired anomalous reflection channel is not 100%.

#### 2. Meta-surfaces with parabolic reflection-phase distributions

Focusing a plane wave to a point image is always fascinating. Conventional lenses are much thicker than wavelength and exhibit certain curved shapes. Several previous works [27, 29, 40] have shown that an *ultra-thin flat* MTM lens exhibiting a parabolic reflection-phase distribution

*l*is the focal length. The key idea is to use the reflection phase gained at the meta-surface to compensate the propagation phase difference for waves radiated from different local positions at the meta-surface [see Fig. 4(a)]. However, previous works [27, 29, 40] only contained simulation and experimental results, without analytical calculations on model systems. In this subsection, we employ the newly developed mode-expansion theory to study the scattering properties of a meta-surface satisfying Eq. (29), as another application of our theory. Still based on the geometry shown in Fig. 1, we assume that the material parameters for the capping layer are given bywhere $v$ is an arbitrary value to keep all parameters positive (here we set it as 30). Figure 4(b) depicts the profiles of ${\epsilon}_{M}^{y}(x)$ as well as $\Phi (x)$ for the system under study, where we have assumed

_{focus}*l*= 5

_{focus}*λ*and took a super periodicity

*L*= 10

*λ*to truncate the otherwise divergent parameter profile. We studied the scattering properties of such a meta-surface, and Fig. 4(c) depicts the calculated $|\rho ({k}_{x}^{r}){|}^{2}$ spectrum. We note that the spectrum is symmetrical for $\pm $

*k*and does not show a delta-like peak at some particular

_{x}*k*position, which is different from the cases studied in last subsection. This is reasonable since we do not expect the reflected beam to be a plane wave. However, it is difficult to see the focusing effect from the $|\rho ({k}_{x}^{r}){|}^{2}$ spectrum. Therefore, we computed the real-space field distribution from the $|\rho ({k}_{x}^{r}){|}^{2}$ spectrum based on Eq. (4). Figure 4(d) presents the calculated result (with incoming plane wave deducted), from which we can clearly identify a focal point at the pre-decided position. The

_{x}**E**field at the focus is enhanced roughly 6 times with respect to that of the incident wave. However, the focusing effect suffers some distortions, which is due to the finite size of the flat lens (i.e., the super-periodicity

*L*here), as already discussed in [29].

## 4. Comparisons with the local response model

Although the theory developed in last section is rigorous, it looks quite complicated and how it is connected with previously established theories (say, the LRM) is unclear. In this section, we show that our theory can recover the LRM [19, 21–26, 30, 31, 41] under several approximations. However, we note that the LRM inevitably face energy non-conservation problems in non-specular reflection cases, while our full theory (without taking approximations) always yields correct (energy-conserved) results. For simplicity, we study the impedance-matched meta-surfaces with ${\epsilon}_{M}^{y}={\mu}_{M}^{x}=f(x)$ throughout this section.

The *first approximation* is to set *k _{z}* /

*k*

_{0}= 1 in the second equation of Eq. (17). The physics and limitation of this approximation will be discussed later. Under this approximation, Eq. (17) can be rewritten as,

*G*functions, we get the solution for ${C}^{+}({q}_{z})$ as,

*k*) and the inner eigen wave-function specified by the perpendicular wave-vector

_{x}*q*. With Eq. (39), we can finally rewrite Eq. (38) as a standard T matrix form,

_{z}Equation (40) clearly shows that the considered problem is a second-order scattering process *under the adopted approximation*. When a plane wave with parallel wave-vector ${k}_{x}^{in}$ strikes on the meta-surface, it first couples into all possible eigenstates inside the MTM with coupling strength ${V}_{{q}_{z},{k}_{x}^{in}}^{}$. These eigenstates propagate inside the MTM (along *z* direction) without interacting each other, and after reflection by the MTM/PEC interface, they propagate along -*z* direction and couple out of the MTM to the out-side plane wave mode ${k}_{x}^{r}$ with coupling coefficient ${V}_{{k}_{x}^{r},{{q}^{\prime}}_{z}}^{*}$. Obviously, Eq. (40) neglected the multiple scattering processes.

Base on Eq. (40), we can obtain an analytical solution of $\rho ({k}_{x}^{r})$ if the *V* matrix (Eq. (39)) is known. However, the eigen wave-function $G({q}_{z},x)$ is difficult to solve analytically from Eq. (7). Fortunately, numerical solutions of $G({q}_{z},x)$ function give us enough hints to “guess” an analytical form. Figure 5 shows the computed $G({q}_{z},x)$ functions (*in a discretized version*) with different *q _{z}*, obtained by the numerical approach described in Sec. 2. The most striking feature of the $G({q}_{z},x)$ function is that it resembles very much as a $\delta $ function in a

*discretized version*, with peak appearing at $x$ which makes the condition

*q*values. Excellent agreement between dashed lines and the peak positions is noted.

_{z}Such an important observation motives us to take the *second approximation*. We assume that the eigen wave function inside the inhomogeneous MTM is given by

*q*and therefore do not talk to each other. Then, the final eigen wave-function associated with a particular

_{z}*q*will be highly localized at the very sub-cell that supports this wave-vector, as shown in Fig. 5.

_{z}Put Eq. (42) into Eq. (38), we finally obtain that

*r*

_{local}(x) = −exp[2

*if*(

*x*)

*k*

_{0}

*d*] is the local reflection coefficient of the system at the position

*x*. An inverse Fourier transform of Eq. (43) gives the distribution of scattered field (measured at the

*z*= 0 plane) as

Equations (43) and (44), derived from our full theory with two approximations, are exactly the same as the LRM widely used for such systems [19, 21–26, 30, 31]. However, we must point out that Eq. (43) and thus the LRM actually suffer severe energy non-conservation problems. We take the *f* (*x*) = $1+\kappa x$ model to explicitly illustrate this point. Put *κ* ≡ *ξ*/2*k*_{0}*d* into Eq. (43), we get that

Equation (45) tells us that the LRM predicts that the reflected wave always carries an additional wave vector *ξ* with amplitude 1. However, energy-conservation law requires the *normal* components of energy fluxes to be conserved after the EM wave is reflected by a flat surface. Let us define a function as

*R*≡ 1 for such a system since there is no loss and transmission here. Unfortunately, the LRM cannot yield the energy-conserved results in many cases, and the total energy of reflected wave can be either larger or smaller than that of the incident beam, depending on the values of

*ξ*and ${k}_{x}^{in}$. As an illustration, we employed both our full theory (without taking the two approximations) and the LRM to calculate the reflection efficiency

*R*for meta-surfaces with different

*ξ*under normal-incidence excitations, and Fig. 6(a) compares the

*R*~

*ξ*relations calculated by two theories. Obviously, the energy non-conservation issue is more severe in large

*ξ*cases for the LRM while our full theory

*always*yields energy-conserved results. As another example, we show in Fig. 6(b) how

*R*depends on the incidence parallel wave-vector ${k}_{x}^{in}$ for a meta-surface with a fixed

*ξ*= 0.4

*k*

_{0}. Again, we found that the LRM cannot yield energy-conserved results in general, and can even yield reflection efficiency exceeding 1 when ${k}_{x}^{in}<-0.2{k}_{0}$. In contrast, our full theory always gives energy-conserved results.

To understand the inherent physics accounting for the failure of LRM, we re-examined the two approximations adopted. For the first one, we found it can be justified only in small-$\xi $ cases, since in such cases for the most relevant channel (the anomalous reflection channel) we have *k _{z}* /

*k*

_{0}= (1 − (

*ξ*/

*k*

_{0})

^{2})

^{1/2}≈1. When

*ξ*is large, such a simplification is no longer valid. In fact, the term

*k*/

_{z}*k*

_{0}represents the impedance mismatch between the scattered and incident waves. When the scattered wave is not along the specular channel, the term

*k*/

_{z}*k*

_{0}generates important local-field corrections which cannot be neglected. However, the LRM assumed that the response of each part of the system is solely dependent on the incident field on that very point, but is independent on the direction of outgoing wave. Obviously, this approximation is too rough to completely neglect the local-field corrections for the non-specular reflections. The second approximation neglected the couplings between adjacent units, which is also questionable in general cases. Therefore, we conclude that the failure of the LRM in certain cases is due to its neglecting the local field corrections and the coupling effects.

## 5. Experimental and simulation verifications

In this section, we design and fabricate realistic gradient meta-surfaces to experimentally verify the theoretical predictions presented in Fig. 2(b). In practice, it is difficult to construct a MTM system representing the model depicted in Fig. 1(a) with *continuous* ${\overleftrightarrow{\epsilon}}_{M}(x),{\overleftrightarrow{\mu}}_{M}(x)$ distributions. Instead, typically the designed/fabricated MTM layers exhibit stepwise *discontinuous* distributions of ${\overleftrightarrow{\epsilon}}_{M}(x),{\overleftrightarrow{\mu}}_{M}(x)$, depending on how many building blocks adopted in one super cell. To model such realistic situations, we truncate the continuous distributions of the $\mu \left(x\right)$ profiles in the [${\epsilon}_{M}^{y}=1,\text{\hspace{0.17em}}{\mu}_{M}^{x}$**]** models for both *ξ* = 0.4*k*_{0} and *ξ* = 0.8*k*_{0} cases to stepwise versions as shown in Fig. 7(a). We then employed the mode-expansion theory to calculate the scattering properties of such systems. Figure 7(b) shows that, for such stepwise models which can better represent the realistic situations, the $|\rho ({k}_{x}^{r}){|}^{2}$ spectra remain essentially the same as those of their continuous counterparts, indicating that such models still work very well to achieve the desired anomalous reflection effect.

We can therefore design realistic samples according to the stepwise *μ*(*x*) profiles shown in Fig. 7(a). As already discussed in Sec. 3, there are multiple ways to realize such model practically, and here we choose one of them. In [6, 33, 34], it had already been proved that the HIS can be very well represented by a double-layer model with a thin homogeneous magnetic MTM layer put on a PEC (e.g., homogeneous version of Fig. 1(a)). The physics is that near-field interaction between the top metallic layer and the bottom PEC ground plane in a HIS can generate a magnetic resonance with anti-parallel currents induced on two layers [33]. Both propagating-wave and surface-wave properties of a HIS can be *accurately* reproduced by calculations based on such an effective-medium model [33], demonstrating the validity of the model. Thus, we can use the HIS as a building block to design our gradient meta-surfaces according to the *μ*(*x*) profiles depicted in Fig. 7(a).

A building block is shown in the inset to Fig. 8, which is a sandwich system consisting of a metallic “H” and a metallic ground plane, separated by a dielectric spacer (with *ε _{r}* = 3.9) of thickness

*d*. We note that although the building block is essentially the same as that adopted in meta-surfaces designed previously [20], here we have to carefully re-perform the designs since the present meta-surface works for TE-polarized incidence wave rather than for TM case considered in [20].

We first employed finite-difference-time-domain (FDTD) simulations [43] to study how the *μ _{eff}* parameter of a building block is “tuned” by varying its geometrical parameters. By changing the central bar length

*L*

_{1}of the “H”, the magnetic resonance frequency of a building block can be efficiently changed, so that the

*μ*parameter can be dramatically modified at a given frequency. Choosing the working frequency as 15 GHz, we performed FDTD simulations to retrieve

_{eff}*μ*parameters of such structures, and depicted the retrieved

_{eff}*μ*parameters as functions of

_{eff}*L*

_{1}in Fig. 8. Such calculations greatly facilitate us to select appropriate building blocks to construct the desired stepwise

*μ*(

*x*) profiles as shown in Fig. 7, for both

*ξ*= 0.4

*k*

_{0}and 0.8

*k*

_{0}meta-surfaces [44]. We then fabricated the two samples and a picture of the

*ξ*= 0.4

*k*

_{0}sample is shown in Fig. 9(a).

We performed microwave experiments to characterize the functionalities of the fabricated samples. As schematically shown in Fig. 9(b), we illuminated these meta-surfaces by normally incident TE-polarized (with $\overrightarrow{E}\left|\right|\widehat{y}$) microwaves with a double-ridged horn antenna [45], and then measured the far-field (FF) scattering patterns using another identical double-ridged horn antenna. Both emitting and receiving horn antennas were connected to a vector-field analyzer (Agilent E8362C). The measured signals were normalized against a reference single, which was obtained through replacing the meta-surface by a metal plate of the same size. Figures 9(c) and 9(d) depict the normalized scattering patterns for the samples with *ξ* = 0.4*k*_{0} and *ξ* = 0.8*k*_{0}, respectively. The experimental results are in excellent agreement with FDTD simulations on realistic structures.

It is difficult to make direct comparisons between experimental/simulation results and the mode-expansion model calculations, since the former are obtained with finite-sized samples while the latter are with infinite systems. Nevertheless, meaningful comparisons can still be made in terms of reflection angle and reflection efficiency. We can easily identify from Figs. 9(c) and 9(d) that the peaks of two scattering patterns appear at 22.5° and 52.5°, respectively. Using the formula ${k}_{x}^{r}={k}_{0}\mathrm{sin}{\theta}_{r}$ to retrieve the parallel **k** vectors of the reflected beams, we find that the measured ${k}_{x}^{r}$ are about 0.38*k*_{0} and 0.79*k*_{0}, respectively, which are in good agreement with the rigorous mode expansion calculations in Fig. 7(b). Meanwhile, both measured and simulated radiation patterns exhibit only a *single* main peak at the anomalous reflection angle in both *ξ* = 0.4*k*_{0} and *ξ* = 0.8*k*_{0} cases, which are again in good agreement with the *ρ*(*k _{x}*) spectra calculated by the mode-expansion theory (Fig. 7(b)).

We also studied the scattering patterns of the fabricated meta-surfaces under oblique-angle excitations. Figures 10(a) and 10(b) present the normalized FF patterns of the two meta-surfaces under external illuminations at different incident angles. Again, the measured patterns are in excellent agreement with FDTD simulations, both showing that the reflected beams have been redirected to non-specular channels with very high efficiencies. Two points are worthy being emphasized. First, we note that in the case of Fig. 10(b) the reflection beam and the incident one appear at the same side with respect to the surface normal, indicating that the reflection is the so-called “negative” reflection. Second, in the measured scattering patterns, the normal (specular) reflection modes are nearly completely suppressed, which are also consistent with the theoretical calculations depicted in Fig. 2. This later property is an important advantage of present system over previous structures suffering the multi-mode conversion shortcomings [19, 21], and explains why the present system can have such high conversion efficiencies (nearly ~100%) for the anomalous reflection.

Finally, we systematically measured the scattering patterns for two samples under illuminations of input waves with incident angles varying within the whole angle region allowed. Figure 10(c) depicts the obtained ${k}_{x}^{r}$ for the reflected beam versus ${k}_{x}^{in}$ which is the parallel **k** vector of the incident wave for the two samples, which are again in excellent agreement with corresponding FDTD simulation results. We note that all measured/simulated data fall into two separate lines defined by Eq. (27) - the general Snell’s law, which are in turn, agreeing perfectly with the mode-expansion results recorded in Fig. 3(b).

## 6. Conclusions

In summary, based on a new model for inhomogeneous meta-surfaces, we established a general theoretical framework to study the scattering properties of such systems, and applied it to two particular examples. Our theory recovers previously established local response model under several simplifications, but the full theory has overcome the energy non-conservation problems encountered by previous theory. We designed and fabricated realistic structures according to theoretical calculations, and performed microwave experiments and full-wave simulations to verify the key theoretical predictions of our theory.

## Acknowledgments

This work was supported by NSFC (60990321, 11174055), the Program of Shanghai Subject Chief Scientist (12XD1400700) and MOE of China (B06011). QH acknowledges financial supports from NSFC (11204040) and China Postdoctoral Science Foundation.

## References and links

**1. **R. A. Shelby, D. R. Smith, and S. Schultz, “Experimental verification of a negative index of refraction,” Science **292**(5514), 77–79 (2001). [CrossRef] [PubMed]

**2. **J. B. Pendry, “Negative refraction makes a perfect lens,” Phys. Rev. Lett. **85**(18), 3966–3969 (2000). [CrossRef] [PubMed]

**3. **N. Fang, H. Lee, C. Sun, and X. Zhang, “Sub-diffraction-limited optical imaging with a silver superlens,” Science **308**(5721), 534–537 (2005). [CrossRef] [PubMed]

**4. **D. O. S. Melville and R. J. Blaikie, “Super-resolution imaging through a planar silver layer,” Opt. Express **13**(6), 2127–2134 (2005). [CrossRef] [PubMed]

**5. **S. A. Ramakrishna, J. B. Pendry, M. C. K. Wiltshire, and W. J. Stewart, “Imaging the near field,” J. Mod. Opt. **50**(9), 1419–1430 (2003).

**6. **J. M. Hao, Y. Yuan, L. X. Ran, T. Jiang, J. A. Kong, C. T. Chan, and L. Zhou, “Manipulating electromagnetic wave polarizations by anisotropic metamaterials,” Phys. Rev. Lett. **99**(6), 063908 (2007). [CrossRef] [PubMed]

**7. **U. Leonhardt, “Optical conformal mapping,” Science **312**(5781), 1777–1780 (2006). [CrossRef] [PubMed]

**8. **J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,” Science **312**(5781), 1780–1782 (2006). [CrossRef] [PubMed]

**9. **D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial electromagnetic cloak at microwave frequencies,” Science **314**(5801), 977–980 (2006). [CrossRef] [PubMed]

**10. **W. S. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, “Optical cloaking with metamaterials,” Nat. Photonics **1**(4), 224–227 (2007). [CrossRef]

**11. **Y. Lai, J. Ng, H. Y. Chen, D. Z. Han, J. J. Xiao, Z. Q. Zhang, and C. T. Chan, “Illusion optics: the optical transformation of an object into another object,” Phys. Rev. Lett. **102**(25), 253902 (2009). [CrossRef] [PubMed]

**12. **Y. Lai, H. Y. Chen, Z. Q. Zhang, and C. T. Chan, “Complementary media invisibility cloak that cloaks objects at a distance outside the cloaking shell,” Phys. Rev. Lett. **102**(9), 093901 (2009). [CrossRef] [PubMed]

**13. **H. Chen, C. T. Chan, and P. Sheng, “Transformation optics and metamaterials,” Nat. Mater. **9**(5), 387–396 (2010). [CrossRef] [PubMed]

**14. **U. Levy, M. Abashin, K. Ikeda, A. Krishnamoorthy, J. Cunningham, and Y. Fainman, “Inhomogenous dielectric metamaterials with space-variant polarizability,” Phys. Rev. Lett. **98**(24), 243901 (2007). [CrossRef] [PubMed]

**15. **A. O. Pinchuk and G. C. Schatz, “Metamaterials with gradient negative index of refraction,” J. Opt. Soc. Am. A **24**(10), A39–A44 (2007). [CrossRef] [PubMed]

**16. **O. Paul, B. Reinhard, B. Krolla, R. Beigang, and M. Rahm, “Gradient index metamaterial based on slot elements,” Appl. Phys. Lett. **96**(24), 241110 (2010). [CrossRef]

**17. **N. Kundtz and D. R. Smith, “Extreme-angle broadband metamaterial lens,” Nat. Mater. **9**(2), 129–132 (2010). [CrossRef] [PubMed]

**18. **B. Vasić, G. Isić, R. Gajić, and K. Hingerl, “Controlling electromagnetic fields with graded photonic crystals in metamaterial regime,” Opt. Express **18**(19), 20321–20333 (2010). [CrossRef] [PubMed]

**19. **N. F. 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 **334**(6054), 333–337 (2011). [CrossRef] [PubMed]

**20. **S. L. Sun, Q. He, S. Y. Xiao, Q. Xu, X. Li, and L. Zhou, “Gradient-index meta-surfaces as a bridge linking propagating waves and surface waves,” Nat. Mater. **11**(5), 426–431 (2012). [CrossRef] [PubMed]

**21. **X. Ni, N. K. Emani, A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Broadband light bending with plasmonic nanoantennas,” Science **335**(6067), 427 (2012). [CrossRef] [PubMed]

**22. **S. Sun, K.-Y. Yang, C.-M. Wang, T.-K. Juan, W. T. Chen, C. Y. Liao, Q. He, S. Y. Xiao, W.-T. Kung, G.-Y. Guo, L. Zhou, and D.-P. Tsai, “High-efficiency broadband anomalous reflection by gradient meta-surfaces,” Nano Lett. **12**(12), 6223–6229 (2012). [CrossRef] [PubMed]

**23. **P. Genevet, N. Yu, F. Aieta, J. Lin, M. A. Kats, R. Blanchard, M. O. Scully, Z. Gaburro, and F. Capasso, “Ultra-thin plasmonic optical vortex plate based on phase discontinuities,” Appl. Phys. Lett. **100**(1), 013101 (2012). [CrossRef]

**24. **F. Aieta, P. Genevet, N. Yu, M. A. Kats, Z. Gaburro, and F. Capasso, “Out-of-plane reflection and refraction of light by anisotropic optical antenna metasurfaces with phase discontinuities,” Nano Lett. **12**(3), 1702–1706 (2012). [CrossRef] [PubMed]

**25. **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 Lett. **12**(12), 6328–6333 (2012). [CrossRef] [PubMed]

**26. **L. Huang, X. Chen, H. Mühlenbernd, G. Li, B. Bai, Q. Tan, G. Jin, T. Zentgraf, and S. Zhang, “Dispersionless phase discontinuities for controlling light propagation,” Nano Lett. **12**(11), 5750–5755 (2012). [CrossRef] [PubMed]

**27. **F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, and F. Capasso, “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” Nano Lett. **12**(9), 4932–4936 (2012). [CrossRef] [PubMed]

**28. **X. Chen, L. Huang, H. Mühlenbernd, G. Li, B. Bai, Q. Tan, G. Jin, C. W. Qiu, S. Zhang, and T. Zentgraf, “Dual-polarity plasmonic metalens for visible light,” Nat. Commun. **3**, 1198 (2012). [CrossRef] [PubMed]

**29. **X. Li, S. Y. Xiao, B. G. Cai, Q. He, T. J. Cui, and L. Zhou, “Flat metasurfaces to focus electromagnetic waves in reflection geometry,” Opt. Lett. **37**(23), 4940–4942 (2012). [CrossRef] [PubMed]

**30. **D. Berry, R. Malech, and W. Kennedy, “The reflectarray antenna,” IEEE Trans. Antennas Propag. **11**(6), 645–651 (1963). [CrossRef]

**31. **S. Larouche and D. R. Smith, “Reconciliation of generalized refraction with diffraction theory,” Opt. Lett. **37**(12), 2391–2393 (2012). [CrossRef] [PubMed]

**32. **F. L. Zhang, Q. Zhao, L. Kang, J. Zhou, and D. Lippens, “Experimental verification of isotropic and polarization properties of high permittivity-based metamaterial,” Phys. Rev. B **80**(19), 195119 (2009). [CrossRef]

**33. **J. M. Hao, L. Zhou, and C. T. Chan, “An effective-medium model for high-impedance surfaces,” Appl. Phys. A Mater. Sci. Process. **87**(2), 281–284 (2007). [CrossRef]

**34. **D. Sievenpiper, L. Zhang, R. Broas, N. G. Alexopolous, and E. Yablonovitch, “High-impedance electromagnetic surfaces with a forbidden frequency band,” IEEE Trans. Microwave Theory Tech. **47**(11), 2059–2074 (1999). [CrossRef]

35. These reflection channels could also be understood as the Floquet modes diffracted by our super-periodic system.

**36. **In our computational approach, we have to set the number of sub-cells divided identical to the number plane waves adopted in region (both are 2*N* + 1), in order to ensure that the number of restraints equals to that of variables.

**37. **For two boundary indexes, we have the following off-diagonal matrix elements${H}_{1,2N+1}={\mu}_{M,1}^{x}\gamma $, ${H}_{2N+1,1}={\mu}_{M,2N+1}^{x}\gamma $ according to the periodic boundary condition.

**38. **Since the super-cell length *L* is very large, the distribution of those discretized *k _{x}^{r,n}* is almost continuous. Thus, in what follows, we use

*ρ*(

*k*) to represent

_{x}^{r}*ρ*(

*k*) for simplicity.

_{x}^{r,n}**39. **C. Qu, S. Y. Xiao, S. L. Sun, Q. He, and L. Zhou, “A theoretical study on the conversion efficiencies of gradient meta-surfaces,” Europhys. Lett. **101**(5), 54002 (2013). [CrossRef]

**40. **A. Pors, M. G. Nielsen, R. L. Eriksen, and S. I. Bozhevolnyi, “Broadband focusing flat mirrors based on plasmonic gradient metasurfaces,” Nano Lett. **13**(2), 829–834 (2013). [CrossRef] [PubMed]

**41. **M. Albooyeh, D. Morits, and C. R. Simovski, “Electromagnetic characterization of substrated metasurfaces,” Metamaterials **5**(4), 178–205 (2011). [CrossRef]

**42. **P. Sheng, “Wave scattering formalism,” in *Introduction to Wave Scattering, Localization and Macroscopic Phenomena*, R. Hull, R. M. Osgood, eds. (Springer, 2006).

**43. **EastFDTD v2.0 Beta, DONGJUN Science and Technology Co., China.

**44. **For the *ξ* = 0.4*k*_{0} sample, a super cell contains 10 pairs of “H” (altogether 20 ones) in one supercell, with *L*_{1} values of those 10 pairs set as 1.3 mm, 2.68 mm, 2.98 mm, 3.14 mm, 3.24 mm, 3.36 mm, 3.48 mm, 3.66 mm, 4.08 mm, and 5.7 mm. For the *ξ* = 0.8*k*_{0} sample, a super cell contains 10 “H” in one super cell with *L*_{1} parameters the same as the case of *ξ* = 0.4*k*_{0}.

**45. **The gain of the employed double-ridged horn antenna is about 14dB~15dB in this frequency region.