Method for Electromagnetic Property Extraction of Sublayers in Metal-Backed Inhomogeneous Metamaterials

A scattering (S-) parameter method has been proposed for electromagnetic property extraction of a target layer within metal-backed inhomogeneous metamaterial (MM) structures. It relies on recursive S-parameters for two different polarizations (parallel and perpendicular). Two different algorithms depending on the value of the incidence angles were proposed to add flexibility to our method. The algorithms were first validated by a metal-backed quasi-one-port method applicable only for one-layer (homogeneous) samples. Then, they were applied to extract electromagnetic properties of the split-ring-resonator-wire MM slab of a metal-backed two-layer inhomogeneous structure and the Omega MM slab of a metal-backed four-layer inhomogeneous structure, and extracted properties by our algorithms were compared with those retrieved by different methods applicable for homogeneous samples only. The accuracy of our method was also examined when there was some noise in S-parameters, when there was a misalignment in incidence angles, and a when large value of iteration number was used.


I. INTRODUCTION
Eelectromagnetic material characterization is an active and important subject of research field and can be directly associated with numerous applications including medical treatments, bioelectromagnetics, agriculture, food engineering, and civil engineering [1], [2]. Electromagnetic wave propagation, reflection and transmission properties of materials can be tailored (e.g., electromagnetic cloak [3] and superfocusing [4]) and accurately be predicted if materials are fully characterized. This characterization can be implemented at high frequencies by measuring complex scattering (S-) parameters and then extracting electromagnetic properties (refractive index n, normalized wave impedance z, relative complex permittivity ε r , and relative complex permeability µ r ) of materials under test using these S-parameters.
There are many microwave methods, which could be grouped into resonant type or non-resonant type, to achieve The associate editor coordinating the review of this manuscript and approving it for publication was Abhishek K Jha . such characterization [1]. Selection of a suitable method for a given problem highly depends on some criteria or restrictions such as measurement accuracy, ease of calibration, suitability of the domain (time-domain or frequency-domain), destructive versus non-destructive, invasive versus noninvasive, broadband or narrowband characterization, and accessibility of sample surface (one-side reflection measurement or two-side reflection and transmission measurements) [1], [2], [5]. For example, if the measurement accuracy is the primary concern at narrowband frequencies, then the conventional waveguide cavity resonant method can be applied for the accurate characterization of low-loss test specimen over a narrow band [6]. On contrary, if the material will be characterized over a broadband frequency using a relatively accurate measurement method, then non-resonant methods are the first choice [7]- [49], [52]- [57]. Besides, frequency-domain [7]- [51] or time-domain [52]- [57] measurements could be carried to characterize the material under test. Timedomain measurements are based on the measurement of wave travel and attenuation of a short-duration pulse incident on VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the material. However, these measurements still require expensive measuring equipment such as short-duration pulse and oscilloscope. On the other hand, frequency-domain measurements can be implemented using reflection and/or transmission properties of the material for the electromagnetic signal over a band. Frequency-domain measurements necessitate relatively inexpensive measurement instrument such as vector network analyzer (VNA) operating at microwave frequencies.
Many of the above mentioned methods are either based on direct application of boundary conditions or state-transitionmatrix. They use reflection and transmission measurements of the material under test in frequency-domain or timedomain. While reflection properties give information mainly of surface characteristics, transmission properties provide full information on the whole volume of the material since electromagnetic signals either in time-domain or frequencydomain must propagate through the material [34]. However, the methods using transmission measurements (in addition to reflection measurements) [7]- [34], [50]- [52] become inappropriate for applications where only one-side of the material is accessible (the sample is terminated in a highly reflection surface or metal) such as reflective metasurface applications (a 2D form of 3D metamaterials) [49], ground-penetratingradar applications, or target discrimination in stealth applications [58]. These circumstances necessitate the development of reflection-only methods [35]- [49], [53]- [56] for analysis and characterization of electromagnetic properties of materials, because wave cannot propagate through (reflected back from) the material with highly reflective termination. However, most of these reflection-only methods [35]- [49], [53] are not applicable for electromagnetic characterization of inhomogeneous structures. The methods [54]- [56] can be used for electromagnetic characterization of inhomogeneous multilayered structures. Nonetheless, the methods [55], [56] are developed for temporal analysis of reflection properties of multilayered structures using the concept of subregions and thus are not directly applicable for reconstruction of their electromagnetic properties. On the other hand, the method [54] can be used to reconstruct the electromagnetic properties of multilayered dielectric stacks. Nevertheless, this method is applicable in time-domain measurements which may require expensive measuring instruments. In the research paper we propose a frequency-domain extraction algorithm for electromagnetic characterization of a single layer in multilayered structures with short-circuit termination using recursive S-parameters with priori knowledge of electromagnetic parameters of other layers.
The organization of the remainder of the manuscript is as follows. First, expressions of recursive S-parameters for parallel and perpendicular polarizations were derived in Section II. Next, in Section III we give the details of our extraction algorithms for retrieval of electromagnetic properties of inhomogeneous structures. Then, numerical and simulation analyses were performed to validate our proposed extraction algorithms and evaluate their performance against noise is S-parameters, iteration error, and angle misalignment in Section IV. Finally, all results and remarks related to our method were recapitulated in Section V. Fig. 1 illustrates the schematic view of an inhomogeneous structure which is assumed to be composed of m homogeneous layers and terminated by a short-circuit termination. Each layer has a refractive index n s and a wave impedance Z s where s = 0, 1, · · · , m. It is assumed that a uniform plane wave is incident obliquely from layer 0 (air). In the following two subsections, we will derive underlying expressions for two different polarizations. In our analysis, we will consider the time dependence in the form exp(+jωt).

A. PERPENDICULAR POLARIZATION
For this polarization, electric and magnetic fields in the layer s can be expressed bȳ where H i s⊥ = E i s⊥ /η s , H r s⊥ = E r s⊥ /η s , and Here, the superscripts 'i' and 'r' denote quantities for incident and reflected waves in each respective layer; k 0 = ω/c is the free-space wavenumber; ω = 2πf is the angular frequency; f is the linear frequency; c is the velocity of light; θ i s and θ r s are the angles of incident and reflected waves of the sth layer; and η s is the intrinsic impedance of the sth layer; Enforcing the boundary conditions for the y component ofĒ s and x component ofH s over each interface at z = d t , we derive the following recursive relation between S-parameters at each interface for the incidence angle θ i 0 151706 VOLUME 8, 2020 (see Fig. 1) for the perpendicular polarization where Here, S 11⊥ is the recursive reflection S-parameter defined at z = d 1 + d 2 + · · · + d s−1 for the incidence angle θ i 0 for the perpendicular polarization [31], [33] and denotes the overall reflection coefficient of a structure composed by the layer s + 1 to the short-circuit termination; (s,θ i 0 ) ⊥ is the (interfacial) reflection coefficient right at the interface formed by the layer s − 1 and the layer s for the incidence angle θ i 0 for the perpendicular polarization; and k zs is the wavenumber in z direction within the layer s.

B. PARALLEL POLARIZATION
For this polarization, electric and magnetic fields in the layer s can be written as where E i s = H i s η s and E r s = H r s η s . Applying the boundary conditions for the y component ofH s and x component ofĒ s over each interface at z = d t , we obtain the following recursive relation between S-parameters at each interface for the incidence angle θ i 0 for the parallel polarization where S (m,θ i 0 ) 11 where T (s,θ i 0 ) is defined in (6). We note from (4)-(12) that when the incidence angle θ i 0 is zero, then Z  where When there is only one-layer (m = 1), then (13) further reduces to which is a well-known expression in the literature [49], which validates our forward problem.

III. EXTRACTION PROCESS
For extracting the electromagnetic properties of each layer in Fig. 1, measurement/simulation of complex S (0,θ i 0 ) 11 and S (0,θ i 0 ) 11⊥ for a given inhomogeneous multilayer structure with m layers itself is not sufficient to fully extract (fully characterize) the electromagnetic properties of this multilayer structure. The following extraction procedure with the flowchart in Fig. 2 can be applied for electromagnetic characterization of any target layer u of the inhomogeneous structure with a priori knowledge of parameters of other layers. Since the procedures for the extraction of electromagnetic properties for both polarizations are the same, without loss of generality only the perpendicular polarization will be considered in the explanation of the extraction procedure. 1) Electromagnetic properties of all layers except for the target layer u are initially supplemented by one-layer material characterization techniques such as [15] when θ i 0 = 0 or [9], [10] when θ i 0 = 0. It is noted that this step is implemented only once to initialize our extraction procedure.
2) Set the simulated/measured S-parameter for the first incidence angle θ i 01 to S (0,θ i 01 ) 11⊥ and iterate (4) for u ≥ 2. Repeat this VOLUME 8, 2020 process for the second incidence angle θ i 02 . It is noted that either θ i 01 or θ i 02 can be equal to zero (see the step 4) [A similar strategy can be followed for finding S (u−1,θ i 01 ) 11 and S (u−1,θ i 02 ) 11 from (10)].
3) If u = m, re-express (4) in the following form from the steps 2 and 3, we determine η 0 u in terms of n u for normal incidence as where (u) The correct sign for η 0 u can be ascertained by enforcing the passivity principle; that is, e{η 0 u } ≥ 0 where e{ } means the real part of ' '. Next, substitute η 0 u from (17) into the following objective function where Here, | | denotes the magnitude of the complex quantity ' '. As a figure of merit, n u value resulting in minimum F obj (n u , θ i 02 ) is considered as the optimized solution. After compute n u using a suitable 2-D numerical method from the objective function F obj (n u , θ i 02 ) in (20). More information about the applied numerical method will be given in Section IV. Once n u is computed and η 0 u is evaluated from (17), the relative complex permittivity ε ru and relative complex permebility µ ru of the uth layer can be finally evaluated from A similar algorithm can be exercised for the parallel polarization using S (u−1,0) 11 and S (u+1,0) 11 after following the steps 2 and 3. 5) When θ i 01 = 0 and θ i 02 = 0, minimize the following objective function can be utilized for simultaneous computation of η u and n u where and Here, T (21) and (22). As a figure of merit, η u and n u values resulting in minimum G obj (η u , n u , θ i 01 , θ i 02 ) are considered as the optimized solution. After computing η u and n u from the objective function G obj (η u , n u , θ i 01 , θ i 02 ), ε ru and µ ru of the uth layer can be determined using (25). 6) Run the steps 2 through 5 to determine the electromagnetic properties η u , n u , ε ru , and µ ru of the uth layer till the iteration number n iter reaches the maximum iteration number (N max ). It is instructive to discuss the advantages and disadvantages of the proposed method. The first advantage of the proposed method is that it does not require transmission S-parameter for electromagnetic property extraction as compared with the methods [7]- [34], [50]- [52]. Second advantage of the proposed method is that it is applicable for retrieval of electromagnetic properties of a target layer within an inhomogeneous structure contrary to the transmission-reflection methods [7]- [10], [12]- [15], [17]- [19] and the reflectiononly methods [35]- [49] which are applicable for one-layer (homogeneous) structures. It should be specifically mentioned here that this flexibility of our method distinguishes it from the method in [49], which extracts electromagnetic properties of a metal-backed homogeneous (one-layer) medium using oblique and normal incidence S-parameters, in that it could be applied for electromagnetic property extraction of any target layer within a metal-backed inhomogeneous structure using oblique and normal incidence S-parameters. In other words, it can be stated that our proposed method could be considered as a generalization of the method [49] for electromagnetic property extraction of any target layer within a metal-backed inhomogeneous structure. The third advantage of the proposed method is that it has two different algorithms to increase its flexibility. While the first algorithm uses the objective function F obj (n u , θ i 02 ) in (20) to determine ε ru and µ ru of the target layer u, the second algorithm tries to minimize the objective function G obj (η u , n u , θ i 01 , θ i 02 ) in (26) for the same goal. Because the first algorithm searches for complex n u for given S-parameters and the uniquely determined η 0 u from (17), it is more efficient in terms of computational complexity than the second algorithm which optimizes simultaneously complex n u and η u for given S-parameters. On the other hand, both algorithms require prior knowledge of the electromagnetic properties of layers except for the target layer u within an inhomogeneous structure for electromagnetic characterization of the target layer u. Besides, as its another disadvantage, solutions for η u and n u optimized from the objective functions F obj (n u , θ i 02 ) in (20) and G obj (η u , n u , θ i 01 , θ i 02 ) in (26) by our algorithms might not be the true solutions, if the initial guesses for electromagnetic properties of all layers except for the target layer u are not appropriately given. This is due to the presence of complex- (20) and (26) which produce multiple (non-unique) solutions for η u and n u in the inversion process [21], [25]. This problem is especially a big issue for a target layer with a substantial thickness and/or a larger ε r and/or µ r .

A. VALIDATION AND ACCURACY ANALYSIS
We performed some numerical analysis to validate our method for a metal-backed homogeneous (one-layer) scenario before applying it for extracting electromagnetic properties of any target layer within a metal-backed inhomogeneous structure. Therefore, the method in [49] applicable for a metal-backed homogeneous (one-layer) medium was selected as a benchmark method to compare the accuracy of our proposed method. Toward this end, we calculated S (0,θ i 0 ) 11⊥ and S (0,θ i 0 ) 11 of a metal-backed homogeneous conventional material with ε 1 = ε 0 (3.32−i0.021), µ 1 = µ 0 , and d 1 = 1 mm for f = 10 GHz. Here, θ i 01 = 0, and θ i 02 angle is varied from 5 to 20 degrees with one degree increment. In application of our proposed method based on minimization of the objective function F obj (n u , θ i 02 ) in (20), we applied the 'fmincon' function of MATLAB with the 'active-set' algorithm for 1 ≤ ε r ≤ 7 and 0 ≤ ε r ≤ 0.2 where ε r = ε r −iε r . For example, Figs. 3(a) and 3(b) illustrate the extracted real and imaginary parts of the relative permittivity (ε r1 = ε 1 /ε 0 ) of this material over θ i 02 for parallel and perpendicular polarizations. We note from Figs. 3(a) and 3(b) that while our proposed method extracts accurate ε r1 for all incidence θ i 02 angles, the method [49] retrieves ε r1 values diverging from its prescribed value 3.32−i0.021 especially for larger θ i 02 angle values. For example, the retrieved real part of ε r1 by the method [49] is 18% less than its assumed 3.32 value for θ i 02 = 20 degrees. In addition, the performance of the method [49] decreases further for the perpendicular polarization since impedance variation for this polarization is greater than that for the parallel polarization for the considered analysis, which is not shown for conciseness.

B. SIMULATION RESULTS
We considered the structures shown in Figs. 4(a) and 4(b) to apply our method for extraction of electromagnetic properties of inhomogeneous structures. The commercial 3D electromagnetic simulation program -CST Microwave Studio-was used to obtain S-parameters. Simulation details are given as follows. We used the frequency-domain solver with adaptive mesh refinement feature to simulate oblique incidence S-parameters for perpendicular and parallel polarizations at the waveguide port located at z = 0. While unit cell boundary conditions were implemented in xz and yz planes, open boundary condition was set right at the port plane (xy plane). PEC material (E t = 0) was located right at the end of each structure at z = d 1 + d 2 + · · · + d m . For comparison and validation of our proposed method, we set another waveguide port right at z = d 1 + d 2 + · · · + d m to get simulated transmission S-parameters for the application of the method [15]. The frequency band in our analysis was arranged between 6 and 20 GHz.
The first inhomogeneous structure (m = 2) consists of a metamaterial (MM) unit cell consisting of split-ring resonators (SRRs) with wire and FR4 material (d 2 = 2.0 mm, ε r,FR4 = 4.4(1 − j0.02), and µ r,FR4 = 1). The geometry of the MM cell is identical to that considered in our previous study [33] except that the slit axis of SRRs is rotated by 90 degrees. The SRR MM slab will exhibit electric coupling when the electric field is normal to the slit axis (y axis in Figs. 4(a) and 4(b)) and magnetic coupling when the magnetic field is normal to the SRR plane (yz plane in Fig. 4(a) and 4(b)) [9], [59], [61]. The MM cell in cubical form with side d 1 = 2.5 mm has SRRs over the front face of the FR4 substrate with thickness 0.25 mm. The SRR has the outer ring of 2.2 mm, the separation distance of 0.15 mm, the linewidth of 0.2 mm, and the gap of 0.15 mm. Besides, the wire with 0.2 mm linewidth extends in y direction at the back of the FR4 substrate. Metallic parts of SRRs and wire are assumed to be made by the copper material with thickness of 35 µm.
We then applied our first algorithm (θ i 01 = 0 and θ i 02 = 0) using the objective function F obj (n u , θ i 02 ) in (20) and our second algorithm (θ i 01 = 0 and θ i 02 = 0) using the objective function G obj (η u , n u , θ i 01 , θ i 02 ) in (26) to extract ε r1 , µ r1 , z 1 , and n 1 of the SRR MM slab of the two-layer inhomogeneous structure for parallel and perpendicular polarizations. In implementation of our method, we applied the 'fmincon' function of MATLAB with the 'active-set' algorithm using initial electromagnetic properties obtained from the methods [9], [10], [15]. Figs. 5(a)-5(d), Figs. 6(a)-6(d), and Figs. 7(a)-7(d) illustrate, respectively, the real and imaginary parts of these parameters for the perpendicular and parallel polarizations using N max = 3 when θ i 01 = 0 • and θ i 02 = 10 • and when θ i 01 = 5 • and θ i 02 = 10 • . In order to validate our extraction procedures, we also extracted ε r1 , µ r1 , z 1 , and n 1 of the SRR MM slab immersed in air for perpendicular and parallel polarizations using the method [15] for θ i 01 = 5 • and for the normal incidence case [9], [10]. Besides, we also computed the reconstructed S 11 of the two-layer inhomogeneous structure by using the extracted electromagnetic parameters by our method for both polarizations to evaluate whether it approaches its simulated S 11 . Toward this end, we applied the transfer matrix method (TMM) using the expressions (5) and (6) in [60]. Figs. 5(e), 5(f), 6(e), 6(f), 7(e), and 7(f) illustrate the reconstructed magnitudes and phases of S   We note the following points from Figs. 5(a)-7(f). First, extracted electromagnetic properties by our both algorithms for both polarizations follow the physical constraint -passivity principle; e.g., e{z 1 } ≥ 0 and m{n 1 } ≤ 0 for the e i ωt time-reference. Second, extracted electromagnetic properties by our method and the methods [9], [10], [15] are in good agreement over the whole frequency band for both polarizations. Third, reconstructed magnitudes and phases of S in Figs. 7(e) and 7(f) of the whole structure obtained from retrieved electromagnetic properties by our method are in good agreement with those of the simulated ones. Finally, it is seen that while the SRR-Wire MM slab in Fig. 5(b) demonstrate magnetic resonance behavior for perpendicular polarization due to the presence of magnetic field in x direction (normal to the SRR-Wire plane), the same slab is in resonance-free state for the parallel polarization (no magnetic field normal to the SRR-Wire plane and no electric field normal to the slit axis) [61].
The second inhomogeneous structure (m = 4) has the same SRR-Wire mm slab (d 1 = 2.5 mm), the same FR4 material (d 2 = 2.0 mm), an Omega ( ) mm slab, and a Teflon   1 − j0.0002), and µ r,Tef = 1). The geometry of the Omega mm slab resembles to that of the Omega mm slab (in cubical formd 3 = 2.5 mm) in [25], [29] except for 90 • rotation (without any wire at the back of the substrate). The Omega mm slab will exhibit electric coupling when the electric field is parallel with its legs (z axis in Fig. 4(b)) and magnetic coupling when the magnetic field is normal to the Omega plane (yz plane in Fig. 4(b)) [62]. The -shape has the radius r = 1.1 mm, the width w = 0.1 mm, and the tail l = 0.8 mm, and is located over the front face of the FR4 substrate (1.6 mm thick). A 35 µm thick copper is used to obtain the pattern. We applied our proposed techniques to extract the electromagnetic properties of the mm slab of this inhomogeneous structure. In implementation of our method, we applied the 'fmincon' function of MATLAB with the 'active-set' algorithm using initial electromagnetic properties obtained from the methods [9], [10], [15]. For example, Figs. 8(a)-8(d) and 9(a)-9(d) illustrate the extracted ε r3 , µ r3 , z 3 , and n 3 of this slab for perpendicular and parallel polarizations for N max = 3 when θ i 01 = 0 • and θ i 02 = 10 • using our first algorithm (F obj (n u , θ i 02 ) in (20)). Since our second algorithm (G obj (η u , n u , θ i 01 , θ i 02 ) in (26)) when θ i 01 = 5 • and θ i 02 = 10 • produces similar results for both polarizations, they are not demonstrated here for conciseness. For cross checking the performance of our proposed method, we also extracted ε r3 , µ r3 , z 3 , and n 3 of the Omega mm slab immersed in air for perpendicular and parallel polarizations using the method [15] for θ i 01 = 5 • and for the normal incidence case [9], [10]. Additionally, we obtained the reconstructed magnitudes and phases of S-parameters of this inhomogeneous structure for perpendicular polarization (θ i 01 = 0 • ) and for parallel polarization (θ i 02 = 5 • ) whether the retrieved properties are accurate, and then compared these parameters with the simulated ones, as shown in Figs. 8(e), 8(f), 9(e), and 9(f).
We note the following points from Figs. 8(a)-9(f). First, extracted electromagnetic properties by our method for both polarizations generally follow the physical constraint -passivity principle; e.g., e{z 3 } ≥ 0 and m{n 3 } ≤ 0 for the e iωt time-reference. However, extracted m{ε r3 } ≥ 0 in Fig. 8(a) over 14-16 GHz is nonphysical and attributed to spatial dispersion effects closely associated with the discreteness of conducting elements [11], [18], [25]. Second, extracted electromagnetic properties by our method, the method [15], and the methods [9], [10] are in good harmony over the whole frequency band for both polarizations. Third, reconstructed magnitudes and phases of S in Figs. 9(e) and 9(f) of the whole structure obtained from retrieved electromagnetic properties by our method are in good agreement with those of the simulated ones. Nonetheless, reconstructed |S Fig. 8(e) of the four-layer inhomogeneous structure is slightly divergent from its simulated one as compared with the reconstructed Fig. 5(e) of the two-layer inhomogeneous structure and its simulated one. We suspect that the main reason for such a discrepancy is the coupling effect between resonating SRR-Wire and Omega MM slabs which effect the electromagnetic response of individual MM slabs as well as of the whole structure [31], [33]. Finally, the SRR-Wire MM slab and the Omega MM slab behave almost indistinguishable for the parallel polarization in terms of their extracted electromagnetic properties in Figs. 6 and 9 due to non-existence of any electric or magnetic resonance phenomenon for this polarization (electric (magnetic) field is normal to (parallel with) the SRR or Omega plane) [61].
It is instructive to give some important information about the computation time to evaluate the performance of the proposed algorithms in terms of computational cost. Optimizations of the objective function F obj (n u , θ i 02 ) in (20) (G obj (η u , n u , θ i 01 , θ i 02 ) in (26)) lasted less than 1.5 (2.0) seconds for the two-layer inhomogeneous structure and 1.6 (2.2) seconds for the four-layer inhomogeneous structure to obtain the dependencies in Figs. 5-9 over the full frequency range between 6 and 20 GHz with 1001 data. The hardware we used in optimizations had the following properties: Intel(R) Core(TM) i7-5500U CPU with 2.40 GHz clock speed and 12 GB RAM. The computation cost of the algorithms will definitely vary with the hardware used in optimizations.

C. ANALYSES OF NOISE EFFECT, ITERATIVE ERROR, AND MISALIGNMENT
In application of our method, we have used the simulated S-parameters and fixed the value of N max (N max = 3) to extract electromagnetic properties of the SRR-Wire MM slab and Omega MM slab for perpendicular and parallel polarizations in Figs. 5-9. However, in real measurements, the effects of noise present in S-parameters along with any angle misalignment for parallel and perpendicular plarizations VOLUME 8, 2020 should be incorporated into the extraction to monitor their impact on the performance of the proposed method. In this regard, we first examined the effect of noise in S-parameters in the extraction of electromagnetic properties. For instance, Figs. 10(a)-10(d) demonstrate the retrieved ε r1 , µ r1 , z 1 , and n 1 of the SRR-Wire MM slab of the two-layer inhomogeneous structure in Fig. 4(a) for perpendicular polarization (resonance behavior) for θ i 01 = 0 • and N max = 3. In implementation of the noise into our analysis, we included a normally distributed random noise (zero-mean value and standard deviation σ = 0.05) into both real and imaginary parts of the simulated S we also extracted ε r1 , µ r1 , z 1 , and n 1 of the SRR-Wire MM slab by the method [15] (using simulated noise-free S-parameters of only the slab in air when θ i 01 = 10 • ). It is noted from Figs. 10(a)-10(d) that although the presence of noise augments the difference between ε r1 , µ r1 , z 1 , and n 1 properties of the SRR-Wire MM slab extracted by our method and by the method [15], retrieved properties by our method could still follow the same properties extracted by the method [15]. For instance, as observed from Figs. 5(d) and 9(d), e{n 1 } values at f = 10.2 GHz extracted by the method [15] and our method are −4.95 and −4.45, respectively -showing roughly a 10% variation, as compared to 7% variation for the dependencies obtained by the same methods in Fig. 5(d) for noise-free circumstance.
We also carried out some numerical analysis to assess the effect of iterative errors on the electromagnetic properties extracted by our method. In this analysis, we mainly focused on the effect of iteration number on these properties since round-off error and inappropriate initial seeds for electromagnetic properties (extracted by other techniques [9], [10], [15] before running our algorithm) are highly dependent on external factors. In addition to the analysis of the effect of iteration number, we also performed numerical analysis to investigate the effect of angle misalignment on retrieved electromagnetic properties by our method. For these analyses we examined a three-layer (m = 3) inhomogeneous structure composed of a FR4 material (d 1 = 3.0 mm, ), a MM synthesized slab (d 2 = 5.0 mm), and a Teflon material (d 3 = 4.0 mm). For the synthesis of the MM slab, we applied the Lorentz dispersion model [12] with ε r2 and µ r2 as where f e and f m are the electric and magnetic resonance frequencies; f is the operating frequency; δ e and δ m are the electric and magnetic damping frequencies; and F e and F m are the coefficients depending on structure of the material [12]. In our analysis, F e = 0.4, F m = 0.4, f e = 5.0 GHz, f m = 8.0 GHz, δ e = 0.6 GHz, and δ m = 0.5 GHz were used for the frequency range of 1 − 11 GHz, producing a passive material. Using the electromagnetic properties of the FR4 material, the synthesized MM slab, and the Teflon material, the reflection S-parameter of the whole configuration (metal-backed three-layer inhomogeneous structure) was computed for parallel and perpendicular polarizations for θ i 01 = 0 and θ i 02 = 10 • . Then our method was implemented to retrieve ε r2 and µ r2 of the synthesized MM slab. In implementation of our method, we applied the 'fmincon' function of MATLAB with the 'active-set' algorithm using initial electromagnetic properties of the FR4 material, the synthesized MM slab in (35) and (36), and the Teflon material. Figs. 11(a) and 11(b) illustrate the real parts of the extracted ε r2 and µ r2 (shown by circle symbol) for the This circumstance can arise due to the iterative nature of our algorithm just as other iterative methods [31], [33] used for electromagnetic characterization of inhomogeneous structures by S-parameters. As noted in the study [33], as a rule of thumb, implementation of our method can start with a smaller N max value which should be continuously increased up to the point where the difference between extracted properties for the similar N max values is within the limit of pre-assigned accuracy.
As the last numerical analysis, we analyzed the impact of angle misalignment on the performance of our method. Toward this end, we used the same aforementioned threelayer inhomogeneous structure with synthesized MM slab. In this analysis, we deliberately modified the incidence angles θ i 01 and θ i 02 from 0 • and 10 • to 5 • and 15 • , respectively, to interrogate the change in electromagnetic properties of the synthesized MM slab for perpendicular polarization (N max = 3). For instance, Figs. 12(a) and 12(b) show the retrieved real parts of ε r2 and µ r2 of the synthesized MM slab when the incidence angles are modified. It is noticed from Figs. 12(a) and 12(b) that a 5-degrees change in θ i 01 and θ i 02 values does not significantly alter the values of extracted ε r2 and µ r2 by our method. For example, extracted e{ε r2 } by our method with angle misalignment problem at 4.5 GHz is around 2.362 whereas the same quantity with no-angular misalignment problem is 2.289 (around 3% variation). It means that our proposed method is more sensitive to iterative errors than the angle misalignment problem.

D. GUIDELINES FOR FREE-SPACE MEASUREMENTS AND FABRICATION TOLERANCES
There are some key points that should be considered to perform frequency-domain free-space measurements to extract electromagnetic properties of metal-backed inhomogeneous metamaterials by our method [1]. First, because the whole theory in Section II is based upon the assumption of plane wave propagation, MM slabs should be positioned in the farfield region of the transmitting/receiving antenna. The wellestablished criterion r ≥2 D 2 /λ 0 should be exercised for such positioning. Here, r is the distance between antenna and front face of MM slabs, D is the maximum dimension of the antenna (diagonal length of rectangular antennas or diameter of circular antennas), and λ 0 = c/f is the freespace wavelength corresponding to mininmum frequency f . Second, diffraction effects which might appear around sharp corners or edges of inhomogeneous MM slabs should be minimized, since the theoretical analysis in Section II does not consider these effects. There are some options in reducing these effects. For example, MM slabs having larger transverse dimensions (much greater than the antenna aperture) could be fabricated. As another example, antennas with spot focusing facilities could be employed to minimize the foot print of the antennas over the MM slabs. As further example, calibration technique which directly take into account the diffraction effects into measurements could be applied to eliminate undesired diffraction effect. As a final example, microwave absorbers (in the form of pyramid) could be implemented for elimination of this effect, too. Third, unwanted multiple reflections from ground (or reference level) and/or neighbouring objects should as well be minimized, if not totally eliminated. As a simple, yet highly effective, remedy for minimization of this effect, antennas and MM slabs should be elevated to considerable hight from ground (or reference level) and neighbouring objects should be taken away from antennas and MM slabs. As another remedy which is as efficient as the former one, but costly (the time-domain feature of the VNA should be already installed), is to eliminate undesired reflection from ground (or reference level) and neighbouring objects by gating (filtering) some desired signals within a predefined time-domain. Finally, the inaccurately measured incidence angle can also alter the accuracy of the extracted parameters by the proposed method. The accuracy of incidence angle measurement is closely related to the design of the experimental setup. The measurement setup of the study can be referred to an interested reader for how a specially designed bistatic free-space measurements can be implemented. The appropriate minimum incidence angle increment that the free-space measurement setup should have for accurate free-space measurements can be specified only if the range of electromagnetic properties of the inhomogeneous MM structure is partly known or only after sufficient information about electromagnetic response of the MM structure is gained by some preliminary measurements, because electromagnetic waves will interact differently with structures having various electromagnetic behaviors. Then, with this information at hand, by specifying a criterion for percentage change of the electromagnetic properties of any layer (e.g., front layer) of the whole inhomogeneous structure corresponding to incidence angle for a given frequency value (or range), the minimum incidence angle increment should be determined. Our experiences show that the minimum incident angle increment should not be less than 10 • for accurate freespace measurements. Besides, the higher the incidence angle difference between two free-space measurements, the better electromagnetic property extraction is feasible, because S-parameters measured at two largely separated incidence angles will be considerably different to be easily detected.
On the other hand, in addition to important points which might affect, if not properly considered, the accuracy of free-space measurements, fabrication tolerances arising from improper construction of inhomogeneous MM slabs could as well seriously limit the performance of the proposed method if not taken into consideration. In this regard, an uncertainty analysis should be followed to evaluate the most influential parameter before fabrication of inhomogeneous MM structures. Among many others, one of the most critical parameters in the fabrication of MM cells is the resonator length (width for square resonators or radius for circular resonators). Besides, one of the most important parameters in the design of metal-backed inhomogeneous MM structures (or other artificially synthesized inhomogeneous materials or structures) is the possible, although minor, gap between the metal ending and the last layer of these structures and gap between adjacent layers of these structures. Sensitivity and uncertainty analyses could simultaneously or separately be performed to monitor such effects on the extraction procedure of our proposed method. To this end, the well-known differential uncertainty model or a model based on simulated S-parameters can be applied.

V. CONCLUSION
A retrieval method has been proposed to extract electromagnetic properties of a target layer within inhomogeneous structures backed by a metal plate using frequency-domain S-parameters. Toward this end, recursive S-parameters were first derived for both parallel and perpendicular polarizations and next two different algorithms based on the value of the incidence angle were presented to extract these properties. The algorithms were then validated against a recent extraction method for a metal-backed homogeneous (one-layer) sample. After, electromagnetic properties of two different inhomogeneous structures (a two-layer and a four-layer structures) were retrieved by our method and compared with those extracted by other methods, which are only applicable for homogeneous (one-layer) samples suspending in air. From the comparison, the following important points are noted. First, our proposed algorithms extracted accurate electromagnetic properties for each target layer within inhomogeneous structures. Second, the accuracy of our method partly decreases as a consequence of spatial dispersion and coupling effects between nearby resonating MM slabs. In addition, we examined the effects of noise in S-parameter, iteration error, and misalignment of the incidence angle on the extracted electromagnetic properties by our method. From this examination, we noted that our extraction method is relatively resistant to noise in S-parameter which are generally present in many modern S-parameter measuring VNA instruments. We also noted that maximum iteration number N max should be properly selected (starting from a lower value to a higher value beyond which accuracy does not much improve) in order to accurately extract electromagnetic properties by our method. Finally, not as much important as the iteration number, angle misalignment, which might decrease the accuracy of extraction process, should also be decreased so as to retrieve correct electromagnetic properties by our method.