The Bianisotropic Formulation of the Plane Wave Method from Faraday’s and Ampere’s Laws

Robert Gauthier

Department of Electronics, Carleton University, Ottawa, Canada.

DOI: 10.4236/opj.2021.118026                

Abstract

The plane wave numerical technique is recast from Ampere’s and Faraday’s laws for materials that are characterized with a bianisotropic form of the constitutive relations. The populating expressions are provided for the eigenvalue matrix system that can be directly solved for the angular frequencies and field profiles when bianisotropy is included. To demonstrate the computation process and expected state diagrams and field profiles, numerical computation examples are provided for a bianisotropic Bragg Array with central defect. It is shown that the location of the magnetoelectric tensor elements has a significant effect on the eigenstates of an equivalent isotropic (anisotropic) structure. One form of the magnetoelectric tensor (diagonal elements only) leads to the observation of merging states and the formation of exceptional points. The numerical approach presented can be implemented as an add-on to the familiar plane wave numerical technique.

Keywords

Bianisotropic, Magnetoelectric, Plane Wave Method, Exceptional Point

1. Introduction

The plane wave method (PWM) is a numerical solver often employed in physics and engineering [1] [2] [3]. In the field of photonics, PWM is commonly applied to structures such as the Bragg array, photonic crystals, micro-structured fibers… [4] [5] [6] and can provide information on the field profiles of supported states (eigenvectors) and propagation information through the band diagrams (eigenvalues). It is common practice to extend the outer boundary beyond the “unit cell” when structural defects are present such that localized states can be determined [7] [8] [9]. In most of the analysis, the PWM numerical solver is configured from either the electric or magnetic field wave equations and for structures that include only the anisotropic form of the constituent relations relating D to E and B to H. In this paper, the PWM is reformulated from Faraday’s and Ampere’s laws and includes the bianisotropic nature of media where D and B are related to both E and H [10] [11] [12]. In the next section, the theoretical framework leading to three distinct forms of the eigenmatrix system are provided. Section 3 provides structure details of a Bragg array with central defect and serves as a familiar starting point when the magnetoelectric tensors are included in section 4. The Appendices provide the expressions suitable for populating the matrices and key equations for non-degenerate perturbation theory. The six-field component form of the matrix expression can be directly solved as the angular frequency contribution of the terms related to the magnetoelectric tensors is to first order and additional constitutive terms absorbed into the eigenvalue square matrix. The numerical approach presented here for treating bianisotropy thus builds on the very familiar plane wave method and provides the researcher with a numerical tool suitable for the study of many different material and structure configurations.

2. Theoretical Considerations

The analysis of electromagnetic structures usually proceeds with the material properties displaying single field dependence. That is, the relation between electric flux density depends only on the electric field, $\stackrel{\to }{D}=\stackrel{\leftrightarrow }{{\epsilon }_D}\stackrel{\to }{E}$, and the magnetic flux density depends only on the magnetic field strength, $\stackrel{\to }{B}=\stackrel{\leftrightarrow }{{\mu }_B}\stackrel{\to }{H}$. For certain materials (chiral, gyromagnetic, metamaterial) and/or under certain environmental conditions (rotation, linear translation, gravitational field) [13] [14] [15], the electric and magnetic flux densities can show a dependence to both E and H fields. Such a material is called a bi-isotropic material when the relationship is linear and bianisotropic when the material properties display directional dependencies. The constitutive relations in the more general tensor form are:

$\stackrel{\to }{D}=\stackrel{\leftrightarrow }{{\epsilon }_D}\stackrel{\to }{E}+\frac{\stackrel{\leftrightarrow }{{\mu }_D}}{c}\stackrel{\to }{H}$ (1)

$\stackrel{\to }{B}=\stackrel{\leftrightarrow }{{\mu }_B}\stackrel{\to }{H}+\frac{\stackrel{\leftrightarrow }{{\epsilon }_B}}{c}\stackrel{\to }{E}$ (2)

The notation is such that a subscript of D relates tensor parameters for the electric flux density vector, and the subscript B relates tensor parameters for the magnetic flux density vector. A 1/c has been pulled out of the magnetoelectric tensors as is common practice when dealing with bianisotropy [16] [17]. For a computation environment that is charge and current free, time dependence of the fields taken as ${\text{e}}^{-j\omega t}$ and the constitutive relations in (1)-(2) are used, Faraday’s and Ampere’s laws can be expressed as:

$\stackrel{\to }{\nabla }\times \stackrel{\to }{E}=j\omega \left(\stackrel{\leftrightarrow }{{\mu }_B}\stackrel{\to }{H}+\frac{\stackrel{\leftrightarrow }{{\epsilon }_B}}{c}\stackrel{\to }{E}\right)$ (3)

$\stackrel{\to }{\nabla }\times \stackrel{\to }{H}=-j\omega \left(\stackrel{\leftrightarrow }{{\epsilon }_D}\stackrel{\to }{E}+\frac{\stackrel{\leftrightarrow }{{\mu }_D}}{c}\stackrel{\to }{H}\right)$ (4)

The relative permittivity and relative permeability are defined through $\stackrel{\leftrightarrow }{{\epsilon }_D}=$

${\epsilon }_o\stackrel{\leftrightarrow }{{\epsilon }_{rD}}$ and $\stackrel{\leftrightarrow }{{\mu }_B}={\mu }_o\stackrel{\leftrightarrow }{{\mu }_{rB}}$, the free space speed of light is $c=\frac{1}{\sqrt{{\epsilon }_o{\mu }_o}}$ and impedance is $Z_o=\sqrt{\frac{{\mu }_o}{{\epsilon }_o}}$. It is convenient to redefine the magnetic field strength vector using $j{Z}_o\stackrel{\to }{H}=\stackrel{\to }{\Im }$, Complex Scaled (SC) magnetic field. Using these definitions in (3)-(4), similar symbolic forms of Faraday’s and Ampere’s are obtained:

${\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\left[\stackrel{\to }{\nabla }\times -j\frac{\omega }{c}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\stackrel{\to }{E}=\frac{\omega }{c}\stackrel{\to }{\Im }$ (5)

${\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\left[\stackrel{\to }{\nabla }\times +j\frac{\omega }{c}\stackrel{\leftrightarrow }{{\mu }_D}\right]\stackrel{\to }{\Im }=\frac{\omega }{c}\stackrel{\to }{E}$ (6)

The solution of these equations is facilitated using matrix operators. The inverse relative permittivity and permeability tensors can be obtained from the relative permittivity and permeability of the structure expressed in the $\left(x,y,z\right)$ coordinate system and may be non-diagonal matrices. The magnetoelectric constitutive parameters may also be expressed in matrix form with elements that depend on the material properties and possibly local environment conditions.

$\begin{array}{l}{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}={\left[{\epsilon }_r\right]}^{-1}={\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{yx}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{zx}& {\epsilon }_{zy}& {\epsilon }_{zz}\end{array}\right]}^{-1}\\ {\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}={\left[{\mu }_r\right]}^{-1}={\left[\begin{array}{ccc}{\mu }_{xx}& {\mu }_{xy}& {\mu }_{xz}\\ {\mu }_{yx}& {\mu }_{yy}& {\mu }_{yz}\\ {\mu }_{zx}& {\mu }_{zy}& {\mu }_{zz}\end{array}\right]}^{-1}\end{array}$ (7)

$\begin{array}{l}\stackrel{\leftrightarrow }{{\epsilon }_B}=\left[{\epsilon }_B\right]=\left[\begin{array}{ccc}{\epsilon }_{Bxx}& {\epsilon }_{Bxy}& {\epsilon }_{Bxz}\\ {\epsilon }_{Byx}& {\epsilon }_{Byy}& {\epsilon }_{Byz}\\ {\epsilon }_{Bzx}& {\epsilon }_{Bzy}& {\epsilon }_{Bzz}\end{array}\right]\\ \stackrel{\leftrightarrow }{{\mu }_D}=\left[{\mu }_D\right]=\left[\begin{array}{ccc}{\mu }_{Dxx}& {\mu }_{Dxy}& {\mu }_{Dxz}\\ {\mu }_{Dyx}& {\mu }_{Dyy}& {\mu }_{Dyz}\\ {\mu }_{Dzx}& {\mu }_{Dzy}& {\mu }_{Dzz}\end{array}\right]\end{array}$ (8)

In PWM, a vector field, $\stackrel{\to }{A}$, can be expanded as a series of plane waves using reciprocal lattice vectors to generate the basis functions.

$\stackrel{\to }{A}={\displaystyle {\sum }_{qpn}\stackrel{\to }{{\kappa }_A}{\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}}$ (9)

The curl of $\stackrel{\to }{A}$ is organized into matrix block form using the standard definition of the field expansion coefficients for each component of the vector:

$\left[\stackrel{\to }{\nabla }\times \right]\left[\stackrel{\to }{{\kappa }_A}\right]=\left[\begin{array}{ccc}0& A_{\Downarrow }{}_{{}_{\nabla \times }}\left(yz\right)& A_{\Downarrow }{}_{{}_{\nabla \times }}\left(zy\right)\\ A_{\Downarrow }{}_{{}_{\nabla \times }}\left(xz\right)& 0& A_{\Downarrow }{}_{{}_{\nabla \times }}\left(zx\right)\\ A_{\Downarrow }{}_{{}_{\nabla \times }}\left(xy\right)& A_{\Downarrow }{}_{{}_{\nabla \times }}\left(yx\right)& 0\end{array}\right]\left[\begin{array}{c}{\kappa }_{A_x}\\ {\kappa }_{A_y}\\ {\kappa }_{A_z}\end{array}\right]$ (10)

The symbol, $\Downarrow$, indicates that once the derivatives are performed, the expansion coefficients of the vector field are extracted and used to build the column vector on the right side. The matrix operator form of Faradays’ and Ampere’s law is:

$\left(\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]-j\frac{\omega }{c}\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\right)\left[\stackrel{\to }{{\kappa }_E}\right]=\frac{\omega }{c}\left[\stackrel{\to }{{\kappa }_{\Im }}\right]$ (11)

$\left(\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]+j\frac{\omega }{c}\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\right)\left[\stackrel{\to }{{\kappa }_{\Im }}\right]=\frac{\omega }{c}\left[\stackrel{\to }{{\kappa }_E}\right]$ (12)

The expressions (11) and (12) have the same operator form and as such the computation engines required to generate the matrices in (11) can be utilized for (12) with only a substitution of the numerical values for the structure under examination. Expressions (11) and (12) require only the first derivative of the field component (rather than the second derivative using the wave equation) and no derivatives of the material properties (rather than the first derivative using the wave equation). The reduced level of derivatives typically results in a faster convergence when compared to utilizing either wave equation.

2.1. Six-Field Component Matrix Operator Formulation

The six-field component form is a matrix structure in which the column vector is composed from all 6 field components, three from the electric field and three from the CS magnetic field. Other examples of six field combination of Faraday’s and Ampere’s laws can be found in (XX [18] [19] [20] [21] [22] ) and are suitable for bianisotropic materials. The expression is obtained through a direct combination of Equations (11) and (12) into a square matrix.

$\left(\left[\begin{array}{cc}0& \left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\\ \left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]& 0\end{array}\right]+\frac{\omega }{c}\left[\begin{array}{cc}0& j\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\\ -j\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]& 0\end{array}\right]\right)\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]=\frac{\omega }{c}\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]$ (13)

Each element in these matrices is in fact a square array of numbers with an order equal to three times the number of basis functions utilized to series expand the field components. Equation (13) when displayed using matrix operator blocks is anti-diagonal which highlights the coupling between electric and magnetic fields. If the magnetoelectric contributions, $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$, display a frequency independence over the range of interest (as is usually done with the permittivity and permeability in a plane wave analysis), then expression (13) is first order in angular frequency. The solution of Equation (13) directly provides the angular frequencies, as opposed to the square of the angular frequencies when the wave equation is utilized. The field components for both electric and CS magnetic fields are also obtained at solution time.

When $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$ are zero, the leftmost matrix in (13) is referred to as the original matrix system $\left[E{H}_6\right]$ and is dependent only on the permittivity, permeability and geometry of the structure being examined. The matrix contribution from the magnetoelectric properties of (13) is referred to here as the Magneto

Electric matrix $\left[M{E}_6\right]$ and when the $\frac{\omega }{c}$ is included this may be considered as

an angular frequency dependent perturbation matrix being applied to an anisotropic optical structure. See later in text for a form of the perturbation matrix that is frequency independent. Expression (13) can be compactly written as:

$\left(\left[E{H}_6\right]+\frac{\omega }{c}\left[M{E}_6\right]\right)\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]=\frac{\omega }{c}\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]$ (14)

and when cast as an eigenvalue problem the expression is:

${\left[I-M{E}_6\right]}^{-1}\left[E{H}_6\right]\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]=\frac{\omega }{c}\left[\begin{array}{c}\stackrel{\to }{{\kappa }_E}\\ \stackrel{\to }{{\kappa }_{\Im }}\end{array}\right]$ (15)

Should $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$ be frequency dependent and the values are such that the matrix operator $\left[M{E}_6\right]$ is small compared to the curl operator blocks, then a perturbative approach to solving when the bianisotropy is present can be applied. The needed key expressions for a non-degenerate perturbative analysis are provided in appendix C. One numerical disadvantage in using this matrix system (13) is that the order of the matrix is twice as large as would be obtained using PWM formulated using eitherwave equations. However, (13) can be directly solved as an eigenvalue problem, while utilizing the wave equations in a bianisotropic formulation cannot and one must revert to a perturbative approach. The order of the matrix system, originating from Faraday’s and Ampere’s operator form, can be reduced to half that of (13) without loss of solution accuracy as is shown in the single field matrix formulation presented next. However, these forms would require a perturbative approach when matrix elements of the magnetoelectric terms are non-zero.

2.2. Single-Field Component Matrix Operator Formulation

The single field matrix formulation can be cast for solving for either the electric field or the CS magnetic field. The single field formulation for the electric field vector requires the use of Equation (12) and eliminating the CS magnetic field using (11). The resulting expression is:

$\begin{array}{l}(\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\nabla \times \right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]+j\frac{\omega }{c}{(\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]}_{{}_{}^{}}^{}\\ -\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right])+{\left(\frac{\omega }{c}\right)}^2\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right])\left[\stackrel{\to }{{\kappa }_E}\right]={\left(\frac{\omega }{c}\right)}^2\left[\stackrel{\to }{{\kappa }_E}\right]\end{array}$ (16)

This matrix expression is a quadratic of the angular frequency and possibly of higher order if the bianisotropy is frequency dependent. Once the single field matrix system is solved, the missing CS magnetic field components can be obtained using (11). If the bianisotropic contribution is zero, then the matrix system

reduces to $\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[\stackrel{\to }{{\kappa }_E}\right]={\left(\frac{\omega }{c}\right)}^2\left[\stackrel{\to }{{\kappa }_E}\right]$ which can be solved as an

eigenvalue problem. The isotropic/anisotropic matrix operator when solving for the electric field can be defined as:

$\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]=\left[H{E}_3\right]$ (17)

The eigenvalues returned are the square of the angular frequency, and the eigenvectors are the electric field components of the states. The order of the matrix system is half that of matrix expression (13). The remaining terms in (16), not associated with $\left[H{E}_3\right]$ result from the magnetoelectric properties and can be collected in a matrix operator and considered as an angular frequency dependent perturbation to the isotropic/anisotropic structure.

$\begin{array}{c}\left[P_E\left(\omega ,{\omega }^2\right)\right]=j\frac{\omega }{c}\left(\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]-\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\right)\\ +{\left(\frac{\omega }{c}\right)}^2\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\end{array}$ (18)

The bianisotropic matrix system to solve is:

$\left(\left[H{E}_3\right]+\left[P_E\left(\omega ,{\omega }^2\right)\right]\right)\left[\stackrel{\to }{{\kappa }_E}\right]={\left(\frac{\omega }{c}\right)}^2\left[\stackrel{\to }{{\kappa }_E}\right]$ (19)

This expression may be difficult to solve directly due to the different powers of the angular frequency. Iterative or perturbative approaches may be more desirable if the magnetoelectric contribution to the system matrix is small. A similar form to (19) can be obtained for the CS magnetic field. The matrix expression is:

$\left(\left[E{H}_3\right]+\left[P_{\Im }\left(\omega ,{\omega }^2\right)\right]\right)\left[\stackrel{\to }{{\kappa }_{\Im }}\right]={\left(\frac{\omega }{c}\right)}^2\left[\stackrel{\to }{{\kappa }_{\Im }}\right]$ (20)

with

$\left[E{H}_3\right]=\left[\stackrel{\leftrightarrow }{{\mu }_{rB}{}^{-1}}\stackrel{\to }{\nabla }\times \right]\left[\stackrel{\leftrightarrow }{{\epsilon }_{rD}{}^{-1}}\stackrel{\to }{\nabla }\times \right]$ (21)

and

$\begin{array}{c}\left[P_{\Im }\left(\omega ,{\omega }^2\right)\right]=j\frac{\omega }{c}\left(\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]-\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\right)\\ +{\left(\frac{\omega }{c}\right)}^2\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\end{array}$ (22)

The missing electric field components would be obtained using (12). In the Appendices, the matrix generating expressions are provided for all the matrix operators involved in setting up the eigenvalue expressions for (15), (19) and (20).

3. Bragg Plane Array

The effects of incrementally increasing the magnetoelectric contribution to the constituent expressions will be examined using a Bragg array of planes with a missing central plane as the defect region, Figure 1-Left. The spatial parameters of the structure and coordinate orientation are provided in the figure and caption. The relative permittivity and permeability media are first chosen as isotropic. The relative permittivity of the dark regions has ${\epsilon }_r=10.000$ and the light regions has ${\epsilon }_r=1.000$. The relative permeability of the entire structure has ${\mu }_r=1.000$. Anisotropy is introduced in the latter part of the simulations and analysis. The structure is oriented such that it is uniform-infinite along the y and z axis. A cubic super cell with lattice constant, a, of 9 µm is chosen. The grey border in the structure image shows the placement of the unit cell in the $\left(x,y\right)$ plane. The supercell is discretized using 50 points per micron. The constituent parameter series of (A2) utilize $-400\le p_{\epsilon }\le +400,q_{\epsilon }=0,n_{\epsilon }=0$ for each element in the inverse relative permittivity tensor. The structure is first examined with no magnetoelectric contribution such that the band structure and centrally localized states can be determined. When the structure under examination is non-magnetic and uniform-infinite in y andz, as examined here, the non-zero matrix blocks of (13), for $\left[E{H}_6\right]$ are:

$\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& -k_z& k_y\\ k_z& 0& -G_p-k_x\\ -k_y& G_p+k_x& 0\end{array}\right]$ (23)

$\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& -k_z{\kappa }_{{\epsilon }_r}& k_y{\kappa }_{{\epsilon }_r}\\ k_z{\kappa }_{{\epsilon }_r}& 0& -{\kappa }_{{\epsilon }_r}{G}_p-k_x{\kappa }_{{\epsilon }_r}\\ -k_y{\kappa }_{{\epsilon }_r}& {\kappa }_{{\epsilon }_r}{G}_p+k_x{\kappa }_{{\epsilon }_r}& 0\end{array}\right]$ (24)

The field components are series represented using ( $-400\le p_f\le +400$, $q_f=0$, $n_f=0$ ). The band structure can be computed using either expression

Figure 1. Left: Bragg array shown in the (x,y) plane. Dark region has relative permittivity of 10.000, light region is air with relative permittivity of 1.000. Mediums are non-magnetic. Units of measures are microns. Grey area defines border of 9 µm × 9 µm supper cell used in the plane wave expansion. Right: Band structure obtained for the k vector directions as indicated. Each component and component pairing range from 0 to π/a. Grey bands on right indicate band gap locations, arrows indicate defect states.

(15), (19), (20) with the magnetoelectric tensor contribution set to zero. The resulting band structure is shown in Figure 1-right for k-vectors that contain components along the three-coordinate axis $\left(k_x,k_y,k_z\right)$ and along the combined directions of $\left(k_{x,y},k_{y,z},k_{x,y,z}\right)$. K vector components run from 0 to π/a. Band gaps are observed between 1.301 - 2.766, 3.560 - 4.992, 5.750 - 7.346, 8.490 - 8.762, 9.968 - 11.296 and 12.137 - 13.517. For this structure orientation and geometrical parameters, the states are degenerate and bandgap states are observed at 2.282, 3.623, 6.337, 9.035 and 12.261. The states which are TM have ( $E_z$, $H_y$, $H_x=0$ ), while the states which are TE have ( $H_z$, $E_y$, $E_x=0$ ). The $E_z$ field component magnitude for the TM polarized localized states is plotted in Figure 2.

4. Bianisotropy and the Magnetoelectric Tensor

Each element in the magnetoelectric tensors, $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$ Equation (8), must be specified at every discretization point of the supercell. The individual element values may be real, imaginary, or complex and inter-related through constraints and symmetry operations [10]. For the Bragg array structure being examined, the magnetoelectric contributions to the constituent relations are present for only the dark material regions of Figure 1-Left. The light region, air, has no magnetoelectric contribution. The spatial properties of $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$ are significantly

Figure 2. Magnitude of the E_z field component for the defect region localized states identified in Figure 1 that are TM polarized. Plots produce with zero k vector. Plot shown in the (x,y) plane (10 µm, 10 µm). Grey edging shows extent of supercell. Normalized angular frequencies of states plotted: (a) 2.282 (b) 3.623 (c) 6.337 (d) 9.035 (e) 12.261.

different from the inverse relative permittivity, ${\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}$ and inverse relative permeability, ${\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}$. The magnetoelectric permittivity, $\stackrel{\leftrightarrow }{{\epsilon }_B}$, varies between 0 (not 1 as for air) and the magnetoelectric value being considered and the magnetoelectric permeability, $\stackrel{\leftrightarrow }{{\mu }_D}$, also varies between 0 (not 1 for non-magnetic) and the magnetoelectric value being applied. The magnetoelectric properties, $\stackrel{\leftrightarrow }{{\epsilon }_B}$ and $\stackrel{\leftrightarrow }{{\mu }_D}$, are periodic in the x-axis direction. The matrix operators associated from inclusion of the magnetoelectric portion in the constituent relations are composed of a tensor multiplication of a permittivity (permeability) with a magnetoelectric tensor $\stackrel{\leftrightarrow }{{\mu }_D}$ ( $\stackrel{\leftrightarrow }{{\epsilon }_B}$ ). See expressions (11) and (12). The matrix products are obtained for every discretization grid point and these new tensors can be considered as an “effective” material which may be series expanded in the PWM formulation using:

${\epsilon }_{ab}^{■}{\mu }_{Dcd}={\displaystyle {\sum }_{p_D{q}_D{n}_D}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}{\text{e}}^{j\left(G_{p_D}x+G_{q_D}y+G_{n_D}z\right)}}$ (25)

${\mu }_{ab}^{■}{\epsilon }_{Bcd}={\displaystyle {\sum }_{p_B{q}_B{n}_B}{\kappa }_{{\mu }_{ab}{\epsilon }_{cd}}{\text{e}}^{j\left(G_{p_B}x+G_{q_B}y+G_{n_B}z\right)}}$ (26)

The symbol ( $■$ ) is utilized to represent that the tensor element utilized is obtained after the inverse of the relative permittivity or permeability tensors are computed for each grid point of the discretized structure. The generating expressions required to populate the matrix operators when the magnetoelectric properties are included in Appendix B.

The tensor elements created by the series expansion of ${\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}$ results in complex valued expansion coefficients as the values of ${\epsilon }_{rD}{}^{-1}$ are period (between 1 and 10.000) and act as multiplier to an x-axis period $\stackrel{\leftrightarrow }{{\mu }_D}$. The series expansion of the tensor elements related to ${\stackrel{\leftrightarrow }{{\mu }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}$ are also period in the x-axis direction as $\stackrel{\leftrightarrow }{{\epsilon }_B}$ varies between 0 and the applied magnetoelectric value, with ${\stackrel{\leftrightarrow }{{\mu }_{rD}}}^{-1}=1$ as a constant multiplier.

When the magnetoelectric properties are to be included in this theoretical examination, a systematic approach to the selection of the non-zero tensor elements is facilitated by examining the corresponding matrix structure for the two operators, $\left[E{H}_6\right]$ and $\left[M{E}_6\right]$. The non-zero matrix elements in $\left[E{H}_6\right]$ are governed by the components of $\stackrel{\to }{k}$ as explicitly indicated in (23) and (24). The $k_x$ component contributes to matrix block elements at positions (23, 32), the $k_y$ component contributes to matrix block locations (13, 31) and the $k_z$ component contributes to matrix block element locations (12, 21). The placement of the non-zero matrix elements when the magnetoelectric properties are present will follow those of the $\stackrel{\to }{k}$ in $\left[E{H}_6\right]$. Examination of the matrix operators that build up in $\left[E{H}_6\right]$ shows that they have no diagonal elements. The effect of having diagonal element in the magnetoelectric tensors is also examined. Representative examples of the state space obtained when the magnetoelectric properties are incrementally increased are provided for a sampling of the tensor value combinations using the Γ point of the band diagram. This point has zero $\stackrel{\to }{k}$ and the 6-field matrix operator has the block form shown in Figure 3. The operator has non-zero element blocks at (26, 35, 53, 62). Each matrix

Figure 3. Top: Block matrix operator form for [EH₆] at the Γ point for the Bragg array. Block rows and columns 1 and 4 are zero indicating zero x directed field components for the states. Matrix element blocks of interest are located at 26, 35, 53, 62 and couple electric to magnetic field components producing TE and TM polarized states. Lower: Figures show relative size of the real and imaginary contributions to the system matrix.

block has an order 801 ( $-400\le p_f\le +400$ ) determined by the series expansion space for the field components. Since rows and columns blocks 1 and 4 are padded with zeros, the states of the structure with no magnetoelectric presences will have zero $E_x$ and $H_x$ field components.

4.1. Magnetoelectric Tensor Elements at Equivalent “k_y” and “k_z” Positions

The rotational symmetry about the x-axis dictates that the results obtained with the magnetoelectric tensor elements at the “k_y” equivalent positions should be the same as those obtained when the tensor elements are at the equivalent “k_z” positions. The magnetoelectric tensor structure is shown in (27) and is chosen as anti-symmetric such that it follows the anti-symmetric nature of (23) and (24). Real, imaginary, or complex tensor elements are chosen by setting $\left(\alpha ,\beta \right)\iff \left(0\leftrightarrow 1\right)$. The strength of the magnetoelectric contribution is adjusted through the multiplier in the front of each matrix. For the computations, the tensor elements values are as indicated in the figures and the magnetoelectric properties are incremented through the multiplier from 0 to 1 in 0.005 steps. State diagrams are produced using the 6-field matrix operator form.

$\begin{array}{l}\stackrel{\leftrightarrow }{{\epsilon }_B}={\epsilon }_y\left[\begin{array}{ccc}0& 0& \alpha +j\beta \\ 0& 0& 0\\ -\alpha -j\beta & 0& 0\end{array}\right]+{\epsilon }_z\left[\begin{array}{ccc}0& \alpha +j\beta & 0\\ -\alpha -j\beta & 0& 0\\ 0& 0& 0\end{array}\right]\\ \stackrel{\leftrightarrow }{{\mu }_D}={\mu }_y\left[\begin{array}{ccc}0& 0& \alpha +j\beta \\ 0& 0& 0\\ -\alpha -j\beta & 0& 0\end{array}\right]+{\mu }_z\left[\begin{array}{ccc}0& \alpha +j\beta & 0\\ -\alpha -j\beta & 0& 0\\ 0& 0& 0\end{array}\right]\end{array}$ (27)

Figure 4 shows the state diagrams when the magnetoelectric contribution is incrementally increasing for real (left), imaginary (center) and complex (right) tensor elements. Interpretation of the state evolution is facilitated by examining the nature of the matrix system solved when the magnetoelectric properties are present. A perturbation matrix resulting from the presence of the magnetoelectric values, $\left[P_6\right]$, can be computed through the following expression:

$\left[P_6\right]={\left[I-M{E}_6\right]}^{-1}\left[E{H}_6\right]-\left[E{H}_6\right]$ (28)

and represents the deviation to the matrix elements from the Γ point computed in Figure 2. Non-zero matrix block elements for magnetoelectric tensor elements at the “k_y” equivalents are located at positions indicated in the 6 by 6 matrix in (29). The mixing of field components between different states is obtained through the perturbative expression (C6) and provides the column vector on the right-side of (29):

$\begin{array}{l}\left[\begin{array}{cc}\begin{array}{ccc}{E^{\prime }}_x& {E^{\prime }}_y& {E^{\prime }}_z\end{array}& \begin{array}{ccc}{H^{\prime }}_x& {H^{\prime }}_y& {H^{\prime }}_z\end{array}\end{array}\right]\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ P_{21}& 0& 0& 0& 0& P_{26}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& P_{53}& P_{54}& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\\ \begin{array}{c}{H}_x\\ H_y\\ H_z\end{array}\end{array}\right]\\ =\left[\begin{array}{c}\begin{array}{c}0\\ {E^{\prime }}_y\left(P_{21}{E}_x+P_{26}{H}_z\right)\\ 0\end{array}\\ \begin{array}{c}0\\ {H^{\prime }}_y\left(P_{53}{E}_z+P_{54}{H}_x\right)\\ 0\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}0\\ {E^{\prime }}_y\left(P_{26}{H}_z\right)\\ 0\end{array}\\ \begin{array}{c}0\\ {H^{\prime }}_y\left(P_{53}{E}_z\right)\\ 0\end{array}\end{array}\right]\end{array}$ (29)

This matrix system shows that states that are TE in nature mix, through the perturbation with other TE states and that TM states mix with other TM states. The level of mixing is governed by the blocks P₂₆ and P₅₃ as the states have no E_x

Figure 4. State diagrams for the Γ point when the magnetoelectric contribution is incrementally increased with non-zero tensor elements at the “k_y” equivalent matrix block positions. Left—Real. Center—Imaginary. Right—Complex. X-axis parameters referring to (27): $\left({\epsilon }_y,{\mu }_y\right)=\left(0:0.005:1\right)$, $\left({\epsilon }_z,{\mu }_z\right)=\left(0,0\right)$.

component. The magnetoelectric tensors produce P₂₆ and P₅₃ matrix blocks that act as though an effective real “k_y” is produced. When the matrix elements are complex, the elements P₂₆ and P₅₃ are not simply the complex conjugates as the real and imaginary contributions to the permittivity are not the same. This leads to a slightly different state evolution as the imaginary magnetoelectric properties are incremented. A deeper analysis shows that the complex nature of the tensor elements is similar to having an imaginary effective propagation constant, or equivalently, introducing an imaginary component to the relative permittivity. The complex element state diagram is produced with magnetoelectric tensor elements that have both real and imaginary contributions. The net effect is to extend the trace as the magnitude of the tensor elements exceeds those of either real or imaginary by 1.4142 ( $\sqrt{2}$ ). Due to the large number of matrix elements present and the fact that the states are degenerate, a more detailed analysis is difficult to perform and is worthy of a more focused follow up publication.

When the matrix system has non-zero magnetoelectric tensor elements at the equivalent “k_z” locations, the mixing of field components is given in the column vector on the right-side of (30). Comparison of (30) with (29) shows that although the non-zero matrix blocks locations in the two perturbation matrices are different, the net effect is to mix the same field component pairs through equivalent strength perturbation block elements ( $P_{26}=P_{62},P_{53}=P_{35}$ ).

$\left[P_6\right]=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ P_{31}& 0& 0& 0& P_{35}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& P_{62}& 0& P_{64}& 0& 0\end{array}\right]\Rightarrow \left[\begin{array}{c}\begin{array}{c}0\\ 0\\ {E^{\prime }}_z\left(P_{31}{E}_x+P_{35}{H}_y\right)\end{array}\\ \begin{array}{c}0\\ 0\\ {H^{\prime }}_z\left(P_{62}{E}_y+P_{64}{H}_x\right)\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}0\\ 0\\ {E^{\prime }}_z\left(P_{35}{H}_y\right)\end{array}\\ \begin{array}{c}0\\ 0\\ {H^{\prime }}_z\left(P_{62}{E}_y\right)\end{array}\end{array}\right]$ (30)

4.2. Magnetoelectric Tensor Elements at Equivalent “k_x” Positions

When the non-zero magnetoelectric tensor elements occupy matrix positions 23 and 32, they are located at the equivalent positions of the “k_x” contributions of the propagation constant in $\left[E{H}_{\text{6}}\right]$. The tensor form is shown in (31) with $\alpha ,\beta$ and the matrix multipliers as define in Section 4.1.

$\stackrel{\leftrightarrow }{{\epsilon }_B}={\epsilon }_x\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \alpha +j\beta \\ 0& -\alpha -j\beta & 0\end{array}\right],\stackrel{\leftrightarrow }{{\mu }_D}={\mu }_x\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \alpha +j\beta \\ 0& -\alpha -j\beta & 0\end{array}\right]$ (31)

Figure 5 shows the state diagram when the magnetoelectric contribution is incremented for real (left), imaginary (center) and complex (right) tensor elements. The similarity in the effect of introducing “k_x” equivalent terms to the magnetoelectric tensor follows that of increasing the x-axis propagation constant when tensor elements are real. For imaginary and complex tensor elements, the curvature is related to the equivalency of adding an imaginary contribution to the propagation constant (k_x effective) or making the permittivity have an

Figure 5. State diagrams for the Γ point when the magnetoelectric contribution is incrementally increased with non-zero tensor elements at the “k_x” equivalent positions. Left—Real. Center—Imaginary. Right—Complex. X-axis parameters referring to (31): $\left({\epsilon }_x,{\mu }_x\right)=\left(0:0.005:1\right)$.

imaginary contribution. The perturbation matrix and its effect on field coupling is expressed in (32) and indicates that states of similar polarization will couple while states of different polarization do not.

$\left[P_6\right]=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& P_{22}& 0& 0& 0& P_{26}\\ 0& 0& P_{33}& 0& P_{35}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& P_{53}& 0& P_{55}& 0\\ 0& P_{62}& 0& 0& 0& P_{66}\end{array}\right]\Rightarrow \left[\begin{array}{c}\begin{array}{c}0\\ {E^{\prime }}_y\left(P_{22}{E}_y+P_{26}{H}_z\right)\\ {E^{\prime }}_z\left(P_{33}{E}_z+P_{35}{H}_y\right)\end{array}\\ \begin{array}{c}0\\ {H^{\prime }}_y\left(P_{53}{E}_z+P_{55}{H}_y\right)\\ {H^{\prime }}_z\left(P_{62}{E}_y+P_{66}{H}_z\right)\end{array}\end{array}\right]$ (32)

4.3. Diagonal Magnetoelectric Tensor Elements, Isotropic Bragg Array

The diagonal magnetoelectric tensor elements may be real, imaginary, or complex and representative computation examples for each are provided. The magnetoelectric tensor matrix structure is provided in (33) with $\alpha ,\beta$ and the matrix multipliers as previously defined.

$\begin{array}{l}\stackrel{\leftrightarrow }{{\epsilon }_B}={\epsilon }_{\text{diagonal}}\left[\begin{array}{ccc}\alpha +j\beta & 0& 0\\ 0& \alpha +j\beta & 0\\ 0& 0& \alpha +j\beta \end{array}\right]\\ \stackrel{\leftrightarrow }{{\mu }_D}={\mu }_{\text{diagonal}}\left[\begin{array}{ccc}-\alpha -j\beta & 0& 0\\ 0& -\alpha -j\beta & 0\\ 0& 0& -\alpha -j\beta \end{array}\right]\end{array}$ (33)

Diagonal matrix elements in the magnetoelectric tensors provide diagonal matrix blocks in the corresponding operators and results in a perturbation matrix of the form provided in (34). The location of the non-zero blocks in the perturbation matrix indicates that states with similar and opposite polarizations will couple. The opposite polarization state coupling terms are set in bold type in (34).

$\left[P_6\right]=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& P_{23}& 0& 0& P_{26}\\ 0& P_{32}& 0& 0& P_{35}& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& P_{53}& 0& 0& P_{56}\\ 0& P_{62}& 0& 0& P_{65}& 0\end{array}\right]\Rightarrow \left[\begin{array}{c}\begin{array}{c}0\\ \left({E^{\prime }}_y{P}_{23}{E}_z+{E^{\prime }}_y{P}_{26}{H}_z\right)\\ \left({E^{\prime }}_z{P}_{32}{E}_y+{E^{\prime }}_z{P}_{35}{H}_y\right)\end{array}\\ \begin{array}{c}0\\ \left({H^{\prime }}_y{P}_{53}{E}_z+{H^{\prime }}_y{P}_{56}{H}_z\right)\\ \left({H^{\prime }}_z{P}_{62}{E}_y+{H^{\prime }}_z{P}_{65}{H}_y\right)\end{array}\end{array}\right]$ (34)

Figure 6 shows the state diagrams when the magnetoelectric contribution is incremented for real (left), imaginary (center) and complex (right) tensor elements. For real tensor elements, there is a merging of adjacent states as the magnetoelectric properties increase. Close examination shows that the strength of magnetoelectric values required to complete a two state merger is inversely proportional to angular frequency difference of these states. This is expected as the coupling strength is inversely proportional to eigenvalue difference as indicated in (C6). The center figure (imaginary) shows that the states that merged when the magnetoelectric values were real, are now diverging. The divergent nature of the states leads to the criss-cross appearance. Note that in general, the matrix blocks of $\left[P_6\right]$ are complex. Changing the magnetoelectric properties from real to imaginary results in a complex conjugate exchange of the real and imaginary parts of the numerical values making up the matrix blocks in $\left[P_6\right]$. State components that were mixing as real-real or imaginary-imaginary are now mixing as real-imaginary and imaginary-real. The right figure (complex) is produced from tensor elements that contain equal real and imaginary values. Over the range of the plot, the initial band gaps and defect states observed are maintained, and state merging has been suppressed. The nature of the plot is similar to the complex plots produced for the presence of off-diagonal tensor elements shown previously.

Figure 6. State diagrams for the Γ point when the magnetoelectric contribution is incrementally increased with non-zero tensor elements at the diagonal block locations. Left—Real. Center—Imaginary. Right—Complex. X-axis parameters referring to (33): $\left({\epsilon }_{\text{diagonal}},{\mu }_{\text{diagonal}}\right)=\left(0:0.005:1\right)$.

4.4. Diagonal Magnetoelectric Tensor Elements, Anisotropic Bragg Array

The state evolution when diagonal magnetoelectric tensor elements are present is now examined when the Bragg planes are composed of a biaxial anisotropic relative permittivity. The initially isotropic regions with ${\epsilon }_r=10.000$ is replaced with a medium having ${\epsilon }_{xx}=8.000$, ${\epsilon }_{yy}=10.000$, ${\epsilon }_{zz}=12.000$. The anisotropy will remove the degeneracy of TE and TM states but result in a more complex band structure as shown in Figure 7. Two of the band gap states, one TE and one TM, have the field profiles plotted in Figure 8. For the anisotropic structure and orientation, there are no $E_x$ or $H_x$ components.

The magnetoelectric properties are incrementally applied through expression (30) for diagonal elements and the state diagrams produced for real (left), imaginary (center) and complex (right) tensor elements are shown in Figure 9. As expected, the lifting of the degeneracy results in more elaborate state diagrams with the same general behavior observed for the isotropic Bragg array state diagrams of Figure 6. An enlargement of the state diagram in the 5.6 to 6.65 is shown in Figure 10 for real tensor elements.

The evolution of the normalized propagation constant for TE polarized band gap state initially at 6.348 and the TM polarized state it merges with at 6.612 are shown in Figure 10 top left for real and top right for imaginary contributions. The field profiles are plotted in Figure 11 for the 4 points indicated. Before the merger points, there is no imaginary contribution to the propagation constant. Up to the merger point the states couple and the coupling is responsible for the mixing of the two polarizations. See column vector in (34). TE states begin to acquire a TM contribution and vice versa. At the merge point, the two states have the same normalized angular frequency and have identical field profiles,

Figure 7. Band diagram when the structure is made biaxial anisotropic. Several complete band gaps exist. Two band gap states are, initially degenerate when the material is isotropic, are shown no longer degenerate.

Figure 8. Top pair—TE state at normalized frequency 6.348 showing absolute value of the E_y and H_z fields. Bottompair—TM state at normalized frequency 6.087 showing absolute value of the E_z and H_y fields.

Figure 9. State diagrams for the Γ point for a anisotropic Bragg array when the magnetoelectric contribution is incrementally increased with non-zero tensor elements at the diagonal block locations. Left—Real. Center—Imaginary. Right—Complex. X-axis parameters referring to (33): $\left({\epsilon }_{\text{diagonal}},{\mu }_{\text{diagonal}}\right)=\left(0:0.005:1\right)$.

trace point C in field profiles of Figure 11. They form a degenerate pair. After the merger point, opposite sign complex contributions to the normalized angular frequency are present making the point at C acquire the behavior of an exceptional point [23] [24] [25]. The traces in Figure 11 at point D show that the field profiles are still the same even with the imaginary contribution to the angular frequency.

The evolution of the normalized propagation constant for TM polarized band gap state initially at 6.087 and the TE polarized state it merges with at 5.774 are

Figure 10. Real (left) and imaginary (right) contributions to the normalized angular frequency. Top pair—Merger of the localized state at normalized angular frequency 6.348 with the state at 6.612. Bottom pair—Merger of the localized state at normalized angular frequency 6.087 with the state at 5.774. Points A to E have field profiles plotted in Figure 11 and Figure 12. X-axis parameters referring to (33): $\left({\epsilon }_{\text{diagonal}},{\mu }_{\text{diagonal}}\right)=\left(0:0.005:1\right)$.

shown in Figure 10 bottom left for real and bottom right for imaginary contributions. The field profiles for are plotted in Figure 12 for the 5 points indicated. The behavior of these states is similar to those above and an exceptional point (at different alpha location) is observed.

The effect of having imaginary and complex tensor elements on the band gap states of 9 is shown in Figure 13. The attractive nature of adjacent states observed when the tensor elements are real now have a repulsive nature and state merger is no longer observed. The upper band gaps are rapidly closed off and lower band gaps would eventually close provided the tensor elements are sufficiently increased. For complex diagonal tensor elements, previously merging states are also repulsive. A significant difference is that complex diagonal tensor elements can retain the original bandgaps.

5. Conclusion

The plane wave method was recast using Faraday’s and Ampere’s laws into a matrix operator form suitable for analysis of structures that are characterized as bianisotropic. Of the three matrix forms provided, the 6-field eigensystem can be

Figure 11. Traces for the absolute value of the E_z (top two) and E_y (lower two) field component for the states (6.348, 6.612) that merge with increasing magnetoelectric real tensor elements. Inserts show H_z (top inserts) and H_y (lower inserts) for these states. At the exceptional point (C), field profiles for the two states are the same.

directly solved for the eigenstates and angular frequencies in the presence of non-zero magnetoelectric tensor elements. The implementation of this numerical approach provides the researcher with a tool suitable for the study of numerous and different configurations were the bianisotropic properties and structure geometry can be chosen freely. Representative numerical results are presented based on the location of the magnetoelectric tensor elements for a familiar optical structure and serve to illustrate the computation process flow and expected numerical results. Features extracted from non-degenerate perturbation theory are utilized to display the nature of the field component coupling taking place. The diagonal element form of the magnetoelectric tensor is shown to lead to the existence of exceptional points in the state diagrams. Provided the bianisotropic structures for examination can be contained within a periodic super cell, the numerical approach presented is suitable for the examination of material properties such as that may be encountered when dealing with environmental effects such as the Sagnac effect, or material properties such as Chirality

Figure 12. Traces for the absolute value of the E_z (top two) and E_y (lower two) field component for the states (6.087, 5.774) that merge with increasing magnetoelectric real tensor elements. Inserts show H_z (top inserts) and H_y (lower inserts) for these states. At the exceptional point (D), field profiles for the two states are the same.

Figure 13. Segment of the state evolution diagram when the diagonal tensor elements are imaginary (left) and complex (right) for the two merging pairs of Figure 10. States are no longer observed to merge and generate exceptional points.

and metamaterials.

Acknowledgements

The author wishes to acknowledge support from NSERC.

Appendix A: Matrix Populating Expressions for Curl Containing Operators

The expressions provided are used to populate matrices constructed from terms containing the curl operator. Each of the matrix operators $\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]$ and $\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]$ can be populated using the same numerical algorithm. Working from the first of these the matrix written in expanded form is:

$\begin{array}{l}\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\stackrel{\to }{{\kappa }_{\Im }}\\ =\left[\begin{array}{ccc}{\epsilon }_{xy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xz\right)+{\epsilon }_{xz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xy\right)& {\epsilon }_{xx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yz\right)+{\epsilon }_{xz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yx\right)& {\epsilon }_{xx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zy\right)+{\epsilon }_{xy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zx\right)\\ {\epsilon }_{yy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xz\right)+{\epsilon }_{yz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xy\right)& {\epsilon }_{yx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yz\right)+{\epsilon }_{yz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yx\right)& {\epsilon }_{yx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zy\right)+{\epsilon }_{yy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zx\right)\\ {\epsilon }_{zy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xz\right)+{\epsilon }_{zz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xy\right)& {\epsilon }_{zx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yz\right)+{\epsilon }_{zz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yx\right)& {\epsilon }_{zx}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zy\right)+{\epsilon }_{zy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zx\right)\end{array}\right]\stackrel{\to }{{\kappa }_{\Im }}\end{array}$ (A1)

Six distinct populating factors need to be determined, one related to each field component and inverse material property combination. The tensor elements of the inverse relative permittivity and permeability can be series expanded as:

${\epsilon }_{ab}^{■}={\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{ab}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}$ (A2)

${\mu }_{ab}^{■}={\displaystyle {\sum }_{p_{\mu }{q}_{\mu }{n}_{\mu }}{\kappa }_{{\mu }_{ab}}{\text{e}}^{j\left(G_{p_{\mu }}x+G_{q_{\mu }}y+G_{n_{\mu }}z\right)}}$ (A3)

and the field components in series, similar to (9), the derivatives are readily obtained for each term through the following intermediaries:

${\epsilon }_{iy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xz\right)={\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{iy}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_n+k_z\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$ (A4)

${\epsilon }_{iz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(xy\right)=-{\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{i3}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_q+k_y\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$ (A5)

${\epsilon }_{ix}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yz\right)=-{\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{i1}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_n+k_z\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$ (A6)

${\epsilon }_{iz}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(yx\right)={\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{i3}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_p+k_x\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$ (A7)

${\epsilon }_{ix}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zy\right)={\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{i1}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_q+k_y\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$ (A8)

${\epsilon }_{iy}^{■}{\Im }_{\Downarrow }{}_{{}_{\nabla \times }}\left(zx\right)=-{\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{i2}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}{\displaystyle {\sum }_{pqn}j}\left(G_p+k_x\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}$< (A9)

These are then multiplied by the complex conjugate of a single normalized basis function with complex propagation exponential as well and integrated over the unit cell of the structure being examined. Using (A4) as the computation example, the integral to evaluate is:

$\begin{array}{l}\frac{1}{V}{\displaystyle \int {\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }}{\kappa }_{{\epsilon }_{iy}}{\text{e}}^{j\left(G_{p_{\epsilon }}x+G_{q_{\epsilon }}y+G_{n_{\epsilon }}z\right)}}}{\displaystyle {\sum }_{pqn}j}\left(G_n+k_z\right){\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}\\ \times {\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}{\text{e}}^{-j\left(G_{p^{\prime }}x+G_{q^{\prime }}y+G_{n^{\prime }}z\right)}{\text{e}}^{-j\left(k_{x}x+k_{y}y+k_{z}z\right)}\text{d}V\\ =\frac{1}{V}{\displaystyle \int {\displaystyle {\sum }_{p_{\epsilon }{q}_{\epsilon }{n}_{\epsilon }pqn}j{\kappa }_{{\epsilon }_{iy}}\left(G_n+k_z\right){\text{e}}^{j\left(\left(G_{p_{\epsilon }}+G_p-G_{p^{\prime }}\right)x+\left(G_{q_{\epsilon }}+G_q-G_{q^{\prime }}\right)y+\left(G_{n_{\epsilon }}+G_n-G_{n^{\prime }}\right)z\right)}\text{d}V}}\\ =j{\kappa }_{{\epsilon }_{iy}}\left(G_n+k_z\right)\delta \left(\begin{array}{c}{p}_{\epsilon }+p,p^{\prime }\\ q_{\epsilon }+q,q^{\prime }\\ n_{\epsilon }+n,n^{\prime }\end{array}\right)\end{array}$ (A10)

Dropping the summation, delta functions and representing each matrix element block expression by the appropriate ${\kappa }_{{\epsilon }_{im}}\left(G_{pqn}+k_i\right)$ factor, the matrix operator can be expressed in a compact form as

$\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}{\kappa }_{{\epsilon }_{xy}}\left(G_n+k_z\right)-{\kappa }_{{\epsilon }_{xz}}\left(G_q+k_y\right)& -{\kappa }_{{\epsilon }_{xx}}\left(G_n+k_z\right)+{\kappa }_{{\epsilon }_{xz}}\left(G_p+k_x\right)& {\kappa }_{{\epsilon }_{xx}}\left(G_q+k_y\right)-{\kappa }_{{\epsilon }_{xy}}\left(G_p+k_x\right)\\ {\kappa }_{{\epsilon }_{yy}}\left(G_n+k_z\right)-{\kappa }_{{\epsilon }_{yz}}\left(G_q+k_y\right)& -{\kappa }_{{\epsilon }_{yx}}\left(G_n+k_z\right)+{\kappa }_{{\epsilon }_{yz}}\left(G_p+k_x\right)& {\kappa }_{{\epsilon }_{yx}}\left(G_q+k_y\right)-{\kappa }_{{\epsilon }_{yy}}\left(G_p+k_x\right)\\ {\kappa }_{{\epsilon }_{zy}}\left(G_n+k_z\right)-{\kappa }_{{\epsilon }_{zz}}\left(G_q+k_y\right)& -{\kappa }_{{\epsilon }_{zx}}\left(G_n+k_z\right)+{\kappa }_{{\epsilon }_{zz}}\left(G_p+k_x\right)& {\kappa }_{{\epsilon }_{zx}}\left(G_q+k_y\right)-{\kappa }_{{\epsilon }_{zy}}\left(G_p+k_x\right)\end{array}\right]$ (A11)

Each matrix element in (A11) is in fact a square matrix with the order equal to the number of basis functions utilized in the series expansion of the field components. The matrix (A11) can be segmented into several matrices that are added together. The maximum value for $k_x=k_{x,\max }$, $k_y=k_{y,\max }$, $k_z=k_{z,\max }$ are utilized. These can be individually scaled by $\left(\Omega ,\Psi ,\Phi \right)$ to generate any of the k-vectors required when computing the band diagram. The expanded matrix form is:

$\begin{array}{l}\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}{\kappa }_{{\epsilon }_{xy}}{G}_n-{\kappa }_{{\epsilon }_{xz}}{G}_q& -{\kappa }_{{\epsilon }_{xx}}{G}_n+{\kappa }_{{\epsilon }_{xz}}{G}_p& {\kappa }_{{\epsilon }_{xx}}{G}_q-{\kappa }_{{\epsilon }_{xy}}{G}_p\\ {\kappa }_{{\epsilon }_{yy}}{G}_n-{\kappa }_{{\epsilon }_{yz}}{G}_q& -{\kappa }_{{\epsilon }_{yx}}{G}_n+{\kappa }_{{\epsilon }_{yz}}{G}_p& {\kappa }_{{\epsilon }_{yx}}{G}_q-{\kappa }_{{\epsilon }_{yy}}{G}_p\\ {\kappa }_{{\epsilon }_{zy}}{G}_n-{\kappa }_{{\epsilon }_{zz}}{G}_q& -{\kappa }_{{\epsilon }_{zx}}{G}_n+{\kappa }_{{\epsilon }_{zz}}{G}_p& {\kappa }_{{\epsilon }_{zx}}{G}_q-{\kappa }_{{\epsilon }_{zy}}{G}_p\end{array}\right]\\ +j\Omega k_{x,\max }\left[\begin{array}{ccc}0& {\kappa }_{{\epsilon }_{xz}}& -{\kappa }_{{\epsilon }_{xy}}\\ 0& {\kappa }_{{\epsilon }_{yz}}& -{\kappa }_{{\epsilon }_{yy}}\\ 0& {\kappa }_{{\epsilon }_{zz}}& -{\kappa }_{{\epsilon }_{zy}}\end{array}\right]+j\Psi k_{y,\max }\left[\begin{array}{ccc}-{\kappa }_{{\epsilon }_{xz}}& 0& {\kappa }_{{\epsilon }_{xx}}\\ -{\kappa }_{{\epsilon }_{yz}}& 0& {\kappa }_{{\epsilon }_{yx}}\\ -{\kappa }_{{\epsilon }_{zz}}& 0& {\kappa }_{{\epsilon }_{zx}}\end{array}\right]\\ +j\Phi k_{z,\max }\left[\begin{array}{ccc}{\kappa }_{{\epsilon }_{xy}}& -{\kappa }_{{\epsilon }_{xx}}& 0\\ {\kappa }_{{\epsilon }_{yy}}& -{\kappa }_{{\epsilon }_{yx}}& 0\\ {\kappa }_{{\epsilon }_{zy}}& -{\kappa }_{{\epsilon }_{zx}}& 0\end{array}\right]\end{array}$ (A12)

A compact form of (A12) is:

$\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]=j{\left[G\right]}_{\epsilon }+j\Omega {\left[k_{x,\max }\right]}_{\epsilon }+j\Psi {\left[k_{y,\max }\right]}_{\epsilon }+j\Phi {\left[k_{z,\max }\right]}_{\epsilon }$ (A13)

and duplicated for the twin operator, $\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]$ :

$\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]=j{\left[G\right]}_{\mu }+j\Omega {\left[k_{x,\max }\right]}_{\mu }+j\Psi {\left[k_{y,\max }\right]}_{\mu }+j\Phi {\left[k_{z,\max }\right]}_{\mu }$ (A14)

Once the material is specified, the material expansion coefficient space for the series (A2), (A3) and (13) can be determined by orthogonal integration over a unit cell. The matrix operators $\left[\bar{G}\right],\left[\overline{k_{x,\max }}\right],\left[\overline{k_{y,\max }}\right],\left[\overline{k_{z,\max }}\right]$ can be populated using the maximum component values value of the propagation constant. The position in k-space, also the point in the band diagram is determined through the propagation constant component scaling factors $\left(\Omega ,\Psi ,\Phi \right)$.

For the example of computations presented here, the materials involved are non-magnetic. The structure is oriented with planes contained in $\left(y,z\right)$ plane of the coordinate system. The expansion spaces $\left\{G_q\right\}=\left\{0\right\},\left\{G_n\right\}=\left\{0\right\}$. The $\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]$ matrix operator simplifies to:

$\begin{array}{l}\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -G_p\\ 0& G_p& 0\end{array}\right]+j\Omega k_{x,\max }\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right]\\ +j\Psi k_{y,\max }\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right]+j\Phi k_{z,\max }\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\end{array}$ (A15)

The matrix form of $\left[\stackrel{\leftrightarrow }{{\epsilon }_{rD}{}^{-1}}\stackrel{\to }{\nabla }\times \right]$ when the electric properties are anisotropic has the form:

$\begin{array}{l}\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -{\kappa }_{{\epsilon }_{yy}}{G}_p\\ 0& {\kappa }_{{\epsilon }_{zz}}{G}_p& 0\end{array}\right]+j\Omega k_{x,\max }\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -{\kappa }_{{\epsilon }_{yy}}\\ 0& {\kappa }_{{\epsilon }_{zz}}& 0\end{array}\right]\\ +j\Psi k_{y,\max }\left[\begin{array}{ccc}0& 0& {\kappa }_{{\epsilon }_{xx}}\\ 0& 0& 0\\ -{\kappa }_{{\epsilon }_{zz}}& 0& 0\end{array}\right]+j\Phi k_{z,\max }\left[\begin{array}{ccc}0& -{\kappa }_{{\epsilon }_{xx}}& 0\\ {\kappa }_{{\epsilon }_{yy}}& 0& 0\\ 0& 0& 0\end{array}\right]\end{array}$ (A16)

Forms (A15) and (A16) are useful when examining the effect that the magnetoelectric material properties have on the states and band structure. The recollected form of the operators, for the calculation example presented here are:

$\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& -k_z& k_y\\ k_z& 0& -G_p-k_x\\ -k_y& G_p+k_x& 0\end{array}\right]$ (A17)

$\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\to }{\nabla }\times \right]\Rightarrow j\left[\begin{array}{ccc}0& -k_z{\kappa }_{{\epsilon }_{xx}}& k_y{\kappa }_{{\epsilon }_{xx}}\\ k_z{\kappa }_{{\epsilon }_{yy}}& 0& -{\kappa }_{{\epsilon }_{yy}}\left(G_p+k_x\right)\\ -k_y{\kappa }_{{\epsilon }_{zz}}& {\kappa }_{{\epsilon }_{zz}}\left(G_p+k_x\right)& 0\end{array}\right]$ (A18)

Appendix B: Matrix Populating Expressions for Bianisotropic Containing Operators

The expressions provided are used to populate matrices constructed from terms containing the magnetoelectric contributions. Each of the matrix operators $j\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]$ and $-j\left[{\stackrel{\leftrightarrow }{{\mu }_{rB}}}^{-1}\stackrel{\leftrightarrow }{{\epsilon }_B}\right]$ can be populated using the same numerical approach. Working from the first of these the matrix written in expanded form is:

$j\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\stackrel{\to }{{\kappa }_{\Im }}\Rightarrow j{\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{yx}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{zx}& {\epsilon }_{zy}& {\epsilon }_{zz}\end{array}\right]}^{-1}\left[\begin{array}{ccc}{\mu }_{Dxx}& {\mu }_{Dxy}& {\mu }_{Dxz}\\ {\mu }_{Dyx}& {\mu }_{Dyy}& {\mu }_{Dyz}\\ {\mu }_{Dzx}& {\mu }_{Dzy}& {\mu }_{Dxx}\end{array}\right]\stackrel{\to }{{\kappa }_{\Im }}$ (B1)

$j\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\stackrel{\to }{{\kappa }_{\Im }}\Rightarrow j\left[\begin{array}{ccc}{\epsilon }_{xx}^{■}{\mu }_{Dxx}+{\epsilon }_{xy}^{■}{\mu }_{Dyx}+{\epsilon }_{xz}^{■}{\mu }_{Dzx}& {\epsilon }_{xx}^{■}{\mu }_{Dxy}+{\epsilon }_{xy}^{■}{\mu }_{Dyy}+{\epsilon }_{xz}^{■}{\mu }_{Dzy}& {\epsilon }_{xx}^{■}{\mu }_{Dxz}+{\epsilon }_{xy}^{■}{\mu }_{Dyz}+{\epsilon }_{xz}^{■}{\mu }_{Dzz}\\ {\epsilon }_{yx}^{■}{\mu }_{Dxx}+{\epsilon }_{yy}^{■}{\mu }_{Dyx}+{\epsilon }_{yz}^{■}{\mu }_{Dzx}& {\epsilon }_{yx}^{■}{\mu }_{Dxy}+{\epsilon }_{yy}^{■}{\mu }_{Dyy}+{\epsilon }_{yz}^{■}{\mu }_{Dzy}& {\epsilon }_{yx}^{■}{\mu }_{Dxz}+{\epsilon }_{yy}^{■}{\mu }_{Dyz}+{\epsilon }_{yz}^{■}{\mu }_{Dzz}\\ {\epsilon }_{zx}^{■}{\mu }_{Dxx}+{\epsilon }_{zy}^v{\mu }_{Dyx}+{\epsilon }_{zz}^{■}{\mu }_{Dzx}& {\epsilon }_{zx}^{■}{\mu }_{Dxy}+{\epsilon }_{zy}^{■}{\mu }_{Dyy}+{\epsilon }_{zz}^{■}{\mu }_{Dzy}& {\epsilon }_{zx}^{■}{\mu }_{Dxz}+{\epsilon }_{zy}^{■}{\mu }_{Dyz}+{\epsilon }_{zz}^{■}{\mu }_{Dzz}\end{array}\right]\stackrel{\to }{{\kappa }_{\Im }}$ (B2)

Note that the complex scaled magnetic field has been included in the term being developed. The inclusion is required when matrix populating expressions are being generated. In general, each of the 9 tensor elements is composed of a sum of two terms having the same form. At each discretization grid point, the permittivity value is multiplied with the permeability value and the product decomposed as though it was the effective material value at that grid point. The approach avoids the need to decompose the permittivity and then decompose the permeability and then multiply the two series together. Either approach should provide the same result. The series expansion for any of these terms is:

$\begin{array}{l}j{\epsilon }_{ab}^{■}{\mu }_{Dcd}\left[{\kappa }_f\right]\\ =j{\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}{\text{e}}^{j\left(G_{p_{\epsilon \mu }}x+G_{q_{\epsilon \mu }}y+G_{n_{\epsilon \mu }}z\right)}}{\displaystyle {\sum }_{pqn}{\kappa }_{pqn}{\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{j\left(k_{x}x+k_{y}y+k_{z}z\right)}}\end{array}$ (B3)

These are then multiplied by the complex conjugate of a single normalized basis function with complex propagation exponential as well and integrated over the unit cell of the structure being examined. Using (B3) as the computation example, the integral to evaluate is (propagation exponential factors are pre-cancelled):

$\frac{j}{V}{\displaystyle \int {\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}{\text{e}}^{j\left(G_{p_{\epsilon \mu }}x+G_{q_{\epsilon \mu }}y+G_{n_{\epsilon \mu }}z\right)}}}{\displaystyle {\sum }_{pqn}{\kappa }_{pqn}{\text{e}}^{j\left(G_{p}x+G_{q}y+G_{n}z\right)}{\text{e}}^{-j\left(G_{p^{\prime }}x+G_{q^{\prime }}y+G_{n^{\prime }}z\right)}\text{d}V}$ (B4)

which simplifies to:

$\begin{array}{l}\frac{j}{V}{\displaystyle \int {\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }pqn}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}{\kappa }_{pqn}{\text{e}}^{j\left(\left\{G_{p_{\epsilon \mu }}+G_p-G_{p^{\prime }}\right\}x+\left\{G_{q_{\epsilon \mu }}+G_q-G_{q^{\prime }}\right\}y+\left\{G_{n_{\epsilon \mu }}+G_n-G_{n^{\prime }}\right\}z\right)}\text{d}V}}\\ =j{\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }pqn}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}{\kappa }_{pqn}\delta \left(\begin{array}{c}{p}_{\epsilon \mu }+p,p^{\prime }\\ q_{\epsilon \mu }+q,q^{\prime }\\ n_{\epsilon \mu }+n,n^{\prime }\end{array}\right)}\end{array}$ (B5)

This expression can be recast such that the SC column vector expansion coefficients are external:

$j{\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}\delta \left(\begin{array}{c}{p}_{\epsilon \mu }+p,p^{\prime }\\ q_{\epsilon \mu }+q,q^{\prime }\\ n_{\epsilon \mu }+n,n^{\prime }\end{array}\right)\left[{\kappa }_f\right]}$ (B6)

In the matrix element populating process a considerable amount of bookkeeping is required to generate meaningful matrix operators. The row being populated is specified by the field indices, $p^{\prime }{q}^{\prime }{n}^{\prime }$. The column being populated is specified by the field indices, $pqn$. The matrix block being populated is specified by the field vector component. The matrix element numerical value is generated using the summation in (B6) simplified to a single summation since specifying the indices for the permittivity results in a single value of the bianisotropic expansion indices satisfying the non-zero value of the delta function:

$\left[Row/Column\right]\left[p^{\prime }{q}^{\prime }{n}^{\prime }/pqn\right]=j{\displaystyle {\sum }_{p_{\epsilon \mu }{q}_{\epsilon \mu }{n}_{\epsilon \mu }}{\kappa }_{{\epsilon }_{ab}{\mu }_{cd}}\delta \left(\begin{array}{c}{p}_{\epsilon \mu }+p,p^{\prime }\\ q_{\epsilon \mu }+q,q^{\prime }\\ n_{\epsilon \mu }+n,n^{\prime }\end{array}\right)}$ (B7)

A considerable amount of matrix populating simplification can take place if the inverse permittivity matrix is diagonal. It may still be anisotropic with unequal diagonal elements. In this case, the simplified matrix is:

$j\left[{\stackrel{\leftrightarrow }{{\epsilon }_{rD}}}^{-1}\stackrel{\leftrightarrow }{{\mu }_D}\right]\stackrel{\to }{{\kappa }_{\Im }}\Rightarrow j\left[\begin{array}{ccc}{\epsilon }_{xx}^{-1}{\mu }_{Dxx}& {\epsilon }_{xx}^{-1}{\mu }_{Dxy}& {\epsilon }_{xx}^{-1}{\mu }_{Dxz}\\ {\epsilon }_{yy}^{-1}{\mu }_{Dyx}& {\epsilon }_{yy}^{-1}{\mu }_{Dyy}& {\epsilon }_{yy}^{-1}{\mu }_{Dyz}\\ {\epsilon }_{zz}^{-1}{\mu }_{Dzx}& {\epsilon }_{zz}^{-1}{\mu }_{Dzy}& {\epsilon }_{zz}^{-1}{\mu }_{Dzz}\end{array}\right]\stackrel{\to }{{\kappa }_{\Im }}$ (B8)

Further simplifications to the matrix operator are possible once the zero elements of the bianisotropic permeability are specified as is the case in the computation examples provided in the main text.

Appendix C: Non-Degenerate Perturbation Theory Expressions

The 6 field matrix systems of (15) is structured such that the magnetoelectric contributions can be extracted and appear as a perturbation, expression (28), to a structure displaying only anisotropic properties. The unperturbed matrix system may be written as:

$\left[E{H}_6\right]|{\varphi }_{oi}\rangle ={\lambda }_{oi}|{\varphi }_{oi}\rangle$ (C1)

with unperturbed eigenvalues, ${\lambda }_{oi}$, and eigenvectors $|{\varphi }_{oi}\rangle$. When the perturbation is included the matrix expression becomes:

$\left\{\left[E{H}_6\right]+\left[P_6\right]\right\}|{\varphi }_i\rangle ={\lambda }_i|{\varphi }_i\rangle$ (C2)

with new eigenvalues, ${\lambda }_i$, and eigenvectors, $|{\varphi }_i\rangle$. If the perturbations are small and the states are non-degenerate, the eigenvalues will be very close to the original values, offset by a small correction:

${\lambda }_i={\lambda }_{oi}+\partial {\lambda }_i$ (C3)

The eigenstates of the perturbed system will also very closely match those of the unperturbed states and may be expressed as a slight correction:

$|{\varphi }_i\rangle =|{\varphi }_{oi}\rangle +|\partial {\varphi }_i\rangle$ (C4)

The full matrix system eigenvalues can be computed from, $\left[P_6\right]$ defined in (28):

${\lambda }_i={\lambda }_{oi}+\langle {\varphi }_{oj}|\left[P_6\right]|{\varphi }_{oi}\rangle$ (C5)

and the full matrix system eigenvectors can be computed from:

$|{\varphi }_i\rangle =|{\varphi }_{oi}\rangle +{\displaystyle {\sum }_j\frac{\langle {\varphi }_{oj}|\left[P_6\right]|{\varphi }_{oi}\rangle }{{\lambda }_{oi}-{\lambda }_{oj}}|{\varphi }_{oj}\rangle }$ (C6)

If the states are degenerate, then degenerate perturbation theory must be applied using a linear combination of the degenerate states. The general result is that the perturbation will lift the degeneracy. The perturbed states are highly determined through coupling elements of the form in the non-degenerate perturbation analysis. The presence and location of the non-zero matrix elements in $\left[P_6\right]$ will determine the type of field component coupling that may occur. Such an approach is useful in interpreting the results of incrementally applied magnetoelectric values in the constituent expressions for the materials of a photonic structure.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Farjallah, H., Hassiba, H. and Ghozlen, M.H.B. (2016) Application of the Plane Wave Expansion and Stiffness Matrix Methods to Study Transmission Properties and Guided Mode of Phononic Plates. Wave Motion, 64, 68-78. https://doi.org/10.1016/j.wavemoti.2016.03.001
[2]	Zhou, Y., Chen, H. and Zhou, A. (2019) Plane Wave Methods for Quantum Eigenvalue Problems of Incommensurate Systems. Journal of Computational Physics, 384, 99-113. https://doi.org/10.1016/j.jcp.2019.02.003
[3]	Pan, Y., Dai, X., de Gioronocoli, S., Gong, X., Rignanese, G. and Zhou, A. (2017) A Parallel Orbital-Updated Based Plane-Wave Basis Method for Electronic Structure Calculations. Journal of Computational Physics, 348, 482-492. https://doi.org/10.1016/j.jcp.2017.07.033
[4]	Joannopoulos, J.D., Meade, R.D. and Winn, J.N. (1995) Photonic Crystals Modeling the Flow of Light. Princeton University Press, Princeton.
[5]	Johnson, S.G. and Joannopoulos, J.D. (2002) Photonic Crystal the Road from Theory to Practice. Kluwer Academic Publishing, Boston.
[6]	Sukhoivanov, I.A. and Guryev, I.V. (2009) Photonic Crystals Physics and Practical Modeling. Springer, Berlin.
[7]	Liew, S.F., Knitter, S., Xiong, W. and Cao, H. (2015) Photonic Crystals with Topological Defects. Physical Review A, 91, Article ID: 023811.
[8]	Cornelious, S.R., Sahoo, T., Sujatha, A. and Chelliah, P. (2019) Study of Symmetric and Asymmetric Defect Positions in a Photonic Crystal Supercell. Optics Communications, 450, 22-27. https://doi.org/10.1016/j.optcom.2019.05.014
[9]	Zhang, T., Sun, J., Yang, Y. and Li, Z. (2018) Photonic Crystal Filter Based on Defect Mode and Waveguide Mode Symmetry Matching. Optics Communications, 428, 53-56. https://doi.org/10.1016/j.optcom.2018.07.039
[10]	Mackay, T.G. and Lakhtakia, A. (2019) Electromagnetic Anisotropy and Bianisotropy A Field Guide. Second Edition, World Scientific, London.
[11]	Razzaz, F. and Alkanhal, M.A.S. (2018) Resonances in Bianisotropic Layers. IEEE Photonics Journal, 10, Article ID: 6100112. https://doi.org/10.1109/JPHOT.2017.2782742
[12]	Jakoby, B. and Baghai-Wadji, A. (1996) Analysis of Bianisotropic Layered Structures with Laterally Periodic Inhomogeneities—An Eigenvector Formulation. IEEE Transactions on Antennas and Propagation, 44, 615-626. https://doi.org/10.1109/8.496247
[13]	Kong, J.A. (2008) Electromagnetic Wave Theory. EMW Publishing, Cambridge.
[14]	Post, E.J. (1967) Sagnac Effect. Reviews of Modern Physics, 39, 475-493. https://doi.org/10.1103/RevModPhys.39.475
[15]	Lakhtakia, A. and Mackay, T. (2004) Towards Gravitationally Assisted Negative Refraction of Light in Vacuum. Journal of Physics A: Mathematical and General, 37, L505-L510.
[16]	Ivannenko, Y., Gustafsson, M. and Nordebo, S. (2019) Optical Theorems and Physical Bounds on Absorption in Lossy Media. Optics Express, 27, 34323-34341. https://doi.org/10.1364/OE.27.034323
[17]	Xu, X., Quan, B., Gu, C. and Wang, L. (2006) Bianisotropic Response of Microfrabricated Metamaterials in the Terahertz Region. Journal of the Optical Society of America B, 23, 1174-1180. https://doi.org/10.1364/JOSAB.23.001174
[18]	Mulkey, T., Dillies, J. and Durach, M. (2018) Inverse Problem of Quartic Photonics. Optics Letters, 43, 1266-1229. https://doi.org/10.1364/OL.43.001226
[19]	Gangaraj, A.A.H. and Hanson, G.W. (2017) Berry Phase, Berry Connection, and Chern Number for a Continuum Bianisotropic Material from a Classical Electromagnetics Perspective. IEEE Journal on Multiscale and Multiphysics Computational Techniques, 2, 3-17. https://doi.org/10.1109/JMMCT.2017.2654962
[20]	Durach, M., Williamson, R., Laballe, M. and Mulkey, T. (2020) Tri- and Tetra-Hyperbolic Iso-Frequency Topologies Complete Classification of Bi-Anisotropic Materials. Applied Sciences, 10, 763. https://doi.org/10.3390/app10030763
[21]	Durach, M. (2020) Tetra-Hyperbolic and Tri-Hyperbolic Optical Phase in Anisotropic Metamaterials without Magnetoelectric Coupling Due to Hybridization of Plasmonic and Magnetic Block High-K Polaritons. Optics Communications, 476, Article ID: 126349.
[22]	Kong, J.A. (1972) Theorems of Bianisotropic Media. Proceedings of the IEEE, 60, 1036-1046. https://doi.org/10.1109/PROC.1972.8851
[23]	Kang, M., Zhang, T., Zhao, B., Sun, L. and Chen, J. (2021) Chirality of Exceptional Points in Bianisotropic Metasurfaces. Optics Express, 29, 11582-11590. https://doi.org/10.1364/OE.419511
[24]	Heiss, W.D. (2004) Exceptional Points of Non-Hermitian Operators. Journal of Physics A: Mathematical and General, 37, 2455-2464. https://doi.org/10.1088/0305-4470/37/6/034
[25]	Miri, M. and Alu, A. (2019) Exceptional Points in Optics and Photonics. Science, 363, eaar7709. https://doi.org/10.1126/science.aar7709