Ray Optical Scattering from Uniform Reflective Cylindrical Metasurfaces using Surface Susceptibility Tensors

A ray optical methodology based on the uniform theory of diffraction is proposed to model electromagnetic field scattering from curved metasurfaces. The problem addressed is the illumination of a purely reflective uniform cylindrical metasurface by a line source, models the surface with susceptibilities and employs a methodology previously used for cylinders coated in thin dielectric layers [1]. The approach is fundamentally based on a representation of the metasurface using the General Sheet Transition Conditions (GSTCs) which characterizes the surface in terms of susceptibility dyadics. An eigenfunction description of the metasurface problem is derived considering both tangential and normal surface susceptibilities, and used to develop a ray optics (RO) description of the scattered fields; including the specular geometrical optical field, surface diffraction described by creeping waves and a transition region over the shadow boundary. The specification of the fields in the transition region is dependent on the evaluation of the Pekeris caret function integral and the method follows [1]. The proposed RO-GSTC model is then successfully demonstrated for a variety of cases and is independently verified using a rigorous eigenfunction solution (EF-GSTC) and full-wave Integral Equation method (IE-GSTC), over the entire domain from the deep lit to deep shadow.


I. INTRODUCTION
Electromagnetic (EM) metasurfaces (MS) have been introduced as a possible solution to a wide range of applications at many frequencies from microwave to optical [2]- [6]. Applications include; 5G routing [7]- [9], the creation of holograms, cloaking, and possible camouflaging of objects [3], [10], [11]. Metasurfaces are sub-wavelength periodically tiled surfaces forming an array of scattering elements. The surface is thin (sub-wavelength) and transforms the incident radiation into scattered and transmitted waves. This transformation can be controlled by the design of the scattering element (unit-cell) and is capable of creating a large array of magnitude, phase and polarization responses for specified incident fields.
Although the detailed structure of the MS is sub-wavelength (thickness and periodicity) applications generally require an electrically large surface, placed in an extensive environment. A good example of this is the use of MS to route 5G signals within a workplace, in which surfaces, which extend over 100 of wavelengths (λ), could be strategically placed inside Scott Stewart, T. J. Smy and S. Gupta, are with the Department of Electronics, Carleton University, Ottawa, Ontario, Canada. Email: scott.stewart@carleton.ca. a building, steering, splitting and spreading the signal carriers for optimal signal propagation and covering blind spots.
This design of MS systems therefore has a multi-scale nature, which poses challenges. The detailed design of the unit-cell to achieve a particular transformation such as beam forming will need to be done in a full-wave simulator assuming periodic boundary conditions. However, such approaches are very limited, by computational factors, for the simulation of surfaces comprised of many unit-cells or when placed in a large environment.
A very productive approach to simulation on a larger scale is the introduction of the Generalized Sheet Transition Conditions (GSTC) to model field scattering off the surface. In such an approach the surface is assumed to be zero thickness and is represented by four homogenized surface susceptibility dyadics (χ) [12]- [14]. For a MS that is deeply sub-wavelength in thickness and unit-cell size, the GSTCs provide a rigorous macroscopic field solution for a given incident field, allowing for the determination of the propagating reflected and transmitted fields. Full wave unit-cell simulations can be used to extracted the susceptibilities needed to correctly model the surface [15], [16].
The GSTCs can be used for analytical solutions and have been incorporated into volumetric techniques using a Yee cell descritization or finite elements [16]- [21]. The analytical techniques are limited by various fundamental assumptions and simplifications and the volumetric numerical techniques become cumbersome if the domain is more then a few hundreds of wavelengths in size. The incorporation of the GSTCs into a surface integral approach to model scattering problems [11], [22], [23] enables a larger scale length to be addressed, but again as the domain size is increased this will become challenging and impractical.
To simulate MS systems accurately at very large scale lengths, an approximate method of analysis needs to be used. Traditionally the use of geometrical optics (GO) and ray tracing (RT) enhanced by a theory of diffraction have been used for many decades to solve similar scattering problems. The incorporation of a uniform theory of diffraction (UTD) creates a ray optics (RO) infrastructure and has enabled the accurate prediction of scattering in a wide variety of situations; including, radio propagation modeling [24], analysis of EM structures such as impenetrable gratings and diffusers [25], conventional curved surfaces [26], and diffraction gratings with slowly varying parameters [27].
Very limited work has been done to develop ray optical methods for metasurfaces [28]- [30]. Work on the use of RT and UTD for planar modulated MS is presented in [28] and [29]. In this work, a UTD formulation is developed for scattering off a non-uniform planar finite MS represented by a locally periodic variation in the susceptibilities. The work uses a high frequency approximation to incorporate planar and edge effects into a RO model, comparing the resultant fields to a full wave simulator based on integral equations. The non-uniformity of the MS allows for modeling complex field transformations. An obvious, and a needed, extension of this work is to model curved surfaces, where electrically thin metasurface architectures may be envisioned to conform to non-planar surfaces for greater wave control. The introduction of a curved surface into a UTD framework has a long history [31]- [34] and has been primarily concerned with analysis of the "canonical" problem of a cylindrical surface illuminated by either a line source (cylindrical waves) or a plane wave. The introduction of a curved surface considerably complicates the RO picture of the fields due to surface diffraction creating field components in the shadowed region. Initially this was analysed by Keller who characterized these field components as "creeping" waves -these waves hug the surface, shedding field rays as they propagate [31], [35]. The rays that represent the creeping waves can be derived from a rigorous solution to the cylindrical problem using a summation of eigenfunctions. This work was then expanded upon by Pathak, and many others [32], [33], to correctly model the transition between the lit region and shadowed region.
Indeed, a key factor in the development of this work is a requirement for a smooth transition from the lit region, where GO optics provides an accurate solution, to the deep shadow region, where a solution based on creeping waves is appropriate. The analysis of this transition region is complex and requires the development of a surface reflection characterized by an appropriate version of the Pekeris caret function -a Fock-Type integral. The nature of this integral is dependent on the boundary condition present at the cylindrical surface and approaches have been developed for canonical cases of perfectly conducting, scalar impedance, coated and dielectric surfaces [32], [36], [37].
In this work, to initiate the development of an RO framework for cylindrical metasurfaces, we propose for the first time, a UTD based RO approach to model a uniform purely reflective MS cylinder in 2D. The uniform cylindrical metasurface is described using both tangential and normal surface susceptibilities and conditions are provided for perfectly reflecting surfaces. Unlike a perfect electrical conductor surface (PEC) which depends on only the tangential field components, a fully reflective metasurface provides a very rich scattering behaviour. Such surfaces have been successfully modelled using the GSTCs with both tangetial and normal surface susceptibility components, and prove to be accurate when compared to volumetric unit cell simulations of the same surface [38].
The applicability of the homogenous surface susceptibility description of the metasurface for curved topographies has been studied in [39], where the GSTCs have proven to ac-curately capture the scattering behaviour of a physical curved metasurface. In addition to the novelty of developing the RO-GSTC model for the cylindrical metasurface, the proposed work naturally extends the previous works on UTD focusing on specific canonical structures only to generalized boundary conditions via the GSTCs.
The paper is structured as follows. In Sec. II the canonical problem will be described. First the GSTC formulation for the cylinder and the conditions for a purely reflective MS are provided. The appropriate eigenfunction (EF) solution for the MS cylinder is then developed. This analytical solution will be used to develop the RO approach and can be used as a base-line result for comparison. The final part of this section will present two alternative Greens' functions for the problem under analysis. The first representation of the Greens' function is based on an integral form in a complex plane. The second is based on a summation of eigenfunctions and is essentially the EF solution for the line source. These Greens' functions will be used to develop the UTD formulation for the transition region. Section III provides the UTD approach to the problem, describing the GO fields present in the lit region, the creeping waves in the deep shadowed region and developing the transition equations in terms of a Pekeris caret function. The approach to dealing with this integral for the MS is presented in Sec. IV. The integral is solved using a heuristic approach presented in [1] which results in an efficient and a pragmatic solution. Finally, in Sec. V, RO-GSTC simulations of three different cases are presented and compared with baseline results to verify their accuracy. The first two cases are for TM z (E z , H φ , H ρ ) polarization and the third for a TE z (H z , E φ , E ρ ) polarization. For the TM z example two surfaces are used; the first with only tangential χ components and the second with a normal component included. The TE z surface has both tangential and normal components present. In this section the RO results are verified by comparison to both the original EF solution and a boundary element method (BEM) simulator.

II. PROBLEM STATEMENT
The problem under consideration is shown in Fig. 1. A metasurface is configured into a circular cylinder of radius b centered around the z−axis. A line source, which produces either a TM z or TE z polarized electromagnetic field, is oriented parallel to the cylinder axis and located at ρ = (ρ > b, φ = π). This electromagnetic field interacts with the cylinder and produces a scattered field. This field is determined by the boundary conditions of the cylinder, as prescribed by the surface susceptibilities, and the position of the source.
A few simplifying assumptions will be made to make this problem more tractable. Firstly, the MS is assumed to be purely reflective throughout this work, i.e. fields inside the cylinder are identically zero. Transmissive cylinders on the other hand, such as the simpler case of a dielectric cylinder, experience more complicated field effects such as internal refraction and interior creeping waves due to fields being present in the interior of the cylinders [37]. Replicating these effects in a UTD formulation is involved and will be left for future consideration. Secondly, a metasurface may also introduce a polarization rotation to the scattered fields, which results in both TM z and TE z polarizations being present even with an incident field of a single polarization. This case, while straightforward to develop, is tedious to analyze due to the couplings between the two polarizations [40]. This complexity also involves finding appropriate conditions on the surface susceptibility tensors to ensure perfect reflection at all angles from the metasurface. Consequently, only the case of a non-polarizing metasurface is considered here for simplicity and to maintain the focus on the proposed method. The subsequent development is therefore limited to purely reflective, non-polarizing and a uniform cylindrical metasurface with spatially constant surface susceptibilities.

A. Generalized Sheet Transition Conditions (GSTCs)
When considered as a zero thickness structure, metasurfaces create a spatial discontinuity in both the phase and amplitude of an electromagnetic field. This discontinuity is described by the General Sheet Transition Conditions (GSTCs) and relates the difference and average of the fields across the metasurface via the following relations [13] [12]: where are the differences between the electric and magnetic fields on the positive (exterior) side and the negative (interior) side of the cylinder, respectively, and the normal vector (n) points from the interior of the cylinder to the exterior. In addition to this, and µ are the permittivity and permeability of the surrounding medium, and ω is the angular frequency of the fields. The operator ∇ (·) is the evaluation of the gradient in the contoured tangential plane of the surface and (·) ⊥ indicates the normal component of the field. The differences in the fields are related to the electric and magnetic dipolar surface polarization densities, P and M respectively, which are given by where are the average fields across the surface. In general, the surface susceptibilities χ ab , with a, b ∈ {e,m}, are full dyadic tensors, which may introduce an arbitrary polarization rotation to the scattered fields. For the case of a circular cylindrical metasurface, the problem may be conveniently formulated in a cylindrical co-ordinate system so that the surface's normal can be defined asn =ρ for the entire arc of the surface. Therefore in cylindrical coordinates, the surface susceptibilities take the form There are four sets of tensors, which when combined together gives a total of 36 terms. The number of unique terms can be reduced depending on the properties of the analyzed surface, such as reciprocity and energy conservation [41], [42]. In this paper, we are dealing specifically with metasurfaces that do not induce a polarization rotation in the scattered fields which, when combined with reciprocity, simplifies the tensors to [16] This gives a total of eight possible non-zero surface susceptibility terms which are complex in general for lossy surface.

B. Purely Reflective Metasurfaces
As in this paper we will be dealing with purely reflective MSs, we will now develop the conditions on the susceptibilities which ensure zero transmission across the surface for all angles of incidence. For a non-polarizing metasurface cylinder, the transmission and reflection coefficients, S respectively, for a plane wave incident on the exterior of the cylinder at an angle of incidence of θ are given by, [16] S In order to synthesize a non-transmissive surface, the transmission coefficient must be zero at all angles and for both polarizations. This can be done by careful removal of the angular dependence and then a simple rearrangement of (4). From the above equations, it can be seen that this occurs for two different cases: 1) When ζ {TE,TM} = 0 which requires, χ zz ee = χ zz mm = χ ρρ ee = χ ρρ mm = 0 and, This leaves only the tangential susceptibilities χ φφ ee , χ φφ mm as free and arbitrary. 2) When χ φφ ee = χ φφ mm = 0 which again requires, k and leaves χ ρρ ee , χ zz ee , χ ρρ mm and χ zz mm free and arbitrary. One can also synthesize a zero-transmission metasurface by applying ζ {TE,TM} → ∞; because the resulting reflection profiles as a function of θ are considered trivial, this case is not explored in this paper.

C. Eigenfunction Solution
The development of the ray optics formulation for the canonical problem starts from a rigorous solution of Maxwell's equations for the problem and the definition of an appropriate Greens function. In this section we will develop an eigenfunction (EF) solution for the MS cylinder. This solution helps formulate the problem for an RO analysis and can be used as a test of the RO implementation as it is a completely rigorous solution of the problem containing all scattered field components.
Any electromagnetic field must be able to be represented as a sum of eigenfunctions that satisfies the given boundary conditions of the problem. For the case of a two dimensional field (using e jωt time convention), the spatial dependence of the eigenfunction u(ρ, φ) satisfies the following partial differential wave equation: where k is the wavenumber of the field in the given region. By applying separation of variables, the radial and angular dependence can be separated. The angular dependence must be periodic such that u(ρ, φ) = u(ρ, φ + 2π). The radial dependence is characterized as a sum of either the Bessel functions of the first and second kind, J n (x) and Y n (x) respectively, or of the Hankel functions of the first and second kind (otherwise known as the Bessel functions of the third kind), H n (x) and H (2) n (x) respectively, defined as where n ∈ Z is the order of the function. While J n (x) and Y n (x) are most useful for standing waves, H n (x) and H (1) n (x) are used to describe radially outgoing and incoming waves, respectively. As such, a linear combination of two Bessel functions of any kind can be used to describe a given order's radial dependence. Because of the constant radius b of the cylinder under analysis, the eigenfunction for a given order can thus be written as where a n , b n are constant coefficients. There are two separate regions that are under consideration: 1 the case of the fields inside (ρ < b) of the cylinder and the fields outside (ρ > b) of the cylinder. The total field in either case can be represented as a sum of all of the eigenfunctions, with different coefficients depending on the boundary conditions. For the interior fields, the fields must be resolvable at ρ = 0; because H (2) n (kρ) becomes singular as ρ approaches 0, the eigenfunction of the interior fields must only consist of J n (kρ) terms. Meanwhile, the exterior scattered fields have to follow the radiation condition and go towards zero as ρ → ∞, which is satisfied by H (2) n (kρ). Therefore, the total field in the exterior region can in general be represented as a sum of J n (kρ) (exclusively for incident fields) and H (2) n (kρ) terms (for both incident and scattered fields).
Because there is no cross-coupling of the TE z and TM z fields for the considered metasurfaces, the effects of the two polarizations can be considered separately. For the case of the TM z polarization, the electric field in the z−direction in each region is given as a sum over all possible eigenfunctions, given by: where E + z /E − z is the field of the exterior/interior region. Note, for the case of purely reflective MS and an exterior source C n = 0 by construction.
To complete the eigenfunction solution for the case of a metasurface, the GSTCs have to be applied as the boundary conditions at ρ = b. As metasurfaces relate both the electric and magnetic fields across the surface, expressions for the magnetic field in both regions are also required. Applying Maxwell's equations to (10a) and (10b), these fields are given by: where the total magnetic field is given by H ± =ρH ± ρ +φH ± φ and the presence of a prime indicates a derivative with respect to the argument of the function. These expressions can be substituted into (1) and (2) and, after some rearranging, the relationships between the total field coefficients are given by: where the auxiliary variables Υ are given by For future use in the development of the ray optics, it is helpful to rewrite the B n coefficient in terms of a linear operator Q which operates on the Bessel functions J n and H (2) n [32], [35]. By rearranging the numerator and denominator of (12a), this operator for the TM z polarized field of the metasurface takes the form, These expressions give a compact final form for the total electric field outside of the cylinder, which is of the form In (15) the total field outside of the metasurface is calculated by a set of weighted Bessel and Hankel functions. Because H (2) n (kρ) represents a radially outgoing wave, these terms can be assumed to encompass the scattered fields from the metasurface cylinder. As a result, the J n (kρ) terms represent only the incident field interacting with the cylinder. If a known incident field interacting with the metasurface is decomposed into a series of J n (kρ) terms, it is possible to extract the relative weights of the unknown coefficients A n , which allows the calculation of the total/scattered fields using the eigenfunction method. For the case of a line source excitation at This solution ρ < x 0 is appropriate as the line source is exterior to the cylinder. This gives after some manipulation, It should be noted that this process can be repeated for the case of an incident TE z polarized field. For this instance, the q e (n) variable in (13) can be replaced by where the auxiliary variables β are given by

D. Green's Function Representation
An alternative representation of the fields scattering from a cylinder is through the use of the two-dimensional Green's function [1]. For the same canonical problem defined previously, consider the case of the cylinder excited by a line source of electric (TM z ) and/or magnetic (TE z ) currents parallel to the cylinder's axis at a location Q (ρ > b, φ = 0). Instead of performing a discrete sum of Bessel functions of order n, the Green's function is found by performing a contour integral over the continuous, complex order ν [43], [44]. The contour path is an infinite semicircle in the lower half of the ν-plane.
The authors in [1] use this formulation to produce a form of the Green's function that is the sum of characteristic radial Greens' function and a residue sum of the poles present in the lower half of the complex plane, resulting in an expression that is rapidly convergent for large values of kb but valid for any cylinder radius: where the notation ρ = (ρ, φ) indicates a radial and angular position, φ l = φ + 2lπ, and ν p are defined by d(ν p ) = 0. The above Green's function is suitable for all values of |φ−φ | ≤ π.
The function g ρ (ρ, ρ ; γ) is the radial characteristic Green's function, which is given as where ν = √ γ, ρ < stands for the smaller of the two radial terms ρ and ρ , and ρ > is the larger of the two terms. The prime on the sum over l denotes the omission of the l = 0 term in the summation. Due to the exponential dependence on the quantity |φ − φ l |, the sum over l converges very quickly. The summation over ν p represents the contributions of all of the poles ν p , and is dependent on the boundary conditions of the given cylinder through d(ν). This sum captures the creeping wave contributions. For large cylinders, this contribution is dominated by a few poles, which reduces the computational complexity.
A second form of the two-dimensional Green's function is given in terms of angular eigenfunctions [1]. These have a form similar to the total field eigenfunction solution of (15), and is expressed as: and It can be seen that G 0 (ρ, ρ ) is the free-space Green's function for the incident field and G s (ρ, ρ ) represents the scattered field from the cylinder. Both of the presented Green's functions (19) and (21)   on the surface and are equally valid for the problem being discussed. However, the two representations are useful as they represent the scattering system in two slightly different ways. Later on, these two separate representations will be used together in order to describe the complicated field behavior in the transition region around the shadow boundary.

III. UNIFORM THEORY OF DIFFRACTION (UTD) FIELDS
In ray optics the field from a set of scattering elements interacting with an incident field is not calculated using a rigorous full-wave description of the system. Instead, in the high frequency limit (k → ∞) the wavefronts of the fields scattering from the objects can be described using a ray description. By summing up the contributions of the discrete rays that converge to some position P , the total field can be obtained.
For the case of 2D simulations of planar metasurfaces, which was covered in-depth in [29], two types of critical points contribute to the scattered fields. Critical points of the first kind lie in the interior of the surface and are associated with reflected and transmitted specular fields. For metasurfaces with uniform surface susceptibilities, the locations of these scattering points on the surface satisfy Snell's law of reflection, and can thus be found using Fermat's principle [45]. Critical points of the second kind lay on the edges of the surface and give rise to edge-diffracted fields. These two types of ray interactions provide a uniform description of the field across the shadow boundary for the planar surface.
For a fully reflective cylinder there are four types of UTD rays: incident rays generated by the source, and scattered rays of three types -1) surface reflected, 2) those that generate the shadow; and 3) surface diffraction rays. Fig. 2 shows the general ray diagram of a set of rays interacting with a fully reflective cylinder. The incident ray and the shadow ray are straight-forward to determine: the incident ray can be determined from the source description and the distance s i to a point of interest; and the shadow ray is simply determined by extending the ray intersecting the surface by a distance s h and negating the field to cancel the incident field for regions This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2022.3184540 behind the surface [28] as would be the case for the points ρ s and ρ d .
The characteristics of the two scattered rays (reflected and surface diffracted) are calculated differently for the various locations of the field point, being determined by their relation to the incident field and the shadow boundary created by the cylinder. For Region I and III, the field is sufficiently far away from the shadow boundary and can be calculated using the standard General Theory of Diffraction (GTD) which includes the creeping waves generated by surface diffraction [35]. Therefore the total field in Region I is evaluated using the incident and reflected rays of geometrical optics, in addition to contributions from the creeping waves that occur when an incident ray interacts with the cylinder at a point of tangency (Q 2 in the figure). Because the metasurface cylinder is impenetrable, the shadow region of Region III does not receive contributions directly from the GO rays and the field is composed of creeping wave contributions only.
Region II has the most complicated behavior; if the field point is sufficiently close to the shadow boundary, the description of the fields using GTD does not converge towards the analytical solution, nor are the fields uniform across the shadow boundary. More rigorous analysis is required in order to provide an accurate description of the fields in this transition region. The remainder of this section goes through the equations required to calculate the total field in each of the regions discussed above. It should be noted that the formulas developed in the section are for fields far away from the surface; those close to the cylinder will not necessarily be accurate. Moreover, the GSTC description of the surface is only applicable for describing macroscopic (and not microscopic) fields. As we are interested in a large scale analysis of the metasurface scattered fields, the high frequency approximations are thus justified.

A. Specular Geometrical Optics Fields
For the situation under consideration the incident field is generated by a line source with a distribution prescribed by (16). The ray is characterized by the propagating radial field (E i (ρ)) and the distance travelled (s i ) to a point of interest: typically an intersection with the surface or the direct illumination of a point in the lit region. 2 An incident ray intersecting with the fully reflective metasurface cylinder at a point Q r (ρ c ) within the lit region gives rise to specular scattering components E s,a , with a ∈ {h, r}, whose scattering directions are given byk a . The case of a = h creates the shadow field for the cylinder, the scattering directionk h is the same as that of the incident ray's, and the contribution from this ray is equal in magnitude but opposite in sign to the incident field at the point of the detector, or in other words 2 A more general ray formulation of the incident field could be constructed by obtaining the field over a source surface and constructing a set of rays described by an intensity, direction and curvature.
where ρ i is the radius of curvature of the incident ray at the point of the surface, and s h is the distance from the critical point to the detector at ρ h = ρ c + s hki . This field cancels out the incident field both within the cylinder as well as behind the cylinder, creating the shadow region.
If sufficiently away from the shadow boundaries, the case of reflection (a = r) can also be described using GO rays. The scattering direction of the reflected specular rays is determined by Snell's law of reflection, which mathematically is evaluated using The field amplitude corresponding to the specularly reflected ray is evaluated as, where R(ρ c ) is the reflection coefficient, and ρ r is the radius of curvature of the reflected wavefront. The reflection coefficient is dependent on the metasurface parameters, the polarization of the incident ray, as well as the angle of incidence between the surface and the incident ray, the latter of which is given by θ i = − cos −1 (n ·k s ). The calculation of R has been discussed in various literature, and was studied in a ray optics framework in [29]. For planar surfaces, except when a modulation of the local surface parameters occurs [28], the term ρ r would be the same as that of the incident ray's. However, the addition of a smooth curved topography causes a shift in the location of the reflected ray's caustic, resulting in ρ r taking the form [44], [46] For the case of a general metasurface there is also a transmitted scattered field component: this would take the form of scattered rays that would travel in the interior of the cylinder, repeatedly scattering off of the concave interior and emitting further scattered fields outside. While the GO component could be modeled in a straightforward fashion, the contributions of other scattering effects such as the surface diffracted field is more complicated. Therefore the metasurfaces that will be studied here are limited to the case of non-transmissive metasurfaces for simplicity.

B. Creeping Waves
When an incident ray interacts with a convex surface at a point of tangency Q 1 , or in other words with an incident angle of |θ i | = π/2 (see Fig. 2), some of the energy of the incident ray excites a surface ray that travels along the arc of the surface. Originally described by Keller in [31], [35], these surface rays emit a portion of their energy while following the convex surface in the form of diffracted rays, which travel tangentially from where they leave the surface, at Q 2 . If the surface that the excited ray is traveling on is a closed convex surface the ray will continue to loop around the entirety of the surface, shedding rays with progressively lower energy until it is considered negligible.
These "creeping waves" have multiple modes that can be excited, each of which contain a different proportion of the energy from the initial incident ray and each with a different azimuthal propagation constant. The waves are captured in the Greens function representation given by (19) with the creeping wave modes characterized by poles in a plane described by the complex Hankel order ν. Each mode corresponds to a root ν p of the characteristic equation of the boundary conditions, d(ν), and the contribution of a single ray for a given mode is given by [1]: where t is the arc length between Q 1 and Q 2 that the ray has travelled on the surface, s d is the distance from where the ray leaves the surface to the field point ρ d , and D p is the diffraction coefficient, which is given by Due to the cylinder under analysis being a uniform circular cylinder, the principle of reciprocity states that the diffraction coefficients at the points of incidence and emission are the same. The total field from the creeping waves resulting from a single point of diffraction is the sum of all modal contributions and those rays with multiple encirclements of the surface. Mathematically, this is given by For most applications, the dominant mode of the creeping waves is sufficient, as the relative values of D p decreases the higher the mode is. Similarly, the e −jνpt/b term causes significant attenuation of the surface ray as it travels around the cylinder, which makes large values of m corresponding to m full encirclements produce a negligible contribution.
The approximations in the literature used to derive the above creeping wave formulation do introduce some limitations. This ray picture of the surface diffracted fields is most valid for plane-wave incidence and far-zone scattering [1], [31]. When the field point is close to the shadow boundary, however, this formulation breaks down; for this situation, a more sophisticated evaluation of the fields is required.

C. Transition Region
In the previous sections the Geometrical Optics solutions to the lit and shadow regions were discussed. The GTD solution, which encompasses the specular reflected and creeping wave rays, is valid for most regions, however it breaks down when sufficiently close to the shadow boundary, otherwise known as the transition region. This region, which spans an angular range on either side of the shadow boundary on the order of (2/kb) 1/3 [32], requires a more rigorous evaluation of the fields.
The asymptotic high-frequency analysis of the canonical problem in this transition region has been studied in several papers, most notably by Pathak in [32], who outlines the process for calculating these fields for a plane wave scattered by a perfectly conducting circular cylinder. This process has since been extended for other boundary conditions, such as impedance boundaries [1] and dielectric coated cylinders [36], as well as for non-circular cylinders with incident fields that are not plane waves [45].
There are two separate regions in the transition region: the lit region and the shadow region. The evaluation of the fields in either transition region differ from one another, but provide a uniform transition of the fields across the shadow boundary. The requisite equations are given below.
1) Lit Region: For a field in the lit transition region, the surface diffracted field takes the form of a reflected field. The UTD expression for this reflected field at an observation point ρ l is given by where the reflection point Q r on the cylinder is determined using the law of reflection in the same manner as the specular GO fields. The reflected ray caustic ρ r and the reflected path distance s r are also calculated in the same way. The main difference between the GO rays and fields in the lit transition region is the generalized reflection coefficient, which is given by [32], [45] where F (x) is the UTD transition function [32], which is used commonly in standard edge diffraction [28] and whose argument is given by X p = 2kL p cos 2 θ i , where L p = (s r s i )/(s r + s i ) is the distance parameter, s i = ρ i is the distance from the presumed cylindrical source to the diffraction point, is the form of the Fock parameter ξ for the lit transition region, and m = (kb/2) 1/3 is the curvature parameter. Because m is always positive for convex surfaces, the parameter ξ = ξ p is always negative in the lit transition region. Both X p and ξ p are larger when away from the shadow boundary and are zero at the shadow boundary. Equations (29) and (30) are a generalization of (24) which was formulated in terms of a local reflectivity of the surface. In the deep lit region for sufficiently large cylinders (29) and (30) will converge to the GO result predicted by the direct use of (24). The use of (24) when appropriate can be of practical use as the computational cost will be significantly lower than using the generalized formulation.
The functionP g (ξ) is the Pekeris caret function 3 , otherwise known as the Fock-Type integral function, which is denoted as such for notational consistency. The function is parameterized by the boundary conditions of the cylinder. Evaluation of this function must be handled with care and will be discussed specifically in the subsequent section.
2) Shadow Region: The surface diffracted field for a field point ρ d in the shadow transition region is given by The field point ρ d , the incident ray diffraction point Q 1 and the distance s d are defined in the same way as for the creeping waves; the only change in these parameters is that ρ d can be any point within the shadow transition region, even on the shadow boundary itself, provided that it is some distance away from the surface. The UTD surface diffraction coefficient T g is given by [32], [45] is the distance parameter, s i is the distance from the source to Q 1 , and ξ d = mt b is the Fock parameter associated with the shadow zone. The parameters t and s d are defined in the same way as in the creeping waves. Similar, but opposite, to the lit transition region, the terms X d and ξ d are always positive in the shadow region and are zero at the surface boundary.
To summarize, the process of calculating the ray contributions anywhere around the metasurface is depicted in Fig. 3 for a number of illustrative situations: 1) Lit or Lit Transition: If the detector, i.e. a field observation point, is in the lit or the lit transition region, the total field is calculated as the sum of the incident ray (16) at the detector, the reflected ray resulting from the critical point through (24) or (29), and the creeping wave contribution from the diffraction point on the opposite side through (26). 2) Shadow Transition: If within a shadow transition region the total field is the sum of the surface diffracted field from the closest diffraction point plus the creeping wave contribution from the other diffraction point, using (31) and (26), respectively. 3) Deep Shadow: If within the shadow boundary, but not within a transition region, the total field is the sum of the creeping wave contributions (26) from both diffraction points. This procedure can be simply extended to find the fields in any 2D region and construct the scattered fields in a desired region of space.

IV. TRANSITION INTEGRAL
A. Evaluation of the Fock-Type integral functionP g (ξ) The remaining task to complete the RO formulation is to evaluate the Fock-Type integral functionP g (ξ) which depends on the boundary condition of the convex surface as well as the input argument ξ, which was shown in the previous section to be related to the relative proximity of the field point to the shadow boundary. If the characteristic operator of the surface is given by (13), then the Fock-Type integral function is written as [32] where q e,m (τ ) = mq e,m (ν). The functions V (τ ) and W 2 (τ ) are the Fock-type Airy functions, and the transition integral p * (ξ), which factors out the singular term found at ξ = 0, is commonly used. The relationship between the order ν and the integration variable τ is through a change of variables ν = kb+mτ , which is required to relate the Bessel functions to the Fock-type Airy functions under certain approximations [32]. A direct evaluation of this integral is numerically difficult. While the integrand tends towards zero as τ → ∞, for negative τ values the integrand oscillates with a variable amplitude, and these oscillations only increase in frequency as τ → −∞ with no reduction in average amplitude. A number of approaches to evaluate this integral for a variety of surfaces have been reported. One of the more consistent methods is to deform the contour of the integral in the complex τ plane to perform the integral over discrete intervals, which is used for different types of impedance cylinders [47], [48]. This method is most accurate with surfaces that have a constant q e,m value, which is not the case for metasurfaces, and using the method for non-constant values of q e,m results in erroneous calculations of the fields. These errors become apparent in the lit region, where the fields calculated using (29) do not properly converge to the specular fields for large negative values of ξ p .
More rigorous approaches of deforming the contour path of the integral do exist [49], but are tedious to implement for a general metasurface. Without a detailed analysis of the poles of the surface integral, and an appropriate deformation of the contour integral path to account for these poles, it is challenging to compute the integral directly. Considering these complexities, a simpler and more heuristic approach to calculate the transition integral will be used instead, which will be shown to work well for the metasurfaces considered.
B. Heuristic Approach to e −jπ/4 p * (ξ, q) The approach found in [1] to calculate the transition integral is based on using a Green's function representation of the system. The two alternative representations of the Green's function, (19) and (21), are used together to obtain The above Green's function contains all the geometricaloptical, transition, and dominant diffraction effects [1], except for the creeping waves which have been removed, and is dependent on the cylinder radius b, the metasurface's surface susceptibilities, and the pair of field and source reference positions, ρ and ρ respectively. In the lit transition region the field without the creeping waves present is also prescribed by (29) and in the shadow transition region by (31). Using (34) the expressions (29) and (31) can be inverted to solve for the transition integrals, respectively: 35a) is used for field points in the lit region of the transition region and ξ is determined by the angle of incidence θ i . On the other hand, in the shadow transition region (35b) is used and ξ is related to the arc length t.

A. Numerical Setup
In order to test the accuracy of the formulation for the surface diffracted fields, a UTD simulation (i.e. RO-GSTC) was developed in MATLAB and compared to an equivalent 2D Boundary Element Method (BEM) simulation based on an Integral Equations Solver (i.e. BEM-GSTC) [10], [11]. The UTD results are also compared to the eigenfunction solution (i.e. EF-GSTC) given by (15), which represents the actual rigorous field solution of the problem.
Each simulation has three primary components: the surface itself, one or more sources, and multiple field detectors. The metasurface S of constant radius b is placed in the ρ − φ plane centered at the origin, and its surface susceptibilities χ ab are known and uniform across the surface. Cylindrical waves are radiated from a source Tx at (ρ , φ = π), and interact with the surface; to simulate the effect of a plane wave incidence, the line source is moved far away from the origin (ρ → ∞) such that the cylindrical waves appear planar at the surface. In order to capture the scattered field as a function of angle, several numerical field detectors are placed evenly around the metasurface at locations {ρ, φ} with an angular separation of ∆φ between detectors. Based on each detector's position relative to the shadow boundaries (see Fig. 3), the evaluation of the field at the detector is calculated using a summation of the scattering and diffraction effects.
To analyze the scattered fields from the metasurface cylinder, the detectors are first divided into those that are directly illuminated and those that are shadowed by the cylinder. For each detector in the lit region, the specular reflection point ρ c is calculated by applying Fermat's principal of least time. The transition integral is calculated based on the approach outlined in Sec. IV-B for the desired metasurface and polarization being analyzed. For the synthesis, a source reference position ] is a set of equally spaced discrete angular points. Once the transition integral p * (ξ) is found for a range of ξ values, the surface diffracted field for a given diffraction point can be calculated by first determining ξ based on its position relative to the shadow boundary and then using linear interpolation to find p * (ξ) utilizing the numerically calculated values.
For points in the transition or shadow region the two diffraction points on the cylinder are found by determining the incident points of tangency in relation to the source. The scattering points of tangency for each detector is then calculated, which are used to determine the ξ parameter for both diffraction points for each ray. The total field is then calculated for each detector based on their positions relative to the two shadow boundaries. Fig. 4 presents the transition integral as a function of ξ for the three example surfaces considered below to illustrate the proposed method. Two surfaces are illuminated by a TM z polarization; the first surface has only tangential susceptibility components (TM z, ), however, for the second a normal component is present (TM z,⊥ ). The third surface is for a TE z polarization and uses both normal and tangential susceptibility components. Details on the surface descriptions are found in the next section.
Shown in the first plot is the magnitude of p * (ξ) for the three surfaces with the transition integral synthesized using ρ ref = 2b and ρ ref = 3b. As can be seen the procedure outlined in Sec. IV-B provides a clean and smooth transition from the lit to shadow region. In the second plot (Fig. 4b) the variation in p * (ξ) is presented for three different cylinder radii for the TM z,⊥ surface. Again we see a smooth transition from then lit to the shadow region. These two plots show the need to characterize p * (ξ) for a specific surface configuration defined by the susceptibilities and curvature. The final plot shows the independence of the method to the choice of the synthesis parameters ρ ref and ρ ref by plotting p * (ξ) for three different configurations. This last result is important as it shows that the transition integral is insensitive to the nature of the incident field. The overall accuracy and applicability of the method will be addressed in the following sections.

B. Results
To illustrate the proposed formulation and address its accuracy, several examples will be presented. In all cases, the simulations are performed in free-space at a frequency of f = 60 GHz with a metasurface radius of b = 0.03 m, i.e., approximately 6λ or ≈ 38λ in circumference. To validate the computed field solution, the results are compared to two different methods: the first method is the eigenfunction solution of the fields, which was presented in Sec. II-C. This method is a rigorous analytical method only limited by the number of expansion terms taken and accurately captures all of the scattering effects from the metasurface interacting with a given incident field, and is thus assumed to be a base-line solution.
The second method that is used for comparison is a well established 2D BEM-GSTC numerical tool [10], [11], which rigorously solves for the scattered fields from arbitrarily shaped conformal metasurfaces with a specified surface susceptibility profile. This simulation tool has been verified against both analytical approaches and full-wave simulations of detailed unit-cell level simulations from commercial simulators [23], as well as for surfaces with modulated parameters that affect the curvature of the scattered fields [29]. The primary source of error in the BEM simulations is the level of surface discretization and a surface mesh density of 30 divisions per wavelength was found to be sufficient for the fields external to the cylinder for proper convergence. Figure 5 shows the scattering profile of the three metasurfaces that will be considered in this section, as a function of the angle of incidence of the illumination. These scattering profiles are calculated using the synthesized susceptibilities with (4a) and (4c). As seen from the amplitude and phase profiles, the scattering behavior of each metasurface varies with the angle of incidence of an interacting ray, and differs considerably between the examples.
We will now proceed with a detailed comparison for the three surfaces.  1) TM z Case using Tangential χ: The first example is the case of a reflective metasurface whose susceptibilities feature only tangential components and is case 1 in Sec. II-B. For this case we set χ zz ee = χ zz mm = χ ρρ ee = χ ρρ mm = 0 and the bianisotropic terms are prescribed to be χ zφ em = −χ φz me = 2j/k 0 . We have two free and arbitrary parameters χ zz ee and χ φφ mm .
Let us start with a metasurface that is characterized by a finite transmission and reflection at normal incidence. The surface susceptibilities of a planar surface with a transmission coefficient of T and a reflection coefficient of R at normal incidence can be synthesized using the method from [15], and To make the metasurface fully reflective at all angles, the remaining surface susceptibility terms were chosen following Sec. II-B. The values for the constructing scattering coefficients were chosen arbitrarily as T = 0.5j and R = 0.5 √ j to provide a surface with a rich reflective response (See Fig. 5). These synthesized values provide realistic values of the surface susceptibility terms assuming plane-wave incidence, and thus the reflectance is expected to be different when the surface is conformed to a cylinder and excited with an arbitrary incident field. Figure 6 shows the computed total fields for this case in the presence of a line source at the exterior of the surface. Shown here are the results for the RO and EF simulations, where the RO results are presented as two separate sets of data: 1) the contribution of only the GO-reflected, incident, and shadow rays, and 2) the surface diffracted rays in addition to the previous contributions. These results are normalized with respect to the incident field at the leading edge of the metasurface cylinder (ρ = (b, π)).
In Fig. 6a) which presents the RO fields over a 2D region, we can see the importance of the surface diffracted fields, which causes fields to diffract into the shadow region behind the cylinder and provides a smooth transition between the lit and shadow regions. Figure 6b) then compares both of these fields to the EF solution on a circle of radius 3b. Here it is seen that the GO fields alone are sufficient to model the scattering behavior of the cylinder when in the deep lit region. However when close enough to the shadow boundary the GO solution alone is not enough, and the complete RO solution matches the expected results. When looking in the deep shadow region, the interference pattern that is caused by the two diffraction points of the cylinder is the same for the RO and EF simulations, which further confirms the formulation's accuracy. The same simulation is then compared to the BEM, and the results are compiled in Fig. 7. To illustrate the effect of the MS compared to more traditional boundary conditions, an EF simulation for a PEC surface with the same cylindrical dimensions is also shown [43]. In Fig. 7c) we can see that the BEM simulation, which is formulated using a different methodology, matches the scattered fields of the RO simulation very well exterior to the cylinder. There is a bit of a discrepancy in the fields within the deep shadow region, however this occurs below -40dB and can be considered acceptable with the current discretization of the BEM. As can be further seen, the BEM solution for the interior of the cylinder is not precisely zero, which it is by construction for both the RO and EF formulations. This is due the quantization of the surface currents and was found to reduce for finer surface elements.
There is also a small discrepancy between the EF solution and the RO solution in the deep shadow region right behind the metasurface (Fig. 7 a and b). This is due to the presented UTD solution not including the formulation for fields close to the surface of the cylinder.
2) TM z Case using Normal χ: The second example utilizes a metasurface that has normal susceptibilities in addition to the tangential susceptibilities (case 2 in Sec. II-B). For this surface we have χ zφ em = −χ φz me = −2j/k 0 and we choose the tangential components χ zz ee and χ φφ mm by setting T = 0.5j and R = 0.5 √ j. In addition we have two other free susceptibilities to specify and we set χ ρρ ee = −10χ φφ ee , and χ ρρ mm = 3χ φφ mm ; all other susceptibility parameters are 0. This produces the reflectivity response shown in Fig. 5. Figure 8 shows the simulation results of the metasurface  cylinder interacting with a TM z polarized plane wave travelling in the +x direction. Figure 8a) presents 2D field depictions for the RO and BEM simulations showing a close match in the lit and transition regions. In Fig. 8b), we can see that the RO simulation matches very well with the EF solution for the entire angular extent 2b away from the surface of the cylinder. In the inset, which gives a better look at the interference pattern in the shadow region, the BEM simulation diverges slightly from the RO/EF solution. This can be seen more clearly in Fig. 8a), where a small leakage field passes through the, in principle, non-transmissive surface, causing a resonance in the interior and subsequently distorting the region behind the metasurface. As before this effect can be reduced by increasing the discretization of the BEM surface. It should be noted that this example has used a planewave excitation (which is very different than the excitation used to characterize the transition integral p * (ξ)) and confirms the robustness of the method described in Sec. IV-B with respect to the incident field.
3) TE z Case using Tangential and Normal χ: The final example is the same metasurface as the second example, instead excited with a TE z polarized line source at the exterior of the surface. Fig. 9a) shows the simulation results comparing the 2D fields of the RO simulation with the BEM. Similar to the previous examples, the RO matches the BEM quite well visually in the lit and transition regions, whereas there is some discrepancy when observing the diffraction pattern behind the cylinder. This is further emphasized in Fig. 9b); here the RO simulation results match the EF solution very well for the entire arc of the circle, including the shadow region. However, the BEM results diverge more apparently for this example in the shadow region, where the calculated fields are small (less than -20 dB) compared to the fields in the lit region. qaw

VI. CONCLUSION
A ray optical methodology based on the uniform theory of diffraction has been proposed to model electromagnetic field scattering from uniform curved metasurfaces with both normal and tangential surface susceptiblities. This work may be seen as an extension of the previous work in [29] on the creation of ray optical infrastructure for metasurface systems. That work only dealt with planar surfaces and consequently This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2022.3184540 only specular reflection and edge diffraction. The work in this paper develops a methodology for the incorporation of curved surfaces and the resultant need to incorporate "creeping wave" diffraction and transition effects across the shadow boundary. The paper addresses the issue of illumination of a purely reflective uniform cylindrical MS illuminated by a line source, dealing with both polarizations, but not a MS that produces cross-polarized reflections, for simplicity. This canonical problem has a long history in UTD development for a wide variety of surfaces, and the extension to a MS provided here exploits a methodology previously used for cylinders coated in thin dielectric layers [1].
The formulation has been fundamentally based on a representation of the MS using the GSTCs which characterizes the surface in terms of susceptibilities. After establishing the conditions for a fully reflective surface the characteristic eigenfunction description of the MS problem has been developed. This rigorous analytical solution to the problem is then used to develop a RO-GSTC description of the scattered fields; including the specular GO field, surface diffraction described by creeping waves and a transition region over the shadow boundary. The specification of the fields in the transition region is dependent on the evaluation of the Pekeris caret function integral -a traditionally challenging problem. The approach taken in the paper follows [1] using a pragmatic method which allows for a smooth transition across the shadow boundary.
It should be noted that unlike the EF solution, the proposed RO-GSTC method can be integrated inside a more sophisticated RT framework where multiple scattering objects may be present. Traditional RO tracing has been applied extensively to electrically large scattering problems, and commercial software developed, as there was a clear need for methods able to handle problems that were too intensive for full-wave simulations [50]. In this paper a well established BEM fullwave simulation method was also used as a base-line for comparison of the predicted fields. However, the use of the BEM method for large problems results in the solution of a very large sparse matrix and even 2D simulations would be restricted to problems not significantly larger than presented here -primarily due to memory limitations. On the other hand, traditional ray tracing has been used successfully for much larger problems [51] and the methodology presented here would be appropriate for these scale lengths.
In conclusion the approach presented here is found to work well for the reflective MS considered here. The effectiveness of the RO approach is confirmed in the final section where three surfaces are used to demonstrate that the RO-GSTC fields very closely match the analytical solution over the entire domain of deep lit to deep shadow. The fields are also compared to an integral based numerical approach and a good match is found. The addition of surface diffraction effects to the RO-GSTC framework for a fully reflective metasurface will provide a basis for the continued development of ray optics approaches incorporating metasurfaces. The need for very large scale simulation of metasurface systems within complex environments motivates the extension of this work to nonuniform curved surfaces and transparent metasurfaces, which will be a topic of future investigation.