Complex Bianisotropy Effect on the Propagation Constant of a Shielded Multilayered Coplanar Waveguide Using Improved Full Generalized Exponential Matrix Technique

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.


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, 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][42][43][44][45][46][47][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][45][46][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][52][53][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 ) 2 and (α/κ 0 ) 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. Bianisotropic materials, in their general form, are characterized by the following constitutive relations [30,41,45].
where ψ stands for [ ] . 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: α and β are the Fourier variables corresponding to the space domain wavenumbers x κ and y  Bianisotropic materials, in their general form, are characterized by the following constitutive relations [30,41,45].
The tensors of the relative permittivity [ε], relative permeability [µ] and magneto-electric elements [ξ] and [η] are represented in the Cartesian coordinate system as follows: where ψ stands for 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: α and β are the Fourier variables corresponding to the space domain wavenumbers κ x and κ y , with and where with where κ 0 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: with and 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: where a j 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.

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 κ x = 0 and κ y = 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 a i (i = 0, 1, 2, 3) in Equation (10) are determined by solving the Vandermode linear algebraic system: : are eigenvalues of [P] which correspond to propagating waves [56] defined by: The coefficients α i (i = 0, 1, 2, 3) are given in terms of the matrix [P], explicit expressions may be found in [40].

Derivation of the Initial Value Expression of the Effective 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 κ x = 0 and κ y = 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: In this case, the normalized eigenvalues of matrix [P], with respect to κ 0 , are: 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.

. Uniaxial anisotropic case
In this example, we consider the following special case of biaxial anisotropy the derived matrix [P] is: where the four normalized eigenvalues are found to be and c. Diagonal bianisotropy case As an example, we take the following case: where [I] is a 3 × 3 identity matrix. The derived [P] is: and the four normalized eigenvalues of [P] are: The initial value depends on the constitutive parameters ε r , µ r , ξ xx and η xx .

d. Gyrotropic bianisotropy case
Considering the following medium case: The derived corresponding [P] is: which gives as solutions: and 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: and Case2: 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.

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 j x and j y . A matrix of the Green's tensor elements G ij (α n , β) for the CPW structure is achieved. It is arranged in the following system of equations with 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]: where and 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 −12 .

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. 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. 20   In this study, we consider a suspended three-layer CPW structure implanted on a complex bianisotropic dielectric material (Figure 1 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 magnetoelectric elements, whether they are real, imaginary, positive or negative directly affect the phase 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, d1 = 4.5 mm, d2 = 1 mm, d3 = 4.5 mm, w = 1 mm, s = 1 mm, ε r = 2.53, µ r = 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.
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 ) 2 (Figure 3a), with reciprocal effect (curves superposition for ξ = ±a √ ε 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.      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 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 ) 2 is almost identical for both cases ξ yy,2 and ξ zz,2 (Figure 4c,e). isotropic [41], ii= (xx, yy, zz)

Effect of Gyrotropic Bianisotropy
For the gyrotropic elements, five achiral cases are considered: 1.

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 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

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 ) 2 are similar (reversed cases with similar effects), however, (α/κ 0 ) 2 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. 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.