**Complex Bianisotropy E**ff**ect on the Propagation Constant of a Shielded Multilayered Coplanar Waveguide Using Improved Full Generalized Exponential Matrix Technique**

**Djamel Sayad 1, Chemseddine Zebiri 2,3, Issa Elfergani 4,\*, Jonathan Rodriguez 5, Hasan Abobaker 6, Atta Ullah 3, Raed Abd-Alhameed 3, Ifiok Otung <sup>5</sup> and Fatiha Benabdelaziz <sup>7</sup>**


Received: 23 December 2019; Accepted: 28 January 2020; Published: 2 February 2020

**Abstract:** A theoretical study of the electromagnetic propagation in a complex medium suspended multilayer coplanar waveguide (CPW) is presented. The study is based on the generalized exponential matrix technique (GEMT) combined with Galerkin's spectral method of moments applied to a CPW printed on a bianisotropic medium. The analytical formulation is based on a Full-GEMT, a method that avoids usual procedures of heavy and tedious mathematical expressions in the development of calculations and uses matrix-based mathematical expressions instead. These particularities are exploited to develop a mathematical model for the characterization of wave propagation in a three-layer shielded suspended CPW structure. This study is based on the development of mathematical formulations in full compact matrix-based expressions resulting in Green's functions in a matrix form. The implemented method incorporates a new accelerating procedure developed in the GEMT which provides an initial value used to speed up searching for the exact solution in the principal computation code. This helped us to obtain accurate solutions with tolerable computing time. Good agreements have been achieved with the literature in terms of accuracy and rapid convergence. The results for different cases of bianisotropy have been investigated, and particularly, the effect on the dispersion characteristics is presented and compared with the isotropic case.

**Keywords:** chiral; Tellegen; multilayer CPW structure; dispersion characteristics; full-GEMT; Muller's method; complex propagation constant; acceleration procedure

#### **1. Introduction**

To establish strong foundations for the development of modern and mass-market applications in the field of telecommunications, microwave designers have to develop further efficient devices, which aim to meet the specific needs of modern telecommunication systems, particularly the 5G

technology. In addition to providing far better levels of reliability and performances by offering new services, the new technology should be fully consistent with traditional services including 2G, 3G, 4G, Wi-Fi, and other relevant wireless systems, which requires suitable and high performance microwave devices in terms of miniaturization and ultra-wideband characteristics that are proven to be the most challenging issues of all time.

An important class of existing microwave devices is that exploiting the particular properties of bianisotropic media [1] for the development of special and innovative devices that may respond to the needs of modern technologies [2]. In this class of promising materials, we may mention, for example, non-reciprocal, gyrotropic, ferrites, chirals, metamaterials, metasurfaces, etc. [1,3]. Over the last three decades, the electromagnetism of bianisotropic media have gained a great deal of interest from scientists and researchers within the frame of artificial media with new and exciting properties [4]. However, practical exploitation of these phenomena did not develop on a large scale; many physics and engineering problems needed to be solved. Recently, as the science of materials has tremendously advanced, the concept of bianisotropic media has substantially reemerged as a field of importance in microwaves and optics technology [5–7].

The particular properties of bianisotropic media arise from a coupling between the electric and magnetic fields that can be explicitly described by general constitutive relations. Due to their diversity, they have found many potential applications from microwaves to optical frequencies such as polarization transformers, directional couplers, antenna and transmission line substrates, antenna radomes, radar systems, chiral waveguides and others [8–13].

The electromagnetic properties of bianisotropic media should be analyzed to perceive their exotic characteristics. Several studies have been conducted to characterize the electromagnetic behavior of bianisotropic media [14–18], ferrites [19], metamaterials [20], chiral [21,22], nonreciprocal [23] for simple and complex dielectric based microwave planar structures using numerical and analytical methods [8,14,17,20,24–35]. In [8] and [35], the method of lines is used to analyze planar transmission lines with conductor losses and to analyze integrated optical waveguide structures, respectively. In [17] and [20], the transmission line matrix (TLM) method is used for modeling dispersive chiral media and the analysis of dispersion in metamaterials. In [25], a fast computation of planar microstrip lines using the generalized equivalent circuit method of moments is presented. The finite difference technique and the iterated moment method are utilized for the analysis of dielectric [33] and optical waveguides [34]. In [27], a complex image method based on genetic algorithm (GA) is proposed to calculate the Green's functions of a coplanar waveguide structure. Anisotropic based multilayers, microstrip and waveguide structures are treated in [28,29,31,32]. Recently, in [14] and [24], Karma et al. studied microstrip transmission lines with anisotropic and uniaxial anisotropic substrates using the discrete mode matching method.

To extract the effective constitutive parameters of bianisotropic materials, various techniques such as stepwise method, S-parameters method, resonator method, coaxial probe method, free-space characterization method, rectangular waveguide measurements and recursive algorithms have been employed [36–39].

By knowing the intrinsic physical properties of complex media, designers can predict the response of microwave components for the development of inventive devices. However, the complexity of mathematical modeling of bianisotropic media is a real challenge in the characterization of microwave components. This has recently become a central area of research in microwaves and optics.

This paper introduces a mathematical modeling of complex media characterized by full 3 × 3 bianisotropic tensors of permittivity, permeability and magneto-electric parameters. The objective of this work is to sufficiently develop predictive mathematical models to judiciously characterize the propagation of electromagnetic waves in a suspended shielded bianisotropic three-layer CPW structure using a Full-GEM technique [28,40,41]. Three primary considerations have been exploited to accomplish new and promising results. These include: the consideration of the most general reciprocal and non-reciprocal chiral and achiral (Tellegen) complex media using the Full-GEMT, estimation of the

effective permittivity constant to be used as an initial value to search for the exact solution, and the use of Muller's method for the extraction of the complex solution of the associated propagation constant in both chiral and achiral media.

The efficient spectral Galerkin-based method of moments (SGMoM) is extensively used to analyze microwave planar structures [18,22,25,27,29,30,41–48]. In our recent work [42], we presented an analytical modeling of a shielded microstrip line based on an anisotropic medium with full 3 × 3 permittivity and permeability tensors. To reduce calculations, some conditions on the permittivity and permeability constants were considered ensuring the decoupling of the TE and TM modes. On the other hand, the computational time drastically increases with the required accuracy because of the slow convergence of integrals and series summations making this approach time intensive [43]. For accelerated convergence and efficient computation, several techniques have been used [44–47]. The time intensive part in the SGMoM is the evaluation of the matrix elements and the determinant calculation, as the matrix size is large in most cases. In [48], the SGMoM calculation of the propagation constant (β) for a shielded microstrip line is accelerated using asymptotic expansions for the Bessel's and the Green's functions with the aid of super convergent series in the approximation of the summation of the leading terms. In [49], for the multilayered shielded microstrip analysis, series summation calculation are accelerated using the Levin's transformation in the spectral method. In [28], the Green's function-based volume integral equation computation is accelerated using the fast Fourier transform technique. In [44,45], the impedance matrix elements based on Sommerfeld-type infinite double integral Green functions evaluation is accelerated by converting the infinite double integral of the impedance-matrix elements into a finite one-dimensional integral by using the asymptotic Green's functions and triangular basis functions with edge condition. In [47], integral-equation formulation in the SGMoM is accelerated by extracting suitable half-space parts of the kernels which leads to an exponentially decaying integrand functions. The integrals of the extracted parts are expressed as combinations of proper integrals and fast converging improper integrals. In [50], an efficient quasi-static analysis is presented, which can be used for speeding up full-wave SGMoM computations as well. Accelerated versions of full-wave spectral domain approach are also reported in [51–54].

In this paper, we propose a novel approach for the numerical acceleration of the SGMoM for the analysis of bianisotropic medium-based microstrip structures by accelerating and fixing problems of the convergence of the series summation in the elements of the Galerkin's matrix based on Green's functions.

The herein considered complex medium is bianisotropic with non-zero magneto-electric tensors (ξ - 0 and η - 0 at the same time). Our previous study, presented in [41], did not treat the case of bianisotropic media with both magneto-electric tensors; only one tensor was considered non zero. The relative simplicity of this case of medium does not require a more complex resolution technique or longer calculation time. Note that this technique failed to provide accurate solutions for general complex bianisotropic media due to the round-off errors and the highly oscillatory fields behavior. A new procedure is implemented to improve the technique and expand it to support the general case of bianisotropy for a CPW structure. Due to the complexity of the considered bianisotropic medium, the resolution method has required a more efficient technique to overcome the drawbacks in terms of non-convergence or considerable calculation time for a tolerable accuracy. An improvement is achieved by introducing an intermediate calculation procedure based on the GEMT to retrieve an approximated initial value of the relative effective permittivity of the three layer-structure as a function of ε*r*, μ*r*, ξ and η: the bianisotropic layer constitutive parameters. This value is used for searching for the exact solutions of the normalized complex constant of propagation ((β/κ0) <sup>2</sup> and (α/κ0) 2).

#### **2. Exponential Matrix Technique Formulation**

The general CPW geometry and the appropriate coordinate system with the z-axis as the direction of propagation are shown in Figure 1. The considered structure is based on a complex bianisotropic medium (region 1) characterized by full 3 × 3-magneto-electric tensors expressing the cross coupling between electric and magnetic fields.

**Figure 1.** Geometry of the shielded suspended 3-layer CPW structure.

Bianisotropic materials, in their general form, are characterized by the following constitutive relations [30,41,45].

$$\begin{array}{lcl}\overrightarrow{D} &=& \varepsilon \omicron [\varepsilon] \overrightarrow{E} + \sqrt{\varepsilon \omicron \mu\_0} [\xi] \overrightarrow{H} \\ \xrightarrow{\longrightarrow} & \mu\_0 [\mu] \overrightarrow{H} + \sqrt{\varepsilon \underline{\mu\_0}} [\eta] \overrightarrow{E} \end{array} \tag{1}$$

The tensors of the relative permittivity [ε], relative permeability [μ] and magneto-electric elements [ξ] and [η] are represented in the Cartesian coordinate system as follows:

$$
\begin{array}{rcl}
\psi &=& \begin{bmatrix}
\psi\_{xx} & \psi\_{xy} & \psi\_{xz} \\
\psi\_{yx} & \psi\_{yy} & \psi\_{yz} \\
\psi\_{zx} & \psi\_{zy} & \psi\_{zz}
\end{bmatrix}
\end{array}
\tag{2}
$$

where ψ stands for [ε], [μ], [ξ], or [η].

Starting from Maxwell's equations and using the GEMT in the spectral domain, we come to four coupled first-order differential equations for the transverse electromagnetic field components as functions of their derivatives [40,41] given in the Fourier domain:

$$\frac{\partial \left[ \widetilde{f}^{(i)}(a,\beta,z) \right]}{\partial z} = \left[ P^{(i)} \right]\_{4 \times 4} \left[ \widetilde{f}^{(i)}(a,\beta,z) \right]\_{\prime} \tag{3}$$

α and β are the Fourier variables corresponding to the space domain wavenumbers κ*<sup>x</sup>* and κ*y*, with

$$\begin{bmatrix} \overleftarrow{f}^{(i)}(\alpha,\beta,z) \end{bmatrix} = \begin{bmatrix} \overleftarrow{E}\_{\text{x}}^{(i)}(\alpha,\beta,z) \\ \overleftarrow{E}\_{\text{y}}^{(i)}(\alpha,\beta,z) \\ \overleftarrow{H}\_{\text{x}}^{(i)}(\alpha,\beta,z) \\ \overleftarrow{H}\_{\text{y}}^{(i)}(\alpha,\beta,z) \end{bmatrix}' \tag{4}$$

and

$$\left[ \left[ P \right]\_{4 \times 4} = \; \not\text{\&} \left\{ \left[ \mathbf{M} \right] + \left[ \mathbf{N} \right] \left[ \mathbf{Q} \right] \left[ \mathbf{R} \right] \right\} \right] \tag{5}$$

where

$$\begin{aligned} \text{[M]} \quad = \begin{bmatrix} -\eta\_{yx} & -\eta\_{yy} & -\eta\_{\ell}\mu\_{yx} & -\eta\_{\ell}\mu\_{yy} \\ \eta\_{xx} & \eta\_{xy} & \eta\_{\ell}\mu\_{xx} & \eta\_{\ell}\mu\_{xy} \\ \chi\_{0}\varepsilon\_{yx} & \chi\_{0}\varepsilon\_{yy} & \xi\_{yx} & \xi\_{yy} \\ -\chi\_{0}\varepsilon\_{xx} & -\chi\_{0}\varepsilon\_{xy} & -\xi\_{xx} & -\xi\_{xy} \end{bmatrix} = \begin{bmatrix} [\eta\_{T}] & \imath\_{0}[\mu\_{T}] \\ -\varchi\_{0}[\varepsilon\_{T}] & -[\xi\_{T}] \end{bmatrix} \end{aligned} \tag{6a}$$

$$\begin{array}{rcl} \text{[N]} &=& \begin{bmatrix} -\left(\eta\_{yz} + \mathbf{x}\_{x}^{n}\right) & -\eta\_{0}\mu\_{yz} \\ \left(\eta\_{xz} - \mathbf{x}\_{y}^{n}\right) & \mathbf{0}\mu\_{xz} \\ \mathbf{Y}\_{0}\varepsilon\_{yz} & \left(\underline{\boldsymbol{\zeta}}\_{yz} - \mathbf{x}\_{x}^{n}\right) \\ -\mathbf{Y}\_{0}\varepsilon\_{xz} & -\left(\underline{\boldsymbol{\zeta}}\_{xz} + \mathbf{x}\_{y}^{n}\right) \end{bmatrix} \end{array} \tag{6b}$$

$$[\mathbf{Q}] = \frac{-1}{\left(\varepsilon\_{zz}\mu\_{zz} - \underline{\zeta}\_{zz}\eta\_{zz}\right)} \begin{bmatrix} \, \_0\mu\_{zz} & \, \_{\overline{\zeta}\_{zz}} \\ \, \_0\eta\_{zz} & \, \_1\chi\_{\overline{\zeta}\_{zz}} \end{bmatrix} \tag{6c}$$

$$\begin{array}{rcl} \text{[R]} &=& \begin{bmatrix} \chi\_{0\varepsilon\_{xx}} & \chi\_{0\varepsilon\_{yz}} & \left(\xi\_{xx} - \kappa\_y^n\right) & \left(\xi\_{yz} + \kappa\_x^n\right) \\ -\left(\eta\_{xx} + \kappa\_y^n\right) & -\left(\eta\_{yz} - \kappa\_x^n\right) & -\eta\_{0\varepsilon\_{xz}} & -\eta\_{0}\mu\_{yz} \end{bmatrix} \end{array} \tag{6d}$$

$$
\kappa\_x^n = \frac{\kappa\_x}{\kappa\_0},\tag{6e}
$$

$$
\kappa\_y^n = \frac{\kappa\_y}{\kappa\_0},
\tag{6f}
$$

$$
\rho\_0 = \frac{1}{Y\_0} = \sqrt{\frac{\mu\_0}{\varepsilon\_0}},\tag{6g}
$$

with

$$
\begin{bmatrix} \eta \mathbf{r} \end{bmatrix} = \begin{bmatrix} -\eta\_{\text{yx}} & -\eta\_{\text{yy}} \\ \eta\_{\text{xx}} & \eta\_{\text{xy}} \end{bmatrix} \text{ [ $\varepsilon\text{r}$ ] } = \begin{bmatrix} -\varepsilon\_{\text{yx}} & -\varepsilon\_{\text{yy}} \\ \varepsilon\_{\text{xx}} & \varepsilon\_{\text{xy}} \end{bmatrix} \text{ [ $\xi\text{r}$ ] } = \begin{bmatrix} -\xi\_{\text{yx}} & -\xi\_{\text{yy}} \\ \xi\_{\text{xx}} & \xi\_{\text{xy}} \end{bmatrix} \text{ [ $\mu\text{r}$ ] } = \begin{bmatrix} -\mu\_{\text{yx}} & -\mu\_{\text{yy}} \\ \mu\_{\text{xx}} & \mu\_{\text{xy}} \end{bmatrix} \text{ (\"\epsilon\text{h}\") } = \begin{bmatrix} \varepsilon\_{\text{x}} & -\mu\_{\text{xy}} \\ \mu\_{\text{xx}} & \mu\_{\text{xy}} \end{bmatrix} \text{ (\"\epsilon\text{h}\") } = \begin{bmatrix} \varepsilon\_{\text{x}} & -\mu\_{\text{xy}} \\ \varepsilon\_{\text{xx}} & \mu\_{\text{xy}} \end{bmatrix} \text{ (\"\epsilon\text{h}\") } = \begin{bmatrix} \varepsilon\_{\text{x}} & -\mu\_{\text{xy}} \\ \mu\_{\text{xx}} & \mu\_{\text{xy}} \end{bmatrix} \text{ (\"\epsilon\text{h}\") } = \begin{bmatrix} \varepsilon\_{\text{x}} & -\mu\_{\text{xy}} \\ \mu\_{\text{xx}} & \mu\_{\text{xy}} \end{bmatrix} \text{ (\"\epsilon\text{h}\") } = \begin{bmatrix} \varepsilon\_{\text{x}} & -\mu\_{\text{xy}} \\ \mu\_{\text{xx}} & \mu\_{\text{xy}} \end$$

where κ<sup>0</sup> is the free space wavenumber and ω is the angular frequency.

This study is essentially based on the development of mathematical formulations in compact matrix-based forms; this is deemed as a promising approach, since it avoids excessive and complex calculation developments. This can dramatically reduce the complexity of wave propagation modeling in complex media.

The matrix [*P*] (Equation (5)) is the first foundation for this technique associated with the studied bianisotropic-medium based CPW structure. Its elements are given as functions of the constitutive tensors elements. In previous works [30,45], for particular cases of media calculations were explicitly developed, which is not obvious with heavy mathematical calculations case that characterize the herein studied complex bianisotropic structure.

Equation (3) admits a general solution of the form:

$$
\begin{bmatrix} E\_x(z) \\ E\_y(z) \\ H\_x(z) \\ H\_y(z) \end{bmatrix} = \begin{bmatrix} T(\kappa\_x, \kappa\_y, z) \\ E\_y(0) \\ H\_x(0) \\ H\_y(0) \end{bmatrix} \tag{7}
$$

with

$$T(\kappa\_{x\_{\prime}}\kappa\_{y\_{\prime}}z) = \exp([P]\cdot z),\tag{8}$$

and

$$
\widetilde{f}^{(i)}(\alpha,\beta,z(i)) = \, \_T (\kappa\_{x\_\prime} \kappa\_{y\_\prime} z) \widetilde{f}^{(i)}(\alpha,\beta,z(i-1)) . \tag{9}
$$

The 4 × 4 transfer matrix *T* κ*x*, κ*y*; *z* is calculated in the formulation of the GEMT by means of the Cayley Hamilton theorem for the determination of the complex function roots [40]. It is expressed in the following polynomial form:

$$T(z) \ = \ a\_0[I] + a\_1[P] + a\_2[P]^2 + a\_3[P]^3,\tag{10}$$

where *aj* are scalar expansion coefficients, determined by solving the Vandermode linear algebraic system [40], and [*I*] is a 4 × 4 identity matrix. The transfer matrix is easily obtained by multiplying the different transfer matrices related to the different layers of the structure, this constitutes the main advantage of this new technique. By imposing the appropriate boundary conditions between the heterogeneous medium layers, the appropriate Green's tensor which models the CPW structure is derived. Details can be found in [41]. This technique exhibits a compact matrix form with the advantage of being easily inserted in the calculation code.

#### *2.1. Implementation of the Acceleration Procedure*

To overcome the drawbacks of the resolution method used in [41], when applied for a general complex bianisotropic medium, the resolution method has to be improved in terms of convergence and accuracy for a tolerable computing time. A new procedure is introduced in the calculation technique. This latter is based on the GEMT technique (detailed calculations can be found in [40]) used to retrieve an approximated value for the effective relative permittivity of the whole inhomogeneous structure to be used as an initial value to search for the exact solution of the propagation constant. This value is evaluated for each frequency point by extracting the eigenvalues of matrix [*P*] by resolving Equation (12) for κ*<sup>x</sup>* = 0 and κ*<sup>y</sup>* = 0. By applying this procedure, analytical expressions of the effective relative permittivity maybe obtained as functions of the constitutive parameters of the bianisotropic layer ε*r*, μ*r*, ξ and η. The application of this procedure allows the acceleration of the Matlab® [55] calculation code and provides a better solution accuracy.

The expansion coefficients *ai* (*i* = 0, 1, 2, 3) in Equation (10) are determined by solving the Vandermode linear algebraic system:

$$
\begin{bmatrix} 1 & \lambda\_0 & \lambda\_0^2 & \lambda\_0^3 \\ 1 & \lambda\_1 & \lambda\_1^2 & \lambda\_1^3 \\ 1 & \lambda\_2 & \lambda\_2^2 & \lambda\_2^3 \\ 1 & \lambda\_3 & \lambda\_3^2 & \lambda\_3^3 \end{bmatrix} \begin{bmatrix} a\_0(z) \\ a\_1(z) \\ a\_2(z) \\ a\_3(z) \end{bmatrix} = \begin{bmatrix} \exp(\lambda\_0 z) \\ \exp(\lambda\_1 z) \\ \exp(\lambda\_2 z) \\ \exp(\lambda\_3 z) \end{bmatrix} \tag{11}
$$

λ*<sup>i</sup>* (*i* = 0, 1, 2, 3): are eigenvalues of [*P*] which correspond to propagating waves [56] defined by:

$$\det(\lambda[I] - [P]) = \lambda^4 + a\mathbf{1}\lambda^3 + a\mathbf{2}\lambda^2 + a\mathbf{3}\lambda + a\mathbf{4} = \mathbf{0}.\tag{12}$$

The coefficients α*<sup>i</sup>* (*i* = 0, 1, 2, 3) are given in terms of the matrix [*P*], explicit expressions may be found in [40].

#### *2.2. Derivation of the Initial Value Expression of the E*ff*ective Relative Permittivity*

#### *a. Isotropic case*

The initial effective permittivity value expression is calculated in terms of the medium constitutive parameters ε*r*, μ*r*, ξ and η using the total transfer matrix (Equation (10)) for κ*<sup>x</sup>* = 0 and κ*<sup>y</sup>* = 0, which permits the extraction of an approximated value. As an example of calculations, we present the derived analytical expressions of the initial effective permittivity of the isotropic and some bianisotropic cases. For the isotropic case, the derived matrix [*P*] is given by:

$$\begin{aligned} \mathbf{[}P\] &= \begin{bmatrix} 0 & 0 & 0 & -\mu\_r Z\_0 \kappa\_0 \\ 0 & 0 & \mu\_r Z\_0 \kappa\_0 & 0 \\ 0 & \varepsilon\_r \kappa\_0 / Z\_0 & 0 & 0 \\ -\varepsilon\_r \kappa\_0 / Z\_0 & 0 & 0 & 0 \end{bmatrix} . \end{aligned} \tag{13a}$$

In this case, the normalized eigenvalues of matrix [*P*], with respect to κ0, are:

$$
\lambda\_{n0} = \sqrt{\varepsilon\_l \mu\_r}, \; \lambda\_{n1} = -\sqrt{\varepsilon\_l \mu\_r}, \; \lambda\_{n2} = \; \sqrt{\varepsilon\_l \mu\_r} \text{ and } \lambda\_{n3} = -\sqrt{\varepsilon\_l \mu\_r}; \tag{13b}
$$

functions of the relative permittivity and permeability of the isotropic medium. For the fundamental propagating mode, the numerical maximal value is taken as the initial value.

$$
\varepsilon\_{ref0} = \lambda\_{n0}^2 = \varepsilon\_r \mu\_r \tag{13c}
$$

#### *b. Uniaxial anisotropic case*

*Electronics* **2020**, *9*, 243

In this example, we consider the following special case of biaxial anisotropy

$$\begin{aligned} \begin{bmatrix} \varepsilon \end{bmatrix} &= \begin{bmatrix} \varepsilon\_x & 0 & 0\\ 0 & \varepsilon\_y & 0\\ 0 & 0 & \varepsilon\_z \end{bmatrix} \begin{bmatrix} \mu \end{bmatrix} = \begin{bmatrix} \mu\_x & 0 & 0\\ 0 & \mu\_y & 0\\ 0 & 0 & \mu\_z \end{bmatrix} \begin{bmatrix} \xi \end{bmatrix} = \begin{array}{c} 0, \ [\eta] \end{array} = \begin{array}{c} \tag{14a} \end{array} \tag{14a} \end{aligned} \tag{14b}$$

the derived matrix [*P*] is:

$$
\begin{bmatrix}
[P\] \ &= \begin{bmatrix}
0 & 0 & 0 & -\mu\_y Z\_0 \kappa\_0 \\
0 & 0 & \mu\_x Z\_0 \kappa\_0 & 0 \\
0 & \varepsilon\_y \kappa\_0 / Z\_0 & 0 & 0 \\
\end{bmatrix}
\tag{14b}
$$

where the four normalized eigenvalues are found to be

$$\begin{array}{rcl}\lambda\_{n0} &=& \sqrt{\frac{\mu\_{\pi}}{\mu\_{\pi}} \{\varepsilon\_{y}\mu\_{z} - 1\}} \; \lambda\_{n1} =& -\sqrt{\frac{\mu\_{\pi}}{\mu\_{\pi}} \{\varepsilon\_{y}\mu\_{z} - 1\}}\\\lambda\_{n2} &=& \sqrt{\frac{\varepsilon\_{z}}{\varepsilon\_{z}} \{\varepsilon\_{z}\mu\_{y} - 1\}} \; \lambda\_{n3} =& -\sqrt{\frac{\varepsilon\_{z}}{\varepsilon\_{z}} \{\varepsilon\_{z}\mu\_{y} - 1\}} \end{array} \tag{14c}$$

and

$$\varepsilon\_{\rm rff}|\_{0} = \max \{ \lambda\_{\rm ni}^{2} \} = \max \left( \sqrt{\frac{\mu\_{\rm x}}{\mu\_{z}} (\varepsilon\_{\rm y} \mu\_{z} - 1)}, \sqrt{\frac{\varepsilon\_{\rm x}}{\varepsilon\_{z}} (\varepsilon\_{z} \mu\_{\rm y} - 1)} \right) \tag{14d}$$

*c. Diagonal bianisotropy case*

As an example, we take the following case:

$$\begin{aligned} \left[\boldsymbol{\varepsilon}\right] &= \left. \varepsilon\_{\mathrm{I}} \left[\boldsymbol{l}\right], \left[\boldsymbol{\mu}\right] = \left. \mu\_{\mathrm{I}} \left[\boldsymbol{l}\right], \left[\boldsymbol{\xi}\right] \right] = \left[\begin{array}{ccc} \boldsymbol{\xi}\_{\mathrm{xx}} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right], \left[\boldsymbol{\eta}\right] &= \left[\begin{array}{ccc} \eta\_{\mathrm{xx}} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{array}\right]. \end{aligned} \tag{15a}$$

where [*I*] is a 3 × 3 identity matrix. The derived [*P*] is:

$$\begin{bmatrix} \left[P\right] \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & -\mu\_r \mathbf{Z}\_0 \mathbf{x}\_0 \\ \eta\_{\text{xx}} \kappa\_0 & 0 & \mu\_r \mathbf{Z}\_0 \kappa\_0 & 0 \\ 0 & \varepsilon\_r \kappa\_0 / Z\_0 & 0 & 0 \\ -\varepsilon\_r \kappa\_0 / Z\_0 & 0 & -\xi\_{\text{xx}} \kappa\_0 & 0 \end{bmatrix},\tag{15b}$$

and the four normalized eigenvalues of [*P*] are:

$$\begin{array}{rcl}\lambda\_{n0} &=& \sqrt{\varepsilon\_r \mu\_r + \sqrt{\varepsilon\_r \mu\_r \xi\_{\text{xx}} \eta\_{\text{xx}}}} \\ \lambda\_{n2} &=& \sqrt{\varepsilon\_r \mu\_r - \sqrt{\varepsilon\_r \mu\_r \xi\_{\text{xx}} \eta\_{\text{xx}}}} \\ \end{array} \\ \lambda\_{n1} &=& -\sqrt{\varepsilon\_r \mu\_r + \sqrt{\varepsilon\_r \mu\_r \xi\_{\text{xx}} \eta\_{\text{xx}}}} \\ \lambda\_{n3} &=& -\sqrt{\varepsilon\_r \mu\_r - \sqrt{\varepsilon\_r \mu\_r \xi\_{\text{xx}} \eta\_{\text{xx}}}} \\ \end{pmatrix} \tag{15c}$$

The initial value depends on the constitutive parameters ε*r*, μ*r*, ξ*xx* and η*xx*.

$$
\varepsilon\_{\pi f f 0} = \max \{ \lambda\_{\text{ri}}^2 \} = \varepsilon\_r \mu\_r + \sqrt{\varepsilon\_r \mu\_r \xi\_{\text{xx}} \eta\_{\text{xx}}} \tag{15d}
$$

#### *d. Gyrotropic bianisotropy case*

Considering the following medium case:

$$\begin{aligned} \left[ \begin{smallmatrix} \varepsilon \end{smallmatrix} \right] &= \left. \varepsilon\_{\mathbb{I}} [\mathrm{I}] , \ [\mu] = \mu\_{\mathbb{I}} [\mathrm{I}] , \ [\xi] &= \left[ \begin{array}{ccc} 0 & 0 & \xi\_{\mathrm{xz}} \\ 0 & 0 & 0 \\ \xi\_{\mathrm{xz}} & 0 & 0 \end{array} \right] , \ [\eta] &= \left[ \begin{array}{ccc} 0 & 0 & \eta\_{\mathrm{xz}} \\ 0 & 0 & 0 \\ \eta\_{\mathrm{xz}} & 0 & 0 \end{array} \right] . \end{aligned} \tag{16a}$$

The derived corresponding [*P*] is:

$$\begin{array}{c} \begin{bmatrix} \begin{array}{cc} 0 & 0 & 0 & -\mu\_{r}Z\_{0}\kappa\_{0} \\ 0 & 0 & -\frac{\kappa\_{l}\overline{\upsilon}\_{\mathsf{x}}}{\varepsilon\_{r}}(\varepsilon\_{r}\mu\_{r}-\xi\_{zz}\eta\_{zz}) & 0 \\ 0 & \varepsilon\_{r}\kappa\_{0}/Z\_{0} & 0 & 0 \\ -\frac{\kappa\_{0}}{\mu\_{r}Z\_{0}}(\varepsilon\_{r}\mu\_{r}-\xi\_{zz}\eta\_{zz}) & 0 & 0 \end{array} \end{array} \end{array} \tag{16b}$$

which gives as solutions:

$$\begin{array}{rcl}\lambda\_{n0} &=& \sqrt{\varepsilon\_r \mu\_r - \underline{\xi}\_{zx} \eta\_{zx}} & \lambda\_{n1} &=& -\sqrt{\varepsilon\_r \mu\_r - \underline{\xi}\_{zx} \eta\_{zx}} \\ \lambda\_{n2} &=& \sqrt{\varepsilon\_r \mu\_r - \underline{\xi}\_{zx} \eta\_{zx}} & \lambda\_{n3} &=& -\sqrt{\varepsilon\_r \mu\_r - \underline{\xi}\_{zx} \eta\_{zx}} \end{array} \tag{16c}$$

and

$$\varepsilon\_{nff0} = \max \{ \lambda\_{ni}^2 \} = \max \{ \varepsilon\_I \mu\_{\varGamma} - \underline{\zeta}\_{\mathfrak{X}\mathfrak{X}} \eta\_{\varXi \urcorner \mathfrak{x} \urcorner} \varepsilon\_I \mu\_{\varGamma} - \underline{\zeta}\_{\mathfrak{X}\mathfrak{X}} \eta\_{\varXi \urcorner} \}. \tag{16d}$$

It is medium-case dependent.

Notice that if ξ*xz*, ξ*zx*, η*xz* and η*zx* are taken so that ξ*xz*η*zx* ξ*zx*η*xz*, two different solutions are obtained corresponding to bifurcating modes [21]. This may constitute an independent issue, which is outside the scope of this work. This shows the efficiency of the procedure, without which, the calculation code would diverge to adjacent solutions or give spurious ones [57], since the electromagnetic fields of bianisotropic media are highly oscillatory [56]. The result presented by Equation (15d) shows that the use of this approach has not only made it feasible to get the optimal initial value, in some cases, but also to more accurately infer the appropriate approximation, mainly in the presence of the bifurcating modes phenomenon, in which two neighboring modes are excited.

In order to show the benefits of using the new procedure, two examples of the studied cases of bianisotropy are considered. Case1: <sup>ξ</sup>*zx*,1 <sup>=</sup> <sup>ξ</sup>*xz*,1 <sup>=</sup> <sup>−</sup>0.75 <sup>√</sup>ε*r*, <sup>η</sup>*zx*,1 <sup>=</sup> <sup>−</sup>ξ*xz*,1 and <sup>η</sup>*zx*,1 <sup>=</sup> <sup>η</sup>*xz*,1, and Case2: η*xz*,8 = −ξ*xz*,8 = −*j* <sup>√</sup>ε*r*, <sup>ξ</sup>*zx*,8 <sup>=</sup> <sup>−</sup>ξ*xz*,8 and <sup>η</sup>*zx*,8 <sup>=</sup> <sup>η</sup>*xz*,8. By the introduction of the new procedure, the computing time for Case1, for a frequency point, is reduced by about 33% from 1.937 s to 1.302 s. In Case2, with the aid of the procedure, we get a solution in 1.442 s while without the procedure, the technique failed to find a solution and gave a spurious value instead, after a long execution time. This is due to the oscillating behavior of the series summations of the manipulated complex Galerkin's matrix.

#### **3. Method of Solution**

By applying the boundary conditions, the expressions of the electric and magnetic tangential components are evaluated at the interface air-dielectric in terms of the tangential current densities on the strips *jx* and *jy*. A matrix of the Green's tensor elements *Gij*(α*n*, β) for the CPW structure is achieved. It is arranged in the following system of equations

$$
\begin{bmatrix}
\overline{\dot{j}\_x} \\
\overline{\dot{j}\_y}
\end{bmatrix} = \frac{1}{\Delta\_\mathcal{G}} \begin{bmatrix}
\text{G}\_{22}(\boldsymbol{a}\_{\boldsymbol{n}}\boldsymbol{\beta}) & -\text{G}\_{12}(\boldsymbol{a}\_{\boldsymbol{n}}\boldsymbol{\beta}) \\
\end{bmatrix} \begin{bmatrix}
\overline{\mathcal{E}}\_x \\
\overline{\mathcal{E}}\_y
\end{bmatrix} \tag{17a}
$$

with

$$
\Delta\_{\mathbb{G}} = \mathcal{G}\_{11}\mathcal{G}\_{22} - \mathcal{G}\_{12}\mathcal{G}\_{21}.\tag{17b}
$$

and α*n*: the discrete Fourier transform variable with *n* the Fourier number of terms *n* = 1, 2, 3, ... , *N*.

For the resolution of the problem, the SGMoM is applied, the spectral electric field components are expanded in terms of trigonometric basis function sets [58,59]. A homogeneous system of linear equations arranged in a compact matrix form is derived [58]:

$$\begin{bmatrix} \begin{bmatrix} M(\beta) \end{bmatrix} \begin{bmatrix} a\_1 \ a\_2 \ \dots \ a\_M \ b\_1 \ b\_2 \ \dots \ b\_N \end{bmatrix}^T = \begin{bmatrix} 0 \end{bmatrix},\tag{18}$$

*Electronics* **2020**, *9*, 243

where

$$[\mathcal{M}(\beta)]\;=\begin{bmatrix}\mathcal{M}^{1,1}\_{q',p}\;\mathcal{M}^{1,2}\_{q',q}\\\mathcal{M}^{2,1}\_{p',p}\;\mathcal{M}^{1,1}\_{p',q}\end{bmatrix}\;\tag{19}$$

and

$$M\_{q'p}^{1,1}(\beta) = \sum\_{n} \frac{1}{\Delta\_{\mathcal{G}}} G\_{22}(\alpha\_{n'}\beta) \,\, \widetilde{l\_{x,p}} \,\widetilde{l\_{y,q'}^{\*}} \,\, \tag{20a}$$

$$M\_{q',q}^{1,2}(\beta) = -\sum\_{n} \frac{1}{\Delta\_G} G\_{12}(a\_{n'}\beta) \,\, \widetilde{f}\_{y,q}\widetilde{J}\_{y,q'}^\*;\tag{20b}$$

$$M\_{p',p}^{2,1}(\beta) = -\sum\_{\mathbf{n}} \frac{1}{\Delta\_{\mathbf{G}}} G\_{12}(a\_{\mathbf{n}}, \beta) \, \overleftarrow{J}\_{\mathbf{x},p} \overleftarrow{J}\_{\mathbf{x},p'}^\*;\tag{20c}$$

$$\mathcal{M}\_{p',q}^{2,2}(\beta) = \sum\_{\mathfrak{n}} \frac{1}{\Delta\_{\mathbb{G}}} G\_{11}(a\_{\mathfrak{n},\mathfrak{l}}\beta) \,\,\widetilde{f}\_{\mathfrak{l},\mathfrak{l}}\widetilde{f}\_{\mathfrak{x},p'}^{\*}.\tag{20d}$$

The system admits nontrivial solutions when *det* [*M*(β)] = 0 [41,45,58,60], from which the frequency dependent propagation constant can be determined. For lossy media, a complex constant solution is expected.

Using the new technique, original results for the dispersion characteristics of complex bianisotropic chiral and achiral media are obtained through the ratio (β/κ0) 2, discussed and compared with the isotropic case (ξ = η = 0) using the technique in [41].

Due to the great number of possible medium cases, we have restricted our analysis to highlighting the main results of achiral media that have been less addressed in the literature. Accurate solutions of the determinant roots in Equation (18) are obtained within a tolerance of 10<sup>−</sup>12.

#### **4. Results and Discussions**

In order to validate our calculations and test the efficiency of the proposed method, three magnetic anisotropic cases have initially been considered. Numerical results have been computed and compared with available literature [61,62] (Figure 2) and good agreements are observed. On the other hand, a rapid convergence has been achieved with a reduced Fourier number (N = 500) and basis functions (K = 8) against (N = 3000) used by Khodja et al. [61] for the same number of basis functions.

**Figure 2.** (β/κ0) <sup>2</sup> for the dominant mode of a shielded CPW with magnetic anisotropy, (2*a* = 3.556 mm, 2*w* = 2*s* = 0.7112 mm, *d*1 = 2.8448 mm, *d*2 = 0.7112 mm, *d*3 = 3.556 mm and ε*r* = 3).

In this study, we consider a suspended three-layer CPW structure implanted on a complex bianisotropic dielectric material (Figure 1) with the following geometrical dimensions *a* = 10 mm, *d*1 = 4.5 mm, *d*2 = 1 mm, *d*3 = 4.5 mm, *w* = 1 mm, *s* = 1 mm, ε*<sup>r</sup>* = 2.53, μ*<sup>r</sup>* = 1. Different sub-figures

are differentiated by the included legends where only the sign of the constitutive element changes respectively to the previous case in the same figure.

In order to examine the effect of the magneto-electric parameters on the dispersion characteristics, we first start with the examination of the axial bianisotropy effect. The magneto-electric elements, whether they are real, imaginary, positive or negative directly affect the phase constant as well as the attenuation factor.

The obtained results, treat two principal cases of diagonal bianisotropic medium: achiral and chiral. In each case, the magneto-electric pair (ξ*ij*, η*ij*) is considered non-zero, so that the new original results of the achiral medium case can be validated and compared with the chiral case. In addition, in this parametric study, we examined the effects of the gyrotropic elements of the magneto-electric tensors on the complex propagation coefficient. Results are grouped in figures according to the constitutive parameters effects.

#### *4.1. E*ff*ect of Diagonal Bianisotropy*

For diagonal bianisotropy three cases of magneto-electric elements are considered:


LVRWURSLF>@ [\\ H <sup>U</sup> [\\ H

 <sup>U</sup> [\\ H

 <sup>U</sup> [\\ H

where (*i* = *x*, *y*, *z*) and *a* = (−1, −0.75, −0.5, −0.25, 0.25, 0.5, 0.75, 1). Two main cases are distinguished: achiral (ξ*ii* = *a* <sup>√</sup>ε*r*) and chiral (ξ*ii* <sup>=</sup> *ja* <sup>√</sup>ε*r*).

In Figure 3, the effect of element ξ*yy*, for achiral and chiral media cases is presented. An identical effect is observed on the ratio (β/κ0) <sup>2</sup> (Figure 3a), with reciprocal effect (curves superposition for ξ = ±*a* <sup>√</sup>ε*r*), for both chiral and achiral cases. These media cases show low losses for achirality (case (i)) with η*yy* = ξ*yy* and almost zero losses for chirality with η*yy* = −ξ*yy* (case (ii)) (Figure 3b). It can be concluded that the effect of a reciprocal chirality is almost the same as for an achiral medium with equal magneto-electric elements. In these cases (*ii* = *yy*) propagating modes are excited in the guiding structure, however, for achiral with η*ii* = ξ*ii* and chiral with η*ii* = −ξ*ii*, no solutions are obtained for *ii* = *xx* and *ii* = *zz*, hence, the medium does not support any propagating modes.

> U

[\\ *M* H <sup>U</sup> [\\ *M* H <sup>U</sup> [\\ *M* H <sup>U</sup> [\\ *M* H U

**Figure 3.** Effect of diagonal bianisotropy ((i) ξ*yy* real with η*yy* = ξ*yy* and (ii) ξ*yy* imaginary with η*yy* = −ξ*yy*) on (**a**): (β/κ0) <sup>2</sup> and (**b**): (α/κ0) 2.

Figure 4 illustrates the effect of diagonal bianisotropic media for the case achiral with η*ii* = −ξ*ii*. Unlike the previous case, the fundamental propagating mode is excited for all diagonal magneto-electric

LVRWURSLF>@*LL [[\\]]*

elements (*ii* = *xx*, *yy*, *zz*). However, each of the elements has its own effect. According to Figure 4, for (*ii* = *xx*), a non-reciprocity for the achiral case (Figure 4a) is observed on (β/κ0) 2. Higher losses are observed for *a* ≥ 0.5 (Figure 4b). The effect of element ξ*yy*,2 is presented in Figure 4c and d. These cases are non-reciprocal and exhibit relatively lower losses. The effect on (β/κ0) <sup>2</sup> is almost identical for both cases ξ*yy*,2 and ξ*zz*,2 (Figure 4c,e).

**Figure 4.** Effect of diagonal bianisotropy elements, ξ*xx* on: (**a**) (β/κ0) <sup>2</sup> and (**b**) (α/κ0) 2; ξ*yy* on: (**c**) (β/κ0) 2 and (**d**) (α/κ0) 2; ξ*zz* on: (**e**) (β/κ0) <sup>2</sup> and (**f**) (α/κ0) <sup>2</sup> with <sup>η</sup>*xx* <sup>=</sup> <sup>−</sup>ξ*xx*.

#### *4.2. E*ff*ect of Gyrotropic Bianisotropy*

For the gyrotropic elements, five achiral cases are considered:

*Electronics* **2020**, *9*, 243


In Figure 5a–d, five cases were grouped, each differs from the other by a single change in sign of the magneto-electric element. It can be seen that the combination of the constitutive parameters shows non-reciprocity and a distinct effect on (β/κ0) <sup>2</sup> and (α/κ0) <sup>2</sup> parameters. In Figure 5a,b, for ! ! !ξ*xy*,1 ! ! ! and ! ! !ξ*xy*,2 ! ! ! <sup>≤</sup> 0.25 <sup>√</sup>ε*r*, we observe a weak effect on the ratio (β/κ0) <sup>2</sup> compared to the isotropic case with lower losses for ξ*xy*,2 (case (ii)). The sign change between η*xy*,1 = η*yx*,1 (case (i)) and η*xy*,2 = −η*yx*,2 (case (ii)) keeps the same effect on (β/κ0) <sup>2</sup> (Figure 5a). Only the ratio (α/κ0) <sup>2</sup> is affected with the appearance of non-reciprocity (Figure 5b). As shown in Figure 5c,d, for the combinations ξ*xy*,3, ξ*xy*,4 and ξ*xy*,5, a sign change has only a weak effect on (β/κ0) <sup>2</sup> and (α/κ0) 2.

**Figure 5.** Effect of gyrotropic achiral bianisotropy for different ξ*xy*, ξ*yx*, η*xy* and η*yx* combinations on: (**a**) (β/κ0) 2, (**b**) (α/κ0) 2, (**c**) (β/κ0) <sup>2</sup> and (**d**) (α/κ0) 2.

In Figure 6 are presented results of the element ξ*xz* combinations. The remarkable effect is that the ratio (β/κ0) <sup>2</sup> is almost constant with respect to frequency and decreases with ξ*xz*,1 to reach the unity for ! ! !ξ*xz*,1 ! ! ! <sup>=</sup> 0.75 <sup>√</sup>ε*<sup>r</sup>* and zero for ! ! !ξ*xz*,1 ! ! ! <sup>=</sup> <sup>√</sup>ε*<sup>r</sup>* (Figure 6a), all with negligible losses (Figure 6b,d). The elements, ξ*xz*,2 and ξ*xz*,3 have the same effect on both (β/κ0) <sup>2</sup> and (α/κ0) <sup>2</sup> (Figure 6c,d). In this case, only 3 combinations: ξ*xz*,1, ξ*xz*,2 and ξ*xz*,3 support propagating modes, with the appearance of the non-reciprocal effect. For the two other cases ξ*xz*,4 and ξ*xz*,5, no solutions are obtained. The (β/κ0) 2 and (α/κ0) <sup>2</sup> variations of both cases (ii) and (iii) are completely different from case (i). The exchange of sign between ξ*xz*,3 = −ξ*zx*,3 and η*xz*,2 = −η*zx*,2 preserved the same effect ((ii) and (iii) curves are superposed), as shown in Figure 6c,d.

**Figure 6.** Effect of gyrotropic achiral bianisotropy for different ξ*xz*, ξ*zx*, η*xz* and η*zx* combinations on: (**a**) (β/κ0) 2, (**b**) (α/κ0) 2, (**c**) (β/κ0) <sup>2</sup> and (**d**) (α/κ0) 2.

Figure 7a,b presents the effect of sign exchange between the magneto-electric element of achiral medium cases. The positive sign ξ*yz*,2 = η*yz*,2 and ξ*yz*,3 = η*yz*,3 reveals a non-reciprocity on (β/κ0) 2, while it shows a significant effect on (α/κ0) <sup>2</sup> and reduced losses compared to isotropic case.

**Figure 7.** Effect of gyrotropic achiral bianisotropy for different ξ*yz*, ξ*zy*, η*yz* and η*zy* combinations on: (**a**) (β/κ0) <sup>2</sup> and (**b**) (α/κ0) 2.

#### **5. Conclusions**

In this work, an analytical modeling of a three-layer CPW structure implanted on a complex medium is presented. This study is based on the Full-GEMT developed in a matrix form for the characterization of the bianisotropic-substrate CPW structure. This resulted in compact matrix form expressions of the Green's tensor. The implemented resolution method includes a new accelerating procedure developed in the GEMT that contributed to accomplishing accurate solutions with improved computing time. The computing time, for one frequency point, has been reduced by 33% for some calculation cases and more for others. A satisfactory calculation convergence is achieved with a significantly reduced Fourier number compared to literature. The presented technique can dramatically reduce the complexity of wave propagation characterization, in highly complex media, in terms of mathematical modelling and computational solution method.

According to our calculations, cases of complex media have been achieved, where the ratio (β/κ0) 2 is close to unity such as for cases ξ*xz*,1, ξ*xz*,2. and ξ*xz*,3. This characteristic is vital for the realization of media with permittivity close to unity for better use in the design of radiating antenna structures.

For cases ξ*xy*,1 and ξ*xy*,2, where there is only one element that changes sign between one case and another, the results of (β/κ0) <sup>2</sup> are similar (reversed cases with similar effects), however, (α/κ0) <sup>2</sup> presents a different variation either in form or in magnitude.

It is noted that for the case of achiral medium when ξ = η, ξ*ij* = ξ*ji* and η*ij* = −η*ji* (ξ*xy*,2), the medium is reciprocal and the effect is well reversed (non-reciprocal) when ξ*ij* = −ξ*ji* (ξ*xy*,3). For the cases ξ*xy*,3 and ξ*xy*,5, one had to have inverted cases, whereas, it is not the case. This can be explained by the properties imposed by the geometry of the studied structure.

An original result that should be taken into consideration is the losses, which show changes with each sign change between the magneto-electric elements. On the other hand, it is found that losses in achiral media are of the same magnitude as those in non-reciprocal chiral media, which must be taken into account. Furthermore, it is worth noting that the transverse elements ξ*xz* are the most influential on the phase constant in bianisotropic a CPW structure, and this may be attributed to the geometry of the studied structure.

Investigation of some achiral media cases has shown new results such as the notion of achiral media with a relative permittivity approaching unity with reduced losses. This new finding could serve as a valuable concept from which designers of unusual synthetic materials may benefit to enhance the media intrinsic properties for future innovative applications use.

Finally, the technique discussed in this paper may be extended to deal with propagation of bifurcated modes and enclosed multilayer microstrip structures.

**Author Contributions:** Design and concept, D.S. and C.Z.; methodology, C.Z. and I.E.; investigation, I.O., A.U.; resources, I.E. and J.R.; writing—original draft preparation, C.Z. and D.S.; writing—review and editing, I.E., A.U., H.A. and R.A.-A.; supervision, R.A.-A.; project administration, J.R., formal analysis H.A. and F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement H2020-MSCA-ITN-2016 SECRET-722424. This work is also funded by the FCT/MEC through national funds and when applicable co-financed by the ERDF, under the PT2020 Partnership Agreement under the UID/EEA/50008/2019 project.

**Acknowledgments:** This work is supported by the European Union's Horizon 2020 Research and Innovation program under grant agreement H2020-MSCA-ITN-2016-SECRET-722424. This work is also funded by the FCT/MEC through national funds and when applicable co-financed by the ERDF, under the PT2020 Partnership Agreement under the UID/EEA/50008/2019 project.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
