Open AccessArticle

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

Praveen Kalarickel Ramakrishnan

^†

and

Mirco Raffetto

^*,†

Department of Electrical, Electronic, Telecommunications Engineering and Naval Architecture, University of Genoa, Via Opera Pia 11a, I–16145 Genoa, Italy

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Electronics 2020, 9(7), 1065; https://doi.org/10.3390/electronics9071065

Submission received: 29 May 2020 / Revised: 19 June 2020 / Accepted: 24 June 2020 / Published: 29 June 2020

(This article belongs to the Special Issue Electromagnetic Scattering and Its Applications: From Low Frequencies to Photonics)

Download

Browse Figures

Versions Notes

Abstract

A recently developed theory is applied to deduce the well posedness and the finite 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 verified, and the ranges of the constitutive parameters of the media involved for which the finite 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.

Keywords:

electromagnetic scattering; time-harmonic electromagnetic fields; bianisotropic media; metamaterials; variational formulation; well posedness; finite element method; convergence of the approximation

1. 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 objects, where the effect of motion induced bianisotropy. However, the theoretical results derived were applicable to a much wider range of bianisotropic materials and metamaterials.

In this paper, we exploit the recently developed theory to obtain novel results. We consider examples of practical problems discussed in the open literature for which, to the best of our knowledge, none of the previous theories were able to guarantee the reliability of the numerical results. In particular, we study the electromagnetic scattering from plasmonic gratings, which exhibits bianisotropy [1], and from bianisotropic obstacles in waveguides [10,11]. We demonstrate the application of the theory in [9] to derive the conditions on the constitutive parameters of these problems that guarantee the reliability of the results. The numerical solutions of the problems are calculated under such conditions, which, owing to the reliability assured by the theory, can be used as references for other numerical solvers. As far as we are aware, no other tool is available that allows obtaining benchmark solutions for these problems.

The paper is organized as follows. In Section 2, the problem is defined and the theory is summarized to guide the reader in its application. Section 3 contains the main results of the paper, where the application of the theory is demonstrated and reliable solutions are obtained. Section 4 provides the conclusions.

2. Mathematical Description of the Problem

In this paper, we are interested in electromagnetic problems that involve bianisotropic media under time-harmonic excitation, which were studied in [9]. While the full details of the problem definition and results are available in the reference, here we provide a summary of the main points in order to ease the understanding of the present developments.

The problem is formulated in an open, bounded, and connected domain

Ω \in R^{3}

, which has a Lipschitz continuous stationary boundary denoted by

Γ

. To take into account electromagnetic problems involving inhomogeneous materials, we assume that

Ω

can be decomposed into m subdomains (see HD3 of [9]) denoted

Ω_{i}

i \in I = {1, \dots, m}

The time-harmonic sources imply that all the resulting fields are in turn time-harmonic and the assumed factor

e^{j ω t}

is ubiquitous and is suppressed. The media involved in the problem are linear and time-invariant and are considered to satisfy the following constitutive relations:

\{\begin{array}{l} D = (1 / c_{0}) P E + L B & in Ω, \\ H = M E + c_{0} Q B & in Ω . \end{array}

(1)

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

U = H_{L^{2}, Γ} (curl, Ω) = {v \in H (curl, Ω) | v \times n \in L_{t}^{2} (Γ)},

(2)

where [12] (p. 48):

L_{t}^{2} (Γ) = {v \in {(L^{2} (Γ))}^{3} | v \cdot n = 0 almost everywhere on Γ} .

(3)

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} \in {(L^{2} (Ω))}^{3}

and the known term

f_{R} \in L_{t}^{2} (Γ)

, involved in admittance boundary condition, find

E \in U

such that:

a (E, v) = l (v) \forall v \in U,

(4)

where:

\begin{matrix} a (u, v) = c_{0} {(Q curl u, curl v)}_{0, Ω} - \frac{ω^{2}}{c_{0}} {(P u, v)}_{0, Ω} - j ω {(M u, curl v)}_{0, Ω} \\ - j ω {(L curl u, v)}_{0, Ω} + j ω {(Y (n \times u \times n), n \times v \times n)}_{0, Γ}, \end{matrix}

(5)

and:

l (v) = - j ω {(J_{e}, v)}_{0, Ω} - c_{0} {(Q J_{m}, curl v)}_{0, Ω} + j ω {(L J_{m}, v)}_{0, Ω} - j ω {(f_{R}, n \times v \times n)}_{0, Γ} .

(6)

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 \in 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 \in I_{a}

when they are anisotropic (that is, with

L = M = 0

) and

i \in I_{b}

when they are bianisotropic. Any matrix A with complex entries can be split into

A = A_{s} - j A_{s s}

with

A_{s} = \frac{A + A^{*}}{2}

and

A_{s s} = \frac{A^{*} - A}{2 j}

. Some of the hypotheses will be stated using the alternative form of constitutive relations defined by:

\{\begin{cases} E = κ D + χ B & in Ω, \\ H = γ D + ν B & in Ω . \end{cases}

(7)

where the constitutive matrices are given by [13]

κ = c_{0} P^{- 1}

χ = - c_{0} P^{- 1} L

γ = c_{0} M P^{- 1}

, and

ν = c_{0} (Q - M P^{- 1} L)

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.

Hypothesis 1 (H1).

\exists \exists C_{κ, d} > 0, C_{ν, d} > 0 : | d e t e r m i n a n t (κ) | \geq C_{κ, d}, | d e t e r m i n a n t (ν) | \geq C_{ν, d}, \forall x \in {\bar{Ω}}_{i}, \forall i \in I

Hypothesis 2 (H2).

l_{1, 3}^{T} κ^{- 1} l_{1, 3} \neq 0, l_{1, 3}^{T} ν^{- 1} l_{1, 3} \neq 0 \forall l_{1, 3} \in R^{3}, l_{1, 3} \neq 0, \forall x \in {\bar{Ω}}_{i}, \forall i \in I_{a}

Hypothesis 3 (H3).

\exists \exists C_{κ, r} > 0, C_{ν, r} > 0 : | l_{1, 3, n}^{T} κ^{- 1} l_{1, 3, n} | \geq C_{κ, r}, | l_{1, 3, n}^{T} ν^{- 1} l_{1, 3, n} | \geq C_{ν, r} \forall l_{1, 3, n} \in R^{3} : {‖ l_{1, 3, n} ‖}_{2} = 1, \forall x \in {\bar{Ω}}_{i}, \forall i \in I_{b}

Hypothesis 4 (H4).

\exists \exists C_{κ, s} > 0, C_{ν, s} > 0

(\sum_{i, j = 1}^{3} | κ_{i j} |) - min_{i = 1, 2, 3} | κ_{i i} | \leq C_{κ, s} \forall x \in {\bar{Ω}}_{k}, \forall k \in I_{b},

(8)

(\sum_{i, j = 1}^{3} | ν_{i j} |) - min_{i = 1, 2, 3} | ν_{i i} | \leq C_{ν, s} \forall x \in {\bar{Ω}}_{k}, \forall k \in I_{b},

(9)

and κ, χ, γ, and ν satisfy:

\frac{4 ((\sum_{i, j = 1}^{3} | γ_{i j} |) - {min}_{i = 1, 2, 3} | γ_{i i} |) ((\sum_{i, j = 1}^{3} | χ_{i j} |) - {min}_{i = 1, 2, 3} | χ_{i i} |)}{(- C_{κ, s} + \sqrt{C_{κ, s}^{2} + 4 C_{κ, d} C_{κ, r}}) (- C_{ν, s} + \sqrt{C_{ν, s}^{2} + 4 C_{ν, d} C_{ν, r}})} < 1

(10)

\forall x \in {\bar{Ω}}_{k}, \forall k \in I_{b}

Among the above hypotheses, H2 needs to hold only in the subdomains

Ω_{i}

i \in I_{a}

. H1 can be verified separately in any subdomain

Ω_{i}

i \in I

, whereas H3 and H4 are to be defined and verified only locally on any subdomain

Ω_{i}

i \in 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.

$\exists C_{L} > 0$ : ${| (L curl u, v)}_{0, Ω} | \leq C_{L} {‖ curl u ‖}_{0, Ω} {‖ v ‖}_{0, Ω}$ for all $u \in H (curl, Ω)$ and $v \in {(L^{2} (Ω))}^{3}$ ,
$\exists C_{M} > 0$ : ${| (M u, curl v)}_{0, Ω} | \leq C_{M} {‖ u ‖}_{0, Ω} {‖ curl v ‖}_{0, Ω}$ for all $u \in {(L^{2} (Ω))}^{3}$ and $v \in H (curl, Ω)$ .

Finally, HM13–HM15 of [9] are reported here as H5–H7.

Hypothesis 5 (H5).

We can find

C_{P S} > 0

such that

{| (P u, u)}_{0, Ω} | \geq C_{P S} {‖ u ‖}_{0, Ω}^{2}

for all

u \in {(L^{2} (Ω))}^{3}

Hypothesis 6 (H6).

We can find

C_{Q S} > 0

such that

{| (Q curl u, curl u)}_{0, Ω} | \geq C_{Q S} {‖ curl u ‖}_{0, Ω}^{2}

for all

u \in H (curl, Ω)

Hypothesis 7 (H7).

C_{P S}

C_{Q S}

C_{L}

, and

C_{M}

(i.e., all media involved) are such that

C_{Q S} - \frac{C_{L} C_{M}}{C_{P S}} > 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_{s s}

is uniformly positive definite in

Ω_{e l} \subset Ω

, that is

\exists C_{1} > 0

such that:

\int_{Ω_{e l}} u^{*} P_{s s} u \geq C_{1} \int_{Ω e l} {| u |}^{2} = C_{1} {| | u | |}_{0, Ω_{e l}}^{2} \forall u \in {(L^{2} (Ω))}^{3} .

(11)

whenever

Ω_{e l} = Ω

, we can simply define

C_{P S} = C_{1}

Whenever

Ω_{e l}

is not the whole Ω, suppose that, in the complementary region,

P_{s}

is uniformly positive or negative definite, that is

\exists C_{5} > 0

such that:

|\int_{Ω \ Ω_{e l}} u^{*} P_{s} u| \geq C_{5} {| | u | |}_{0, Ω \ Ω_{e l}}^{2} .

(12)

whenever

Ω_{e l} = \emptyset

, we simply have

C_{P S} = C_{5}

, and we can set:

C_{P S} = min_{i \in I} inf_{x \in Ω_{i}} λ_{m i n} (P_{s}),

(13)

where

λ_{m i n}

denotes the minimum of the magnitudes of the eigenvalues of the Hermitian symmetric matrix

P_{s}

Finally, whenever

Ω_{e l}

is neither the empty set nor the whole domain, under assumptions HM2 and HM3 of [9], condition H5 is satisfied with

C_{P S}

given by:

C_{P S} = \frac{1}{\sqrt{2}} min (\sqrt{(1 - α)} C_{5}, \sqrt{C_{1}^{2} + (1 - \frac{1}{α}) C_{3}^{2}}),

(14)

where

C_{3} > 0

is defined by:

|\int_{Ω_{e l}} u^{*} P_{s} u| \leq C_{3} {‖ u ‖}_{0, Ω_{e l}}^{2}

(15)

and α is such that

1 > α > \frac{C_{3}^{2}}{C_{1}^{2} + C_{3}^{2}} > 0

Analogously, by replacing P with Q in Equations (11), (12), and (15), we define, respectively,

Ω_{m l}

and the constants

C_{2} > 0

C_{4} > 0

, and

C_{6} > 0

and deduce that condition H6 is satisfied if we set:

C_{Q S} = min_{i \in I} inf_{x \in Ω_{i}} λ_{m i n} (Q_{s}),

(16)

whenever

Ω_{m l} = \emptyset

C_{Q S} = C_{2}

whenever

Ω_{m l} = Ω

, or:

C_{Q S} = \frac{1}{\sqrt{2}} min (\sqrt{(1 - α)} C_{6}, \sqrt{C_{2}^{2} + (1 - \frac{1}{α}) C_{4}^{2}}),

(17)

α

being such that

1 > α > \frac{C_{4}^{2}}{C_{2}^{2} + C_{4}^{2}} > 0

, when

Ω_{m l} \neq Ω

and

Ω_{m l} \neq \emptyset

With respect to the H4, if we define:

C_{χ, s} = max_{i \in I_{b}} sup_{x \in Ω_{i}} ((\sum_{i, j = 1}^{3} | χ_{i j} |) - \min_{i = 1, 2, 3} | χ_{i i} |),

(18)

C_{γ, s} = max_{i \in I_{b}} sup_{x \in Ω_{i}} ((\sum_{i, j = 1}^{3} | γ_{i j} |) - \min_{i = 1, 2, 3} | γ_{i i} |),

(19)

the sufficient condition guaranteeing the validity of the inequality in the hypothesis can be expressed as:

K_{u} = \frac{4 C_{χ, s} C_{γ, s}}{(- C_{κ, s} + \sqrt{C_{κ, s}^{2} + 4 C_{κ, d} C_{κ, r}}) (- C_{ν, s} + \sqrt{C_{ν, s}^{2} + 4 C_{ν, d} C_{ν, r}})} < 1 .

(20)

3. 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,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].

3.1. 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:

\{\begin{cases} D = ε E + ξ H, \\ B = ζ E + μ H, \end{cases}

(21)

with:

ε = ε_{0} [\begin{matrix} ε_{x} & 0 & 0 \\ 0 & ε_{y} & 0 \\ 0 & 0 & ε_{z} \end{matrix}], μ = μ_{0} [\begin{matrix} μ_{x} & 0 & 0 \\ 0 & μ_{y} & 0 \\ 0 & 0 & μ_{z} \end{matrix}],

(22)

ξ = \frac{1}{c_{0}} [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & j ξ_{0} & 0 \end{matrix}], ζ = \frac{1}{c_{0}} [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & - j ξ_{0} \\ 0 & 0 & 0 \end{matrix}] .

(23)

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} \subset Ω

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

P = c_{0} ε_{0} [\begin{matrix} ε_{r} & 0 & 0 \\ 0 & ε_{r} & 0 \\ 0 & 0 & ε_{r} - ξ_{0}^{2} \end{matrix}],

(24)

Q = \frac{1}{c_{0} μ_{0}} I_{3},

(25)

L = M^{T} = \frac{j ξ_{0}}{μ_{0} c_{0}} [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \end{matrix}] .

(26)

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 = \frac{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

I m (ε_{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} - j P_{s s}

with:

P_{s} = \frac{P + P^{*}}{2} = c_{0} ε_{0} [\begin{matrix} R e (ε_{r}) & 0 & 0 \\ 0 & R e (ε_{r}) & 0 \\ 0 & 0 & R e (ε_{r}) - ξ_{0}^{2} \end{matrix}],

(27)

and

P_{s s} = \frac{P^{*} - P}{2 j} = - c_{0} ε_{0} I m (ε_{r}) I_{3}

. Hence, we have

Ω_{e l} = Ω_{s}

, the lossy region where

P_{s s}

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

C_{1} = c_{0} ε_{0} | I m (ε_{r}) |

C_{5} = c_{0} ε_{0}

, and

C_{3} = c_{0} ε_{0} max (| R e (ε_{r}) - ξ_{0}^{2} |, | R e (ε_{r}) |)

. 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_{P S}

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

α_{o p t}

, is given by:

α_{o p t} = \frac{C_{5}^{2} - C_{1}^{2} - C_{3}^{2} + \sqrt{{(C_{5}^{2} - C_{1}^{2} - C_{3}^{2})}^{2} + 4 C_{5}^{2} C_{3}^{2}}}{2 C_{5}^{2}} .

(28)

Thus, we may simply write:

C_{P S} = \sqrt{\frac{1 - α_{o p t}}{2}} c_{0} ε_{0} .

(29)

As mentioned in [9], this does not mean that a better value of

C_{P S}

cannot be found. For example, if

R e (ε_{r}) - ξ_{0}^{2} > 0

, then

P_{s}

is uniformly positive definite in

Ω

, and we can find another candidate for

C_{P S}

, namely

C_{7} = c_{0} ε_{0} min (1, | R e (ε_{r}) - ξ_{0}^{2} |)

(see Equation (31) of [9]). In particular, when

R e (ε_{r}) - ξ_{0}^{2} > \frac{1}{\sqrt{2}}

C_{7}

is always going to give a value for

C_{P S}

that is higher than that obtained from the lemma.

In the rest of the subsection, the discussion focuses on the cases with

R e (ε_{r}) < 0

, which gives a non-definite

P_{s}

Ω

. For this case,

C_{3} = c_{0} ε_{0} | R e (ε_{r}) - ξ_{0}^{2} |

, and we can directly use the value of

C_{P S}

in Equation (29). Since the material is assumed to be non-magnetic, the direct application of the definition gives

C_{Q S} = \frac{1}{c_{0} μ_{0}}

, and H6 is valid. Likewise, for the continuity constants

C_{L}

and

C_{M}

, we can choose the value

\frac{| ξ_{0} |}{c_{0} μ_{0}}

. Then, the inequality in H7 becomes

C_{Q S} - \frac{C_{L} C_{M}}{C_{P S}} = \frac{1}{c_{0} μ_{0}} (1 - \frac{ξ_{0}^{2} \sqrt{2}}{\sqrt{1 - α_{o p t}}}) > 0

, which gives:

| ξ_{0} | < {(\frac{1 - α_{o p t}}{2})}^{1 / 4} .

(30)

Since the right-hand side of Equation (30) also depends on

ξ_{0}

due to the presence of

C_{3}

in the expression for

α_{o p t}

, 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

R e (ε_{r})

is varied in the range

(- 5.0, - 1.0)

, whereas

I m (ε_{r})

assumes values in the range

(- 0.5, - 0.1)

. It can be observed that for a fixed value of

I m (ε_{r})

, the range of

ξ_{0}

over which H7 is valid steadily decreases as

| R e (ε_{r}) |

increases. As for the dependence on

I m (ε_{r})

, the corresponding range increases when the medium becomes lossier due to higher

| I m (ε_{r}) |

, as expected.

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

κ = \frac{1}{ε_{0} ε_{r}} [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \frac{ε_{r}}{ε_{r} - ξ_{0}^{2}} \end{matrix}],

(31)

ν = \frac{1}{μ_{0}} [\begin{matrix} 1 & 0 & 0 \\ 0 & \frac{ε_{r}}{ε_{r} - ξ_{0}^{2}} & 0 \\ 0 & 0 & 1 \end{matrix}],

(32)

γ = - χ^{T} = \frac{j ξ_{0} c_{0}}{ε_{r} - ξ_{0}^{2}} [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{matrix}] .

(33)

The constants of interest can be evaluated directly from the definitions. The determinants of

κ

and

ν

are, respectively,

\frac{1}{{(ε_{0} ε_{r})}^{3} (1 - \frac{ξ_{0}^{2}}{ε_{r}})}

and

\frac{1}{{(μ_{0})}^{3} (1 - \frac{ξ_{0}^{2}}{ε_{r}})}

, which immediately give the values of

C_{κ, d}

and

C_{ν, d}

C_{κ, d} = \frac{1}{| ε_{0}^{3} ε_{r}^{3} (1 - \frac{ξ_{0}^{2}}{ε_{r}}) |},

(34)

C_{ν, d} = \frac{1}{| μ_{0}^{3} (1 - \frac{ξ_{0}^{2}}{ε_{r}}) |} .

(35)

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}

C_{κ, r} = | ε_{0} ε_{r} |,

(36)

C_{ν, r} = μ_{0} .

(37)

Using Equations (36) and (37) of [9], we get

C_{κ, s}

and

C_{ν, s}

C_{κ, s} = \frac{2}{| ε_{0} ε_{r} |},

(38)

C_{ν, s} = \frac{2}{μ_{0}} .

(39)

From Equations (18) and (19), we can easily evaluate

C_{γ, s}

and

C_{χ, s}

C_{χ, s} = C_{γ, s} = |\frac{ξ_{0} c_{0}}{ε_{r} - ξ_{0}^{2}}| .

(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

R e (ε_{r})

I m (ε_{r})

. The results show that the range of

ξ_{0}

for which H4 holds true increases with the increase in

| R e (ε_{r}) |

, while it is practically independent of

I m (ε_{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

μ

Inside

Ω_{s}

, the medium is characterized by

ε_{r} = - 1 - j 0.4

μ_{r} = 1

, and

ξ_{0} = - 0.41

. This value is such that the hypotheses required for well posedness and finite element approximability are satisfied. In fact, for the

ε_{r}

considered, condition H4 is valid for

| ξ_{0} | < 0.4393

, and condition H7 is valid for

| ξ_{0} | < 0.4235

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

\sqrt{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.

Figure 4 shows the significance of the bianisotropic effect on the z component of the electric field along a line parallel to the x axis and passing through the center of gravity of the domain. 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.

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

\{\begin{cases} D = ε_{0} ε_{r} I_{3} E - j ξ_{c} I_{3} B, \\ H = - j ξ_{c} I_{3} E + \frac{1}{μ_{0} μ_{r}} I_{3} B . \end{cases}

(41)

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:

P = ε_{0} ε_{r} c_{0} I_{3},

(42)

Q = \frac{1}{μ_{0} μ_{r} c_{0}} I_{3},

(43)

L = M = - j ξ_{c} I_{3} .

(44)

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} \in 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} \geq 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

Ω_{e l} = \emptyset

, a value of

C_{P S}

can be found by using Equation (13). In particular, H5 is satisfied for

C_{P S} = ε_{0} c_{0}

. The hypothesis H6 is also trivially valid with

C_{Q S} = \frac{1}{μ_{0} c_{0}}

. By Equations (32) and (33) of [9],

C_{L} = C_{M} = ξ_{c}

. Then, the inequality in hypothesis H7 becomes

C_{Q S} - \frac{C_{L} C_{M}}{C_{P S}} = c_{0} (ε_{0} - μ_{0} ξ_{c}^{2}) > 0

, which implies:

ξ_{c} < \sqrt{\frac{ε_{0}}{μ_{0}}} = 2.654 \times 10^{- 3} mho .

(45)

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

κ = \frac{1}{ε_{0} ε_{r}} I_{3},

(46)

ν = \frac{ε_{0} ε_{r} + μ_{0} ξ_{c}^{2}}{μ_{0} ε_{0} ε_{r}} I_{3},

(47)

χ = - γ = \frac{j ξ_{c}}{ε_{0} ε_{r}} I_{3} .

(48)

κ

and

ν

are multiples of the identity matrix with eigenvalues

\frac{1}{ε_{0} ε_{r}}

and

(\frac{ε_{0} ε_{r} + μ_{0} ξ_{c}^{2}}{μ_{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}

C_{κ, d} = {(\frac{1}{ε_{0} ε_{r}})}^{3},

(49)

C_{ν, d} = {(\frac{ε_{0} ε_{r} + μ_{0} ξ_{c}^{2}}{μ_{0} ε_{0} ε_{r}})}^{3} .

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

C_{κ, s} = \frac{2}{ε_{0} ε_{r}},

(51)

C_{ν, s} = 2 \frac{ε_{0} ε_{r} + μ_{0} ξ_{c}^{2}}{μ_{0} ε_{0} ε_{r}} .

(52)

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}

C_{κ, r} = ε_{0} ε_{r},

(53)

C_{ν, r} = \frac{μ_{0} ε_{0} ε_{r}}{ε_{0} ε_{r} + μ_{0} ξ_{c}^{2}} .

(54)

From Equations (18) and (19), we get:

C_{χ, s} = C_{γ, s} = \frac{2 ξ_{c}}{ε_{0} ε_{r}} .

(55)

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}

ξ_{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 \times 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 \times 10^{- 4}

mho for

ε_{r} = 1

, increases with

ε_{r}

, and is

1.78 \times 10^{- 3}

mho for

ε_{r} = 10

. The curve crosses the

2.654 \times 10^{- 3}

mho line at around

ε_{r} = 22.3

, and above that value, Equation (45) imposes the stricter limit.

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

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

T E_{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

T E_{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

\frac{1}{2}

mm,

\frac{1}{4}

mm, and

\frac{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.

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.

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.

3.3. 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:

P = ε_{0} c_{0} I_{3},

(56)

Q = \frac{1}{μ_{0} c_{0}} I_{3},

(57)

L = M = j κ_{0} A,

(58)

where A is the matrix given by:

A = [\begin{matrix} 1 & 1 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}],

(59)

and

κ_{0}

is a positive real number.

The hypotheses H5 and H6 are trivially valid with

C_{P S} = ε_{0} c_{0}

and

C_{Q S} = \frac{1}{μ_{0} c_{0}}

. Further,

L^{*} L = M^{*} M = κ_{0}^{2} A^{2}

whose eigenvalues are zero,

κ_{0}^{2}

, and

4 κ_{0}^{2}

. Therefore, by Equations (32) and (33) of [9], we get

C_{L} = C_{M} = 2 κ_{0}

. The condition in hypothesis H7 then becomes

C_{Q S} - \frac{C_{L} C_{M}}{C_{P S}} = \frac{1}{μ_{0} c_{0}} - \frac{4 κ_{0}^{2}}{ε_{0} c_{0}} > 0

, which gives the following limit on

κ_{0}

κ_{0} < \frac{1}{2} \sqrt{\frac{ε_{0}}{μ_{0}}} = 1.327 \times 10^{- 3} mho .

(60)

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

κ = \frac{1}{ε_{0}} I_{3},

(61)

ν = \frac{1}{μ_{0}} I_{3} + \frac{κ_{0}^{2}}{ε_{0}} A^{2},

(62)

χ = - γ = \frac{- j κ_{0}}{ε_{0}} A .

(63)

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_{κ, d} = \frac{1}{ε_{0}^{3}},

(64)

C_{ν, d} = \frac{(ε_{0} + μ_{0} κ_{0}^{2}) (ε_{0} + 4 μ_{0} κ_{0}^{2})}{μ_{0}^{3} ε_{0}^{2}} .

(65)

C_{κ, s}

and

C_{ν, s}

can be directly obtained from their definitions:

C_{κ, s} = \frac{2}{ε_{0}},

(66)

C_{ν, s} = \frac{2 ε_{0} + 8 μ_{0} κ_{0}^{2}}{μ_{0} ε_{0}} .

(67)

By simple application of the definition:

C_{κ, r} = ε_{0}

(68)

and by Equation (41) of [9],

C_{ν, r}

evaluates to the reciprocal of the largest eigenvalue of the real matrix

ν

, which is given by:

C_{ν, r} = \frac{μ_{0} ε_{0}}{ε_{0} + 4 μ_{0} κ_{0}^{2}} .

(69)

Equations (18) and (19) give:

C_{χ, s} = C_{γ, s} = \frac{4 κ_{0}}{ε_{0}} .

(70)

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} \leq 2.72 \times 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 \times 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.

T E_{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,

\frac{1}{2}

mm, and

\frac{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.

We provide the magnitudes and phases of the x component of the electric field obtained from the simulation in Figure 15, Figure 16 and Figure 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.

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

Author Contributions

Both authors have contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

Kraft, M.; Braun, A.; Luo, Y.; Maier, S.A.; Pendry, J.B. Bianisotropy and magnetism in plasmonic gratings. ACS Photonics 2016, 3, 764–769. [Google Scholar] [CrossRef] [Green Version]
Yazdi, M.; Albooyeh, M.; Alaee, R.; Asadchy, V.; Komjani, N.; Rockstuhl, C.; Simovski, C.R.; Tretyakov, S. A bianisotropic metasurface with resonant asymmetric absorption. IEEE Trans. Antennas Propag. 2015, 63, 3004–3015. [Google Scholar] [CrossRef]
Kildishev, A.V.; Borneman, J.D.; Ni, X.; Shalaev, V.M.; Drachev, V.P. Bianisotropic effective parameters of optical metamagnetics and negative-index materials. Proc. IEEE 2011, 99, 1691–1700. [Google Scholar] [CrossRef]
Kriegler, C.E.; Rill, M.S.; Linden, S.; Wegener, M. Bianisotropic photonic metamaterials. IEEE J. Sel. Top. Quantum Electron. 2010, 16, 367–375. [Google Scholar] [CrossRef]
Fernandes, P.; Raffetto, M. Well-posedness and finite element approximability of time-harmonic electromagnetic boundary value problems involving bianisotropic materials and metamaterials. Math. Model. Methods Appl. Sci. 2009, 19, 2299–2335. [Google Scholar] [CrossRef]
Brignone, M.; Raffetto, M. Well posedness and finite element approximability of two-dimensional time-harmonic electromagnetic problems involving non-conducting moving objects with stationary boundaries. ESAIM Math. Model. Numer. Anal. 2015, 49, 1157–1192. [Google Scholar] [CrossRef]
Ioannidis, A.D.; Kristensson, G.; Stratis, I.G. On the well-posedness of the Maxwell system for linear bianisotropic media. SIAM J. Math. Anal. 2012, 44, 2459–2473. [Google Scholar] [CrossRef] [Green Version]
Cocquet, P.; Mazet, P.; Mouysset, V. On the existence and uniqueness of a solution for some frequency-dependent partial differential equations coming from the modeling of metamaterials. SIAM J. Math. Anal. 2012, 44, 3806–3833. [Google Scholar] [CrossRef]
Kalarickel Ramakrishnan, P.; Raffetto, M. Well posedness and finite element approximability of three-dimensional time-harmonic electromagnetic problems involving rotating axisymmetric objects. Symmetry 2020, 12, 218. [Google Scholar] [CrossRef] [Green Version]
Wu, T.X.; Jaggard, D.L. A comprehensive study of discontinuities in chirowaveguides. IEEE Trans. Microw. Theory Tech. 2002, 50, 2320–2330. [Google Scholar] [CrossRef]
Alotto, P.; Codecasa, L. A fit formulation of bianisotropic materials over polyhedral grids. IEEE Trans. Magn. 2014, 50, 349–352. [Google Scholar] [CrossRef]
Monk, P. Finite Element Methods for Maxwell’s Equations; Oxford Science Publications: Oxford, UK, 2003. [Google Scholar]
Fernandes, P.; Ottonello, M.; Raffetto, M. Regularity of time-harmonic electromagnetic fields in the interior of bianisotropic materials and metamaterials. IMA J. Appl. Math. 2014, 79, 54–93. [Google Scholar] [CrossRef]
Ciarlet, P.G.; Lions, J.L. Handbook of Numerical Analysis; Finite Element Methods, Part 1; North-Holland: Amsterdam, The Netherlands, 1991; Volume II. [Google Scholar]

Figure 1. The plot indicates the maximum of

| ξ_{0} |

guaranteeing that condition H7 is satisfied, for scattering problems involving different media considered in [1]. The curves are computed by assuming

ξ_{0} \in R

R e (ε_{r}) \in (- 5.0, - 1.0)

I m (ε_{r}) \in (- 0.5, - 0.1)

, and

μ_{r} = 1

Figure 1. The plot indicates the maximum of

| ξ_{0} |

guaranteeing that condition H7 is satisfied, for scattering problems involving different media considered in [1]. The curves are computed by assuming

ξ_{0} \in R

R e (ε_{r}) \in (- 5.0, - 1.0)

I m (ε_{r}) \in (- 0.5, - 0.1)

, and

μ_{r} = 1

Figure 2. The plots indicate the maximum

| ξ_{0} |

guaranteeing that condition H4 is satisfied, for scattering problems involving different media considered in [1]. The curves are computed by assuming

ξ_{0} \in R

R e (ε_{r}) \in (- 5.0, - 1.0)

I m (ε_{r}) \in (- 0.5, - 0.1)

, and

μ_{r} = 1

Figure 2. The plots indicate the maximum

| ξ_{0} |

guaranteeing that condition H4 is satisfied, for scattering problems involving different media considered in [1]. The curves are computed by assuming

ξ_{0} \in R

R e (ε_{r}) \in (- 5.0, - 1.0)

I m (ε_{r}) \in (- 0.5, - 0.1)

, and

μ_{r} = 1

Figure 3. Stability of the solution for problem involving the medium in [1]. The magnitude of the z component of the electric field is plotted for four different meshes along a line parallel to the y axis and passing through the center of gravity of the domain.

Figure 4. The magnitude of the z component of electric field along a line parallel to the x axis and passing through the center of gravity of the domain, for the problem involving the medium in [1]. The plot for the magnitude of the field

| E_{z} |

obtained in the bianisotropic case using

ξ_{0} = - 0.41

is shown along with the magnitude of the difference between the two solutions

| E_{z} - E_{z, 0} |

, where

E_{z, 0}

is obtained using

ξ_{0} = 0

| E_{z} |

obtained in the bianisotropic case using

ξ_{0} = - 0.41

is shown along with the magnitude of the difference between the two solutions

| E_{z} - E_{z, 0} |

, where

E_{z, 0}

is obtained using

ξ_{0} = 0

Figure 5. The magnitude of the z component of electric field along a line parallel to the z axis and passing through the center of gravity of the domain, for the problem involving the medium in [1]. The plot for the magnitude of the field

| E_{z} |

obtained in the bianisotropic case using

ξ_{0} = - 0.41

is shown along with the magnitude of the difference between the two solutions

| E_{z} - E_{z, 0} |

, where

E_{z, 0}

is obtained using

ξ_{0} = 0

| E_{z} |

obtained in the bianisotropic case using

ξ_{0} = - 0.41

is shown along with the magnitude of the difference between the two solutions

| E_{z} - E_{z, 0} |

, where

E_{z, 0}

is obtained using

ξ_{0} = 0

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

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

Figure 7. The value of

ξ_{c}

below which the hypothesis H4 is satisfied is plotted against

ε_{r}

. The limit of

2.654 \times 10^{- 3}

mho, arising from Equation (45) required to satisfy H7, is also shown.

Figure 7. The value of

ξ_{c}

below which the hypothesis H4 is satisfied is plotted against

ε_{r}

. The limit of

2.654 \times 10^{- 3}

mho, arising from Equation (45) required to satisfy H7, is also shown.

Figure 8. The geometry of a rectangular waveguide partially filled with the chiral media considered in [10].

Figure 9. Stability of the solution for the problem involving the medium in [10]. 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.

Figure 10. The magnitude and phase of the x component of the electric field along a line parallel to the x axis and passing though the center of gravity of the domain for the problem involving the medium in [10]. The plot for the bianisotropic case using

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

Figure 11. The magnitude and phase of the x component of the electric field along a line parallel to the z axis and passing though the center of gravity of the domain for the problem involving the medium in [10]. The plot for the bianisotropic case using

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

Figure 12. The magnitude and phase of the y component of the electric field along a line parallel to the z axis and passing though the center of gravity of the domain for the problem involving the medium in [10]. The plot for the bianisotropic case using

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

Figure 13. The magnitude and phase of the z component of the electric field along a line parallel to the x axis and passing though the center of gravity of the domain for problem involving the medium in [10]. The plot for the bianisotropic case using

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

ξ_{c} = 1.24 \times 10^{- 3}

mho is compared with the solution obtained in the isotropic case using

ξ_{c} = 0

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.

Figure 15. The magnitude and phase of the x component of the electric field along a line parallel to the x axis and passing though the center of gravity of the domain for the problem involving the medium in [11]. The plot for the bianisotropic case using

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

Figure 16. The magnitude and phase of the x component of the electric field along a line parallel to the y axis and passing though the center of gravity of the domain for the problem involving the medium in [11]. The plot for the bianisotropic case using

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

Figure 17. The magnitude and phase of the x component of the electric field along a line parallel to the z axis and passing though the center of gravity of the domain for the problem involving the medium in [11]. The plot for the bianisotropic case using

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

κ_{0} = 2.7 \times 10^{- 4}

mho is compared with the solution obtained in the isotropic case using

κ_{0} = 0

Table 1. Details of the different meshes used.

Type of Mesh	Maximum Diameter of the Mesh (h in nm)	Number of Nodes	Number of Elements	Number of Boundary Faces
Very coarse	200 $\sqrt{3}$	1331	6000	1200
Coarse	100 $\sqrt{3}$	9261	48,000	4800
Fine	50 $\sqrt{3}$	68,921	384,000	19,200
Very fine	25 $\sqrt{3}$	531,441	3,072,000	76,800

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Kalarickel Ramakrishnan, P.; Raffetto, M. Three-Dimensional Time-Harmonic Electromagnetic Scattering Problems from Bianisotropic Materials and Metamaterials: Reference Solutions Provided by Converging Finite Element Approximations. Electronics 2020, 9, 1065. https://doi.org/10.3390/electronics9071065

AMA Style

Kalarickel Ramakrishnan P, Raffetto M. Three-Dimensional Time-Harmonic Electromagnetic Scattering Problems from Bianisotropic Materials and Metamaterials: Reference Solutions Provided by Converging Finite Element Approximations. Electronics. 2020; 9(7):1065. https://doi.org/10.3390/electronics9071065

Chicago/Turabian Style

Kalarickel Ramakrishnan, Praveen, and Mirco Raffetto. 2020. "Three-Dimensional Time-Harmonic Electromagnetic Scattering Problems from Bianisotropic Materials and Metamaterials: Reference Solutions Provided by Converging Finite Element Approximations" Electronics 9, no. 7: 1065. https://doi.org/10.3390/electronics9071065

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

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

Abstract

1. Introduction

2. Mathematical Description of the Problem

3. Results and Discussion

3.1. Scattering from Plasmonic Gratings Behaving as Bianisotropic Metamaterials

3.2. Scattering from Chiral Obstacles in a Waveguide

3.3. Reflection by a Short-Circuited Waveguide Half Filled with Bianisotropic Media

4. Conclusions

Author Contributions

Funding

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI