Time-Domain Meshless Method Simulation of Metasurfaces With Generalized Sheet Transition Conditions

This letter presents a new time-domain formulation of a meshless method based on the generalized sheet transition conditions (GSTCs). This formulation is applicable for modeling a kind of dispersive metasurface and a time-varying metasurface. This is the first time a meshless method is developed for analysis of the metasurfaces using GSTCs. The significance of time-domain meshless methods originates from their capabilities for multiscale and conformal modeling in addition to analysis of problems in a broad frequency-range. Here, for simplicity, the meshless formulation is presented for analysis of 1-D problems consisting of monoisotropic metasurfaces. The extension of this method to 2-D and 3-D problems or bianisotropic metasurfaces is straightforward. The scheme is verified by numerical examples.

which are composed of subwavelength scattering particles. The diverse capabilities of metasurfaces in controlling the magnitude, the phase, and the polarization of the electromagnetic (EM) waves, in addition to being less bulky, less lossy, and their easier fabrication in comparison to the metamaterials have provided a vast range of applications for them in manipulating EM waves [1]. As a metasurface has subwavelength thickness and it is a complex structure, the computational cost of the full-wave analysis is too high [2]. To improve the simulation efficiency of the numerical analysis of a metasurface, this structure can be considered as a zero-thickness sheet, which creates discontinuities in the electric and magnetic fields [3]. The generalized sheet transition conditions (GSTCs) [4] technique is known as an efficient approach, which can precisely model discontinuities of EM fields due to a general metasurface [3].
Some of the first studies, which have applied GSTCs to several conventional numerical methods for analysis of metasurfaces, can be listed as follows: the finite difference frequency-domain Manuscript  in [5], the spectral-domain integral equation in [6], the finitedifference time-domain (FDTD) in [3], [7], [8], and [9], the finite-element method in [10], the method of moment in [11], and the discontinuous Galerkin time-domain method in [12]. In these conventional techniques, the spatial domains of problems are often discretized by predefined meshes/grids. Thus, in analysis of a complicated structure, meshing scheme is time consuming and remeshing process is difficult. In contrast to the conventional techniques, meshless methods provide a scheme for the spatial discretization without the use of a predefined mesh, which is based on a set of nodes scattered within the domain. In meshless methods, the flexibility in distribution of nodes can significantly reduce the computational effort of remeshing process [13]. Their advantages in conformal and multiscale modeling of complex structures over grid-based methods are demonstrated in [14]. By development of new meshless formulations for simulation of the metasurfaces, we can benefit from the unique features of the meshless methods in modeling these structures.
The application of meshless methods for analysis of some classical 2-D artificial periodic structures can be found in some research works, such as [15], [16]. As these structures are modeled by uniform unit cells, the analysis is simplified to one unit cell by applying periodic boundary conditions or periodic shape functions. In this work, for analyzing a general nonperiodic metasurface, we model it by surface susceptibilities to reduce the computational complexity. Then, by applying the GSTCs to the time-domain radial basis function (RBF) meshless method, we analyze a special kind of dispersive and a time-varying dispersive metasurface. In recent years, the RBF meshless method [17], [18], [19], and the vector RBF meshless method [20] have been developed for analysis of dispersive mediums in [21] and [22] and left-handed materials in [23], [24]. To take into account the discontinuities of EM fields in the GSTC-meshless formulations, we have considered virtual electric and magnetic nodes on both sides of a metasurface. By careful application of the shape functions of the RBF meshless method, we can accurately model the effect of discontinuities.

A. Metasurface Synthesis Equations
The frequency-domain GSTCs for a metasurface lying in the yz-plane and perpendicular to the x-axis of a Cartesian coordinate system are given by [1] x where ∇ || =ẑ∂/∂z +ŷ∂/∂y, P is the electric polarization density, and M is the magnetic polarization density. ΔH and ΔE are the differences between EM fields on the right and left sides of the metasurface, . The subscripts "inc," "ref," and "tr" denote the incident, reflected, and transmitted fields, respectively.
Using the susceptibility tensors (i.e., χ κν where κ, ν = e, m, e denotes the electric, and m denotes the magnetic), the electric and magnetic polarization densities can be expressed as where ΔH av and ΔE av are the arithmetic average of the fields on the two sides of the metasurface, By assuming that in (1) and (2) the normal electric and magnetic polarization densities are zero, i.e., M x = P x = 0, substitution of (3) and (4) into (1) and (2) results in bianisotropic GSTC synthesis equations in the frequency-domain as where χ eν = χ yy eν χ yz eν χ zy eν χ zz eν , χ mν = χ zz mν χ zy mν χ yz mν χ yy mν , and ν = e or m.
In this work, for simplicity, we have assumed that the only nonzero susceptibility terms are χ zz ee and χ yy mm . Therefore, the monoisotropic GSTCs in frequency-domain are obtained as

B. Formulation of GSTC-Meshless Method for Time-domain Analysis of a Metasurface
In order to simulation of a metasurface in this work, we use the RBF meshless method, which is acknowledged as a robust meshless method for analysis of EM problems. We will not give a complete introduction to the RBF meshless method (as it can be found in a number of research works, such as [13] and [14]). However, since the derivation of our proposed formulations is not possible without careful understanding of the concept of shape functions, we have to explain the roles of the RBFs and the shape functions in the RBF meshless method.
In the scalar RBF meshless method, in a domain containing N nodes, an unknown electric or magnetic field is interpolated using the scalar RBFs as where a j are the unknown expansion coefficients, 2 is the scalar RBF with Gaussian type and α is the shape parameter that controls the decaying rate of the exponential function, x j = (x j , y j , z j ) is the location of the jth node surrounding the point of interest at x = (x, y, z), and The unknown expansion coefficients should be calculated through solving a matrix equation, which can be obtained by enforcing (9) to pass through all the N nodes of the domain. Finally, by substitution of a j into (9), we obtain shape functions that can be used to approximate the unknown field variable as where ϕ j (x) is the value of the shape function of the jth node at the point of interest at x (the jth node is located at the point x j ) and u j (x j ) is the known field variable associated with the jth node. Also, the derivative of the field variable can be obtained using the derivative of the shape function as where ∂κϕ j (x) is the value of the derivative of the shape function and ∂κu(x) is the derivative of the field variable with respect to each spatial coordinate (κdenotes x, y, or z).
Here, to reduce the computational effort and without loss of generality, we obtain the formulation of GSTC-meshless for a 1-D problem. The generalization of the method to 2-D and 3-D cases is straightforward. To implement the meshless method, a set of electric field nodes (E-nodes) and a set of magnetic field nodes (H-nodes) should be considered for representing the spatial domain of the problem. Because of the coupling between electric fields and magnetic fields, each E-node of the set of the E-nodes should be surrounded by the H-nodes and vice versa [14]. Fig. 1 shows the arrangement of the E-nodes and H-nodes in the computational domain of the 1-D problem with considering propagation in +x direction. To take into account the effect of metasurface discontinuity, we have located a virtual E-node just before the metasurface at x = 0 − and a virtual H-node just after the metasurface x = 0 + . We assume that the background medium is free-space and perform the numerical simulations for transverse magnetic polarization, i.e., H z = E x = E y = 0. We can express the updating equations for all the nonvirtual nodes in the computational domain of the problem, except that for the H y,i+1 and E z,i+2 nodes that are located in the vicinity of the metasurface, by the following conventional equations: where x H and x E denote the locations of the H-node and E-node, respectively. Also, NE and NH are the numbers of the electric and magnetic field nodes. In Section III, we will explain that how we consider the locations (or the numbers) of the nonvirtual and also virtual nodes. For nodes around the metasurface (H y,i+1 and E z,i+2 nodes), we incorporate the virtual nodes into the interpolating term of conventional equations as Then, by applying the GSTCs to (14) and (15), we will replace the values of electric and magnetic fields at the virtual nodes located at x = 0 − and x = 0 + , i.e., E z (x E 0 − ) and H y (x H 0 + ), and obtain the governing equations for analysis of metasurfaces.
Here, we consider a dispersive metasurface with the frequency-dependent susceptibilities in the form of: χ zz ee (ω) = χ yy mm (ω) = κ/jω. When such a metasurface is illuminated by a normally incident wave, it absorbs a part of the wave and transmits the other part of the wave [3]. The degree of the absorption or transmission depends on the value of κ. Substitution of these susceptibilities into (7) and (8) Hence, the values of EM fields at the virtual nodes can be obtained by substitution of the average of electric field that is )/2 (where q = n and q = n + 1) in (16), and substitution of the average of magnetic field that is H q y,av = (H q y (x H i+1 ) + H q y (x H i+2 ))/2 (where q = n − (1/2) and q = n + (1/2)) in (17). Finally, we replace the obtained relations (14) and (15) to derive the update equations for nodes in vicinity of the metasurface as ). It is worth mentioning that although (18) and (19) are obtained for the nodes around the metasurface, removing the metasurface from the medium converts these equations into the conventional relations of (12) and (13). In this case, we have κ = 0 in (18) and (19), and as there is no need to consider virtual nodes, the derivatives of the shape functions of virtual nodes (i.e., ∂xϕ E 0 − and ∂xϕ H 0 + ) will also be omitted from (18) and (19). According to the presented approach to attain (18) and (19) for the nodes that surround the metasurface, we can provide more explanations about the virtual nodes and the role of them in modeling a metasurface. The reason we call the nodes located at x = 0 − and x = 0 + as the virtual nodes is that there is no need to directly obtain the values of EM fields at these points using the conventional updating relations of (12) and (13). Instead, by the use of GSTCs, we apply (16) and (17) for obtaining the field values at these virtual nodes. Thus, considering the virtual nodes helps us to incorporate the GSTCs into the conventional formulations for taking into account the presence of a metasurface in the medium. In addition to the values of EM fields at the virtual nodes, the derivatives of the shape functions of the virtual nodes at the points which are located in vicinity of a metasurface (∂xϕ E 0 − (x H i+1 ) and ∂xϕ H 0 + (x E i+2 )) are necessary for precise simulation of a metasurface by (18) and (19).
For the case of a time-varying metasurface, if we suppose that its susceptibilities are not frequency-dependent, we can calculate the inverse Fourier transform of (7) and (8) by replacing jω with d/dt and obtain time-domain relations as: ΔH y = ε 0 (d(χ zz ee E z,av )/dt) and ΔE z = μ 0 (d(χ yy mm H y,av )/dt). Similar to the previous case, we can obtain the values of EM fields at the virtual nodes, and substitution of them into (14) and (15) leads to the governing equations of the nondispersive time-varying metasurface. However, in general, as a metasurface is constructed based on subwavelength resonators, its susceptibilities in addition to being time-variant are also frequency-dependent [3]. Therefore, in Section III, we will investigate a type of dispersive time-varying metasurface. One can easily obtain the update equations for modeling such a metasurafce according to our presented formulations.

III. NUMERICAL EXAMPLES
Here, to validate the efficiency of our proposed equations, we have considered a 1-D problem that its computational domain is similar to Fig. 1, and we truncated the problem domain using perfectly matched layer absorbing boundary condition. We model two examples of metasurfaces, which are illuminated by a normally incident plane wave of sin(ω 0 t) using GSTC-meshless method.
To discretize the spatial domain, we considered uniform distributions for the nonvirtual EM field nodes. By examination of various uniform nodal spacing, we found that numerical accuracy and efficiency of the computational cost can be guaranteed when the distance between each E-node (H-node) in the set of E-nodes (H-nodes) is equal to d min = λ 0 /30. Also, the distance between each E-node and H-node is equal to d min /2. To specify the locations of virtual nodes, it should be noticed that shape functions of RBF meshless methods will be obtained through solving a matrix equation. Two almost identical entries in the matrix of this equation leads to ill-conditioning and inaccuracy in the results [13]. To avoid this problem, we do not position the virtual E-node (H-node) very close to a nonvirtual E-node (H-node) and consider the minimum distance of d min /5 between a virtual and a nonvirtual node. In such a nodal distribution, the choice of the shape parameter as α = 0.2 ensures the stability and accuracy of the results. Also, we determine the distance between two virtual nodes equal to d min /10. In the conventional RBF meshless method, the maximum size of time-step can be specified by Δt ≤ d min √ εμ [19]. We found that the largest time-step size to obtain stable results in GSTC-meshless formulations is smaller than this value and it is about 0.6d min √ εμ. The reason for this difference is that we have to define virtual nodes for modeling the metasurface, and the distance between a virtual node and the nearest nonvirtual node to it is smaller than d min . According to [1], the volume susceptibility of a metasurface with subwavelength thickness of d can be approximated using the relation: χ 3D ee,mm = χ 2D ee,mm /d. Assuming d = 2d min in the simulations, we can obtain the bulk (or diluted) susceptibility which is dimensionless.

A. Dispersive Metasurface
To investigate the validity of (18) and (19), we consider a dispersive metasurface with frequency-dependent surface susceptibilities: χ zz ee (ω) = χ yy mm (ω) = −2c 0 /(3jω). We know from [3] that the metasurfaces with equal monoisotropic susceptibilities do not reflect the wave, also according to metasurface synthesize equations of (7) and (8), it can be found that this metasurface absorbs half of the wave and transmits the other half. Fig. 2 shows the spatial profile and time variations of the incident and transmitted waves. According to Fig. 2, the results of simulation  using GSTC-meshless scheme agree well with the predicted results by synthesize equations.
This metasurface is synthesized in such a way that it has linear time variations, which lead to generation of new frequency harmonics. Similar to the presented example in [3], the value of β(t) changes periodically between 1 and 3. When β(t) = 1, the metasurface will absorb all the incident wave, and it will transmit the half of the incident wave when β(t) = 3. We have analyzed this problem by the proposed GSTC-meshless method formulations and compared the results with the GSTC-FDTD method's solutions. Fig. 3 shows that the proposed formulation have correctly modeled the behavior of the metasurface in time and frequency domains.

IV. CONCLUSION
In this letter, we have illustrated that the incorporation of GSTCs into the RBF meshless method leads to an efficient formulation for time-domain analysis of a dispersive and a timevarying metasurface. Using proper auxiliary differential equations, we can also extend the RBF meshless method for analysis of a number of dispersive metasurfaces that their susceptibilities are described by Drude or Lorentz dispersion models.