Three-Dimensional Time-Harmonic Electromagnetic Scattering Problems from Bianisotropic Materials and Metamaterials: Reference Solutions Provided by Converging Finite Element Approximations

: A recently developed theory is applied to deduce the well posedness and the ﬁnite element approximability of time-harmonic electromagnetic scattering problems involving bianisotropic media in free-space or inside waveguides. In particular, three example problems are considered of which one deals with scattering from plasmonic gratings that exhibit bianisotropy while the other two deal with bianisotropic obstacles inside waveguides. The hypotheses that guarantee the reliability of the numerical results are veriﬁed, and the ranges of the constitutive parameters of the media involved for which the ﬁnite element solutions are guaranteed to be reliable are deduced. It is shown that, within these ranges, there can be significant bianisotropic effects for the practical media considered as examples. The ensured reliability of the obtained results can make them useful as benchmarks for other numerical approaches. To the best of our knowledge, no other tool can guarantee reliable solutions.


Introduction
Bianisotropic media have important applications in numerous practical problems ranging from the microwave to photonic frequency bands [1][2][3][4]. The electromagnetic problems involving such media admit analytical solutions only in very specialized cases, and numerical simulators are necessary to solve the vast majority of them.
In this context, the reliability of the numerical solvers is of utmost importance, and results guaranteeing the well posedness of the problems and the convergence of the numerical solutions are crucial. Some of the previous papers that addressed this issue were limited in their applications [5][6][7][8]. For instance, the work in [5] made strong assumptions on the losses to guarantee the reliability of the results. The results in [6] were derived for two-dimensional problems involving axially moving cylinders, whereas [7] dealt only with evolution problems inside cavities. As for [8], the constitutive parameters were taken to be smooth, which did not allow applications to radiation and scattering problems.
A considerable generalization of the conditions that allow ensuring the well posedness of the problems and the convergence of the finite element solutions was recently achieved in [9]. A set of non-restrictive hypotheses was shown to guarantee such results for three-dimensional time-harmonic problems involving bianisotropic media. The authors applied the theory to rotating axisymmetric D = (1/c 0 ) P E + L B in Ω, In the above equation, E, B, D, and H are complex valued functions defined in Ω and represent, respectively, the electric field, magnetic induction, electric displacement, and magnetic field, while c 0 is the speed of light in a vacuum. The space where we will seek E and H is [12] (p. 82; see also p. 69): where [12] (p. 48): Based on Maxwell's equations, boundary conditions, and constitutive relations, the following variational formulation of the problem can be deduced [5]: given ω > 0, the electric and magnetic current densities prescribed by the sources J e , J m ∈ (L 2 (Ω)) 3 and the known term f R ∈ L 2 t (Γ), involved in admittance boundary condition, find E ∈ U such that: where: a(u, v) = c 0 Q curl u, curl v 0,Ω − ω 2 c 0 P u, v 0,Ω − jω M u, curl v 0,Ω −jω L curl u, v 0,Ω + jω Y (n × u × n), n × v × n 0,Γ , and: In [9], we derived a set of sufficient conditions that guarantee the well posedness and finite element approximability of the problem. The developed theory was applied to problems involving rotating axisymmetric objects. In this paper, we apply the theory to a wider range of problems involving bianisotropic materials and metamaterials demonstrating the generality of the developments and obtaining interesting new solutions [1,10,11].
We recall important definitions and hypotheses to fix the notations that are required for proceeding with the application of the theory. The subscript i ∈ I identifying the subdomain Ω i may belong to two subsets, I a and I b of I, according to the properties of the media involved: i ∈ I a when they are anisotropic (that is, with L = M = 0) and i ∈ I b when they are bianisotropic. Any matrix A with complex entries can be split into A = A s − jA ss with A s = A+A * 2 and A ss = A * −A 2j . Some of the hypotheses will be stated using the alternative form of constitutive relations defined by: where the constitutive matrices are given by [13] As mentioned in Section 6 of [9], most of the hypotheses are readily satisfied for important practical problems. That leaves us with seven critical hypotheses (HM9-HM15 in [9]) on the media involved in the problem that need to be verified. Among them, HM9-HM12 of [9] are stated in terms of the constitutive relations involving κ, ν, χ, and γ and are restated here as H1-H4.
Among the above hypotheses, H2 needs to hold only in the subdomains Ω i , i ∈ I a . H1 can be verified separately in any subdomain Ω i , i ∈ I, whereas H3 and H4 are to be defined and verified only locally on any subdomain Ω i , i ∈ I b (see Remark 1 of [9]).
The local continuity of the tensors P, Q, L, and M can be assumed in most practical problems, which allows the definition of the following constants.
Hypothesis 7 (H7). C PS , C QS , C L , and C M (i.e., all media involved) are such that C QS − C L C M C PS > 0.
We refer to Section 6 of [9] for some hints about the calculation of the constants involved in the above conditions. In particular, we recall here Lemma 1 of [9] for easy reference, which is helpful in finding the constant involved in H5.

Lemma 1.
Suppose that P ss is uniformly positive definite in Ω el ⊂ Ω, that is ∃C 1 > 0 such that: Whenever Ω el = Ω, we can simply define C PS = C 1 .
Whenever Ω el is not the whole Ω, suppose that, in the complementary region, P s is uniformly positive or negative definite, that is ∃C 5 > 0 such that: Whenever Ω el = ∅, we simply have C PS = C 5 , and we can set: where λ min denotes the minimum of the magnitudes of the eigenvalues of the Hermitian symmetric matrix P s . Finally, whenever Ω el is neither the empty set nor the whole domain, under assumptions HM2 and HM3 of [9], condition H5 is satisfied with C PS given by: where C 3 > 0 is defined by: and α is such that Analogously, by replacing P with Q in Equations (11), (12), and (15) , we define, respectively, Ω ml and the constants C 2 > 0, C 4 > 0, and C 6 > 0 and deduce that condition H6 is satisfied if we set: whenever Ω ml = ∅, C QS = C 2 whenever Ω ml = Ω, or: α being such that 1 > α > > 0, when Ω ml = Ω and Ω ml = ∅.
With respect to the H4, if we define: the sufficient condition guaranteeing the validity of the inequality in the hypothesis can be expressed as:

Results and Discussion
In this section, we apply the theory developed in [9] to several classes of problems that could not be managed with the previous theories [5][6][7][8]. The conditions are established on the constitutive parameters of such problems, under which the well posedness and finite element approximability can be guaranteed. In particular, we apply the theory and obtain solutions for three examples of bianisotropic media, which are found in the open literature. The first one is that introduced in [1] where the authors considered plasmonic gratings, which are represented by an equivalent bianisotropic media. In this case, we study the scattering from a slab of the equivalent medium, which is placed in empty space, in accordance with the setup considered by the authors [1]. Next, the media introduced in [10] and [11] are analyzed. The authors there considered the scattering from the bianisotropic obstacles placed inside hollow waveguides. Although we stick to the original configurations proposed by the authors, our tools can be used to analyze other problems involving these media.
Under the conditions that guarantee that hypotheses H1-H7 are satisfied, the numerical solutions of these problems are computed. They can be used as reference solutions for other approaches and simulators because, on the one hand, for each problem, our theory guarantees the convergence of the sequence of approximations, and on the other hand, we verify that the outcome we present is close to the limit of the sequence by a stability analysis of the numerical solutions.
We use a first order edge element based Galerkin finite element simulator to obtain the solutions as described in Section 5 of [9].

Scattering from Plasmonic Gratings Behaving as Bianisotropic Metamaterials
In [1], the authors considered a plasmonic grating that exhibits bianisotropy at visible wavelengths. The bianisotropic media considered there are of the form: with: Here, ε x , ε y , and ε z are complex valued functions and are set to be equal to a unique value ε r . Moreover, µ x , µ y , and µ z are each set equal to one, and ξ 0 is taken to be real valued. The region occupied by the scatterer may be denoted as Ω s ⊂ Ω. The above form can be converted into the alternative form of constitutive relations [13] involving the P, Q, L, and M matrices defined in Equation (1), and the final result is shown in the following equations: Here, I 3 is the three by three identity matrix. The complementary region Ω \ Ω s is occupied by the empty space, which is characterized by P = c 0 ε 0 I 3 , Q = 1 c 0 µ 0 I 3 , L = M = 0. This problem cannot be managed by the previous theories and, in particular, by the theory in [5], which relied on strong hypothesis about losses.
The bianisotropic medium is lossy with the imaginary part Im(ε r ) < 0, whereas ξ 0 is assumed real here to avoid some longer calculations. Now, Lemma 1 can be applied to verify hypothesis H5.
Inside Ω s , P can be decomposed as P = P s − jP ss with: and P ss = P * −P 2j = −c 0 ε 0 Im(ε r )I 3 . Hence, we have Ω el = Ω s , the lossy region where P ss is uniformly positive definite and the complementary region with the free space where P s is uniformly positive definite. This means that the conditions of Lemma 1 are satisfied, and as a result, H5 is valid.
From the definitions (see Equations (11), (12), and (15)), To find the minimum of the two expressions in Equation (14), we note that, in the valid range, the value of the first expression decreases monotonically with α, whereas that of the second expression increases with it. The highest estimate for C PS is obtained when the two expressions have the same value. The value of α at which this happens can be evaluated by equating the two expressions and finding the positive root of the resulting quadratic equation. This value of α, denoted as α opt , is given by: Thus, we may simply write: As mentioned in [9], this does not mean that a better value of C PS cannot be found. For example, if Re(ε r ) − ξ 2 0 > 0, then P s is uniformly positive definite in Ω, and we can find another candidate for C PS , namely [9]). In particular, when Re(ε r ) − ξ 2 0 > 1 √ 2 , C 7 is always going to give a value for C PS that is higher than that obtained from the lemma.
In the rest of the subsection, the discussion focuses on the cases with Re(ε r ) < 0, which gives a non-definite P s in Ω. For this case, C 3 = c 0 ε 0 |Re(ε r ) − ξ 2 0 |, and we can directly use the value of C PS in Equation (29). Since the material is assumed to be non-magnetic, the direct application of the definition gives C QS = 1 c 0 µ 0 , and H6 is valid. Likewise, for the continuity constants C L and C M , we can choose the value ) > 0, which gives: Since the right-hand side of Equation (30) also depends on ξ 0 due to the presence of C 3 in the expression for α opt , we do not have a closed-form expression on the limit on ξ 0 below which H7 is satisfied. However, a graphical analysis can be done for estimating such a limit on |ξ 0 |, as shown in Figure 1. The value of Re(ε r ) is varied in the range (−5.0, −1.0), whereas Im(ε r ) assumes values in the range (−0.5, −0.1). It can be observed that for a fixed value of Im(ε r ), the range of ξ 0 over which H7 is valid steadily decreases as |Re(ε r )| increases. As for the dependence on Im(ε r ), the corresponding range increases when the medium becomes lossier due to higher |Im(ε r )|, as expected. |ξ 0 | guaranteeing that condition H7 is satisfied, for scattering problems involving different media considered in [1]. The curves are computed by assuming ξ 0 ∈ R, Re(ε r ) ∈ (−5.0, −1.0), Im(ε r ) ∈ (−0.5, −0.1), and µ r = 1.
Since outside the region occupied by the bianisotropic media (Ω \ Ω s ), we just have the empty space, H1 and H2 are trivially satisfied there. By Remark 1 of [9], since H1, H3, and H4 need to hold only locally, now we have to just analyze them inside Ω s occupied by the bianisotropic medium. We consider the alternative form of constitutive relations, which for the medium inside Ω s becomes as in Equations (31) to (33), for examining the validity of H1, H3 and H4 [13]: The constants of interest can be evaluated directly from the definitions. The determinants of κ and ν are, respectively, , which immediately give the values of C κ,d and C ν,d .
The inverses of the diagonal matrices κ and ν are just the diagonal matrices with the reciprocal entries. Applying Equations (40) and (41) of [9] gives the values of C κ,r and C ν,r .
From Equations (18) and (19), we can easily evaluate C γ,s and C χ,s . (40) The hypotheses H1 and H3 are valid due to the existence of the above constants. Using these constants, the value of K u can be calculated from Equation (20). The critical value of |ξ 0 | below which the condition in H4 is satisfied is plotted in Figure 2, with respect to either Re(ε r ) or Im(ε r ). The results show that the range of ξ 0 for which H4 holds true increases with the increase in |Re(ε r )|, while it is practically independent of Im(ε r ).
Let us try to understand the implications of the theory by applying it to the numerical solution of a specific problem involving the medium of interest. We consider the region with the scatterer Ω s to be a cube filled with homogeneous bianisotropic media. The surrounding region is filled with empty space, and the overall domain of numerical investigation, Ω, has a cubic shape as well and is concentric to Ω s . In the following, Ω and Ω s are characterized by sides of length 2 µm and 0.8 µm, respectively. The axes are taken along the sides of the cubic domain Ω, and the excitation is with a plane wave incident along the x axis, with the electric field polarized along the z axis, having a magnitude of 1 V/m and wavelength of 1 µm.
The solutions are obtained with a first order edge element based Galerkin finite element method. The boundary condition is enforced with Y equal to the admittance of a vacuum and with an inhomogeneous term f R , taking into account the incident field.
The domain is discretized uniformly using tetrahedral meshes. The meshing is done by first dividing the domain into small identical cubes, each of which is in turn divided into six tetrahedra. The parameter h denotes the maximum diameter of all the elements of the mesh [14] (p. 131), and in this case, it is simply given by the side of the small cubes times √ 3. To study the stability of the solution, we consider different levels of refinement of meshes ranked in order of h, ranging from "very coarse" to "very fine". For example, the mesh denoted as very coarse is characterized by cubes of sides 200 nm, and the resulting mesh has 1331 nodes, 6000 tetrahedral elements, and 1200 boundary faces. A summary of the information related to the four different refinements of the meshes that were used is given in Table 1. The outcomes related to the stability of the results of the simulations are shown in Figure 3 by plotting the magnitude of the z component of the electric field along a line parallel to the y axis and passing through the center of gravity of the domain. The difference between successive refinements progressively decreases, and the fine and very fine meshes give stable solutions. The well posedness and convergence result that was predicted using the theory guarantee that our solutions are reliable.  Here, E z denotes the solution obtained with ξ 0 = −0.41, and E z,0 is the solution when ξ 0 = 0, while the plot shows the magnitude of the difference |E z − E z,0 | along with |E z |. The magnitude of the difference is as large as 50% of the incident field. Similarly, the results along a line parallel to the z axis and passing through the center of gravity of the domain are shown in Figure 5, and we get a difference of around 30% of the incident field.
These non-negligible effects imply that to get accurate results, it is necessary to consider the bianisotropy of the medium. Hence, the reliability of the finite element solution in the presence of bianisotropy is important for getting good results for these problems. The application of our theory gives the conditions under which we can guarantee such reliability. The solutions obtained can serve as references for other approaches.

Scattering from Chiral Obstacles in a Waveguide
In [10], the authors considered a metallic waveguide, which was hollow except for an obstacle characterized by a chiral medium with the following constitutive relations.
Here, ε r , µ r , and ξ c are strictly positive real quantities. Thus, from Equation (1), we can easily identify P, Q, L, and M, which are given below: In Section 6 of [5], it was shown that this media could not be managed by the theory developed there, independently of the value of ξ c ∈ R and of any other material involved in the model of interest. However, we show that the generality of the recently developed theory in [9] allows us to apply it to obtain the conditions for well posedness and finite element approximability for this kind of problem of practical interest.
Let us analyze the validity of the hypotheses by considering, as did the authors in [10], ε r ≥ 1 and µ r = 1. We can make use of Lemma 1 of [9] to check H5. P s = P and is equal to ε 0 ε r c 0 inside the material and simply ε 0 c 0 outside. Since Ω el = ∅, a value of C PS can be found by using Equation (13). In particular, H5 is satisfied for C PS = ε 0 c 0 . The hypothesis H6 is also trivially valid with C QS = 1 µ 0 c 0 . By Equations (32) and (33) of [9], C L = C M = ξ c . Then, the inequality in hypothesis H7 becomes C QS − C L C M C PS = c 0 (ε 0 − µ 0 ξ 2 c ) > 0, which implies: This is not a small value considering the chiral effects reported in [10]. As the region outside the obstacle is empty space, H2 is trivially satisfied, and by Remark 1 of [9], we need to verify that the hypotheses H1, H3, and H4 hold true locally inside the region occupied by the bianisotropic medium. To do this, the suitable form of constitutive relations is in terms of κ, ν, γ, and χ, which are given by the following [13]: κ and ν are multiples of the identity matrix with eigenvalues 1 ε 0 ε r and ( ε 0 ε r +µ 0 ξ 2 c µ 0 ε 0 ε r ), respectively. The determinants are just the cubes of the eigenvalues, and hence, according to Equations (34) and (35) of [9], we get the values of C κ,d and C ν,d .
(50) C κ,s and C ν,s , by Equations (36) and (37) of [9], are in this case simply twice the eigenvalue of the corresponding diagonal matrix: The inverse of the matrices is also trivial, and Equations (40) and (41) of [9] simply evaluate to the reciprocals of the eigenvalues of κ and ν, respectively giving C κ,r and C ν,r : (54) From Equations (18) and (19), we get: Having shown that the hypotheses H1 and H3 are satisfied, we can use the above constants to calculate K u to verify H4. Figure 6 shows the dependence of K u on ξ c for various values of ε r . As the value of ε r increases, the hypothesis H4 remains valid for higher and higher values of ξ c . Figure 7 shows the plot of the critical value of ξ c below which H4 is satisfied against ε r . The limit of 2.654 × 10 −3 mho, arising from Equation (45) required to satisfy H7, is also shown in the same figure. It is seen that for low values of ε r , the tighter condition arises from the need to satisfy H4. For example, the limiting value is 5.6 × 10 −4 mho for ε r = 1, increases with ε r , and is 1.78 × 10 −3 mho for ε r = 10. The curve crosses the 2.654 × 10 −3 mho line at around ε r = 22.3, and above that value, Equation (45) imposes the stricter limit.  Figure 6. Plot of K u versus ξ c for the bianisotropic medium described in [10]. The plots are shown for various values of ε r . The hypothesis H4 is satisfied for K u < 1.  . The value of ξ c below which the hypothesis H4 is satisfied is plotted against ε r . The limit of 2.654 × 10 −3 mho, arising from Equation (45) required to satisfy H7, is also shown. Now, we consider a specific numerical problem for which the solution is calculated using our finite element simulator. A rectangular waveguide with a discontinuity due to a block of bianisotropic medium is considered as shown in Figure 8. In the simulation, the rectangular waveguide is characterized by a = 23 mm, b = 10 mm and has a length l = 40 mm. The obstacle is a parallelepiped with c = 11 mm, d = 5 mm and a length w = 10 mm. The origin of the axis is at the lower right corner of the near face of the waveguide shown in Figure 8. The obstacle ranges from x = 6 mm to x = 17 mm along the x axis, from y = 0 to y = 5 mm along the y axis, and from z = 15 mm to z = 25 mm along the z axis. The bianisotropic medium making up the obstacle is characterized by ε r = 5 and ξ c = 1.24 × 10 −3 mho. For this medium, K u = 0.98 < 1, and also, Equation (45) is satisfied; hence, all the hypotheses required to guarantee the well posedness and convergence of finite element solutions hold true. The waveguide is excited with TE 10 mode having an amplitude of 1 V/m and a frequency of 9 GHz. The input port is on the x-y plane, and the output port is closed on a matched homogeneous admittance boundary condition for the TE 10 mode in the empty waveguide. The details of the Galerkin finite element solver is the same as before. The tetrahedral meshes are obtained as discussed in the previous subsection by dividing the domain into small cubes, each of which is in turn subdivided into six tetrahedra. The stability of the solution is verified by checking the solutions for three different meshes, which are characterized by small cubes of sides 1 2 mm, 1 4 mm, and 1 6 mm, which are referred to as, respectively, "coarse", "fine", and "very fine" meshes. There are 10,824 nodes, 55,200 elements, and 6200 boundary faces in the coarse mesh, whereas the fine mesh has 79,947 nodes, 441,600 elements, and 24,800 boundary faces, and finally, the very fine mesh has 262,570 nodes, 1,490,400 elements, and 55,800 boundary faces. It was verified that the solutions obtained with these meshes are stable. For example, Figure 9 shows the magnitude of the x component of the electric field along a line parallel to the y axis and passing through the center of gravity of the domain with the different meshes and illustrates the stability of the result. It is noted that the x component of the electric field along this line is zero for the achiral case (ξ c = 0), and there is a difference of more than 30% of the incident field, which is induced by the bianisotropy. Figure 10 shows the result for the magnitude and phase of the x component of the electric field along the line parallel to the x axis and passing through the center of gravity of the domain. We have a difference of 20% of the magnitude of the incident field for the x component of the electric field along this line. Similarly, Figure 11 shows that the bianisotropic effect for the x component of the field along the line parallel to the z axis and passing through the center of gravity of the domain causes a difference in the magnitude of the electric field of around 13% of the incident field.   [10]. The plot for the bianisotropic case using ξ c = 1.24 × 10 −3 mho is compared with the solution obtained in the isotropic case using ξ c = 0. Figure 12 shows the y component of the electric field along a line parallel to the z axis and passing through the center of gravity of the domain. The bianisotropic effect is again not negligible and causes a difference of more than 10% of the incident field. A similar effect is present in the z component as can be seen in the plot along a line parallel to the x axis and passing through the center of gravity of the domain, which is given in Figure 13. We do not show the other figures to save space, but it is noted that the y component of electric field along the lines through the center of the domain and parallel to the x and y axes shows small differences in magnitude between the chiral and achiral cases, being less than five percent of the incident field. The z component on the other hand along the lines parallel to the y and z axes and passing through the center of the domain shows a difference of more than 15% of the incident field.  [10]. The plot for the bianisotropic case using ξ c = 1.24 × 10 −3 mho is compared with the solution obtained in the isotropic case using ξ c = 0.
Together, these results provide a point of reference for other approaches of solving such problems, owing to the reliability of the results provided here, which is guaranteed by the recently developed theory. The previous theory [5] was not able to manage these problems, and our results are therefore novel. In fact, to the best of our knowledge, there are no other approaches that are available to get reliable results for benchmarking. The significant bianisotropic effects demonstrated in the results show the practical importance of the theory for such media.

Reflection by a Short-Circuited Waveguide Half Filled with Bianisotropic Media
Another relevant bianisotropic medium was discussed in [11], where the authors considered a short-circuited rectangular waveguide, half of which was empty and the other half filled with a lossless bianisotropic material characterized by: where A is the matrix given by: and κ 0 is a positive real number. The hypotheses H5 and H6 are trivially valid with C PS = ε 0 c 0 and C QS = 1 µ 0 c 0 . Further, L * L = M * M = κ 2 0 A 2 whose eigenvalues are zero, κ 2 0 , and 4κ 2 0 . Therefore, by Equations (32) and (33) of [9], we get C L = C M = 2κ 0 . The condition in hypothesis H7 then becomes ε 0 c 0 > 0, which gives the following limit on κ 0 : The hypothesis H2 holds true since the anisotropic part is just empty space. H1, H3, and H4 can be studied using the alternative form of constitutive relations, which is characterized by the following matrices [13]: The determinants of κ and ν can be readily calculated, and by using Equations (34) and (35) of [9], we get C κ,d and C ν,d : C κ,s and C ν,s can be directly obtained from their definitions: By simple application of the definition: and by Equation (41) of [9], C ν,r evaluates to the reciprocal of the largest eigenvalue of the real matrix ν, which is given by: (69) Equations (18) and (19) give: The existence of the above constants verifies H1 and H3, and we can use them to calculate K u to check H4. It can be verified that K u is less than one when κ 0 ≤ 2.72 × 10 −4 mho, which is stricter than the limit obtained from Equation (60).
Finally, let us give some numerical solutions for this problem, which can be used as references for other approaches. The cross-section of the waveguide is 2 cm along the x axis and 1 cm along the y axis, and the length of the waveguide is 2 cm, half of which is filled with the bianisotropic medium characterized by κ 0 = 2.7 × 10 −3 mho (see Figure 1 of [11]). The origin is taken on the corner of the open face of the waveguide on the empty side. TE 10 mode is excited in the waveguide with a source of amplitude 1 V/m and a frequency of 12 GHz.
The first order edge element based Galerkin finite element method is used to obtain the solution. The meshing is carried out by dividing the domain into identical cubes, each of which is then subdivided into six tetrahedra. The stability of the result is ensured by evaluating the solutions on three meshes, termed as "coarse", "fine", and "very fine", characterized, respectively, by cubes of sides 1 mm, 1 2 mm, and 1 3 mm. The coarse mesh has 4852 nodes, 24,000 elements, and 3200 boundary faces. The fine mesh is composed of 35,301 nodes, 192,000 elements, and 12,800 boundary faces. The very fine mesh has 115,351 nodes, 648,000 elements, and 28,800 boundary faces. The results obtained with these meshes are stable. For example, Figure 14 shows the stability of the results obtained with the three meshes for the x component of the electric field along a line parallel to the y axis and passing through the center of gravity of the domain.  Figure 14. Stability of the solution for the problem involving the medium in [11]. The magnitude of the x component of the electric field is plotted for three different meshes along a line parallel to the y axis and passing through the center of gravity of the domain.
We provide the magnitudes and phases of the x component of the electric field obtained from the simulation in Figures 15-17. The bianisotropy causes a difference in magnitude of up to 14% of the incident field. The phases are also significantly affected by the bianisotropy of the medium. The figures for other components are not shown to save space, but it is noted that the y component is affected by the bianisotropy, showing a difference of up to about 10% of the incident field. The z component of the field along the line parallel to the x axis passing through the center of gravity of the domain does show a difference of about 10% of the incident field, but the magnitudes along the other directions are close to zero for both values of κ 0 considered.
Since the theory guarantees the reliability of these results, they can be used as references for other solvers and approaches. It is to be noted that the previous theory developed in [5] could not be applied to this medium since the same reasons mentioned there with respect to the medium in [10] are valid. Hence, our results show the generality of the new theory with respect to the application to interesting bianisotropic media. The results demonstrate that our theory can be applied to problems with significant bianisotropic effects. This consideration underlines its practical importance.  [11]. The plot for the bianisotropic case using κ 0 = 2.7 × 10 −4 mho is compared with the solution obtained in the isotropic case using κ 0 = 0.  [11]. The plot for the bianisotropic case using κ 0 = 2.7 × 10 −4 mho is compared with the solution obtained in the isotropic case using κ 0 = 0.  [11]. The plot for the bianisotropic case using κ 0 = 2.7 × 10 −4 mho is compared with the solution obtained in the isotropic case using κ 0 = 0.

Conclusions
In this paper, we discussed the application of a recently developed theory to electromagnetic scattering problems involving bianisotropic materials and metamaterials of practical interest. The range of constitutive parameters over which the reliability of the results was guaranteed was calculated for these problems, and the solutions were obtained. The ensured reliability of the results made them useful as benchmarks for other numerical techniques. To the best of our knowledge, none of the previous tools were able to get such benchmark solutions.