© M. Nitas et al., Published by EDP Sciences, 2024

1 Introduction

Homogenization of composite periodic media has been a very popular topic during the last decades. A great variety of techniques has been proposed for the retrieval of their effective parameters: Numerical, analytical, semi-analytical and also experimental schemes were introduced in many studies, capable of returning the effective parameters of effectively isotropic, anisotoropic and even bianisotropic composite periodic structures. One of the most popular approaches has been the S-parameter technique, which in its simplest form, returns the effective permittivity and permeability of a composite periodic medium by inverting the Fresnel coefficients relations on a homogeneous slab [1,2]. Based on the Nicolson-Ross-Weir (NRW) method [3,4], these techniques were subsequently extended to bianisotropic structures [5–12]. Another significant analytical approach is based on scatterer interaction models with the utilization of dynamic interaction coefficients. These techniques, initially proposed in [13] and subsequently extended to general bianisotropic scatterers, are generalizations of the Clausius-Mossotti equation and assume that the scatterers may be properly modeled via dipole approximation [14–18]. Yet an interesting homogenization technique based on the so-called T-Matrix approach was reported [19]. In addition to these approaches, a retrieval technique for composite periodic media utilizing the information of the dispersion diagrams and the averaging of the field components of the supported electromagnetic modes was recently reported [20,21]. This technique is capable of accurately returning a total of fifteen effective parameters: three electric, three magnetic and nine magnetoelectric.

A key issue in all homogenization techniques is the accuracy of the returned values of the effective parameters, which may be assessed by specific criteria: Although based on rigorous mathematical derivations and utilizing solid electromagnetic theory relations, many of them provide ambiguous results, for both real and imaginary parts. In this work, we attempt to reconstruct the dispersion diagrams of a composite periodic medium (corresponding to every path of the irreducible Brillouin zone), providing, in this way, a means of confirmation of the accuracy of its retrieved homogenized effective parameters. Reconstruction of the dispersion curves is carried out with the utilization of the analytical solutions for a general homogeneous medium possessing arbitrary bianisotropy [22] combined with the retrieved effective parameters, acquired through our formerly proposed fully numerical parameter retrieval technique. To this end, we obtain the complex propagation constants at every frequency for the homogenized composite medium, assuming wave propagation at all three orthogonal directions of incidence. The reconstructed dispersion curves corresponding to all of the supported modes of the simulated composite medium are in excellent agreement with the original ones returned by the numerical solution of the composite’s medium unit cell. Impressively, there is also an excellent match between the original and the reconstructed signs of the returned real and imaginary parts of the complex propagation constants. Furthermore, it is noted that no simplification or other analytical or numerical assumption is made, as all effective parameters are involved exactly as returned by the parameter retrieval technique. The significance of the reported dispersion reconstruction technique may be viewed from yet a main aspect of general importance in the theory of electromagnetic homogenization techniques. The excellent agreement between the reconstructed and original dispersion diagrams clearly confirms the validity of the underlying homogenization technique introduced in [20]. Thus, its accuracy may be quantified in terms of dispersion diagram comparisons, therefore it is also argued that the proposed reconstruction technique may serve as a means for any homogenization technique, to confirm the accuracy of their returned effective parameters.

The present manuscript is organized as follows: In Section 2 we briefly present our parameter retrieval technique, which is based on eigenvalue analysis of the composite periodic media and averaging of the field components of the supported modes. Taking outset in the previously proposed field-flux FEM formulation [23–25], we first illustrate the dispersion diagrams for the three orthogonal incidences for a bianisotropic composite medium inhere composed of Edge-Coupled-Split-Ring Resonators (EC-SRRs). Next, the process of field averaging used to obtained the average field components of the supported modes is illustrated. To this end, alternative integration parts are proposed, in order to derive the final system of equations with the known average fields and the unknown effective material parameters. Solution of this system retrieves the fifteen complex material parameters of the EC-SRR composite medium: three diagonal elements for the permittivity and permability tensors, respectively, as well as nine elements for the magnetoelectric coupling tensor.

In Section 3, we initially extract the analytical relations which characterize the wave propagation in a homogeneous bianisotropic medium. Subsequently, we involve the retrieved effective parameters of the EC-SRR composite medium in these analytical relations and calculate the complex propagation constants for the three orthogonal cases of incidence. In this way, we reconstruct the dispersion diagrams for every supported mode. There, we compare the original and reconstructed versions of the dispersion curves and calculate the relative errors between them. Finally, we present a rigorous explanation about the slope difference between the in-plane incidence lightline mode dispersion curves and the corresponding ones supported by the normal wave incidence.

2 Homogenization of composite periodic media

2.1 Computational tool

We utilize the dual field-flux periodic eigenmode formulation, originally proposed in [20,23], as the computational approach to acquire the dispersion diagrams and the field distributions of the modes supported by the composite medium. We note in passing that all simulations were carried out in the Weak Form of Comsol Multiphysics^TM, whereas the details of the formulations’ derivation may be found in the aforementioned references. In the FEM formulations, the determined field variables are the Bloch-Floquet periodic envelopes of the electric field e, the magnetic field h, the dielectric displacement d and the magnetic induction b. They are related to the corresponding field components E, H, D and B via the Bloch-Floquet’s theorem

$$E=e{e}^{-jk\cdot r}$$(1a)

$$H=h{e}^{-jk\cdot r}$$(1b)

$$D=d{e}^{-jk\cdot r}$$(1c)

$$B=b{e}^{-jk\cdot r},$$(1d)

where k = k $\widehat{\text{k}}$ is the Bloch-Floquet wavevector of prescribed propagation direction $\widehat{\text{k}}$ (e.g. $\widehat{\text{k}}$ = $\widehat{\text{x}}$) and k = β − jα is the unknown complex propagation constant.

For the illustration of our results, we now consider a composite medium consisting of EC-SRRs, periodic in all three dimensions (x, y and z). Its unit cell and geometric details are depicted in Figure 1. The SRRs are located in an air-host and are further made of copper (conductivity σ = 5.88 × 10⁷ S/m), with a thickness of 17 μm. The resulting dispersion diagrams (propagation constants β and attenuation constants α) for all three orthogonal incidences k_x, k_y and k_z − corresponding to the paths ΓX, XM and MR  of the irreducible Brillouin Zone, respectively − are illustrated in Figures 2, 3 and 4, respectively. Two modes are found to appear for every incidence case: For the in-plane cases (k_x and k_y) a lightline mode and an SRR mode appears, whereas in the normal incidence case (k_z) a straight-line mode with modified characteristics is returned, which interacts with the composite medium (i.e., it is not a lightline mode) along with the SRR mode corresponding to the resonant dispersion curve [26]. It is noted that the computational domain was discretized in 154 726 tetrahedral elements, which correspond to 517 446 degrees of freedom. The final eigenvalue matrix problems were solved with a relative tolerance of 10⁻⁹.

Fig. 1 (a) Unit cell of the EC-SRR composite medium. The resonators are arranged on a cubic lattice (the edge of the cube is d = dx = dy = dz = 8.8 mm). (b) Dimensions of the EC-SRR (in mm): The metallic rings are made of copper and have a thickness of 17 μm.

Fig. 2 Dispersion diagram of the EC-SRR composite medium for kx incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). Propagating SRR modes are plotted in blue, evanescent SRR modes with large propagation constant in green and evanescent SRR modes with small propagation constant in red. A frequency bandgap is present between 5.3 and 6 GHz. Lightline modes are illustrated in black.

Fig. 3 Dispersion diagram of the EC-SRR composite medium for ky incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). A frequency bandgap is observed between 5.3 and 6 GHz. Coloring is as in Figure 2.

Fig. 4 Dispersion diagram of the EC-SRR composite medium for kz (normal) incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). A frequency bandgap is observed between 5.7 and 6 GHz. The modes forming the straight line constitute the first supported mode whereas the two-branch line constitutes the resonant second mode. Coloring is as in Figure 2.

2.2 Field averaging process

The modes supported by the EC-SRR medium exhibit very fast spatial fluctuations, particularly in regions very close to the scatterers, at frequencies approaching the bandgap zones. For the homogenization to be possible, certain conditions must be fulfilled. In the present case, the wavelength of the guided mode must be much larger than the maximum length of the unit cell and, moreover, higher-order Bloch-Floquet modes must be absent. In this case, the microscopic fields, as obtained by the numerical solution, may be substituted by the average, macroscopic fields which account for the interaction of the electromagnetic field with the corresponding homogenized medium. Several integration schemes have been introduced, reflecting the proper integration areas of the microscopic fields [20], illustrated in Figure 5: Specifically, the tangentially continuous field components (e and h) are averaged through the integration of their tangential parts on the edges or the medians of the unit cell. For example, the x component of the average electric and magnetic field are averaged as

$$e_{{\text{av}}_{\text{x}}}=\frac{1}{d}{{\displaystyle \int }}_{l=-d/2}^{l=d/2}e\cdot \widehat{x}dl$$(2a)

$$h_{{\text{av}}_{\text{x}}}=\frac{1}{d}{{\displaystyle \int }}_{l=-d/2}^{l=d/2}h\cdot \widehat{x}dl,$$(2b)

where d is the length of the unit cell’s edge. Here, it is highlighted that in this work, integration was carried out on the facets of the unit cell, and it may be seen as a mean line integration of multiple line integrals on every facet. On the other hand, the corresponding normally continuous field components (periodic envelopes of the dielectric displacement d and the magnetic induction b) are averaged through the integration of their normal components on the facets of the unit cell as

$$d_{{\text{av}}_{\text{x}}}=\frac{1}{d^2}{{\displaystyle \iint }}_{S_{\text{x}}}d\cdot \widehat{x}d{S}_{\text{x}}$$(3a)

$$b_{{\text{av}}_{\text{x}}}=\frac{1}{d^2}{{\displaystyle \iint }}_{S_{\text{x}}}b\cdot \widehat{x}d{S}_{\text{x}},$$(3b)

where S_x, S_y and S_z are the surfaces of the facets normal to x, y and z axes, respectively. Here, we indicatively illustrate the obtained average field components corresponding to the SRR mode of the k_x case in Figure 6. It is obvious that all field averaging configurations almost converge to the same averaged values.

Fig. 5 Integration paths on a cubic unit cell for the tangentially continuous field components parallel to x-axis. Line integration is carried out on the four face medians (black color) and the four edges (red color). Facet integration corresponds to the four facets parallel to the x-axis: Two xz-facets and two xy-facets.

Fig. 6 Average electric (a), (b), (c) and magnetic (d), (e), (f) field components corresponding to the SRR mode of kx incidence (see Fig. 5 for the integration paths). Results obtained from integration on the four face medians is plotted in red while integration on the four unit cell’s edges is plotted in blue. Facet integration on the four facets is plotted in black: Facet #1 and facet #2 correspond to xy master and slave facets, respectively, whereas facet #3 and facet #4 indicate xz master and slave facets, respectively).

2.3 The parameter retrieval technique

A homogeneous bianisotropic medium is characterized by the following constitutive equations, that correlate the field quantities with its effective constitutive parameters

$$D={ϵ}_0{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}E+\frac{1}{c_0}{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}H$$(4a)

$$B={\mu }_0{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}H+\frac{1}{c_0}{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}E,$$(4b)

where ${{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}$ and ${{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}$ are the diagonal tensors of relative permittivity and relative permeability, respectively, while ${\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}$ is the electric-magnetic coupling tensor and ${\displaystyle \overline{{\displaystyle \overline{\zeta }}}}$ is the magnetic-electric coupling tensor. In the case of a composite periodic homogenizable medium, the average field quantities of the periodic envelopes and its effective constitutive parameters are correlated as follows [20]

$$d_{\text{av}}={ϵ}_0{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}{e}_{\text{av}}-\frac{j}{c_0}{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\kappa }}}}{h}_{\text{av}}$$(5a)

$$b_{\text{av}}={\mu }_0{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}{h}_{\text{av}}+\frac{j}{c_0}{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\kappa }}}}}^T{e}_{\text{av}},$$(5b)

which, by expanding the vectors and tensors, gives

$$\left[\begin{array}{l}{d}_{{\text{av}}_{\text{x}}}\hfill \\ d_{{\text{av}}_{\text{y}}}\hfill \\ d_{{\text{av}}_{\text{z}}}\hfill \end{array}\right]={ϵ}_0\left[\begin{array}{lll}{ϵ}_{r_{xx}}\hfill & 0\hfill & 0\hfill \\ 0\hfill & {ϵ}_{r_{yy}}\hfill & 0\hfill \\ 0\hfill & 0\hfill & {ϵ}_{r_{zz}}\hfill \end{array}\right]\left[\begin{array}{l}{e}_{{\text{av}}_{\text{x}}}\hfill \\ e_{{\text{av}}_{\text{y}}}\hfill \\ e_{{\text{av}}_{\text{z}}}\hfill \end{array}\right]$$

$$-\frac{j}{c_0}\left[\begin{array}{lll}{\kappa }_{xx}\hfill & {\kappa }_{xy}\hfill & {\kappa }_{xz}\hfill \\ {\kappa }_{yx}\hfill & {\kappa }_{yy}\hfill & {\kappa }_{yz}\hfill \\ {\kappa }_{zx}\hfill & {\kappa }_{zy}\hfill & {\kappa }_{zz}\hfill \end{array}\right]\left[\begin{array}{l}{h}_{{\text{av}}_{\text{x}}}\hfill \\ h_{{\text{av}}_{\text{y}}}\hfill \\ h_{{\text{av}}_{\text{z}}}\hfill \end{array}\right]$$(6a)

$$\left[\begin{array}{l}{b}_{{\text{av}}_{\text{x}}}\hfill \\ b_{{\text{av}}_{\text{y}}}\hfill \\ b_{{\text{av}}_{\text{z}}}\hfill \end{array}\right]={\mu }_0\left[\begin{array}{lll}{\mu }_{r_{xx}}\hfill & 0\hfill & 0\hfill \\ 0\hfill & {\mu }_{r_{yy}}\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\mu }_{r_{zz}}\hfill \end{array}\right]\left[\begin{array}{l}{h}_{{\text{av}}_{\text{x}}}\hfill \\ h_{{\text{av}}_{\text{y}}}\hfill \\ h_{{\text{av}}_{\text{z}}}\hfill \end{array}\right]$$

$$+\frac{j}{c_0}\left[\begin{array}{lll}{\kappa }_{xx}\hfill & {\kappa }_{yx}\hfill & {\kappa }_{zx}\hfill \\ {\kappa }_{xy}\hfill & {\kappa }_{yy}\hfill & {\kappa }_{zy}\hfill \\ {\kappa }_{xz}\hfill & {\kappa }_{yz}\hfill & {\kappa }_{zz}\hfill \end{array}\right]\left[\begin{array}{l}{e}_{{\text{av}}_{\text{x}}}\hfill \\ e_{{\text{av}}_{\text{y}}}\hfill \\ e_{{\text{av}}_{\text{z}}}\hfill \end{array}\right],$$(6b)

where ${\displaystyle \overline{{\displaystyle \stackrel{‾‾}{\kappa }}}}$ is the tensor of magnetoelectric coupling. It is noted that an assumption of Lorentz-reciprocal medium holds, i.e. the composite medium does not include any materials with nonreciprocal properties.

As discussed in Section 2.2, the dispersion diagrams of the EC-SRR medium consist of two supported modes for each of the three incidence cases. For each of these modes there exist six constitutive equations which correlate the average field quantities with its effective constitutive parameters. Consequently, a linear system of the form Ax = B of 36 equations and 15 unknowns is formed. The matrix A consists of the tangentially continuous average field quantities, the vector B contains the corresponding normally continuous average field components, while the unknown effective constitutive parameters are contained in vector x. It is obvious that the proposed technique is capable of returning all fifteen effective electromagnetic parameters, and that it is thus suitable to properly characterize the described composite medium.

2.4 Retrieved parameters for the EC-SRR medium

The overdetermined system of equations Ax = B is solved numerically with an accuracy of 10⁻¹³ for every simulated frequency, utilizing the mean values of the averaged field components on the unit cell’s facets for both tangentially continuous and normally continuous fields. The effective electromagnetic parameters of the investigated EC-SRR medium are shown in Figures 7–10. A resonant behavior is exhibited for the ${\mu }_{r_{zz}}$, ${ϵ}_{r_{yy}}$ and κ_yz effective parameters at the bandgap frequency range (5.2–5.9 GHz), confirming the expected behavior of magnetic dipole formation at microscopic level, normal to the scatterers’ plane. Moreover, electric dipole existence in y-direction and magnetoelectric coupling at directions yz are also obvious, owing to the non-symmetry of the resonators. All other responses, electric, magnetic or magnetoelectric, in any other direction are computationally zero with magnitudes less than 10⁻³.

Fig. 7 Effective material parameters of the investigated EC-SRR medium. First row of subfigures illustrate the real parts of permittivity (${ϵ}_{r_{\text{xx}}}$, ${ϵ}_{r_{\text{yy}}}$ and ${ϵ}_{r_{\text{zz}}}$) and the second row depicts the real parts of permeability (${\mu }_{r_{\text{xx}}}$, ${\mu }_{r_{\text{yy}}}$ and ${\mu }_{r_{\text{zz}}}$). ${\mu }_{r_{\text{zz}}}$ and ${ϵ}_{r_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies.

Fig. 8 Effective material parameters of the investigated EC-SRR medium. First row of subfigures illustrate the imaginary parts of permittivity (${ϵ}_{r_{\text{xx}}}$, ${ϵ}_{r_{\text{yy}}}$ and ${ϵ}_{r_{\text{zz}}}$) and the second row depicts the imaginary parts of permeability (${\mu }_{r_{\text{xx}}}$, ${\mu }_{r_{\text{yy}}}$ and ${\mu }_{r_{\text{zz}}}$). ${\mu }_{r_{\text{zz}}}$ and ${ϵ}_{r_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies.

Fig. 9 Effective material parameters of the investigated EC-SRR medium. The real part of the full tensor of the effective magnetoelectric coupling ${\displaystyle \overline{{\displaystyle \stackrel{‾‾}{\kappa }}}}$ is illustrated. κyz term exhibits resonant behavior around bandgap frequency zone, whereas the remaining parameters are computationally zero.

Fig. 10 Effective material parameters of the investigated EC-SRR medium. The imaginary part of the full tensor of the effective magnetoelectric coupling ${\displaystyle \overline{{\displaystyle \stackrel{‾‾}{\kappa }}}}$ is illustrated. κyz term exhibits resonant behavior around bandgap frequency zone, whereas the remaining parameters are computationally zero.

3 Reconstruction of the dispersion diagrams

3.1 Dispersion relations for a bianisotropic homogeneous medium − analytical approach

Here we present the theoretical framework utilized throughout this work for the acquisition of the analytical relations which describe the propagation of electromagnetic waves inside an arbitrarily bianisotropic homogeneous medium. Based on the work of Graglia et al. [22], we consider a plane electromagnetic wave of the type exp$\left(j\omega t-j{k}_0{\displaystyle \overline{k}}\cdot r\right)$ − where ω is the angular frequency, t is the time, k_o is the free-space wavenumber, ${\displaystyle \overline{k}}=\left({{\displaystyle \overline{k}}}_x,{{\displaystyle \overline{k}}}_y,{{\displaystyle \overline{k}}}_z\right)$ is the normalized vector wavenumber and r is the position vector − propagating inside a medium with arbitrary relative permittivity ${{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}$, relative permeability ${{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}$, electric-magnetic coupling tensor ${\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}$ and magnetic-electric coupling tensor ${\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}$ with constitutive relations given by equations (4). By introducing the singular skew-symmetric matrix

$${\displaystyle \overline{K}}=\left[\begin{array}{lll}0\hfill & -{{\displaystyle \overline{k}}}_{\text{z}}\hfill & {{\displaystyle \overline{k}}}_{\text{y}}\hfill \\ {{\displaystyle \overline{k}}}_{\text{z}}\hfill & 0\hfill & -{{\displaystyle \overline{k}}}_{\text{x}}\hfill \\ -{{\displaystyle \overline{k}}}_{\text{y}}\hfill & {{\displaystyle \overline{k}}}_{\text{x}}\hfill & 0\hfill \end{array}\right],$$(7)

so that ${\displaystyle \overline{K}}\cdot ={\displaystyle \overline{k}}\times$, we may rewrite Maxwell’s curl equations in a source-free bianisotropic medium as

$$\left[\begin{array}{cc}\hfill {{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}\hfill & {\displaystyle \overline{K}}+{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}\hfill \\ {\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}-{\displaystyle \overline{K}}\hfill & \hfill {{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}\hfill \end{array}\right]\left[\begin{array}{l}E\hfill \\ H\hfill \end{array}\right]=\left[\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right],$$(8)

which lead to the following wave equations:

$$\left[\left({\displaystyle \overline{K}}+{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}\right){{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}^{-1}\left({\displaystyle \overline{K}}-{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}\right)+{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}\right]E=0$$(9a)

$$\left[\left({\displaystyle \overline{K}}+{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}\right){{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}^{-1}\left({\displaystyle \overline{K}}-{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}\right)+{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}\right]H=0.$$(9b)

The conditions under which, equations (9) have non-trivial solutions are

$$\text{det}\left[\left({\displaystyle \overline{K}}+{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}\right){{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}^{-1}\left({\displaystyle \overline{K}}-{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}\right)+{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}\right]=0$$(10a)

$$\text{det}\left[\left({\displaystyle \overline{K}}-{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\zeta }}}}\right){{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{ϵ}}}}}_{\text{r}}^{-1}\left({\displaystyle \overline{K}}+{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\xi }}}}\right)+{{\displaystyle \overline{{\displaystyle \stackrel{¯¯}{\mu }}}}}_{\text{r}}\right]=0\text{.}$$(10b)

By solving equations (10) for any component of the wavevector (${{\displaystyle \overline{k}}}_{\text{x}}$, ${{\displaystyle \overline{k}}}_{\text{y}}$ or ${{\displaystyle \overline{k}}}_{\text{z}}$) via Matlab Symbolic Math Toolbox^TM, we obtain the values of the complex propagation constants for a medium of arbitrary bianisotropy. In [24], the validity of these equations was tested towards our alternative bianisotropic FEM eigenvalue schemes proposed there. Theoretical values extracted from equations (10) coincide with the simulated ones: The computational error is at the order of 10⁻⁹.

3.2 Comparison of original and reconstructed dispersion diagrams

We now combine the analytical relations derived in Section 3.1 and the values of the effective material parameters determined in Section 2.4 to reconstruct the dispersion diagrams of the composite EC-SRR periodic medium. In Figure 11 the dispersion diagrams for k_x incidence are shown. An excellent agreement is found between the original and the reconstructed dispersion diagrams, this being so for both the lightline and the SRR mode, in their real and imaginary parts. Exactly the same conclusion is drawn for the case of k_y incidence, see Figure 12. There is an excellent match between the original and the reconstructed values for the real and the imaginary parts, as well. Finally, the same behavior is exhibited for the normal incidence case, shown in Figure 13 for both modes (in this case, the straight line mode and the SRR mode). It is noted that the exact values of the fifteen effective parameters returned from the parameter retrieval technique were utilized, i.e., no value was set to zero a priori, which means that no simplification of the analytical relations was done, but they were used as is. Furthermore, it is very important to point out that the signs of the reconstructed values totally match the signs of the original ones, as they do not exhibit any sign change at all frequencies, including the frequencies that lay inside the bandgaps. The last observation is of great importance for the assessment of the accuracy in parameter retrieval and may serve as a criterion for any parameter retrieval technique in the literature

In Figure 14, we illustrate the percentage difference between the original and the reconstructed real parts of the complex propagation constants for the lightline modes of all three incidences (straight line mode for normal incidence). It is clear that none of the results exceeds 1% at any frequency. More specifically, the aforementioned error remains at very low values up to 6 GHz (less than 0.07% in all cases) and increases with the increase of frequency due to the lower accuracy of the FEM simulation at higher frequencies. To further support our findings, Figure 15 depicts the same difference (in %) for the EC-SRR modes. In this case, the corresponding errors are found not to exceed 4% outside of the bandgap, while larger errors occur close to the resonant frequencies, as the returned propagation constant values exhibit a very large variance at a very short frequency range. We observe that the discrepancies between the original and the reconstructed dispersion curves are at the level of numerical accuracy: Denser mesh simulations (h-refinement) and/ or utilization of higher order elements (p-refinement) will obviously result to more accurate returned eigenvectors (field distributions), leading to more accurate integrations for the averaged fields [27]. Consequently, the minimization of the difference between the original and reconstructed dispersion diagrams is feasible, particularly for higher frequencies and regions nearby bandgaps.

Fig. 11 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for kx incidence. Real and imaginary part of the complex propagation constant γ = β − jα.

Fig. 12 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for ky incidence. Real and imaginary part of the complex propagation constant γ = β − jα.

Fig. 13 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for kz incidence. Real and imaginary part of the complex propagation constant γ = β − jα.

Fig. 14 Percentage difference between the original and reconstructed real parts of the complex propagation constants for the lightline modes of all three incidences.

3.3 Rigorous explanation for the dispersion curve slopes for in-plane and normal incidence modes

The first mode at normal incidence (i.e. the one that its dispersion curve follows a straight line) is not characterized as a lightline mode, due to its alternated characteristic field distributions, which exhibit a rather complex type of wave that illustrates a specific interaction with the resonators [23]. It is definitely not a TEM wave which ignores the medium (as in the case of the in-plane incidence lightline modes). A typical field distribution is illustrated in Figure 16 in terms of the normal electric and a tangential magnetic component at 4 GHz for the investigated EC-SRR. This mode may be excited by a wave polarized as (E_x, H_y, k_z) at the aforementioned frequency.

By a very good approximation, as easily concluded from the non-zero values of its homogenized parameters, the EC-SRR composite medium may be characterized as an omega medium with its effective parameters ${ϵ}_{r_{yy}}$, ${\mu }_{r_{zz}}$ and κ_yz being non-zero and resonant; its ${ϵ}_{r_{xx}}$ parameter almost constant, with values between approximately 1.5 and 1.6; ${ϵ}_{r_{zz}}$, ${\mu }_{r_{xx}}$ and ${\mu }_{r_{yy}}$ equal to 1 and the rest of the magnetoelectric tensor elements equal to zero at the simulated frequency window [28,29]. At this point, we write the analytical relation for the refractive index that “sees” the electromagnetic wave at normal incidence, when polarized as (E_x, H_y, k_z). Solving Maxwell’s curl equations for a TEM wave of this kind yields $n_{\text{eff},z1}=\sqrt{{ϵ}_{r_{xx}}{\mu }_{r_{yy}}}$. Subsequently, the corresponding refractive index for the k_x incidence of the lightline wave (E_z, − H_y, k_x) is easily calculated as $n_{\text{eff},x}=\sqrt{{ϵ}_{r_{zz}}{\mu }_{r_{yy}}}$. We now plot these lines in Figure 17 for the sake of comparison. It is obvious that the k_z incidence dispersion curve has a different slope than the k_x incidence line, as the electromagnetic wave at normal incidence sees a larger refractive index − and a larger propagation constant β − than at in-plane incidence (the same slope is exhibited at k_y incidence, as well). In this way, it is rigorously shown that the difference between these modes is due to the different refractive indices this composite medium exhibits at these incidences, under the aforementioned polarization of incident waves.

Fig. 15 Percentage difference between the original and reconstructed real parts of the complex propagation constants for the SRR modes of all three incidences.

Fig. 16 Field distributions for the mode corresponding to the straight line of the dispersion diagram for kz incidence of the EC-SRR medium at 4 GHz. Ez (V/m) component (left) and Hx (A/m) component (right). There is a sign change on the distributions between upper surface (upper row) and lower surface (lower row), indicating the discontinuity of the components due to the presence of the resonators.

Fig. 17 Dispersion diagrams corresponding to kx (in-plane) incidence lightline mode and kz normal incidence straight line mode. Their values coincide with the ones derived by the analytical relations for wave propagation in an omega medium.

4 Conclusion

We have proposed a technique for the reconstruction of the dispersion diagrams for a general bianisotropic composite medium. We initially utilize our formerly proposed parameter retrieval technique, which is based on eigenmode analysis and field averaging, to extract the effective parameters of an arbitrarily bianisotropic composite medium. Subsequently, we derive the equations for the wave propagation inside a general homogeneous bianisotropic medium. These are used with the retrieved material parameters as input to determine the complex propagation constant of the composite medium, thereby reconstructing the dispersion diagrams for all three orthogonal incidences. The agreement between the original and the reconstructed curves is excellent, proving in this way the accuracy of our parameter retrieval technique, whereas paving the way for other numerical techniques for a means to confirm their own accuracy. Finally, utilization of our proposed technique leads to the explanation of the difference between the slopes of the in-plane lightline modes and the straight line mode present at normal wave incidence.

Funding

This project has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101064186.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contribution statement

Michalis Nitas: Conceptualization (equal); Methodology (equal); Writing ? original draft (equal). Maria Kafesaki: Supervision (supporting). Samel Arslanagic: Supervision (lead).

References

All Figures

	Fig. 1 (a) Unit cell of the EC-SRR composite medium. The resonators are arranged on a cubic lattice (the edge of the cube is d = dx = dy = dz = 8.8 mm). (b) Dimensions of the EC-SRR (in mm): The metallic rings are made of copper and have a thickness of 17 μm.
In the text

Fig. 2 Dispersion diagram of the EC-SRR composite medium for kx incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). Propagating SRR modes are plotted in blue, evanescent SRR modes with large propagation constant in green and evanescent SRR modes with small propagation constant in red. A frequency bandgap is present between 5.3 and 6 GHz. Lightline modes are illustrated in black.

In the text

	Fig. 3 Dispersion diagram of the EC-SRR composite medium for ky incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). A frequency bandgap is observed between 5.3 and 6 GHz. Coloring is as in Figure 2.
In the text

	Fig. 4 Dispersion diagram of the EC-SRR composite medium for kz (normal) incidence: propagation constant β (rad/m) and attenuation constant α (Np/m). A frequency bandgap is observed between 5.7 and 6 GHz. The modes forming the straight line constitute the first supported mode whereas the two-branch line constitutes the resonant second mode. Coloring is as in Figure 2.
In the text

	Fig. 5 Integration paths on a cubic unit cell for the tangentially continuous field components parallel to x-axis. Line integration is carried out on the four face medians (black color) and the four edges (red color). Facet integration corresponds to the four facets parallel to the x-axis: Two xz-facets and two xy-facets.
In the text

Fig. 6 Average electric (a), (b), (c) and magnetic (d), (e), (f) field components corresponding to the SRR mode of kx incidence (see Fig. 5 for the integration paths). Results obtained from integration on the four face medians is plotted in red while integration on the four unit cell’s edges is plotted in blue. Facet integration on the four facets is plotted in black: Facet #1 and facet #2 correspond to xy master and slave facets, respectively, whereas facet #3 and facet #4 indicate xz master and slave facets, respectively).

In the text

Fig. 7 Effective material parameters of the investigated EC-SRR medium. First row of subfigures illustrate the real parts of permittivity (${ϵ}_{r_{\text{xx}}}$, ${ϵ}_{r_{\text{yy}}}$ and ${ϵ}_{r_{\text{zz}}}$) and the second row depicts the real parts of permeability (${\mu }_{r_{\text{xx}}}$, ${\mu }_{r_{\text{yy}}}$ and ${\mu }_{r_{\text{zz}}}$). ${\mu }_{r_{\text{zz}}}$ and ${ϵ}_{r_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies.

In the text

Fig. 8 Effective material parameters of the investigated EC-SRR medium. First row of subfigures illustrate the imaginary parts of permittivity (${ϵ}_{r_{\text{xx}}}$, ${ϵ}_{r_{\text{yy}}}$ and ${ϵ}_{r_{\text{zz}}}$) and the second row depicts the imaginary parts of permeability (${\mu }_{r_{\text{xx}}}$, ${\mu }_{r_{\text{yy}}}$ and ${\mu }_{r_{\text{zz}}}$). ${\mu }_{r_{\text{zz}}}$ and ${ϵ}_{r_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies.

In the text

	Fig. 9 Effective material parameters of the investigated EC-SRR medium. The real part of the full tensor of the effective magnetoelectric coupling ${\displaystyle \overline{{\displaystyle \stackrel{‾‾}{\kappa }}}}$ is illustrated. κyz term exhibits resonant behavior around bandgap frequency zone, whereas the remaining parameters are computationally zero.
In the text

	Fig. 10 Effective material parameters of the investigated EC-SRR medium. The imaginary part of the full tensor of the effective magnetoelectric coupling ${\displaystyle \overline{{\displaystyle \stackrel{‾‾}{\kappa }}}}$ is illustrated. κyz term exhibits resonant behavior around bandgap frequency zone, whereas the remaining parameters are computationally zero.
In the text

	Fig. 11 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for kx incidence. Real and imaginary part of the complex propagation constant γ = β − jα.
In the text

	Fig. 12 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for ky incidence. Real and imaginary part of the complex propagation constant γ = β − jα.
In the text

	Fig. 13 Original and reconstructed dispersion diagrams of the investigated EC-SRR medium for kz incidence. Real and imaginary part of the complex propagation constant γ = β − jα.
In the text

	Fig. 14 Percentage difference between the original and reconstructed real parts of the complex propagation constants for the lightline modes of all three incidences.
In the text

	Fig. 15 Percentage difference between the original and reconstructed real parts of the complex propagation constants for the SRR modes of all three incidences.
In the text

Fig. 16 Field distributions for the mode corresponding to the straight line of the dispersion diagram for kz incidence of the EC-SRR medium at 4 GHz. Ez (V/m) component (left) and Hx (A/m) component (right). There is a sign change on the distributions between upper surface (upper row) and lower surface (lower row), indicating the discontinuity of the components due to the presence of the resonators.

In the text

	Fig. 17 Dispersion diagrams corresponding to kx (in-plane) incidence lightline mode and kz normal incidence straight line mode. Their values coincide with the ones derived by the analytical relations for wave propagation in an omega medium.
In the text