Surface Susceptibility Synthesis of Metasurface Holograms for creating Electromagnetic Illusions

A systematic approach is presented to exploit the rich field transformation capabilities of Electromagnetic (EM) metasurfaces for creating a variety of illusions using the concept of metasurface holograms. A system level approach to metasurface hologram synthesis is presented here, in which the hologram is co-designed with the desired object to be projected. A structured approach for the classification of the creation of EM illusions is proposed for better organization and tractability of the overall synthesis problem. The deliniation is in terms of the initial incident (reference) illumination of the object to be recreated (front/back-lit), the position of illusion (posterior/anterior), and the illumination used to create the illusion (front/back). Therefore the classification is based on the specific relationship between the reference object to be recreated, the observer measuring the object, the orientation and placement of the reference and illumination field, and the desired placement of the metasurface hologram creating a virtual image. In the paper a general design procedure to synthesize metasurface holograms is presented based on Integral Equations (IE) and Generalized Sheet Transition Conditions (GSTCs), where the metasurface hologram is described as zero thickness sheet with tensorial surface susceptibility densities. Several selected configurations are chosen to illustrate various aspects of the hologram creation in 2D, along with a novel numerical technique to artificially reverse propagate the scattered fields, required in the synthesis process. Finally, the impact of the metasurface size and the illumination field strength on the quality of the reconstructed scattered fields is also discussed.

Abstract-A systematic approach is presented to exploit the rich field transformation capabilities of Electromagnetic (EM) metasurfaces for creating a variety of illusions using the concept of metasurface holograms. A system level approach to metasurface hologram synthesis is presented here, in which the hologram is co-designed with the desired object to be projected. A structured approach for the classification of the creation of EM illusions is proposed for better organization and tractability of the overall synthesis problem. The deliniation is in terms of the initial incident (reference) illumination of the object to be recreated (front/back-lit), the position of illusion (posterior/anterior), and the illumination used to create the illusion (front/back). Therefore the classification is based on the specific relationship between the reference object to be recreated, the observer measuring the object, the orientation and placement of the reference and illumination field, and the desired placement of the metasurface hologram creating a virtual image. In the paper a general design procedure to synthesize metasurface holograms is presented based on Integral Equations (IE) and Generalized Sheet Transition Conditions (GSTCs), where the metasurface hologram is described as zero thickness sheet with tensorial surface susceptibility densities. Several selected configurations are chosen to illustrate various aspects of the hologram creation in 2D, along with a novel numerical technique to artificially reversepropagate the scattered fields, required in the synthesis process. Finally, the impact of the metasurface size and the illumination field strength on the quality of the reconstructed scattered fields is also discussed.

I. INTRODUCTION
Electromagnetic (EM) metasurfaces are 2D arrays of subwavelength resonating particles, where control of the spatial distribution and EM properties of the individual particles allows the scattered fields to be engineered with unprecedented control of both reflection and transmission, and with complete polarization control [1], [2]. Metasurfaces have been proposed for achieving a variety of wave transformation functionalities including -cloaking, waveform generation, and lensing [3]- [7], This has resulted in many proposed thin sub-wavelength surface devices across the EM spectrum, with a large number of structural designs and topologies using metals, dielectrics and other exotic materials [8], [9]. There has been also recent This work was supported and funded by the Department of National Defence's Innovation for Defence Excellence and Security (IDEaS) Program.
T. J. Smy, S. A. Stewart, and S. Gupta, are with the Department of Electronics, Carleton University, Ottawa, Ontario, Canada. Email: TomSmy@cunet.carleton.ca development in the application of metasurface concepts to system-level applications, such as next generation wireless networks for 5G/6G, where such surfaces act as smart reflectors in large Radio Frequency environments to achieve intelligent wave propagation [10]- [12]. In such cases, the metasurfaces act as part of an overall system and therefore must be co-designed and integrated with the rest of the system components.
One such application is an Hologram. Holograms are wellknown in optics where the spatial (and possibly temporal) information of an arbitrary object is encoded onto the surface (typically photographic plates) [13], [14]. This is a two step process, where the information about scattering properties of an object of interest is first recorded using a given reference beam and modulated onto a given surface. Once the information is recorded, the encoded surface, when illuminated with a reconstructing beam, projects an illusion of the object. With increasing sophistication of encoding capability, more complex illusions can naturally be created. This process essentially exploits the wave-transformation capability of a surfacemanipulating the reference beam electromagnetically. Metasurfaces therefore naturally represent a powerful platform to create sophisticated EM illusions with complete control over the scattered fields with respect to both complex amplitude and polarization. Consequently, these surfaces can also be used to create metasurface holograms, due to their advanced information encoding capability [15], [16].
In this work, a systematic description of metasurface holograms is presented and rigorous procedures are defined to design and synthesize these metasurfaces for achieving a desired EM illusion. To enable a general treatment of this problem, practical metasurfaces can conveniently be modeled as zero thickness sheets characterized using frequency dependent electromagnetic surface susceptibility tensorsχ(ω) [17]- [20]. The EM fields around the metasurface then can be described using Generalized Sheet Transition Conditions (GSTCs) [21]. The spatial distribution of surface susceptibilities of the metasurfaceχ(r) dictates the scattered (and thus total) fields produced by the metasurface when illuminated by an incident field. Therefore, the key design objective in creating metasurface based illusions is to synthesize the spatially varying surface susceptibilities,χ(r), to project the desired scattered fields at a given design frequency, ω to the observer location, identical to that of a real object.
Many metasurface synthesis and analysis problems using surface susceptibilities have been reported in the literature. In typical frequency-domain metasurface synthesis procedures, arbitrary incident and desired scattered fields may be spec-arXiv:2002.07910v1 [physics.comp-ph] 18 Feb 2020 ified, which in conjunction with the GSTCs, can be used to numerically solve, or optimize, for the required surface susceptibilities. For planar surfaces, surface susceptibilities can be analytically computed, for instance see [18], [22]. When performing a synthesis the fields must be specified at the metasurface location and not anywhere in space. While such methods are general in nature, their efficacy rests on a physically meaningful specification of the EM fields. On the other hand, metasurface analysis typically involves numerical computations where the GSTCs are coupled into bulk Maxwell's equations using a variety of standard numerical techniques based on Finite-Difference and Finite Element methods [23]- [25], and Integral-Equation (IE) techniques [22], [26]- [30]. Given that the field scattering from a metasurface hologram may need to be evaluated for electrically large domains, IE-GSTC methods are computationally efficient choices and, as will be shown later, are well suited for metasurface synthesis for holographic applications.
Given this context, a general methodology of synthesizing and designing metasurface holograms is presented in this work using the IE-GSTC method from a system level perspective. The metasurface hologram design problem is systemically defined for various geometrical relationships between the desired object illusion, a reference beam, the observer location and the illuminating beam. The 2D IE-GSTC based numerical platform is further developed to synthesize metasurface susceptibilities with an integrated approach, where the desired fields, specified anywhere in space and not necessarily at the metasurface, are generated using a system level description and fed-back into the metasurface design. A novel wavepropagation technique (a numerical inverse propagation) is further proposed for accurate determination of the desired scattered fields in certain specific scenarios. This ensures a physically meaningful field specification for arbitrary shaped metasurface designs and avoid ill-posed metasurface synthesis problems. Consequently, the proposed IE-GSTC based numerical platform computes spatially varying surface susceptibilities,χ(r) and is demonstrated to project complex EM illusions in a variety of physical scenarios.
The paper is structured as follows. Sec. II defines the overall problem of metasurface holograms from a functional pointof-view in the context of observer and field illuminations, along with basic properties of the general EM metasurfaces. Sec. III presents the IE-GTSC approach: developing the field equations to be applied to various hologram situations, the discretized form of the proposed IE-GSTC field solver, and the numerical procedure for surface susceptibility synthesis. Sec. IV shows several results demonstrating the metasurface hologram operation. Sec. V discusses some practical aspects of metasurface hologram designs such as the effect of finite-size surfaces and the field illuminations strengths on reconstructing the desired scattered fields. Finally, conclusions are provided in Sec. VI.

II. METASURFACE HOLOGRAMS A. Principle of Creating Illusions
Consider an object of arbitrary shape and material composition (dielectric or metal), subject to an incident refer-ence wave, ψ 0 (r, ω), as shown in Fig. 1(a). 1 The reference wave interacts with the object producing propagating scattered fields, ψ s (r, ω). We now place an observer located to the left of the object measuring the scattered fields ψ s (r 0 , ω) within a certain field-of-view. Assume that the incident fields are propagating from left of the object, so that from the perspective of the observer, it appears as a Front-Lit Object. By measuring and analyzing the incoming left propagating scattered fields, the observer perceives (or detects) the presence and properties of the object, such as its geometrical shape and material characteristics. If, on the other hand, the reference wave illuminates the object from the right side of the domain, the observer measures both the scattered and the reference fields (as both are left propagating), including the shadow produced by an object. This is termed as a Back-Lit Object.
Let us now remove the object, and introduce an artificial surface, i.e. a metasurface hologram, as shown in Fig. 1(b) at r = r m . The metasurface is excited with an illumination field, ψ i (r, ω), which may or may not be the same as the reference field ψ 0 (r, ω). Can this surface be engineered to rigorously recreate the scattered fields produced by the original object ψ s (r 0 , ω) within the field-of-view of the observer? In such a case, the real physical object or a virtual object produced by an artificial metasurface are indistinguishable from the perspective of the observer. This creation of a false perception by the metasurface hologram will be referred to as an Electromagnetic Illusion.
There exist several possibilities in the placement of the metasurface hologram, which while irrelevant from the perspective of the observer, is important from the practical design point of view and impacts how the metasurface may later be synthesized. If the metasurface is located between the object and observer, we refer to it as a Posterior Illusion, otherwise, when the metasurface is behind the illusion, we refer to it as an Anterior Illusion. In case of a posterior illusion, the virtual object is formed behind the metasurface, while in case of anterior, it is formed in front of the metasurface. Another variable of importance is how the metasurface is excited, i.e. the relative location and direction of the Illuminating Field. The illuminating field, in general, could be entirely different from the incident field that was used to determine the desired scattered fields from the real object. If it strikes the surface from the same side as the observer, it is referred to as Front Illumination, else, as Back Illumination.
With this description, the various illusion scenarios may be termed as: Front/Back-Lit Posterior/Anterior Illusions with Front/Back Illumination. The front/back-lit configuration determines which fields the observer detects -either scattered fields or total fields. The Posterior or Anterior position of the metasurface determines how the desired fields are to be constructed to form the virtual image (with a finite physical separation from the surface) and are numerically propagated to the metasurface location. This is important since the metasurface synthesis requires fields infinitesimally close to the surface region. It will be shown later that unlike the posterior   configuration anterior illusions require an unusual inverse field propagation as an initial step before metasurface synthesis can be performed. Finally, the front vs back illumination choice will determine whether a fully reflective metasurface is required or a transmissive one. This has important implications for practical realization of the synthesized metasurface, where compared to reflective ones, the transmission-type metasurface requires both reflection and transmission control.

B. Metasurface Descriptions
In order to generate arbitrarily complex and fully-vectorial scattered fields ψ s (r 0 ) from another equally arbitrary illuminating field ψ i (r), the metasurface hologram located at r = r m = r 0 must be capable of general EM wave transformations with complete wave control. Also, it should be noted, that while the prescribed fields may be arbitrary and complex, they are completely physical, making the metasurface synthesis problem a well-posed physically meaningful problem. This is due to the fact that the desired scattered fields are first computed using physical objects under well-defined incident fields, and the metasurface is illuminated with a physical field description as well.
The general wave transformation capability of physical EM metasurfaces can be described by expressing them as mathematical space discontinuities (zero thickness interfaces), owing to their their sub-wavelength thick nature and characterizing their EM wave interaction using electric and magnetic polarizabilities. To correctly model zero-thickness sheets, [31] introduced Generalized Sheet Transition Conditions (GSTCs) which were later applied to metasurface problems by [21], [32]. GSTCs (frequency-domain) relate the tangential EM fields around the metasurface to its tangential and normal surface polarization response as, wheren is the normal vector to the surface, ∆ψ T = (ψ + − ψ − ) is the difference between the fields across the surface; {M T , P T }, {M n , P n } are the average tangential and normal magnetic and electric surface polarizability densities of the surface. The polarization densities can be seen as a response of the metasurface to the average fields, related through the surface susceptibility densities as [17], where ψ av = (ψ + + ψ − )/2 is the average field at the metasurface, andχ is a general 3 × 3 tensor accounting for various microscopic EM characteristics of the metasurface. Eqs. 1 and 2 together rigorously model the EM interaction with the metasurface, while Eq. 2 captures its field transformation capabilities via 36 variables inside the tensors. Therefore, the metasurface hologram synthesis problem of Fig. 1 can now be defined in terms of the determination of these unknown susceptibility tensors needed to produce the desired field configuration across the surface. This field configuration is, of course, determined by the various incident and scattered EM fields needed to create the illusion.

A. Propagation and Problem Formulation
To elucidate the metasurface hologram synthesis problem, consider, for simplicity, the 2D case of Fig. 2, where all field interactions with the metasurface happen in the x − y plane and there is no field variation along z, and ∂/∂z = 0. Consider a reference simulation first, where the first task is to compute the desired scattered fields from a given object for the specified reference wave, as shown in Fig. 2(a). This is achieved by using the Integral Equation (IE) form of the Maxwell's equations and the creation of field propagators.
The EM fields radiated into free-space from electric and magnetic current sources, {J, K}, can be generally expressed as [33], [34]: where the various field operators are given by: The operators L{·} and R{·} produce a field response at location r due to a distribution of current sources located at r . This can be conveniently expressed in a matrix form as: which can be further simplified as, where, are the current source and radiated field vectors respectively and {·} denotes a matrix transpose. The P(r, r ) propagator matrix thus relates the current sources at r to the EM fields produced by them at an arbitrary location r. Finally, G(r, r ) inside Eq. 3, represents the Green's function, which for a 2D case is given by the 2 nd Hankel function, where J 0 and Y 0 are the Bessel functions of the 1st and 2nd kind and the function represents outwardly propagating radial waves. Consider again the illumination of an object producing scattered fields as in Fig. 2a. These scattered fields can alternatively be produced by a set of equivalent currents C(r o ) on the surface of the object, which are determined by Eq. 4. By specifying the appropriate boundary conditions at the known object surface, these unknown currents can be easily solved using standard IE solvers for a specifed reference field [33], [34].
To obtain the desired scattered fields to be later reconstructed by the metasurface hologram, let us introduce a Horizon Plane at r h , which is always located between the object (r o ) and the observation plane (r 0 ). Its utility will become clear shortly. The scattered fields at this horizon, F s (r h ) can be calculated by forward-propagating the fields using Eq. 4 giving, Thus, F ref. s (r h ) represents the desired reference fields that our metasurface hologram must reconstruct in the absence of the object. This completes the first task.
Next, let us remove the object and introduce a metasurface described in terms of its surface susceptibilitiesχ(r m ), which is excited with an arbitrary illumination field F i (r), as illustrated in Fig. 2(b). Dropping the perpendicular terms and assuming scalar susceptibilities, for simplicity, Eq. 1 can be combined with Eq. 2 to give, where ∆ψ T = E(r m+ ) − E(r m− ), and ψ av = {E(r m+ ) + E(r m− )}/2, are expressed in terms of total fields just before and after the metasurface. Since, the metasurface and the horizon are not in general co-located (i.e. r m = r h ), the horizon fields F ref.
s (r h ) must now be used to determine the average fields around the metasurface, F(r m ). The relationship between them depends on whether a Posterior or an Anterior illusion is desired.

B. Anterior vs Posterior Illusions
Let us take the case of front-lit object and front illumination for explaining the procedures for relating the fields at the Reference Simulation Posterior Illusion Anterior Illusion horizon and the metasurface. 2 For the posterior illusion, the metasurface is located in front of the object, as shown in Fig. 2(b). In this case, a judicious choice is to place the metasurface directly at the horizon, so that r m = r h . This simplifies Eq. 7, so that the total fields around the surface are given by, where zero field is arbitrarily enforced on the right-half of the metasurface 3 . All the fields in Eq. 7 are now known, and the unknown surface susceptibilities can now be easily determined to complete the hologram synthesis. With the known susceptibilities and the illuminating field, the metasurface will correctly reconstruct F ref.
s (r h ) everywhere beyond the horizon, so that the observer (also located beyond the horizon) will perceive an illusion of the original object located at its original position. This completes the hologram field specification for this case and the synthesis of the surface can be performed as described in Sec. III-D.
The relationship between F ref. s (r h ) and F(r m ) is however more complex for the anterior illusion case. See Fig. 2(c), where now the metasurface is located behind the object, i.e. r m = r h . In this case, the horizon fields must be reverse propagated to the metasurface location. This reverse propagation is not the usual physical forward-propagation of the fields, but a purely numerical exercise, where horizon wave-fronts are mathematically propagated back to the metasurface location while maintaining the original flow of EM power towards the observer on the left of the horizon. This technique of reverse-propagation requires a few intermediate steps, before the metasurface is ready for synthesis.
To reverse-propagate the horizon fields, the equivalence principle is invoked. The horizon is first represented as a surface with unknown equivalent currents C(r h ), so that where P(r h± , r h ) are the self-propagator operators, and the fields on the left of the horizon are fixed to zero. In principle, the unknown equivalent horizon currents can be obtained by inverting Eq. 9(a) giving . However, in practice, this is not a robust method and the related matrices are not well-behaved due to an ill-definition of the problem. To improve the robustness of the computation we formulate the problem as a set of unknown surface currents, C(r h ), and scattered fields, F s (r h± ). Then in addition to Eq. 9, the equivalent currents, C(r h ) are enforced to be tangential to the horizon with zero normal components and the tangential field components across the horizon are set equal to the reference fields at the horizon, F ref.
s (r h ). These relations can be expressed conveniently in the matrix form as:  wheren × {·} extracts the tangential components of the argument vector. This solution of this system significantly improves the computation of C(r h ) by ensuring that its normal components are zero and the interface boundary conditions are applied using the tangential fields only. The fields produced by the equivalent sources, C(r h ), must now be reverse-propagated to the metasurface location, r = r m . To achieve this operation, the field propagator of Eq. 4 is modified so that the propagation is reversed in a temporal sense, where the operator P r (r m , r h ) is formulated using an alternative Green's function involving a Henkel's function of the Reconstruct Fields at Observer F(r 0 ) BEM Solver, Eq. 19 [26,30] Metasurface Hologram first kind: This function represents inwardly propagating radial waves and with respect to fields generated by surface currents is, of course, a non-physical time reversed solution to Maxwell's equations. However, in this case it is useful as a mathematical tool. Once the desired fields are reverse-propagated to the metasurface, the total fields around the metasurface can be formed as which when substituted inside the GSTCs of Eq. 7 can now be solved for the unknown surface susceptibilities of the metasurface hologram (see Sec III-D). The complete design flow chart for a general Front/Back-Lit Posterior/Anterior Illusions with Front/Back Illumination problem is illustrated in Fig. 3, summarizing various metasurface synthesis scenarios. This process also shows the final confirmation that the synthesized susceptibilities indeed produces desired fields beyond the horizon at the observer location r 0 in which a BEM-GSTC framework is used to determine the output fields for the specified illuminating field used in the metasurface design [26], [30]. The boundary element method (BEM) is a numerical computational method of solving linear partial differential equations which have been re-formulated as descritized integral equations (IE).

C. Discretization -BEM Framework
The field equations for synthesizing metasurface holograms have so far been expressed in terms of the continuous space variable r. However, for numerical computation, each of these equations must be spatially discretized for incorporation into a BEM solver framework.
The surface position vector r S and any equivalent currents (J and K) on an arbitrary shape (such as the bounding region of a reference object) can be spatially discretized using m meshing elements, such that r → r S and J → J, i.e. r S = r s,1 r s,2 . . . r s,m and J = J 1 J 2 . . . J m , for instance. Similarly, the field propagation operators L{·} and R{·}, can also be discretized. For example, The resulting fields anywhere in space, produced by these discretized current sources, are obtained using the field propagation operators of Eq. 3. If the region where the fields are measured is also discretized, r p = r p,1 r p,2 . . . r p,n , we get which further can be written in compact form as F s (r p , r S ) = P(r p , r S )C(r S ), with C = J K . (14) If the observation fields are on the surface itself (required when solving for equivalent currents of the reference object or during reverse-propagation), then r p = r S , then we have for the two sides of the surface (+/-), Defining a surface field configuration S F = F s S+ F s S− T and a surface propagator P S = P S+ P S− we have, The metasurface GSTCs of Eq. 7 can be explicitly written in terms of total E-and H-fields on the left ({· 1 }) and right ({· 2 }) half of the surface as, where the surface susceptibility terms are expressed using auxiliary variables as, The matrix operatorN T performs the operation of extracting the two tangential fields at the surface (one in the x − y plane and the other with respect toẑ) obtaining E T from E for every surface element. The operatorN T × extracts the total tangent field and then rotates these two fields to implement thê n × {·} T operation on every element. The discretized GSTC matrices can now be expressed compactly as

D. Susceptibility Synthesis
The metasurface surface susceptibility synthesis rests on solving the GSTCs matrix equation of Eq. 16(a), once all the desired scattered fields S s F are known and the illumination fields, S i F , are specified. To extract the unknown surface susceptibilities, let us consider scalar susceptibilities, for simplicity. This extraction can be conveniently performed by rearrangingḠ T F of Eq. 16(b), as Considering that the metasurface susceptibilities are discretized over the surface as χ(r m ) and thus become vectors of localized susceptibilities (X), we can express the right hand side of Eq. 16(a), as where • is the point-wise Hagamard product. Each of the terms, is a column vector of one component of tangent fields (xy and z). If we wish to create a distributed χ vector we can form, where we define a diagonal matrix, This form allows a very convenient expression of RHS of Eq. 16(a), where the susceptiblity matrix term is explicitly extracted as Finally using Eq. 16, we now have the explicit relationship for the spatially varying surface susceptibility matrix as which can be used directly for metasurface synthesis for a given S F .

E. Solution of the Final Configuration
To confirm the synthesis of the MS a full simulation of the surface with appropriate illumination is needed. For a single surface subject to an illumination we use Eq. 15 to determine the surface fields, force the normal components of the currents on the metasurface to be zero and enforce the interface conditions prescribed by Eq. 16a. These equations can be assembled into a final matrix equation, where N DC takes the dot product of the currents for all elements and enforces N DC C = ∅ -setting the normal component of the currents to zero. The surface fields S F have been split into two components: 1) the unknown scattered fields S s F , and 2) the known applied field on the metasurface (reference or the illumination fields, for instance) S i F , so that S F = S s F + S i F . The solution of this equation provides C and S s F and using C, the fields at any point in the simulation domain can be obtained using a propagation matrix.

IV. RESULTS AND DISCUSSION
To illustrate the synthesis procedure for metasurface holograms and the ability to reconstruct desired fields, let us consider a few examples. Given that a large number of illusion configurations are possible, however, only few will be  presented for better readability and tractability, and to highlight various aspects of the proposed synthesis procedures. The simulation setup follows the illustration of Fig. 2, in the x − y plane, and the frequency of operation is chosen to be 60 GHz (λ = 0.005 m) where the observer is fixed at x 0 = −15λ. The surface meshing is set to λ/10 based on proper convergence and the metasurface has a physical extent of 120λ unless otherwise noted.

A. Front-lit Posterior Illusion with Front/Back-illumination
We first synthesize a metasurface hologram that creates an illusion of a rectangular PEC (Perfect Electric Conductor) object. The first step is to determine the scattered fields from a real PEC object when excited with a reference field -a uniform plane-wave -incidenting at an angle of θ in as shown in Fig. 4. The reference total and scattered fields are shown in Fig. 4(a) and (b), respectively. The ideal PEC object creates a shadow region behind it, and produces strong scattered fields through its flat faces and sharp corner diffraction. Next, the Horizon is defined at x = 0, which is an arbitrary choice as long as it is lying to the left of the object. The scattered fields (both amplitude and phase) are now recorded on the horizon, which our metasurface hologram must recreate.
We then remove the object and a metasurface hologram is introduced. For posterior illusion, the metasurface is placed directly at the horizon at x = 0. For this first test, the illumination fields are kept identical to the reference fields. The metasurface susceptibilities are next computed using the Horizon fields computed earlier assuming a planar configuration, and it is found that it is sufficient to use only the electric and magnetic surface susceptibilities, χ ee and χ mm , to recreate the fields. They are shown in Fig. 4(c) as a function of space. Their complex nature suggests that the metasurface must produce a spatially varying transformation of both amplitude and phase.
Then using the standard BEM-GSTC technique of [26], the response of this synthesized metasurface is simulated when excited with an illumination field. The resulting total and scattered fields are shown in Fig. 4(d-e). The scattered fields of Fig. 4(e) produced by the metasurface hologram must be compared with that of the object only in Fig. 4(b). The metasurface hologram successfully recreates the scattered fields on the left of the horizon, while producing zero fields on the other side as imposed in the synthesis steps. The 1D recreated scattered fields at the metasurface location are compared with the desired fields at the Horizon and shown to be perfectly superimposed in Fig. 4(f), as expected. Furthermore, at the observer located at x = x 0 , there is an excellent match between the recreated scattered fields and the original fields of the object. Therefore, from the perspective of the observer, the metasurface hologram is perceived as the original PEC rectangular object -a virtual object image is formed behind the metasurface.
The synthesized surface susceptibilities strongly depend on the illumination fields and the shape of the metasurface. The topography of the metasurface is not limited to planar configurations. For instance, if a curvilinear metasurface is preferred, the scattered fields from the object can also be recreated, however, a different set of surface characteristics are needed. Fig. 4(g) shows one example of the synthesized susceptibilities with its corresponding scattered fields shown in Fig. 4(h). As before, this curvilinear metasurface recreates the desired scattered fields perfectly everywhere to the left of the horizon as shown in Fig. 4(i). However, this time the surface susceptibilities of Fig. 4(g) are significantly different from the ones of a planar metasurface of Fig. 4(c). In particular, the curvilinear metasurface exhibits alternating loss-gain regions, compared to the purely passive susceptibilities of the planar metasurface. This illustrates that the choice of metasurface shape is an important design parameter to be considered in practical hologram designs. Now let us consider back-illumination of the metasurface, in this case, the metasurface hologram is excited from the right side of the surface -a uniform plane wave of normal incidence (this is, of course, a case of illumination field being different from the reference field). Fig. 5(a) shows the synthesized susceptibilities with its corresponding scattered fields shown in Fig. 5(b). The reference object's scattered fields are perfectly recreated left of the horizon, with a virtual object image formed behind the metasurface. The 1D plot comparison of Fig. 5(c) confirms the perfect field creation everywhere including that at the observer. Compared to the front illumination, where the metasurface was transforming the illumination fields in reflection, the back illumination requires a transmission-type metasurface, with zero reflection, i.e. a perfectly matched metasurface. Such a metasurface is typically realized using a Huygens' source configuration, and is generally difficult to implement compared to reflection-type surfaces.
It should be noted that in all these examples, once the metasurface hologram is synthesized for a given illumination -front or back -the illusion is produced perfectly in the entire region of space left of the horizon, irrespective of the location and extent of the observer. However, we should be aware that the complexity of the metasurface design rests on the spatial variation of the surface susceptibilities, which is controlled by the choice of the reference object, metasurface shape and the illumination conditions, thereby producing virtually an infinite number of configuration possibilities. The hologram designer will need to make judicious choices to achieve practically realizable illusion configurations. For instance, in practice a purely passive metasurface may be desired compared to a loss-gain metasurface which requires an actively powered metasurface, as some of these above results demand.

B. Front-lit Anterior Illusion with Front/Back-illumination
We shall now deal with the case of an Anterior illusion, where the metasurface hologram is placed behind the object at x m . Fig, 6(a-b) shows the total and scattered fields generated by a parameterized curvi-linear reference PEC object, which is excited by a Gaussian beam from the bottom left of the object producing a front-lit object. Complex scattered fields are produced along with a shadow behind the object. Since, the object is in between the horizon and the metasurface, the horizon fields at x h must be mathmatically reverse-propagated to the metasurface location, before the susceptibility synthesis can be performed.
This reverse-propagation of the horizon fields is shown in Fig. 6(c). In this step, the object is removed, and following the procedure of Sec. III-B, the fields are reverse propagated to the desired metasurface location at x m towards the right. The fields are naturally zero on the left of the horizon. The reverse propagation creates an artificial field distribution which appears to focus the horizon fields before diverging again at the metasurface location. The resulting fields E bck.
z,s (and other associated fields and components) are now the desired scattered fields that the metasurface hologram must recreate in the absence of the reference object under specified illumination fields.
Next, for specified illumination fields (front illuminated normally incident uniform plane-wave) and the reverse-propagated horizon fields, the metasurface surface susceptibilities are synthesized, as shown in Fig. 6(d). The resulting scattered fields are further shown in Fig. 6(e) along with the fields comparisons at the metasurface and observer location with the reference fields. The scattered fields from the metasurface at the metasurface are perfectly recreated, while an excellent reconstruction of the fields at the observer location is seen with slight ripples.
A similar field reconstruction is observed when the metasurface is re-synthesized for back illumination, where a uniform plane-wave is normally incident from the right, as shown in Fig. 7. As expected, for a matched metasurface, the electric and magnetic surface susceptibilities are balanced emulating a Huygens' source configuration. Moreover, in both the front and back illumination case, the synthesized susceptibilities are found to be purely lossy. Finally, in both cases, the virtual image of the object is formed in front of the metasurface. It should be noted that, while the horizon location is arbitrary, it is limited to the left-most extent of the object, beyond which the illusion is perfectly created.

C. Back-lit Anterior Illusion with front-illumination
The last example of this section is a case of back-lit anterior illusion with front illumination. In this case, the same parametrized object of Fig. 6 and 7, is excited with a reference Gaussian beam from the back of the object. The computed total and scattered fields from the object are shown in Fig. 8(a) and (b). Compared to all the previous front-lit cases, the observer this time measures not only the scattered fields, but also the reference fields. Consequently, the total fields are captured at the horizon. These total fields are then reverse propagated to the desired metasurface location behind the object (anterior illusion), as shown in Fig. 8(c).
Next for the specified illumination fields -a front illuminating Gaussian beam -the metasurface is synthesized, and the resulting susceptibilities are shown in Fig. 8(d). The final fields scattered from the metasurface are shown in Fig. 8(e) along with the field comparisons at the metasurface and the observer location in Fig. 8(f). In this case, the scattered fields from the metasurface must be compared to the total fields of the reference object, i.e. Fig. 8(a), in the region left of the horizon. As expected, a near-perfect reconstruction of the reference fields from the metasurface hologram is observed throughout the observation region.
These near-perfect results may be compared to the previous case of a front-lit object ( Fig. 6 and 7), where slight ripples were observed in the recreated fields. While the object is the same in both cases excited with a Gaussian beam, the two reference fields are quite different due to front vs back lit conditions. This further suggests that the nature of the reference fields to be recreated by the hologram impacts the accuracy of the reconstructed fields.

V. PRACTICAL ASPECTS OF METASURFACE HOLOGRAMS A. Effect of Finite-sized Metasurfaces
So far, the metasurface has been treated as an essentially ideally large surface, separating the two half regions. However, practical metasurfaces are finite-sized in nature, and this finite extent of the metasurface may have important consequences on the quality of scattered field recreation using our holograms.
The key distortions are the presence of secondary diffraction from the surface edges and spilling over of the illumination fields around the metasurface, and reducing the field-of-view of the illusion, in some cases.
Let us take the previous example of back-lit anterior illusion with front illumination case, where the length of the metasurface hologram is changed. Fig. 9(a) shows the case of a metasurface of length 20λ, modeled with a dielectric (non-scattering) boundary on each side. The 2D total fields clearly show the strong penetration of the back illumination towards the left of the horizon. Fig. 9(a) also shows the 2D field comparison at the metasurface and the observer compared to the ideal case of a large metasurface. It is clear that the recreated fields significantly differ from the desired fields, with some resemblance near the center of the observer only. Large field oscillations are present on either side of the metasurface due to the illumination field spilling around the metasurface. This could also be attributed to any surface waves reflected off the edge discontinuities of the surface and forming standingwave type fields on the surface. As the metasurface is made larger to 40λ, the field reconstruction improves as shown in Fig. 9(b), which eventually becomes near-identical to the desired ones for a length of 80λ, as shown in Fig. 9(c). A further improvement is shown in Fig. 8(f) in which a surface of length 120λ was used.

B. Effect of Illumination Field Strength
It is clearly evident throughout all these examples, that the synthesized metasurface susceptibilities strongly depend on the nature and configuration of the illumination fieldshow the metasurface hologram is illuminated. In general, the synthesized surfaces exhibit varied regions of loss-gain characteristics to enable the recreation of the desired fields along with a large variation in the susceptibility values. Such surfaces are practically challenging to implement due to design complexity requiring active surfaces and purely passive surfaces maybe desired. One simple way to influence the passivity of the surface, is to increase the illumination field strength. This    could be seen as a practical way to enforce a passive surface using an external control.
In all the previous examples, the illumination field strength was kept the same as the reference field. Fig. 10 shows the effect of the illumination strength on the synthesized surface susceptibilities of the metasurface hologram using the backlit anterior illusion with front illumination case of Fig. 8. Fig. 10(a) shows the nominal case of unity amplitude, where the synthesized susceptibilities show large peaking values of χ ee and χ mm in several local regions, along with several active regions on the surface where the fields are amplified. When the illumination strength is increased to 1.25, the susceptibility values significantly improve to lower values, with the surface becoming near passive and also producing a better reconstruction of the fields, as shown in Fig. 10(b). With a further increase in the illumination strength to 1.5, the susceptibility values drop to small values, although still exhibiting some local regions of gain. However, this time the desired fields are near-perfectly reconstructed, as seen in Fig. 10(c).
This example illustrates the sensitivity of the metasurface hologram design to illumination field strengths and demonstrates that it is an important parameter to take into account. It should be noted that the increased illumination strength is equivalent to using a lower reference field to compute the desired scattered fields. This disparity in the fields will simply manifest as an illusion of lower brightness under normal illumination conditions.

VI. CONCLUSIONS
A systematic and structured approach has been presented to exploit the rich field transformation capabilities of EM metasurfaces for creating a variety of EM illusions using the concept of metasurface holograms. A holistic approach of metasurface hologram synthesis has been undertaken here from a system point of view, where the desired fields detected by the observer are first recreated and then fed back into the overall metasurface synthesis problem. Considering the complexity of the problem and large number of configuration possibilities, the approach of classifying them in terms of front/back-lit posterior/anterior illusions using front/back illumination has been adopted for better organization and tractability of the overall synthesis problem. These classifications are based on specific relationships between the reference object to be recreated, the observer measuring the object, the orientation and placement of the reference and illumination field, and the desired placement of the metasurface hologram creating a virtual image. Consequently, a general design procedure to synthesize metasurface holograms has been proposed based on the IE-GSTC method and its EM illusion creation capabilities has been demonstrated using several examples. The proposed synthesis framework results in the determination of the spatially varying surface susceptibilities describing the EM properties of the metasurface. The synthesis technique combines the integral form of the Maxwell's equations, and the corresponding field propagators with the GSTCs describing the field interaction with zero thickness metasurfaces. It has next been implemented using the BEM approach using a rigorous and compact matrix formulation which is capable of synthesizing planar as well as curvilinear metasurface holograms and for arbitrary specifications of the reference object. Finally, the impact of the metasurface size and the illumination field strength on the quality of the reconstructed scattered fields has been discussed with respect to the feasibility of practical metasurface holograms.
The number of possible configurations and situations for creating EM illusions using metasurface holograms are virtually unlimited. While the handful of examples presented here have been strategically chosen to highlight and illustrate several of the salient features of the hologram synthesis, the presented framework represents a flexible test-bed to explore a wider variety of illusion scenarios before undertaking practical demonstrations. It further highlights the unprecedented capabilities of EM metasurfaces in achieving very complex wave transformations. While only scalar surface susceptibilities have been employed in this work, the BEM-GSTC framework is easily extendable to fully tensorial descriptions of the surface as was done in [26]. Furthermore, the usage of GSTCs combined with surface susceptibility description of zero-thickness surfaces is also a very efficient tool for modeling complex objects (and not limited to PEC objects as used here). This makes the proposed numerical framework complete for handling arbitrarily complex problems using a common IE-GSTC infrastructure. It should also be finally remarked that, while a vast number of techniques have been proposed for designing optical holograms, the proposed technique represents a rigorous full-vectorial field-based approach for metasurface synthesis, as opposed to methods typically based on paraxial approximations at optical frequencies [13]. Furthermore, compared to the existing works on metasurface synthesis, the system level approach undertaken here where the desired fields are first computed from physical considerations, results in a metasurface synthesis problem that is likely to be well-posed, as opposed to the possibility of an otherwise arbitrary and physically disconnected problem. Therefore, this work provides a set of important tools to metasurface hologram designers for creating a myriad of complex EM illusions throughout the EM spectrum.