The challenge of modeling grounded metasurfaces: a new approach to obtain their constitutive parameters

A new technique to retrieve the constitutive parameters of grounded metamaterials is proposed in this article. Previous methods presented in the literature used to retrieve the constitutive parameters of materials in general and metamaterials, in particular. are not applicable since they rely on both reflectivity and transmittivity data. In grounded metamaterials, there is no transmission and hence, a new method is needed to obtain the metamaterial constitutive parameters (both permittivity and permeability) by means of just reflectivity data. The proposed methodology is based on extracting the aforementioned constitutive parameters using the reflectivity response to two different angles of the incident plane wave. First, the solution to the problem is presented from a mathematical perspective. Then, the methodology is applied to a grounded metasurface and their constitutive parameters are extracted. From these parameters, interesting phenomena that occurs in the metasurface are explained. Finally, the range of validity of the method is also analysed.


Introduction
The interaction between waves and materials is a key aspect to understand the propagation of electromagnetic waves. The equations that rule this interaction are known as Maxwell equations, which are sometimes difficult to solve and tricky to understand. Therefore, equivalent circuits or homogenization models are sometimes used to alleviate the formulation and give physical insight to better understand the phenomena beyond it [1].
It should be noticed that metamaterials are complex structures in a microscopic view, but they can be usually analyzed in a macroscopic view by employing different homogenization techniques. Consequently, they can be described by their effective constitutive parameters. In metamaterials, these parameters do not only depend on their composing materials, but also on the geometrical patterns of their microstructure (microscopic crystallography) [2].
Metamaterials are sometimes designed based on the constitutive parameters requirements for the intended application. In this case, two approximations can be followed. One based on retrieving the constitutive parameters from a metasurface (two-dimensional periodic metamaterial) and the other is based on designing the metasurface structure via the required constitutive parameters. In the former one, different homogenization methods can be used and in the latter one, the transformation optics technique can be adopted [3]. This paper will be focused on the first approximation, assuming that the constitutive parameters of the metasurface are wanted to be retrieved for better characterization of the structure.
According to the possible dispersion of metamaterials, different models have been presented in the literature [2][3][4]. It should be noticed that in metamaterials with a non-negligible electrical size of either their unit-cell's metallizations and/or gap between unit-cells, their spatial dispersion should be considered. Consequently, attending to the size and shape of the metamaterial's unit-cells different models can be used.
When these unit-cells are smaller than the wavelength quasi-static models, which only consider the frequency dispersion, can be employed [5,6]. There are two types of models: dielectric and magnetic ones, which are, respectively, based on analyzing the effect of applying a timeharmonic electric or magnetic field. The application of one or the other field, cause a movement of electrons or a change of their orbital and spin motions in the material, respectively.
On the other hand, when the unit-cells are electrically small, but they exhibit complex shapes, weak spatially dispersive models, which also consider the interaction between elements in the same unit-cell, should be employed. In this case, the local effects dominate and the non-local ones (mutual coupling between unit-cells) can be described by the bi-anisotropy and artificial magnetism (effective permeability). In this case, multipolarization of the unit-cell's metallizations has to be considered and the spatial dispersion of first and second order is transferred into additional constitutive parameters (such as the magnetoelectric coupling parameter). On the other hand, when the interaction between unit-cells is strong and hence, the non-local effect is dominant, strong spatially dispersive models should be employed. Therefore, higher order multipoles have to be taken into account and the constitutive parameters will depend on the wave vector of the incident wave [4].
In most cases, metamaterials' unit-cells are much smaller than the effective wavelength and hence, they can be considered as homogeneous media and described at their macroscopic level. Indeed, in these cases, metamaterials are non-or weak spatially dispersive, so that their constitutive parameters can be computed from their scattering ones (S parameters) [7,8]. Therefore, the inverse scattering method introduced by Smith et al. [9], which is based on the wellknown Nicholson and Ross equations [10], and reformulated by Ziolkowski [11] can be used to retrieve the constitutive parameters of certain metasurfaces. However, this method exhibits a phase ambiguity in the retrieval procedure, which was partially solved by considering the metamaterial as passive [9] or using other theories that ensure causality [12]. Group delay measurements were also used to solve the phase ambiguity [13]. Other contributions compare the constitutive parameters retrieved from simulations with the ones obtained from the model to properly determine these parameters [14,15]. In [16], the Kramers-Kroing relations are used to develop an algorithm for choosing the proper branch of the phase.
Aiming at circumventing the previous ambiguity, Holloway et al. based its formulation on the computation of the metamaterials' electric and magnetic susceptibilities [17]. These susceptibilities can be computed by analyzing the reflection and transmission properties of the metamaterial [18].
For bianisotropic metamaterials, Simovski et al. developed a model for retrieving the constitutive parameters by relating the incident electromagnetic fields to the electric and magnetic moments of the inclusions, through a set of collective polarizability tensors [19,20]. These tensors are also extracted from the reflection and transmission properties of the metamaterial [21].
It can be noticed that all aforementioned methods require to compute both the reflection and transmission properties of the metamaterial. However, when considering a metamaterial grounded by a metallic plate, on which no transmission occurs, the previous models fail. Consequently, these methods cannot be applied to most metasurfaces, since they usually entail a frequency selective surface (FSS) on a grounded dielectric.
To properly characterize metasurfaces and better understand their physical behavior, it is crucial to retrieve their constitutive parameters. Consequently, the aim of this paper is to extract an additional equation which can compensate for the lack of transmission in grounded metasurfaces and hence, allowing to compute both the relative permittivity and permeability constitutive parameters. Moreover, the characterization of a metamaterial by its constitutive parameters enables to relax the computation burden for optimizing the performance of these structures. Indeed, the latter opens a new way to swiftly design different types of metasurfaces such as reflectarray antennas.

Retrieving procedure
The proposed method, as most of the cited ones, is based on the assumption that the metasurface can be homogenized (unit-cell's size much smaller than the effective wavelength). In this case, the behavior of the metasurface is identical to the one of a similar medium characterized by its effective permittivity ( ef f ) and permeability ( ef f ). The proposed method will be also focused on retrieving the constitutive parameters from the scattering ones. As it was mentioned, the difficulty here stems from the lack of transmission, since a grounded metamaterial is considered.
Assuming the metasurface (MTS) in the air and laying in an XY plane (see Fig. 1), and after applying proper boundary conditions on the air-metasurface interface [22], the following S 11 parameter, which is the relationship between the reflected and incident electric field, can be extracted under both transverse electric (TE) and transverse magnetic (TM) polarizations: where E r TE,TM is the reflected electric field under TE or TM polarization. E i TE,TM is the incident electric field under TE or TM polarization. Γ TE,TM is the reflection coefficient in the The reflection coefficient in the interface air-MTS can be expressed as: Having a look at the Eqs. 1-3, one can notice that there are two unknown parameters ( ef f and ef f ) and just one is known ( S 11 ). Therefore, an additional equation to extract the constitutive parameters is needed. In this paper, it is proposed to consider two different incidence angles ( ) to obtain the additional equation required to solve the problem. Therefore, a normal incidence and an oblique one will be regarded.
From Eq. 1 and firstly considering normal incidence ( k z ef f = k ef f ; k z air = k air cos( = 0 • ) = k air , where k air is the effective wave number in air), the reflection coefficient can be written as: Therefore, the wave impedance in the MTS ef f = √ can be expressed as: , it is possible to use the Taylor series to expand the following expressions: 1 , c o n s i d e r i n g k air = k ix air + k iz air = k air sin + k air cos and applying several mathematical operations, one can arrive at the following equation when considering a TE polarized incident wave: Where S 11 o is the S 11 parameter under oblique incidence (θ).
The variables are all known except for k z ef f and k ef f . Therefore, it is possible to express one of the latter variables as a function of the other: From the previous expression, one can notice that k z ef f depends on the angle of the incident wave, which is a reasonable conclusion, whereas k ef f should be angularly independent, since a homogeneous medium has been considered ( k ef f = w √ ef f ef f , where w is the angular frequency).
On the other hand, knowing that the tangential wave Then, ef f can be obtained using Eq. 5 and the effective constitutive parameters can be extracted using the following expressions [11]: Once the calculation procedure is exposed, it can be observed that two solutions are obtained for k ef f , as long as the first branch of the phase is chosen (which is a proper selection when considering electrically thin metamaterials [22]). Therefore, aiming at avoiding ambiguities, two Fig. 1 Electromagnetic fields impinging the metasurface for both TE a and TM b polarizations different procedures can be followed. The former one is to repeat the procedure presented above, using another pair of incidence angles and then comparing the k ef f obtained from both procedures at each frequency and select the ones that are closer. The alternative, which will be the one selected for this paper as it is a slightly faster procedure, is to compute the S 11 parameter, using the extracted constitutive parameters as indicated in Eq. 1 and then, to compare them with the ones obtained from simulation (or measurement). The result that provides better matching will determine the proper k ef f .
For TM polarization a similar procedure can be followed but considering Γ TM . After several computations, one can arrive to the following expression: By solving the previous equation, it is possible to express k z ef f as a function of k ef f . However, in this case a nonlinear relation arises, which has to be solved using numerical techniques. As several solutions can be obtained, the passivity condition ( Im k ef f < 0 ) will be used to select the proper one, giving rise to the physical values of the constitutive parameters.

Method applied to a grounded metasurface absorber
Aiming at validating the introduced method, it will be applied to extract the constitutive parameters of the metasurface absorber (MTA) presented in Fig. 2a, which comprises an FSS on a grounded dielectric, so that no transmittance takes place (see Fig. 2b). This structure has been previously published by the authors of this paper in [23] and it exhibits a resonance frequency around 6 GHz with a periodicity (p) of 16.4 mm and a thickness of 0.49 mm. It should be noticed that the periodicity of the unit-cell is much smaller than the effective wavelength.
As it was previously mentioned, to retrieve the constitutive parameters of the MTA using the proposed method, two incidence angles have to be considered. Therefore, the normal incidence and an oblique incidence angle close to the normal incidence (2°) will be first considered.

Retrieving the constitutive parameter assuming a TE polarized incident wave
Applying the method introduced in the previous section and considering a TE polarized incident plane wave impinging the metasurface, the constitutive parameters (relative permittivity and permeability) of the MTA can be retrieved, which are presented in Fig. 3.
From the results, one can observe that the permittivity exhibits an anti-resonant behavior around the resonance frequency (see Fig. 3a), but the permeability has a resonant behavior (see Fig. 3b). Moreover, the imaginary part of the permittivity exhibits a clear negative peak value, which is directly related to the high absorption on the metasurface. Additionally, in [23] it was suspected that the structure might have a predominant inductive behavior, which can be corroborated here from the retrieved constitutive parameters (see Fig. 3b).
On the other hand, the peak observed around 6.5 GHz on the permittivity figure can be attributed to a discontinuity in the retrieved procedure, since it coincides with a null value of the permeability, which is a parameter used to compute the permittivity as shown in Eq. 11. However, it lies far from the resonance frequency of the structure which is the band of interest here.
Another interesting aspect can be extracted from the retrieved wave number and impedance (see Fig. 4). The imaginary part of the wave number Im k ef f is mainly lower or equal to zero, whereas the real part of wave impedance Re ef f is greater or equal to zero. Consequently, it can be concluded that the retrieved constitutive parameters are physical ones.
Finally, the S 11 R parameter, which is computed through Eq. (1) applying the retrieved constitutive parameters of Fig. 3, is computed (see S 11 R in Fig. 5). Then, this parameter is compared with the S 11 scattering parameter computed through a full wave simulation using HFSS software ( S 11 in Fig. 5) and the results are presented in Fig. 5, concluding that almost identical results are obtained.
Aiming at estimating the differences between the S 11 parameter retrieved from the constitutive parameters and the one from simulation, the mean error power is computed as follows: From the latter computation, an almost negligible error is obtained (close to 0%).

Retrieving the constitutive parameter assuming a TM polarized incident wave
Similarly, as proceeded above, the constitutive parameters of the metasurface are computed, but in this case, it is assumed that a TM polarized incident wave impinges on it.
As it was previously stated, a numerical technique, instead of an analytical one, has to be considered to retrieve these parameters (see Eq. 12). Consequently, the wave number ( k ef f ) is numerically extracted from Eq. (12) and using the passivity condition as mentioned in the previous section. After this computation, the constitutive parameters are extracted using Eqs. (10,11) and the results are presented in Fig. 6. Similar conclusions to the ones extracted in the previous section can be drawn: a, respectively, anti-resonant and resonant behavior for the permittivity and permeability parameters, a predominantly inductive behavior and a high negative peak of the permeability imaginary part, which is attributed to the absorption properties of the metasurface.
Following a similar procedure to the one described in Sect. 3.1 both the retrieved constitutive parameters are used to compute the S 11 R parameter and a full wave simulation is used to compute the S 11 parameter, concluding that almost no errors can be noticed (Fig. 7).
The proposed method, as most of the cited ones, is based on the assumption that the metasurface can be homogenized (unit cell's size much smaller than the effective wavelength). In this case.

Validity of the proposed method
In the previous section, the proposed method has been validated by using a pair of incidence angles that are close to each other. Is there any limitation in the selection of these angles?
Aiming at testing the latter, in this section, instead of considering an oblique incidence angle close to the normal incidence one, an angle of 20° will be regarded. This angle has been randomly chosen just to validate the proposed method.
After applying the proposed method, the effective constitutive parameters are presented in Figs. 8, 9, when considering a TE and TM polarized incident wave, respectively. It can be noticed that large discontinuities can be observed in this case. Indeed, when reconstructing the S 11 parameter (see Figs. 10, 11), slightly higher errors can be noticed as compared with the ones obtained in the previous Section (1.12 and 0.05%, respectively, for the TE and TM polarization cases). Therefore, though similar conclusions regarding the metasurface properties can extract (anti-resonant and resonant behavior of the permittivity and permeability    11 Real (a) and imaginary (b) parts of the S 11 , when considering an oblique incidence angle of 20º for TM polarization parameters), small errors are obtained when considering incidence angles that are close together.

Conclusions
In this paper, the difficulties of retrieving the constitutive parameters when dealing with metasurfaces backed by a ground plane have been introduced. A novel method for retrieving these parameters, which is based on the metasurface behavior under two different incidence angles, has been proposed. The method has been validated by using a previously published metasurface absorber and by showing the expected properties that have been mentioned in that paper, regarding its predominant inductive behavior (anti-resonant and resonant behavior of its permittivity and permeability, respectively). Moreover, the well-known loss phenomena that occurs in metasurface absorbers can be also observed in the computed constitutive parameters. Finally, the S 11 parameter is reconstructed using the retrieved constitutive parameters and compared with the results from simulations, concluding that small errors are obtained when comparing both. Finally, the validity of the method is presented and one can conclude that the closer the considered incidence angles are, the better the method will work. It should be highlighted that some approximations have been adopted (as in other well-known papers) to develop the formulation so that this can be considered when applying the method. Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Declarations
Conflict of interest All authors declare that they have no conflicts of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.