Issue 
EPJ Appl. Metamat.
Volume 11, 2024
Special Issue on ‘Metamaterials for Novel Wave Phenomena: Theory, Design and Application in Microwaves’, edited by Sander Mann and Stefano Vellucci



Article Number  10  
Number of page(s)  9  
DOI  https://doi.org/10.1051/epjam/2024008  
Published online  29 April 2024 
https://doi.org/10.1051/epjam/2024008
Research Article
Dispersion diagram reconstruction of effectively bianisotropic composite periodic media
^{1}
Department of Space Research and Technology, Technical University of Denmark, Bld. 348, Ørsteds Plads, 2800 Kgs. Lyngby, Denmark
^{2}
Institute of Electronic Structure and Laser, Foundation for Research and TechnologyHellas (FORTHIESL), 70013 Heraklion, Crete, Greece
^{3}
Department of Materials Science and Technology, University of Crete, Heraklion, Crete 70013, Greece
^{*} email: minit@dtu.dk
Received:
20
October
2023
Accepted:
27
February
2024
Published online: 29 April 2024
A dispersion diagram reconstruction technique is proposed for arbitrarily bianisotropic composite periodic media, which utilizes a previously introduced parameter retrieval technique based on eigenvalue analysis and field averaging. We initially retrieve the effective electromagnetic parameters of a composite periodic medium consisting of EdgeCoupled SplitRing Resonators (ECSRRs) via this homogenization technique using alternative integration approaches for the averaging of the field components. Subsequently, we derive the analytical framework for the wave propagation in a homogeneous medium of arbitrary bianisotropy and extract the appropriate equations which solve for the complex propagation constants. We then involve the retrieved effective parameters in these equations and reconstruct the dispersion diagrams for all three orthogonal directions, thereby spanning the whole irreducible Brillouin zone. An excellent agreement is found between the original dispersion diagrams and the reconstructed ones; a result which further validates the utilized parameter retrieval technique. The reconstruction technique moreover allows one to interpret the slope differences observed in the dispersion diagrams for inplane and normal incidence modes of the examined composite medium. It may also be used as a tool for the confirmation of the accuracy of other formerly proposed homogenization techniques existing in the literature.
Key words: Bianisotropy / dispersion diagram reconstruction / fieldflux FEM formulations / homogenization / metamaterials
© M. Nitas et al., Published by EDP Sciences, 2024
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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, semianalytical 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 Sparameter 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 NicolsonRossWeir (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 ClausiusMossotti equation and assume that the scatterers may be properly modeled via dipole approximation [14–18]. Yet an interesting homogenization technique based on the socalled TMatrix 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 fieldflux FEM formulation [23–25], we first illustrate the dispersion diagrams for the three orthogonal incidences for a bianisotropic composite medium inhere composed of EdgeCoupledSplitRing Resonators (ECSRRs). 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 ECSRR 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 ECSRR 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 inplane 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 fieldflux 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 BlochFloquet 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 BlochFloquet’s theorem
where k = k $\widehat{\text{k}}$ is the BlochFloquet 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 ECSRRs, 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 airhost and are further made of copper (conductivity σ = 5.88 × 10^{7} 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 inplane cases (k_{x} and k_{y}) a lightline mode and an SRR mode appears, whereas in the normal incidence case (k_{z}) a straightline 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^{−9}.
Fig. 1 (a) Unit cell of the ECSRR composite medium. The resonators are arranged on a cubic lattice (the edge of the cube is d = d_{x} = d_{y} = d_{z} = 8.8 mm). (b) Dimensions of the ECSRR (in mm): The metallic rings are made of copper and have a thickness of 17 μm. 
Fig. 2 Dispersion diagram of the ECSRR composite medium for k_{x} 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 ECSRR composite medium for k_{y} 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 ECSRR composite medium for k_{z} (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 twobranch line constitutes the resonant second mode. Coloring is as in Figure 2. 
2.2 Field averaging process
The modes supported by the ECSRR 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, higherorder BlochFloquet 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 xaxis. 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 xaxis: Two xzfacets and two xyfacets. 
Fig. 6 Average electric (a), (b), (c) and magnetic (d), (e), (f) field components corresponding to the SRR mode of k_{x} 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={\u03f5}_{0}{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}E+\frac{1}{{c}_{0}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}H$$(4a)
$$B={\mu}_{0}{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}H+\frac{1}{{c}_{0}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}E,$$(4b)
where ${{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}$ and ${{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}$ are the diagonal tensors of relative permittivity and relative permeability, respectively, while $\stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}$ is the electricmagnetic coupling tensor and $\stackrel{\u203e}{{\displaystyle \overline{\zeta}}}$ is the magneticelectric 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}}={\u03f5}_{0}{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}{e}_{\text{av}}\frac{j}{{c}_{0}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\kappa}}}}{h}_{\text{av}}$$(5a)
$${b}_{\text{av}}={\mu}_{0}{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}{h}_{\text{av}}+\frac{j}{{c}_{0}}{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\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]={\u03f5}_{0}\left[\begin{array}{lll}{\u03f5}_{{r}_{xx}}\hfill & 0\hfill & 0\hfill \\ 0\hfill & {\u03f5}_{{r}_{yy}}\hfill & 0\hfill \\ 0\hfill & 0\hfill & {\u03f5}_{{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 $\stackrel{\u203e}{{\displaystyle \stackrel{\u203e\u203e}{\kappa}}}$ is the tensor of magnetoelectric coupling. It is noted that an assumption of Lorentzreciprocal 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 ECSRR 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 ECSRR medium
The overdetermined system of equations Ax = B is solved numerically with an accuracy of 10^{−13} 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 ECSRR medium are shown in Figures 7–10. A resonant behavior is exhibited for the ${\mu}_{{r}_{zz}}$, ${\u03f5}_{{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 ydirection and magnetoelectric coupling at directions yz are also obvious, owing to the nonsymmetry of the resonators. All other responses, electric, magnetic or magnetoelectric, in any other direction are computationally zero with magnitudes less than 10^{−3}.
Fig. 7 Effective material parameters of the investigated ECSRR medium. First row of subfigures illustrate the real parts of permittivity (${\u03f5}_{{r}_{\text{xx}}}$, ${\u03f5}_{{r}_{\text{yy}}}$ and ${\u03f5}_{{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 ${\u03f5}_{{r}_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies. 
Fig. 8 Effective material parameters of the investigated ECSRR medium. First row of subfigures illustrate the imaginary parts of permittivity (${\u03f5}_{{r}_{\text{xx}}}$, ${\u03f5}_{{r}_{\text{yy}}}$ and ${\u03f5}_{{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 ${\u03f5}_{{r}_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies. 
Fig. 9 Effective material parameters of the investigated ECSRR medium. The real part of the full tensor of the effective magnetoelectric coupling $\stackrel{\u203e}{{\displaystyle \stackrel{\u203e\u203e}{\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 ECSRR medium. The imaginary part of the full tensor of the effective magnetoelectric coupling $\stackrel{\u203e}{{\displaystyle \stackrel{\u203e\u203e}{\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 tj{k}_{0}{\displaystyle \overline{k}}\cdot r\right)$ − where ω is the angular frequency, t is the time, k_{o} is the freespace wavenumber, $\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 \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}$, relative permeability ${{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}$, electricmagnetic coupling tensor $\stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}$ and magneticelectric coupling tensor $\stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}$ with constitutive relations given by equations (4). By introducing the singular skewsymmetric matrix
$$\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 $\overline{K}}\cdot ={\displaystyle \overline{k}}\times $, we may rewrite Maxwell’s curl equations in a sourcefree bianisotropic medium as
$$\left[\begin{array}{cc}\hfill {{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}\hfill & {\displaystyle \overline{K}}+{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}\hfill \\ {\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}{\displaystyle \overline{K}}\hfill & \hfill {{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\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 \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}\right){{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}^{1}\left({\displaystyle \overline{K}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}\right)+{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}\right]E=0$$(9a)
$$\left[\left({\displaystyle \overline{K}}+{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}\right){{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}^{1}\left({\displaystyle \overline{K}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}\right)+{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}\right]H=0.$$(9b)
The conditions under which, equations (9) have nontrivial solutions are
$$\text{det}\left[\left({\displaystyle \overline{K}}+{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}\right){{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\mu}}}}}_{\text{r}}^{1}\left({\displaystyle \overline{K}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}\right)+{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}\right]=0$$(10a)
$$\text{det}\left[\left({\displaystyle \overline{K}}{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\zeta}}}}\right){{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\u03f5}}}}}_{\text{r}}^{1}\left({\displaystyle \overline{K}}+{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\xi}}}}\right)+{{\displaystyle \stackrel{\u203e}{{\displaystyle \stackrel{\xaf\xaf}{\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^{−9}.
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 ECSRR 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 ECSRR 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 (hrefinement) and/ or utilization of higher order elements (prefinement) 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 ECSRR medium for k_{x} incidence. Real and imaginary part of the complex propagation constant γ = β − jα. 
Fig. 12 Original and reconstructed dispersion diagrams of the investigated ECSRR medium for k_{y} incidence. Real and imaginary part of the complex propagation constant γ = β − jα. 
Fig. 13 Original and reconstructed dispersion diagrams of the investigated ECSRR medium for k_{z} 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 inplane 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 inplane 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 ECSRR. 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 nonzero values of its homogenized parameters, the ECSRR composite medium may be characterized as an omega medium with its effective parameters ${\u03f5}_{{r}_{yy}}$, ${\mu}_{{r}_{zz}}$ and κ_{yz} being nonzero and resonant; its ${\u03f5}_{{r}_{xx}}$ parameter almost constant, with values between approximately 1.5 and 1.6; ${\u03f5}_{{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{{\u03f5}_{{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{{\u03f5}_{{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 inplane 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 k_{z} incidence of the ECSRR medium at 4 GHz. E_{z} (V/m) component (left) and H_{x} (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 k_{x} (inplane) incidence lightline mode and k_{z} 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 inplane 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łodowskaCurie 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
 D. Smith, S. Schultz, P. Markoš, C. Soukoulis, Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients, Phys. Rev. B 65, 195104 (2002) [NASA ADS] [CrossRef] [Google Scholar]
 D. Smith, D. Vier, T. Koschny, C. Soukoulis, Electromagnetic parameter retrieval from inhomogeneous metamaterials, Phys. Rev. E 71, 036617 (2005) [CrossRef] [Google Scholar]
 B.W. Weir, Automatic measurement of complex dielectric constant and permeability at microwave frequencies, Proc. IEEE 62, 33 (1974) [CrossRef] [Google Scholar]
 A.M. Nicolson, G.F. Ross, Measurement of the intrinsic properties of materials by timedomain techniques, IEEE Trans. Instrum. Measur. 19, 377 (1970) [CrossRef] [Google Scholar]
 D. Kwon, D. Werner, A. Kildishev, V. Shalaev, Material parameter retrieval procedure for general biisotropic metamaterials and its application to optical chiral negativeindex metamaterial design, Opt. Express 16, 11822 (2008) [CrossRef] [Google Scholar]
 R. Zhao, T. Koschny, C. Soukoulis, Chiral metamaterials: retrieval of the effective parameters with and without substrate, Opt. Express 18, 14553 (2010) [CrossRef] [Google Scholar]
 X. Chen, B. Wu, J. Kong, T. Grzegorczyk, Retrieval of the effective constitutive parameters of bianisotropic metamaterials, Phys. Rev. E 71, 046610 (2005) [CrossRef] [Google Scholar]
 U. Hasar, G. Buldu, Y. Kaya, G. Ozturk, Determination of effective constitutive parameters of inhomogeneous metamaterials with bianisotropy, IEEE Trans. Microwave Theory Tech. 66, 3734 (2018) [CrossRef] [Google Scholar]
 U. Hasar, J. Barroso, C. Sabah, Y. Kaya, M. Ertugrul, Stepwise technique for accurate and unique retrieval of electromagnetic properties of bianisotropic metamaterials, JOSA B 30, 1058 (2013) [CrossRef] [Google Scholar]
 U. Hasar, M. Bute, Method for retrieval of electromagnetic properties of inhomogeneous reciprocal chiral metamaterials, IEEE Trans. Antennas Propag. 68, 5714 (2020) [CrossRef] [Google Scholar]
 Y. Shi, Z. Li, K. Li, L. Li, C. Liang, A retrieval method of effective electromagnetic parameters for inhomogeneous metamaterials, IEEE Trans. Microw. Theory Tech. 65, 1160 (2017) [CrossRef] [Google Scholar]
 U. Hasar, M. Bute, Parameter retrieval of bianisotropic metamaterials without application of the passivity principle, IEEE Trans. Electromagn. Compatib. 63, 951 (2020) [Google Scholar]
 A. Scher, E. Kuester, Extracting the bulk effective parameters of a metamaterial via the scattering from a single planar array of particles, Metamaterials 3, 44 (2009) [CrossRef] [Google Scholar]
 A. Alù, Firstprinciples homogenization theory for periodic metamaterials, Phys. Rev. B 84, 075153 (2011) [CrossRef] [Google Scholar]
 X. Liu, A. Alù, Generalized retrieval method for metamaterial constitutive parameters based on a physically driven homogenization approach, Phys. Rev. B 87, 235136 (2013) [CrossRef] [Google Scholar]
 T. Karamanos, A. Dimitriadis, Polarizability matrix extraction of a bianisotropic metamaterial from the scattering parameters of normally incident plane waves, Adv. Electromagn. 1, 64 (2012) [CrossRef] [Google Scholar]
 T. Karamanos, A. Dimitriadis, N. Kantartzis, Robust technique for the polarisability matrix retrieval of bianisotropic scatterers via their reflection and transmission coefficients, IET Microw. Antennas Propag. 8, 1398 (2014) [CrossRef] [Google Scholar]
 T. Karamanos, S. Assimonis, A. Dimitriadis, N. Kantartzis, Effective parameter extraction of 3D metamaterial arrays via firstprinciples homogenization theory, Photonics Nanostruct. Fundament. Appl. 12, 291 (2014) [CrossRef] [Google Scholar]
 B. Zerulla, R. Venkitakrishnan, D. Beutel, M. Krstić, C. Holzer, C. Rockstuhl, I. Fernandez‐Corbaton, A T‐matrix based approach to homogenize artificial materials, Adv. Opt. Mater. 11, 2201564 (2023) [Google Scholar]
 M. Nitas, T.V. Yioultsis, Electromagnetic parameter retrieval technique utilizing eigenvalue analysis and field averaging, J. Appl. Phys. 131, 114902 (2022) [CrossRef] [Google Scholar]
 M. Nitas, M. Kafesaki, Arslanagić, Metasurface characterization based on eigenmode analysis and averaging of electromagnetic fields, J. Appl. Phys. 134, 12 (2023) [CrossRef] [Google Scholar]
 R.D. Graglia, P.L. Uslenghi, R.E. Zich, Dispersion relation for bianisotropic materials and its symmetry properties, IEEE Trans. Antennas Propag. 39, 83 (1991) [CrossRef] [Google Scholar]
 M. Nitas, T.V. Yioultsis, Characterization of edgecoupled broadsidecoupled and complementary splitring resonator periodic media based on numerical solutions of eigenvalue problems, IEEE Trans. Microw. Theory Tech. 69, 5259 (2021) [CrossRef] [Google Scholar]
 M. Nitas, V. Salonikios, T.V. Yioultsis, Alternative finite element eigenvalue formulations for the simulation of arbitrarily bianisotropic media, IEEE Trans. Microw. Theory Tech. 71, 570 (2023) [CrossRef] [Google Scholar]
 V. Salonikios, M. Nitas, S. Raptis, T.V. Yioultsis, Computational analysis of graphenebased periodic structures via a threedimensional fieldflux eigenmode Finite Element formulation, Prog. Electromagn. Res. 92, 157 (2020) [CrossRef] [Google Scholar]
 M. Nitas, V. Salonikios, C.S. Antonopoulos, T.V. Yioultsis, Numerical calculation of dispersion diagrams and field distributions of waves in 3D periodic splitring resonator media, IEEE Trans. Magn. 55, 1 (2019) [CrossRef] [Google Scholar]
 J.M. Jin, The Finite Element Method in Electromagnetics (John Wiley & Sons, 2015) [Google Scholar]
 T. Mackay, A. Lakhtakia, Electromagnetic Anisotropy and Bianisotropy: A Field Guide (World Scientific, 2010) [Google Scholar]
 A. Serdiukov et al., Electromagnetics of Bianisotropic Materials − Theory and Applications (Gordon and Breach Science Publishers, 2001) [Google Scholar]
Cite this article as: Michalis Nitas, Maria Kafesaki, Samel Arslanagić, Dispersion diagram reconstruction of effectively bianisotropic composite periodic media, EPJ Appl. Metamat. 11, 10 (2024)
All Figures
Fig. 1 (a) Unit cell of the ECSRR composite medium. The resonators are arranged on a cubic lattice (the edge of the cube is d = d_{x} = d_{y} = d_{z} = 8.8 mm). (b) Dimensions of the ECSRR (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 ECSRR composite medium for k_{x} 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 ECSRR composite medium for k_{y} 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 ECSRR composite medium for k_{z} (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 twobranch 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 xaxis. 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 xaxis: Two xzfacets and two xyfacets. 

In the text 
Fig. 6 Average electric (a), (b), (c) and magnetic (d), (e), (f) field components corresponding to the SRR mode of k_{x} 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 ECSRR medium. First row of subfigures illustrate the real parts of permittivity (${\u03f5}_{{r}_{\text{xx}}}$, ${\u03f5}_{{r}_{\text{yy}}}$ and ${\u03f5}_{{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 ${\u03f5}_{{r}_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies. 

In the text 
Fig. 8 Effective material parameters of the investigated ECSRR medium. First row of subfigures illustrate the imaginary parts of permittivity (${\u03f5}_{{r}_{\text{xx}}}$, ${\u03f5}_{{r}_{\text{yy}}}$ and ${\u03f5}_{{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 ${\u03f5}_{{r}_{\text{yy}}}$ exhibit resonant behavior around bandgap frequencies. 

In the text 
Fig. 9 Effective material parameters of the investigated ECSRR medium. The real part of the full tensor of the effective magnetoelectric coupling $\stackrel{\u203e}{{\displaystyle \stackrel{\u203e\u203e}{\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 ECSRR medium. The imaginary part of the full tensor of the effective magnetoelectric coupling $\stackrel{\u203e}{{\displaystyle \stackrel{\u203e\u203e}{\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 ECSRR medium for k_{x} incidence. Real and imaginary part of the complex propagation constant γ = β − jα. 

In the text 
Fig. 12 Original and reconstructed dispersion diagrams of the investigated ECSRR medium for k_{y} incidence. Real and imaginary part of the complex propagation constant γ = β − jα. 

In the text 
Fig. 13 Original and reconstructed dispersion diagrams of the investigated ECSRR medium for k_{z} 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 k_{z} incidence of the ECSRR medium at 4 GHz. E_{z} (V/m) component (left) and H_{x} (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 k_{x} (inplane) incidence lightline mode and k_{z} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.