Uniform Ray Description of Physical Optics Scattering by Finite Locally Periodic Metasurfaces

—This paper presents a uniform ray description of electromagnetic wave scattering by locally periodic metasurfaces of polygonal shape. The model is derived by asymptotically evaluating the critical-point contributions of a physical optics scattering integral. It is valid for metasurfaces whose bulk scattering coefﬁcients are periodic functions of a phase parameter which, in turn, is a continuous and smooth function of surface coordinates. The scattered ﬁeld is expressed in terms of reﬂected, transmitted and diffracted rays that do not generally obey conventional geometrical constraints (i.e., Snell’s law and the Keller cone). An iterative technique is presented to determine the locations of critical points on one or multiple interacting metasurfaces. Numerical results demonstrating the utility and accuracy of the asymptotic physical optics model are also provided.


I. INTRODUCTION
E NGINEERED electromagnetic surfaces (EESs) are a class of printed metasurfaces designed specifically to enhance wireless signal propagation (for example, to improve coverage or reduce interference) in built-up environments [1], [2]. They are fabricated by printing appropriately designed conductive, dielectric, or even ferroelectric ink patterns on substrates such as plastic coatings, window glass, ceramic tile, or drywall. Similar to reflectarrays and transmitarrays [3], these patterns consist of discrete, subwavelength unit cells whose scattering properties can be designed to control the spatial properties of the reflected and transmitted wavefronts.
Many wavefront-shaping EESs employ locally periodic unitcell patterns, meaning that their scattering properties are modulated periodically over space [4], or at least approximately so over any sufficiently small subsurface. An example of a strictly periodic EES is the uniform grating [5], capable of deflecting incident waves in directions that do not obey Snell's law [6]. Examples of locally periodic EESs include reflective diffusers [7], non-uniform diffraction gratings that can be used to disperse wireless signal power, and planar lenses [8] capable of focusing.
Ray tracing [9] is the method of choice for simulating propagation characteristics in complex, electrically large environments such as office buildings and street canyons. To perform ray-based propagation prediction for complex environments equipped with metasurfaces, it is necessary to model the scattering behavior of these surfaces in terms of ray optics.
To the best of the author's knowledge, ray-optical models of metasurface scattering are currently not available, even though a workaround solution has been used with some success for impenetrable gratings and diffusers [10]. The latter approach involves substituting the metasurface with an "equivalent" geometrical structure with similar overall scattering characteristics. Unfortunately, this workaround is costly in terms of simulation time, does not always produce accurate results, and cannot readily be extended to penetrable metasurface designs [7].
The ray-optical metasurface scattering model presented in this paper is superior to the approach in [10] with respect to accuracy and computational complexity, and is applicable to any type of locally periodic metasurface. Its derivation is based in part on the approach taken by Albani et al. to model scattering by curved surfaces [11], but has been appropriately adapted to the related, but different problem of scattering by geometrically planar, space-modulated metasurfaces. The analysis begins by homogenizing [12] the reflection and transmission coefficients associated with individual unit cells, thus allowing the scattered field to be approximated in terms of a physical optics (PO) surface integral. The high-frequency approximation is then used to evaluate this integral as the discrete sum of its critical-point contributions [13], which represent specularly reflected and transmitted, as well as edgeand corner-diffracted rays, depending on where on the surface the critical points are located.
Other relevant prior work includes that of Borovikov et al. on diffraction gratings with slowly varying parameters [14], and work by Chou on scattering from periodic array structures with a truncation boundary [15]. The main contribution of the present work is that it is applicable to a significantly broader class of metasurfaces, including designs with nonuniform periodicity parameters (such as diffusers and lenses) and arbitrary polygonal shapes, while properly accounting for all edge effects, including corner diffraction. This paper also proposes an iterative method for determining the critical (ray interaction) points on one or more mutually interacting metasurfaces -a necessary ray-tracing step. Finally, it aims to provide insight into the relationships between the metasurface parameters (specifically, the gradient and Hessian of its phase parameter) and the properties of its ray-optical scattered fields (i.e., deflection angle and wavefront curvature). This paper is organized as follows. Section II provides a definition and analysis of the scattering problem, and derives mathematical expressions for the ray-optical scattered fields. Section III presents a method for determining the locations of the critical points. Section IV presents numerical results illustrating the utility and accuracy of the proposed model. Two metasurface designs are considered: an impenetrable (reflective) diffuser and a binary Fresnel-zone plate lens (FZPL) [16], the latter being an example of a penetrable metasurface supporting multiple spatial modes. Section V contains concluding remarks, and the appendices are provided as a selfcontained reference to the mathematical details necessary for the reader to implement the equations in Section II-F. This paper largely follows the notation used in [11], including the use of boldface symbols for vectors and dyadic tensors in 3-D space; for example, a andĀ. Hat notation is used for unit vectors; for example,â. Vectors and matrices in 2-D surface space are denoted by roman symbols; for example, a andĀ. An e jωt time-dependence is assumed and suppressed.

II. ANALYSIS
Let S denote a finite, locally periodic metasurface illuminated by an incident field, as illustrated in Fig. 1. S is assumed to be infinitely thin (geometrically), planar, and bounded by straight edges. Its illuminated (positive) and shadow (negative) sides are denoted by S + and S − , respectively. The halfspace on the positive side of S will be referred to as Region I, and the halfspace on the negative side as Region II. The surface normal,n, is chosen to always point into Region I and is, therefore, dependent on the incident field. In the following, the optional argument notation ·(q) indicates that a quantity is evaluated at the surface point where q 0 ,û 1 andû 2 are an arbitrarily chosen origin and a pair of basis vectors in the plane of S.

A. High-Frequency Approximation
Assuming that the dimensions of S are much larger than the wavelength, λ , its scattering properties can be well approximated by high-frequency asymptotic techniques that describe the incident and scattered fields in terms of rays with associated geometrical optics (GO) fields [17]. Under the high-frequency approximation, all ray interactions are local phenomena that depend on the properties of the incident field and the scattering object at the point of ray incidence only.
In the neighborhood of an arbitrary reference point q on S, any incident GO field E i can be approximated by [11] where k is the free-space wavenumber and ϕ i is a wavefront phase function that may be expanded up to second order as Here,k i denotes the incident wave direction at q , and is a curvature dyadic that contains the principal radii of curvature, ρ 1,i and ρ 2,i , and the corresponding principal directions, x 1,i andx 2,i , of the incident wavefront. As with any GO field, the incident magnetic field near q is related to the electric field by the local plane-wave approximation [11], [17] where η is the intrinsic impedance of free space.

B. Physical-Optics Approximation
Starting from a specialization of the Stratton-Chu equation for the electric field outside a bounded, source-free region of space [18], the radiative scattered field due to S can be written as [11], [19], [20] E s (p) = jk C r × J m + ηr ×r × J e e − jkr 4πr dS, where p is an observation point, r = rr is the vector from the integration point to p, and J m and J e are magnetic and electric surface currents on the closed surface C = S + ∪S − (see Fig. 1). These currents are related to the tangential components of the surface fields by J ± m = ∓n × E ± (7) where the superscripts + and − indicate on which side of S the currents and surface fields are evaluated. Under the PO approximation, the surface fields at any point of S are approximated by those that would be present on an infinite, uniformly periodic metasurface that is locally matched to S. While this approximation is clearly invalid near edges and corners, and can be poor if the periodicity varies too rapidly, it generally yields good results for the overall scattered field provided that the surface is electrically large [11] and its periodicity varies slowly.
The (electric) surface field is thus approximated by E i + E r on S + and by E t on S − , where and are the reflected and transmitted surface fields for a locally matching, infinite periodic metasurface.Γ andT are dyadic reflection and transmission coefficients that characterize the localized bulk scattering behavior of S, meaning that the microscopic responses of individual unit cells and their inner structures are homogenized into a (locally periodic) bulk material description. These coefficients generally depend on the propagation direction of the incident field, but this is not made explicit in the notation. If the metasurface is locally periodic, E r and E t can be locally expanded into GO and evanescent surface fields, as will be discussed below.
The subwavelength location-dependence ofΓ andT, in particular their phase characteristics, can be engineered so as to shape, e.g., deflect and/or (de)focus, the reflected and transmitted wavefronts. These coefficients can be directly measured, or computed with the aid of full-wave simulation tools [12], [21]. In this analysis they, or rather, the derivative parameters introduced below, are simply assumed to be computable on demand for any arbitrary surface location and incidence direction, and for fixed λ .

C. Local Periodicity
As they are approximately periodic in the vicinity of q ,Γ andT can be expressed by the Fourier series where m is an integer index referred to as the (spatial) mode number, and ψ denotes a phase parameter. The dependence of Γ andT on ψ is periodic with a constant period equal to one wavelength, λ ; ψ, in turn, is a continuous and smooth function of q, and describes the local subwavelength periodicity 1 of the metasurface, which may vary gradually over S. Assuming that the scale of periodicity variation is much longer than λ , the Fourier coefficients may be approximated as being stationary over subdomains of S whose dimensions are on the order of several wavelengths or less. In practice, the magnitudes of these coefficients often become negligibly small for sufficiently large |m|, so thatΓ andT are accurately described by the coefficients of truncated Fourier series. The variation of ψ near q can be approximated by the quadratic expansion ψ(q) ψ(q ) + g ψ (q ) · (q − q ) 1 It is conceptually straightforward to extend the analysis to (locally) biperiodic surface profiles, which are dependent on two phase parameters with different location dependencies, cf. [ where (14) denote the first-and second-order spatial derivatives of ψ, also referred to as the gradient and Hessian, respectively [11], [13]. The latter can be expressed in terms of a two-element vector and a 2 ×2 matrix, respectively, referenced to the surface basis vectors. Substituting (10a) and (10b) into (9a) and (9b), it is clear that the reflected and transmitted surface fields can be locally expanded into spectra of plane-wave spatial harmonics, in accordance with Floquet's theorem for periodic surfaces [22], [23]. The corresponding wave vectors, denoted by k with a = {r,t}, are slowly-varying functions of q. In order for the incoming and outgoing waves to be phase-matched, these vectors must satisfy whereĪ denotes the unit dyad. Their normal components are determined from the square roots of where care must be taken to select the correct roots for the reflected and transmitted fields. If the right-hand side of (15) does not exceed one,k are real-valued propagation direction vectors whose normal components must be chosen such that they point into Regions I and II, respectively. Modes for which the roots of (16) are imaginary-valued are evanescent and do not contribute to the radiative scattered field. However, they must be taken into account when solving boundary conditions, as in Section IV-B for example. The signs of imaginary normal components ofk (m) r andk (m) t must be chosen such that the amplitudes of the corresponding fields are bounded (i.e., decay exponentially) in Regions I and II, respectively.

D. Summary of Metasurface Properties
To set up the ray-optical metasurface scattering model presented in the remainder of this section, the following parameters are assumed to be available: 1) an origin, q 0 ; 2) a pair of surface basis vectors,û 1 andû 2 ; 3) an ordered set of corner points (vertices); 4) surface profiles of ψ, g ψ andH ψ ; 5) surface profiles ofΓ (m) andT (m) for mode numbers up to a fixed mode truncation order M, i.e., |m| ≤ M. The surface profiles are understood to support the entire area of S. It is convenient to define the surface parameters for discrete points on a 2-D grid aligned withû 1 andû 2 . The grid spacing can typically be much larger than the wavelength, but must be small enough to enable accurate interpolation. Because ψ, g ψ andH ψ are interrelated, it is sufficient to only define, for example, the phase gradient profile and the phase value at a reference point, and precompute the other two profiles by means of numerical differentiation and integration. Finally, instead of specifyingΓ (m) andT (m) directly, it is advantageous to define the scattering properties of S in terms of intrinsic parameters from which the former can be computed on demand, for any given incidence angle. This is illustrated in Section IV, which presents simulation results for metasurfaces characterized by locally periodic sheet impedance profiles. A method for computing the reflected and transmitted GO surface fields from locally periodic surface susceptibilities is presented in [24].

E. Evaluation of Scattering Integral
Using the results from Section II-C, the surface fields in (7) and (8) can be written as in which E Similar to (5), the magnetic surface fields are obtained via for a = {r,t}, which is valid for both propagating and evanescent (locally) plane waves [25]. Through substitution and some algebraic manipulation, (6) can now be rewritten as a vector sum of distinct scattering components with a ∈ {i, r,t}. Here, in which [11] G(k,r) = (n ·r)Ī −nr + (Ī −rr) · (k ·n)Ī −kn .
The summations in (20) are over non-evanescent modes only. The preceding equations show that the scattered field consists of a component due to the incident field, and multiple spatial modes due to the periodically varying reflection and transmission coefficients. As will be shown next, the former cancels out the incident field in the shadow zone behind S, while the latter are responsible for the backscattered field and the remainder of the forward-scattered field.

F. Critical-Point Contributions
In the high-frequency limit, i.e., for k → ∞, (21) is dominated by discrete contributions from a set of critical points, also referred to as stationary-phase points, on S, resulting in a ray description of the scattered field [11]. The scattering contribution from any one of these critical points, denoted by q c , can be written as where is the observation point, as previously, and s andk s are the observation distance and scattering direction. Furthermore, S c is a subdomain of S spanning several wavelengths around q c but becoming infinitesimally small as k → ∞, and ϕ s is the phase function of the scattered wavefront originating at q.
Noting that rr = sk s at q = q c , ϕ s can be approximated by Substituting (12), (3) and (26), (24) becomes The phase function in the integrand of (28) is given by In the following, we consider critical points of the first, second, and third kinds, which correspond with specular, edgediffracted, and corner-diffracted ray contributions, respectively [11], [13]. These critical points must be evaluated separately for all m. The reader is referred to the appendices for mathematical details on how to evaluate the integral in (28) for each kind of critical point.  (29) is stationary with respect to movement alongû 1 andû 2 , such that b (m) ·û 1,2 = 0 (32) andû 1,2 ·C (m) ·û 1,2 = 0. The latter inequality implies that specular points must be spatially isolated, as will generally be the case unless the observation point lies at a focal point, or S is uniformly periodic and ρ 1,i , ρ 2,i and s are infinite. If r . In the latter case, the productF (0) i · E i vanishes. For k s =k i , on the other hand, the critical-point contribution (28), if it exists, becomes equal to This scattered field component is equal in magnitude to the incident field at p but has the opposite sign, therefore cancelling it out for observation points in Region II for which a stationary-phase point exists in the interior of S, i.e., in the optical shadow region behind the metasurface. If a ∈ {r,t}, condition (32) can be thought of as a generalized version of Snell's laws of reflection and transmission [6]. For m = 0, or if g ψ is uniformly zero, the reflection angle equals the incidence angle and, because S is infinitely thin, the transmitted ray continues along the path of the incident ray. For |m| > 0, a non-zero phase gradient "breaks" the conventional law by deflecting the reflected and transmitted rays into different angles and/or planes depending on m. Eq. (32) can be used to determinek s for a ray-optical field incident at a given surface point or, as discussed in Section III, to determine the reflection and/or transmission point(s) given an observation point.
For any givenk i and m, (32) has two solutions fork s , denoted byk (m) r andk (m) t ; as discussed in Section II-C, these are considered valid only if (Ī −nn) ·k i − m g ψ < 1. The productF (m) a · E i is non-zero for only one of these (valid) solutions, such that the specularly reflected and transmitted fields vanish in Regions II and I, respectively. Where they do not vanish, these fields can be asymptotically evaluated as in whichs = s − mψ(q c ), and ρ 1,r = ρ 1,t and ρ 2,r = ρ 2,t are the radii of curvature of the reflected and transmitted wavefronts. For m = 0, these radii are generally different from those of the incident field. As shown in Appendix A, they can be determined from an eigenanalysis involving the projection of Q i − mH ψ (q c ) on S.
2) Edge-Diffracted Fields: Edges are defined by two adjacent corner points, an edge-parallel vectorê pointing from one to the other, and a transverse vectort perpendicular toê and n, pointing into the interior of S. For a critical point of the second kind to exist on a given edge, ϕ (m) must be stationary with respect to movement along the edge, such that b (m) ·ê = 0 (35) andê ·C (m) ·ê = 0. Similarly as before, this condition can be used to determine the half-angle of the Keller cone [6, Sec. 13.3.3] of scattering directionsk s originating at a given incidence point, or to determine the diffraction point(s) given an observation point. Referring to Appendix B for details, the scattered field from any of these diffraction points can be asymptotically evaluated as where the subscript e indicates edge diffraction, and ρ e,a is the radius of curvature of the diffracted field in the plane of the edge, which can be determined by projectingQ contains the transition function F(·) used in the uniform theory of diffraction (UTD). This function can be computed using an approximation in [26]. Its argument is given by with 3) Corner-Diffracted Fields: Critical points of the third kind exist at all corner points of S, regardless of the local behavior of ϕ (m) . Corners are defined by a corner point and two edge-parallel vectors,ê 1 andê 2 , pointing to the previous and next corner points. The scattering contributions from these diffraction points can be asymptotically evaluated as where the subscript c indicates corner diffraction, and is a corner diffraction coefficient [11]. In (42), T (·, ·, ·) denotes the UTD vertex transition function with arguments and This function can be approximated using the algorithm in [27]. Further details are provided in Appendix C. The preceding asymptotic PO analysis, in and of itself, does not reconstruct rays undergoing multiple successive interactions with S. Examples of such rays include multiply diffracted rays representative of coupling between different edges of the same surface [6, Sec. 13.3.7], and so-called complex surface rays excited and radiated back into free space at different surface points [28]. While such ray paths are physically valid and can even be the dominant mechanism in some applications [29], [30], they are believed to be of secondary importance for EESs, the main application considered in this paper.

III. DETERMINATION OF CRITICAL POINTS
Unlike conventional planar surfaces, locally periodic metasurfaces can contain multiple reflection or transmission points for given source and observation points, and their locations can generally not be determined simply using image theory. These locations can instead be found by applying Fermat's principle [6], which states that the phase lag, or optical length, of any ray path must be a local extremum with respect to displacements of its interaction points. To evaluate the optical length, denoted herein bys, of a ray path segment originating at a metasurface, one must correct its geometrical length, s, by the additional, location-dependent phase shift introduced by the metasurface, i.e.,s = s − mψ, m being the mode number as previously.
The following approach is valid for determining critical points of the first kind on consecutive metasurfaces. This procedure must be followed separately for all non-zero mode numbers; the image method [9] can be used for m = 0. It can easily be adapted to arbitrary sequences of ray interactions, including edge diffractions, at conventional surfaces as well as metasurfaces. A brief discussion of how to adapt the approach to critical points of the second kind is provided at the end of this section.
Referring to Fig. 2, the problem at hand is to determine the interaction points for which the overall optical length of the ray path is stationary, i.e., for all n: where ∆ n denotes a scalar displacement of candidate interaction point n in the directiond n . The geometrical length of ray segment n after displacing interaction points n and n + 1 is s n = s n 1 + 2(k n · ∆)/s n + ∆ 2 /s 2 For a ray interaction in the surface interior, as in Fig. 2, considering that the interaction point has two degrees of freedom. For small enough displacements, the additional optical path length induced by the metasurface at interaction point n is equal to −m multiplied by ψ n ψ n + (g ψ n ·d n )∆ n + 1 2 (d n ·H ψ n ·d n )∆ 2 n .
The critical points are determined by iteratively updating an initial vector of candidate points by a vector of scalar displacements, two for each (interior) point. This vector is determined by simultaneously solving (45) for all interaction points, resulting in a banded system of linear equations (of bandwidth 3). For a ray interacting in the interiors of two consecutive metasurfaces, for example, the displacement vector is found by solving the matrix equation resulting from setting to zero. Here, all quantities other than the displacement values are determined by applying (51a)-(51e) to the candidate points found in the previous iteration. The algorithm has converged when the elements of the right-most vector of (52) are sufficiently small: the condition that a p n + b p n = 0 for n = {1, 2} and p = {1, 2} requires that all incident and outgoing ray directions satisfy (32).
Being a Newton-search algorithm, the above procedure converges quadratically if the initial vector of candidate points is sufficiently close to the solution [31]. Suitable initial points can be found as a by-product of a prior shooting and bouncing ray (SBR) stage [9] or, in case of a single, isolated metasurface, by a brute-force search over a point grid on the surface [7].
To determine the location of a critical point of the second kind on an edge of the metasurface, the vectord n in (48) is made identical to the edge-parallel vector,ê n , and only a single linear equation is added to the system. It is not difficult to show that the corresponding matrix elements (51a)-(51e) can be found by substitutingê n for bothû 1,n andû 2,n .

IV. SIMULATION EXAMPLES
This section presents simulation results for two examples of locally periodic metasurfaces, and provides a comparison of the asymptotic PO model proposed herein to direct numerical integration of the scattering integral (21) over a uniform point grid with λ /10 spacing, using the trapezoid rule. The latter method is computationally expensive, but does not invoke the high-frequency approximation.

A. Simulation Setup
The simulation setup is similar for both metasurfaces. Each is centered at the origin of a Cartesian coordinate system, lies in the plane defined by z = 0, and is X × Y = 100λ × 100λ in size. Each is illuminated by a plane wave that is normally incident, along the negative z-direction, and whose electric field is polarized along the y-direction (TM y ). The observation points lie in the y = 0 plane. The surface parameter profiles for the asymptotic PO model are defined on a 10λ point grid, and parameter values in between grid points are computed by means of bilinear interpolation.
The reflection and transmission properties at each grid point are governed by the isotropic sheet impedance boundary conditions [32], which can be written in the form Here, Z e and Z m denote the electric and magnetic sheet impedances, assumed to vary as a locally periodic function of ψ. They are related to the scalar reflection and transmission coefficients for normal incidence, Γ and T , by These equations are used to synthesize the impedance profile required to achieve a given desired reflection and/or transmission behavior.Γ (m) andT (m) , |m| ≤ M, the modal reflection and transmission coefficients for arbitrary incidence directions, are computed fromk i , g ψ , and the spatial harmonics of Z e and Z m at the ray incidence point, using the method presented next. The mode truncation order, M, is set to 5.

B. Computation ofΓ (m) andT (m)
After substitution of (17)- (19) and a transformation to the spatial frequency domain (i.e., multiplication by e − jknψ followed by integration over one period of ψ), the boundary conditions (53a) and (53b) can be written as which must be satisfied for p = {1, 2}. Here, with q = {1, 2, 3}, are the unknown quantities of interest; can be computed a priori;û 1 andû 2 are the surface basis vectors, as before, andû 3 is identical ton. Furthermore, δ n denotes the Kronecker delta function, and Z (n) e and Z (n) m are the nth spatial harmonics of Z e (ψ) and Z m (ψ); see also [24].
Eqs. (55a) and (55b) form a system of linear equations that can be solved for Γ (m) q and T (m) q , |m| ≤ M, by also requiring that which follows from the fact that all field components must be perpendicular to their corresponding propagation directions. If desired, the full dyadicsΓ (m) andT (m) can be constructed from the solutions for two orthogonal polarizations of E i .

1) Diffuser:
The first type of metasurface considered is a reflective diffuser designed to disperse a normally incident wave into a continuous range of horizontal deflection angles [7]. In order to achieve this behavior, the periodicity of the surface profile is chosen such that for |x| ≤ X/2 and |y| ≤ Y /2. In (59), α 1 and α 2 are the desired deflection angles for rays normally incident at x = ±X/2, respectively. The sheet impedance profile is modulated so as to force T in (54a) and (54b) to be uniformly zero, and Γ to be equal to for 0 ≤ ψ < λ , thus supporting a single mode (m = 1). Γ (1) is set to 0.9, which corresponds to an approximately 1-dB reflection loss. Figs. 3 and 4 show simulated distributions of the field scattered by a 30-60 • diffuser (α 1 = 30 • and α 2 = 60 • ). They show the expected behavior for normal wave incidence, consisting of a deflected and divergent reflected field that is uniform (i.e., continuous) across the 30 • and 60 • reflection shadow boundaries associated with the surface edges parallel to the y-axis. The scattered field in the area behind the surface as seen from its illuminated side, visible in Fig. 3, is the shadow field that partially cancels the incident field and thus attenuates the total field in the shadow region. Fig. 4 indicates excellent agreement between asymptotic PO and numerical integration, confirming the validity of rayoptical modeling of scattering by the class of finite, locally periodic metasurfaces considered in this work. This figure also shows the individual contributions of specular, edgediffracted and corner-diffracted fields to the total scattered field. While the specular field dominates the main lobe of the scattering pattern, edge diffraction is the most significant contribution everywhere else; corner diffraction is negligible in this example.
2) Fresnel-Zone Plate Lens: The second example is for a simplified variant of the binary Fresnel-zone plate lens (FZPL) described in [16], and illustrates the capability of asymptotic PO to simulate metasurfaces that give rise to concave wavefronts and support transmitted as well as reflected modes. Binary FZPLs consist of patterns of alternatingly transparent and opaque rings or, in the present example, strips that pass the odd, and reflect the even Fresnel zones on a surface (e.g., a window) between the source of the incident wave and the desired focal point or, in this example, focal line. The periodicity of the surface considered here is defined by for |x| ≤ X/2 and |y| ≤ Y /2; x f and z f are the (x, z) coordinates of the focal line. The impedance profile is modulated such that Γ = 0 and T = 0.9 for 0 ≤ ψ < λ /2, and Γ = 0.9 and T = 0 for λ /2 ≤ ψ < λ .
Figs. 5 and 6 show simulated scattered field intensities for a binary FZPL with a focal line at (x f , z f ) = (350, −350)λ . Similar to the previous example, the scattered field around φ s = 0 • is a shadow field that partially cancels the incident field in the shadow region. The scattered field around φ s = 45 • is the m = 1 transmitted mode, and exhibits the narrow width and high field intensity expected of a focused beam. The lobe around φ s = −45 • is the m = −1 transmitted mode; it and other undesired modes can be suppressed by more sophisticated, non-binary FZPL designs employing dielectric strips, for example, thereby improving aperture efficiency [16]. Fig. 6 again shows excellent agreement between the asymptotic PO model and numerical integration, with the exception of a small region around the focal line at φ s = 45 • , where the former method fails to accurately predict the gain of the focused beam. This failure is due to an inherent limitation of ray-optical methods at and near so-called caustics [6], or observation points where the divergence factors of the specular and edge-diffracted fields possess singularities. Finally, as in the previous example, the relative contribution of corner diffraction to the total scattered field is small. V. CONCLUSIONS This paper has presented a uniform ray description of scattering by locally periodic metasurfaces with polygonal shapes. This broad class of metasurfaces includes many of the EES designs that have recently received increased attention due to their potential to enhance propagation conditions in fifth-generation mobile communications networks operating at millimeter-wave frequency bands. Ray-based models of wireless signal propagation are uniquely suitable for predicting and optimizing wireless performance in electrically large environments, such as the urban areas envisioned to be served by such networks. The work presented herein extends the applicability of these models to large environments outfitted with EESs. Work is currently under way to incorporate the methods presented in this paper in commercial ray-tracing software marketed to wireless engineers and researchers.
Numerical results presented in this paper show excellent agreement between the new ray-based approach and numerical integration of the PO integral. It is important, however, that its accuracy also be compared to methods that do not invoke the PO approximation but, instead, rigorously enforce a given set of metasurface boundary conditions such as, for example, the generalized sheet transition conditions (GSTCs) [32]. This is the topic of an accompanying paper [24], in which the accuracy of the new approach is verified, with generally very good results, against an Integral Equation based Boundary Element Method (BEM) solver applied to a variety of electrically large, 2D scattering geometries. As discussed in [24], full-wave analysis of large 3D geometries similar to those considered in Section IV is beyond the capability of currently available computational hardware, even without considering the interaction of the metasurface with its potentially much larger propagation environment. This reinforces the need for computationally efficient simulation methods, even if they are approximate and therefore sacrifice some of the accuracy of full-wave solvers.
Possible directions for future research include the incorporation of scattering contributions due to rays undergoing multiple successive interactions with the metasurface, such as the complex surface rays [28] that are of special interest for metasurface antenna applications [30], and extensions to more complex scattering geometries involving (collections of) conformal and compound metasurfaces. Additional validations against full-wave simulations or measurements are desirable in order to gain increased confidence in the ray-optical approach.

APPENDIX A ASYMPTOTIC EVALUATION OF SPECULAR FIELDS
Consider the integral over the planar surface S, in which A(q) is a slowly-varying amplitude function, is a quadratic phase function, andC is of the form whereC is independent of s, and s > 0. The surface point q c is a critical point of the first kind if it lies in the interior of S and the integrand phase is stationary with respect to movement around q c , such that the linear phase term vanishes. If this is the case, I(k) can be asymptotically approximated for k → ∞ by [ where the 2 × 2 matrix denotes the projection ofC on the plane spanned by the surface basis vectors of S,û 1 andû 2 . Furthermore, σ is an integer quantity defined such that σ = 2 if the eigenvalues ofC are both positive, σ = −2 if they are both negative, and σ = 0 if they have opposite signs. Using (64),C can be written as is a projection matrix that premultiplies a surface point q in (u 1 , u 2 ) coordinates such that the 2-D matrix-vector product Θ q is its projection in a local frame spanned byx 1,2 . The basis vectorsx 1,2 may be chosen arbitrarily, as long as they are orthogonal to each other and the scattering directionk s . Recognizing that the determinant of a matrix product of square matrices equals the product of their determinants, and that detΘ = detΘ T =k s ·n (69) for any valid choice ofx 1,2 , the determinant ofC can be written as detC = (k s ·n) 2 In this expression, ρ 1 and ρ 2 denote the radii of curvature of the scattered wavefront, which are determined from the eigenvalues of the curvature matrix E is a matrix whose columns, which represent the eigenvectors of the curvature matrix, yield the principal directions of the scattered wavefront with respect to the arbitrarily chosen basis vectorsx 1,2 [35], [36]. Eq. (65) can now be evaluated as in which the so-called divergence factor is computed as ρ 1 ρ 2 (ρ 1 + s)(ρ 2 + s) = ρ 1 ρ 2 (ρ 1 + s)(ρ 2 + s) 1 2 e jνπ/2 (73) ν ∈ {0, 1, 2} being the number of negative eigenvalues ofC. This number is equal to the number of singular points, or caustics, crossed at s = −ρ 1 and/or s = −ρ 2 , when moving from s = 0 to the observation distance s > 0 [17]. Such crossings can only occur if at least one of the two radii of curvature is negative.

APPENDIX B ASYMPTOTIC EVALUATION OF EDGE-DIFFRACTED FIELDS
Consider again the surface integral given by (62)-(64). The surface point q c is a critical point of the second kind if it lies on an edge of S and the integrand phase is stationary with respect to movement along the edge, such that b ·ê = 0, whereê is parallel to the edge. If this is the case, I(k) can be asymptotically approximated as [13,Eq. (36)] with Here,t is the transverse vector defined in Section II-F2, F(·) is the UTD Fresnel transition function [13], and the squareroot function is defined such that it returns negative imaginary values for negative real arguments. Eq. (74) is uniform in the sense that it causes the total (specular plus edge-diffracted) scattered field to be continuous across the shadow boundary, where critical points of the first and second kinds cannot be treated as being isolated from each other [34].
Using (64), the productê ·C ·ê, which is the second-order derivative of the phase function along the edge, can be written asê where β is the diffraction angle defined by (36), and ρ e = sin 2 β e ·C ·ê (77) is the radius of curvature of the diffracted field in the plane of the edge. The radius of curvature in the transverse plane is zero, characteristic of cylindrical waves. Eq. (74) can now be evaluated as

I(k)
√ 2πs e jπ/4 F(X e ) k 3/2 (b ·t) sin β A(q c )e − jkϕ(q c ) ρ e ρ e + s , where the divergence factor is computed as ρ e ρ e + s = ρ e ρ e + s Here,ê 1 andê 2 are parallel to the two edges constituting the corner; each points from the corner point along its edge. T (·) denotes the UTD vertex transition function [13] and, as previously, the branch of the square-root function is chosen such that it returns negative imaginary values for negative real arguments. As a result, x 1 and x 2 can be both real or both imaginary (and w is real), or x 1 (x 2 ) can be real and x 2 (x 1 ) and w are both imaginary [27]. Note that (81) is independent of s, which makes the corner-diffracted fields proportional to 1/s, characteristic of spherical waves, for which both radii of curvature are zero.