Characterizing Scalar Metasurfaces Using Time-Domain Reflectometry

Two efficient methodologies for the determination of electromagnetic (EM) constitutive properties of scalar metasurfaces are introduced and discussed. In contrast to the available methods, and in line with the recent increasing interest in time-domain (TD) analyses of metasurfaces, we show that the material parameters of a scalar metasurface can be readily achieved directly in the TD merely from the EM reflected pulse shape. The two methodologies are based on an analytical TD reflectometry (TDR) approach and a modern stochastic optimization technique. A number of illustrative numerical examples demonstrating the validity and properties of the proposed techniques are presented.


I. INTRODUCTION
Time-domain (TD) reflectometry (TDR) is an efficient nondestructive testing methodology that is based on the detection and subsequent interpretation of a pulsed electromagnetic (EM) field signal reflected by a device under test [1]. It has found a wealth of applications in remote characterization of faulty electric transmission-line systems [2], ultra-wide-band antennas [3], radar targets [4], optical fibers [5] and many others. The present paper aims at proposing new applications of TDR to the determination of EM constitutive properties of a class of metasurfaces.
Generally, metasurfaces can be viewed as purposefully created thin screens that hold great promise for designing novel high-performance antennas, absorbers [6] and lenses [7] (see [8], for other relevant references). Their design can be based on analytical models (e.g. [9]), or more generally, on dedicated numerical techniques (see [10]- [13], for example) that incorporate cross-layer transition conditions [14], [15]. The latter approach has been recently pursued in Ref. [16], where the conjugate gradient minimization is applied to achieve the desired source distributions. While being robust and efficient, this approach is limited to harmonic fields, and thus to linear and time-invariant metasurfaces. The same limitation applies to Refs. [17], [18], where the frequencydomain (FD) method of moments is applied to synthesize metasurface holograms.
In order to exploit application potentialities offered by space-time metamaterials [19], [20], both direct and inverse modeling procedures have to be formulated in the TD [21]. An initial contribution to this effort is presented in this article, where we introduce two efficient methodologies capable of extracting material properties of a scalar metasurface from the pulse reflected by the metasurface. The first TDR approach is analytic and is inspired by the TDR experiment designed for the detection, localization and characterization of faults in power lines [22]. The second methodology relies on a stochastic optimization approach, the feasibility of which has been demonstated in [23]. The (forward) mathematical model employed in both inversion approaches is based on the TD saltus-type conditions applying to thin screens with combined magneto-dielectric properties [15].

II. PROBLEM DEFINITION
The problem under consideration is shown in Fig. 1. Position in the examined configuration in the 3-D space R 3 is specified by coordinates {x, y, z} with respect to an orthogonal Cartesian coordinate system with the origin O and unit vectors {i x , i y , i z } forming the standard base. We shall analyze the pulsed EM field scattered by the metasurface, which is located in a homogenous isotropic and loss-free medium described by the electric permittivity ϵ 0 and magnetic permeability µ 0 . We assume that the layer occupies is the surface of the layer in the xy plane and δ its thickness.
Indeed, the electric permittivity or/and magnetic permeability of a metasurface in the real-frequency domain may have an imaginary part. The present analysis, however, is carried out entirely in the time domain, where losses manifest themselves through additional time-convolution operators in the EM constitutive relations. Metasurfaces are frequently designed by arranging a relatively large number of small scatterers into a 2-D regular pattern of vanishing thickness. To simplify the modeling of their complex structure, the EM scattering by such surfaces is commonly analyzed using bulk material parameters. This model, in line with the broad definition given in Ref. [20], is also adopted in the present work. In particular, for the sake of simplicity, the metasurface is assumed to be described by the (homogenized) relative electric permittivity and relative magnetic permeability, which are the desired material parameters to be extracted. The (mostly negligible) effect of electric conductivity or/and linear magnetic hysteresis losses can be incorporated through the model introduced in Ref. [15].

III. SOLUTION METHODOLOGIES
In this section we shall describe two approaches to achieving the relative electric permittivity ϵ r and relative magnetic permeability µ r of the scalar metasurface (see Fig. 1). The employed scattering model relies heavily on the results introduced in [15].

A. ANALYTICAL APPROACH
The inverse material characterization methodology described in this section is based on the idea hinted at in [22] regarding a TD reflectometric scheme for characterizing faults on a transmission line. Combining the ideas presented in Refs. [15] and [22], the desired material parameters can be precisely determined with minimal computational complexity.
Without any loss of generality we assume that the examined layer is irradiated by an impulsive y-independent, TEpolarized EM plane-wave defined as: where p 0 = sin (β)/c 0 and γ 0 = cos (β)/c 0 are the slowness parameters in the x and z direction, β is the angle of incidence and c 0 is the corresponding EM wave speed. The TD wave reflected from the scalar metasurface represented by [15,Eq. (10)] can be with the aid of the Laplace transform written aŝ where s is the Laplace-transform parameter (= complex frequency). Furthermore, Ψ and Ω represent the influence of the desired parameters ϵ r and µ r via [15,Eqs. (11,12)] where Y 0 = ϵ 0 /µ 0 denotes the wave admittance and δ is the thickness of the layer. Adopting further the TDR methodology presented in Ref. [22], the plane-wave signature of the incident EM field is described by the exponential pulse: where e m is the pulse amplitude and α denotes the pulse decay coefficient. It is noted that the decay coefficient corresponds to 1/t w , where t w is the pulse time width. The sdomain counterpart of Eq. (5) then immediately follows aŝ Combining Eq. (2) with (6), transforming the result to the TD, we find at once where t ′ = t − p 0 x − γ 0 z and H(t) denotes the Heaviside unit-step function.

Calculation Procedure
The closed-form TD expression (7) specifying the pulse reflected against a scalar metasurface will be next used to determine its constitutive material properties. To that end, we shall conduct two experiments considering two distinct exponential-pulse excitations that differ in their pulse decay coefficients. For the sake of clarity, we denote these experiments and its corresponding α-coefficients as follows Consequently, Eq. (5) can be used to represent the two pulses in experiments A and B as e i (t|α A,B ) = e m exp (−α A,B t)H(t). Upon carrying out these experiments, we get two reflected pulses from which we can extract their peaks and the corresponding instants. The first step is to calculate the echograms E r y (x, z, t ′ |α A ) for experiment A and E r y (x, z, t ′ |α B ) for experiment B. For the sake of brevity, we shall further drop the (unit) amplitudes of the incident pulses and consider the normalized reflected Consequently, using Eq. (7) we obtain Next, taking the time differentiation of Eq. (10) we find instants t r;αA and t r;αB at which Pursuing this approach, we end up with the following (independent) equations In practice, the time derivative is calculated using a suitable finite-difference approximation.

Illustrative example
In the following example, the probed metasurface is excited by the uniform plane waves with the unit amplitudes e m = 1 V/m with the angle of incidence β = π/4 rad. The interrogation pulses, shown in Fig. 2, are further described by their decay coefficients α A = 1 ps −1 and α B = 0.5 ps −1 . For this instance, the layer under consideration is characterized by thickness δ = 50 µm and ϵ r = 18, µ r = 6. Figure 3 represents Eq. (10) for the selected parameters of the incident pulses and metasurface. Figure 4 shows a graphical representation of the solutions of Eq. (11) for the selected parameters of the incident pulses and the metasurface. The corresponding peaks, indicated in Fig. 3, are denoted by E p;αA = E r y (x, z, t r;αA |α A ), E p;αB = E r y (x, z, t r;αB |α B ). Substituting the data extracted from the echograms, (t r;αA,B , E p;αA,B ), in Eq. (10), we end up with a system of two (non-linear) Eqs. (12a), (12b) with two uknowns Ψ and Ω. Its solution finally yields via Eqs. (3) and (4) the desired material parameters ϵ r and µ r . (12b)

B. GLOBAL OPTIMIZATION APPROACH
A method for characterization of thin-sheet metasurfaces was introduced in [23]. The method is based on the cooperation of an arbitrary stochastic optimizer with a relatively new VOLUME 4, 2016 model to the TD EM fields in the vicinity of a metasurface with combined magneto-dielectric properties. This TD solution is used as the forward solver that evaluates the candidate solutions u proposed by the global optimization algorithm that solves the single-objective optimization problem that can be formulated as follows : Here, N is the total number of discrete time samples n, r is the position vector where the electric field is observed. Symbols E comp and E true denote the "computed (optimized)" and "true (measured)" observed electric field, respectively. Candidate solutions u are vectors u = {ϵ r , µ r } that can be located at the decision space Γ. Generally, the optimization problem defined by (13) can be solved by any single-objective optimization algorithm, in general.However, the comparative study [23] proved that Particle Swarm Optimization (PSO) algorithm [24] can solve that problem in the most effective way from the set of algorithms containing four other state-of-the-art algorithms namely: Genetic Algorithm [25], Differential Evolution [26], Invasive Weed Optimization [27], and Covariance Matrix Adaptation-Evolutionary Strategy [28]. Therefore, only the PSO algorithm is used in this study.
PSO is the representative of the so-called swarm intelligence algorithms. The set of particles (decision space vectors u) co-operates to search for the position with the best value of the objective function. All the particles move in the decision space with a different velocity. The velocity of every particle is updated based on a combination of three procedures: 1) the particle is forced to continue in the random inertial movement, 2) the particle is attracted to its personal best position (position visited by the particle with the so-far best value of the objective function), and 3) the particle is attracted to the global best position (the best position from all personal best positions in the swarm). Therefore, the particle uses cognitive learning (it benefits from its own experience) and social learning (it benefits from the knowledge of the whole swarm) as well. The balance between the exploration and exploitation of the algorithm is made by the balance of weights for the individual procedures 1-3). For more details about the principles of PSO, the reader is refered to [24]. The MATLAB implementation of the PSO algorithm in the software package FOPS [29] is used in this study.
The trade-off between all three "forces" mentioned above can be set by user-defined controlling parameters. The inertia weight w forces particles to explore the decision space Γ, the cognitive learning factor c 1 supports the local exploitation of the area near the personal best, and the social learning factor c 2 favors the exploitation of the area near the global bets. All the controlling parameters of the PSO algorithm used in this study are summarized in Tab. 1.
All the results presented below are based on statistical data collected for 100 repetitions of every optimization run.

IV. NUMERICAL EXAMPLES A. ANALYTICAL APPROACH
In this subsection, we present the results of the analytical approach for a variety of selected parameters β (= the angle of incidence) or δ (= the layer's thickness). Furthermore, we will perform the calculation for two distinct pairs of interrogation pulses illustrated in Fig. 5. Thus, both experiments A and B will be carried out twice (see Eqs. (8) and (9)).  For the analytical model to apply, the metasurface is assumed to be very thin with respect to the spatial support of the excitation pulse, that is, c 0 /α >> δ. Therefore, the pulse time width of the excitation pulses is to be chosen with respect to this condition.
To find out the influence of β and δ on the desired parameters, we are considering two mutually independent cases. Tables 2  and 3 contain the values of α-coefficients, β and δ. In both cases, the EM constitutive properties of the layer ϵ r = 50 and µ r = 3 remain constant.
Based on the parameters shown in Tabs. 2 and 3, the graphical outputs (Figs. 7 to 10) of the obtained echograms E r y (x, z, t ′ |α A,B ) and its derivative ∂ t E r y (x, z, t ′ |α A,B ), for both experiments, are below. For the sake of clarity, the zerocrossing regions are zoomed and incorporated in the figures as their insets. From these regions, we can read the key values t r;αA,B and E p;αA,B that are next used to extract the desired material parameters.
It is interesting to observe that in contrast to the incident pulse shapes (see Fig. 5), the reflected ones take both positive and negative values. This fully complies with the enforced cross-layer conditions (see Ref. [15]) and the pertaining EM field equations.
Tables 4 to 7 summarise the instants t r;αA,B , normalized to the value of pulse time width t w;αA,B , and corresponding peaks E p;αA,B for both cases.
It is seen from Eq. (2) that if Ψ = Ω, the reflected field response is identically zero and the metasurface is thus electromagnetically transparent. Equivalently, this condition can be written as where β 0 is the angle of incidence at which the metasurface behaves as transparent. For the considered TE polarized EM wave, the Eq. (14) can be satisfied when µ r ≥ ϵ r . Figure 6 illustrates the angles β 0 satisfying Eq. (14) for three different values µ r = {75; 100; 150} and constant ϵ r = 50. Apparently, in this circumstance, the analytical TDR approach fails, which calls for an alternative approach. It will be next demonstrated that a stochastic optimization approach helps to bypass the limitation.

B. GLOBAL OPTIMIZATION APPROACH
The where D is the dimension of the decision space vector u (D = 2 in our case). The superscript true denotes the searched (optimal) values: u true 1 = ϵ r = 50.0, and u true 2 = µ r = 3.0 in our case.
First, we show the influence of the incidence angle β on the accuracy of the proposed characterization method. The DER values for different values of β from range 0 ≤ β ≤ 3π/8 are plotted in Fig. 11. The accuracy slightly decreases with increasing β. The figure compares the accuracy of the method for different number of objective function evaluations VOLUME 4, 2016 (OFE) determined by number of agents and iterations: OFE = N A × N I . While for N A = N I = 30 we get the result with an unacceptable error of order 10 1 , all other combinations provide results with DER < 1.0. It is interesting, that curves for combinations N A = 50, N I = 100 and N A = 100, N I = 50 almost overlap. This implies that it does not matter whether we invest the available OFE rather to more agents or iterations.
Next, we investigate how the DER metric depends on the metasurface thickness δ. This knowledge can be of great importance as the parameters of the characterization method (namely the width t w of the EM wave pulse) can be adjusted according to the thickness of the sample under test. Figure 13 shows the DER metric against values δ from interval 0.1 ≤ δ ≤ 20.0 (in µm). Overall, the accuracy of the characterization gets better with growing δ. This effect can be explained by the variable shape of the objective function that depends on the mutual relation between the layer's EM constitutive parameters and its thickness. There is an obvious discontinuity in DER for all combinations of the number of agents and iterations near value δ ≈ 8.0 µm. To further investigate the mutual influence of the t w and δ parameters, we perform the parametric study over them.
The resulting average values of the DER metric for combinations of values t w and δ are shown in Fig. 12. Please note that the colormap is scaled logarithmically in the figure. The figure clearly shows a slice where the method achieves the error below the order of 10 −4 .
The major advantage in using the global optimization algorithm approach is that it works reliably for all combinations of the searched relative permittivity and permeability values. This can be evidenced by the results shown in Fig. 14. Here, a contour plot with a logarithmic scale of DER metric values for different combinations of metasurface parameters ϵ r and µ r from the whole space Γ is shown. It proves that a reasonably low value of the error (under 10 −2 ) can be achieved for any pair ϵ r -µ r and the error falls below 10 −5 for a significant part of Γ where µ r > 2.5. The orange dashed line in Fig. 14 denotes the region of transparency, where the metasurface is not visible for the used EM wave [15]. The methodology based on the global optimization, however, yields the correct material parameters even if the condition of transparency is met.

V. CONCLUSION
We have presented analytical TDR and stochasticoptimization approaches to characterizing the constitutive properties of a scalar metasurface directly in the TD. It has been shown that the analytical TDR approach allows to obtain the material properties exactly in a computationally effortless manner. On the other hand, the analytical methodology fails in the region of TD EM transparency. In this circumstance, the stochastic optimization approach lends itself to remedy the issue. This has been demonstrated on a number of illustrative numerical examples demonstrating the high accuracy and robustness of the optimization approach.