SINGULAR VALUE DECOMPOSITIONS FOR SINGLE-CURL OPERATORS IN THREE-DIMENSIONAL MAXWELL’S
EQUATIONS FOR COMPLEX MEDIA∗
RUEY-LIN CHERN†, HAN-EN HSIEH‡, TSUNG-MING HUANG§, WEN-WEI LIN¶,
AND WEICHUNG WANG†
Abstract. This article focuses on solving the generalized eigenvalue problems (GEP) arising in
the source-free Maxwell equation with magnetoelectric coupling effects that models three-dimensional complex media. The goal is to compute the smallest positive eigenvalues, and the main challenge is that the coefficient matrix in the discrete Maxwell equation is indefinite and degenerate. To overcome this difficulty, we derive a singular value decomposition (SVD) of the discrete single-curl operator and then explicitly express the basis of the invariant subspace corresponding to the nonzero eigenvalues of the GEP. Consequently, we reduce the GEP to a null space free standard eigenvalue problem (NFSEP) that contains only the nonzero (complex) eigenvalues of the GEP and can be solved by the shift-and-invert Arnoldi method without being disturbed by the null space. Furthermore, the basis of the eigendecomposition is chosen carefully so that we can apply fast Fourier transformation (FFT-) based matrix vector multiplication to solve the embedded linear systems efficiently by an iterative method. For chiral and pseudochiral complex media, which are of great interest in magnetoelectric applications, the NFSEP can be further transformed to a null space free GEP whose coefficient matrices are Hermitian and Hermitian positive definite (HHPD-NFGEP). This HHPD-NFGEP can be solved by using the invert Lanczos method without shifting. Furthermore, the embedded linear system can be solved efficiently by using the conjugate gradient method without preconditioning and the FFT-based matrix vector multiplications. Numerical results are presented to demonstrate the efficiency of the proposed methods.
Key words. singular value decomposition, null space free method, discrete single-curl operator,
the Maxwell equations, chiral medium, pseudochiral medium
AMS subject classifications. 65F15, 65T50, 15A18, 15A23 DOI. 10.1137/140958748
1. Introduction. Understanding the eigenstructure of the discrete single-curl
operator∇× is key to developing efficient numerical simulations for complex materials
that are modeled by the Maxwell equations. Bi-isotropic and bianisotropic materials are two important classes of complex materials [31]. They have consistently been the subject of intensive studies in physical properties and applications. For example, bianisotropic media are a special type of materials whose properties are characterized by the magnetoelectric as well as the permittivity and permeability tensors [18, 26]. Due to the strong modulation of the wave that arises from the magnetoelectric cou-plings, counterintuitive features such as negative refraction and backward waves may
∗Received by the editors February 26, 2014; accepted for publication (in revised form) by
D. O’Leary November 4, 2014; published electronically February 26, 2015. This research was par-tially supported by the National Science Council, the National Center for Theoretical Sciences, the Taida Institute for Mathematical Sciences, and the Chiao-Da ST Yau Center in Taiwan.
http://www.siam.org/journals/simax/36-1/95874.html
†Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan (chernrl@ntu.
edu.tw, [email protected]).
‡Department of Mathematics, National Taiwan University, Taipei 106, Taiwan (D99221002@ntu.
edu.tw).
§Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan (min@
ntnu.edu.tw).
¶Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan
203
appear in bianisotropic media [1, 5, 22, 29]. From the mathematical point of view, the distinctive feature associated with bianisotropic media is the single-curl operator in addition to the double-curl operator in the wave equation, which essentially changes the characters of the eigenwaves.
Mathematically, the propagation of electromagnetic waves in bi-isotropic and bianisotropic materials is modeled by the three-dimensional (3D) frequency domain source-free Maxwell equations [31] with a set of constitutive relations. In particular, we have
∇ × E = iωB, ∇ · B = 0, (1.1a)
∇ × H = −iωD, ∇ · D = 0, (1.1b)
whereω represents the frequency and ε represents the permittivity. E and H are the
electric and magnetic fields, respectively. Based on the Bloch theorem [17], we aim
to find the Bloch eigenfunctionsE and H satisfying the quasi-periodic conditions
E(x + a) =ei2πk·aE(x), H(x + a) =ei2πk·aH(x)
(1.2)
for = 1, 2, 3 [23]. Here, a1, a2, and a3 are the lattice translation vectors. In this
paper, we consider the simple cubic lattice vectors a=ae, where e is theth unit
vector in R3 and a is a lattice constant. Without loss of generality, we set a = 1
throughout this paper. Note that all the techniques developed here can be applied to face-centered cubic lattice media. The Bloch wave vector in the first Brillouin zone is
denoted as 2πk [15]. B and D satisfy the constitutive relations
(1.3) B = μH + ζE and D = εE + ξH,
whereμ represents the permeability, and ξ and ζ are magnetoelectric parameters [26,
p. 26], [31, p. 44]. Note thatε, μ, ξ, and ζ are 3-by-3 matrices in various forms for
describing different types of materials.
The Maxwell equations (1.1) can be rewritten as the following quadratic
eigen-value problems (QEP), which are separate wave equations in terms ofE and H:
∇ × μ−1∇ × E − iω∇ ×μ−1ζE− ξμ−1∇ × E− ω2ε − ξμ−1ζE = 0;
(1.4a)
∇ × ε−1∇ × H − iωζε−1∇ × H − ∇ ×ε−1ξH− ω2μ − ζε−1ξH = 0.
(1.4b)
In the one-dimensional case, we can apply the quasi-periodic conditions (1.2) to (1.4)
and then explicitly define the relations between k and ω [2, 3, 4, 6, 19]. For higher
dimensions, however, solving (1.2) efficiently remains an open question. We illustrate the difficulty by the following example. An explicit eigendecomposition of the discrete
double-curl operator∇×∇× is derived in [11]. Applying this eigendecomposition and
assumingμ = 1 and ζ = ξ = 0, we can explicitly derive the invariant subspace of all
nonzero eigenvalues corresponding to the (discrete) eigenvalue problem (1.4a). Based on the invariant subspace, efficient numerical methods can be developed to solve (1.4a). However, it is not possible to apply this technique to solve the QEPs (1.4)
with ζ = 0 and ξ = 0 due to the following difficulties: (i) The eigendecomposition
of the discrete double-curl operator in [11] cannot be applied to solving the QEP directly because (1.4) contains both double- and single-curl operators. (ii) In general, the double- and single-curl operators in (1.4) cannot be diagonalized simultaneously. Furthermore, should we find the eigendecomposition of the single-curl operator, this decomposition cannot be applied to solve the QEP directly because the single-curl
operator terms in (1.4), e.g.,∇×(μ−1ζE) and ξμ−1∇×E, are coupled with other terms
such asμ−1ζ and ξμ−1. (iii) It is difficult to find the invariant subspace corresponding
to the nonzero eigenvalues in the QEPs.
While solving (1.4) is not recommended, we focus instead on the original Maxwell equations (1.1) and rewrite it as a coupled generalized eigenvalue problem (GEP)
∇× 0 0 ∇× E H = iω ζ μ −ε −ξ E H . (1.5)
For the two-dimensional photonic band structure, the electromagnetic transfer matrix method [8] is applied to the coupled system, similar to (1.5). For the 3D case, to the best of our knowledge no method has yet been proposed to solve the GEP (1.5) efficiently.
We make the following contributions to solve the discrete 3D GEP based on Yee’s finite difference discretization scheme [32]:
• We first derive the singular value decomposition (SVD) of the discrete
single-curl operator∇× in (1.5).
• Using the SVD, we explore an explicit form of the basis for the invariant subspace corresponding to the nonzero eigenvalues of the GEP. Applying this basis, the GEP can be reduced to a null space free standard eigenvalue problem (NFSEP). In this eigenvalue problem, the zero eigenvalues of the GEP are deflated so that the null space does not degrade the computational efficiency.
• We show that all eigenvalues ω of the GEP are real provided the permittivity, permeability, and magnetoelectric parameters satisfy particular assumptions. These assumptions are applicable to a couple of important classes of complex media.
• Under the same assumptions, we can reformulate the NFSEP as a null space
free GEPBrx = ω−1Arx, where Bris Hermitian andAris Hermitian positive
definite. We demonstrate that this problem can be solved efficiently using the generalized Lanczos method algorithmically and numerically.
This paper is outlined as follows. In section 2, we derive the SVD of the discrete single-curl operator. In section 3, by applying the SVD, we derive a null space free eigenvalue problem by deflating the zero eigenvalues and keeping the nonzero eigen-value unchanged. In section 4, we discuss how to improve the solution performance while simulating two important types of complex media. In section 5, we demonstrate numerical results to validate the correctness of proposed schemes and to measure the performance of the schemes. Finally, we present our conclusions in section 6.
Throughout this paper, we use the superscripts and ∗ to denote the transpose
and the conjugate transpose of a matrix, respectively. For the matrix operations,
we let⊗ be the Kronecker product of two matrices. The imaginary number √−1 is
written as i, and the identity matrix of dimensionn is written as In.
2. SVD of the discrete single-curl operator. In this section, we derive an
explicit expression of the SVD of the discrete single-curl operator. Using this SVD, an efficient null space free method to solve the target eigenvalue problem (1.5) is developed in section 3.
We start from the derivation of the matrix representation of the discrete
single-curl operator. By using Yee’s scheme [32], the discrete single-single-curl operators ∇ × E
and∇ × H with a=ae, = 1, 2, 3, can be represented in the matrix form CE and C∗H [13, 14], respectively. Here, C = ⎡ ⎣ C03 −C03 −CC21 −C2 C1 0 ⎤ ⎦ ∈ C3n×3n (2.1) with C1 =δx−1(In3⊗ In2⊗ Ka1,n1)∈ Cn×n, (2.2a) C2 =δy−1(In3⊗ Ka2,n2⊗ In1)∈ Cn×n, (2.2b) C3 =δz−1(Ka3,n3⊗ In2⊗ In1)∈ Cn×n, (2.2c) and Ka,m= ⎡ ⎢ ⎢ ⎢ ⎣ −1 1 . .. . .. −1 1 ei2πk·a −1 ⎤ ⎥ ⎥ ⎥ ⎦∈ Cm×m. (2.3)
We use n1, n2, and n3 to denote the numbers of grid points in the x, y, and z
directions, respectively, and we define n = n1n2n3. We use δx, δy, and δz to denote
the associated mesh lengths along thex, y, and z axial directions, respectively.
It is worth noting that the divergence free conditions (1.1) and the quasi-periodicity conditions (1.2) are satisfied in Yee’s scheme [14, sections 5 and 3]. In addition to Yee’s scheme, there are other approaches to deal with the quasi-periodicity condi-tions. First, the curl operators can be replaced by the shifted operator, so that the transformation of the solution leads to a periodic solution [9, 24, 34]. Second, we can use periodic basis functions in finite element methods to resolve the quasi-periodicity conditions explicitly [25].
It is well known that the eigendecompositions ofC∗C and CC∗are closely related
to the SVD of C. Therefore, we start the derivation from the eigendecompositions
of C∗C and CC∗. First, we introduce the notation to be used later. Defineθm,i =
i2πi m , θa,m=i2πk·am , Da,m= diag 1, eθa,m, . . . , e(m−1)θa,m, (2.4)
um,i = 1 eθm,i . . . e(m−1)θm,i
fori = 0, . . . , m − 1 and
Um =
um,0 . . . um,m−1 ∈ Cm×m,
(2.5)
Λa,m = diageθm,0+θa,m− 1 . . . eθm,m−1+θa,m− 1.
(2.6)
By the definition ofKa,m in (2.3), it can be verified that
Ka,m(Da,mUm) = (Da,mUm) Λa,m.
(2.7) Denote
T = √1n(Da3,n3⊗ Da2,n2⊗ Da1,n1) (Un3⊗ Un2⊗ Un1), (2.8)
where Da,ni and Uni, i = 1, 2, 3, are given in (2.4) and (2.5), respectively. It is
straightforward to check thatT is unitary. The following theorem, which yields the
eigendecomposition ofC, for = 1, 2, 3, can be proved with (2.7) and the property
(A1⊗ A2⊗ A3)(B1⊗ B2⊗ B3) = (A1B1)⊗ (A2B2)⊗ (A3B3).
(2.9)
Theorem 2.1. C1, C2, and C3 can be diagonalized by the unitary matrix T in
the forms C1T = δ−1x T (In3⊗ In2⊗ Λa1,n1)≡ T Λ1, (2.10a) C2T = δ−1y T (In3⊗ Λa2,n2⊗ In1)≡ T Λ2, (2.10b) C3T = δ−1z T (Λa3,n3⊗ In2⊗ In1)≡ T Λ3. (2.10c)
We now define two intermediate matrices Λp and Λq that are used in the
eigen-decompositions ofC∗C and CC∗.
Lemma 2.2. Let Λ1, Λ2, and Λ3 be given in (2.10). Define
Λq = Λ∗1Λ1+ Λ∗2Λ2+ Λ∗3Λ3, Λp= ⎡ ⎣βΛΛ13− αΛ− Λ23 αΛ2− βΛ1 ⎤ ⎦ (2.11)
with α, β = 0. Assume that the vector k = (k1, k2, k3) in (1.2) is nonzero with
0 ≤ k1, k2, k3 ≤ 12. Then, Λq is positive definite, and Λp is of full column rank,
provided that αδx= βδy andδz= βδy.
Proof. See the appendix.
Using the definitions ofC and Λq in (2.1) and (2.11), respectively, and the
eigen-decompositions ofCin Theorem 2.1, the null spaces ofC∗C and CC∗ are derived as
follows.
Theorem 2.3. Assume k = (k1, k2, k3)= 0 with 0 ≤ k1, k2, k3≤1
2. Define Q0= (I3⊗ T ) ⎡ ⎣ΛΛ12 Λ3 ⎤ ⎦ Λ−1/2 q ≡ (I3⊗ T ) Π0, P0= (I3⊗ T ) Π0. (2.12)
Then, Q0 andP0 form orthogonal bases of the null spaces of C∗C and CC∗,
respec-tively.
Next, we apply the techniques developed in [11] to form the orthogonal bases for
the range spaces ofC∗C and CC∗. First, by considering the full column rank matrix
T1= [αT, βT, T] with nonzeroα and β, and taking the orthogonal projection
ofT1with respect to Q0and P0, respectively, we have
Q1 = (I − Q0Q∗0)T1Λ∗pΛpΛ−1q −1/2 = (I3⊗ T ) ⎡ ⎣(αΛ2− βΛ1)Λ ∗ 2− (Λ1− αΛ3)Λ∗3 (βΛ3− Λ2)Λ∗3− (αΛ2− βΛ1)Λ∗1 (Λ1− αΛ3)Λ∗1− (βΛ3− Λ2)Λ∗2 ⎤ ⎦ Λ∗ pΛpΛ−1q −1/2 ≡ (I3⊗ T ) Π1, (2.13a) P1 = (I − P0P0∗)T1Λ∗pΛpΛ−1q −1/2 = (I3⊗ T ) Π1. (2.13b)
Then, Q1 and P1 are orthogonal, and (C∗C)Q1 = Q1Λq and (CC∗)P1 = P1Λq.
Second, to form the remaining part of the orthogonal basis for the range spaces of
C∗C and CC∗, we apply the discrete curl and dual-curl operators onT
1, respectively.
That is, we premultiplyC∗ andC by T1 to obtain
Q2 =C∗T1Λ∗pΛp−1/2= (I3⊗ T ) ⎡ ⎣βΛ ∗ 3− Λ∗2 Λ∗1− αΛ∗3 αΛ∗ 2− βΛ∗1 ⎤ ⎦ Λ∗ pΛp −1/2 ≡ (I3⊗ T ) Π2, (2.14a) P2 =CT1Λ∗pΛp−1/2 = (I3⊗ T )−Π2. (2.14b)
It holds that (C∗C)Q2 = Q2Λq and (CC∗)P2 = P2Λq. From (2.12) to (2.14), we
define
Q ≡Q1 Q2 Q0= (I3⊗ T )Π1 Π2 Π0,
(2.15a)
P ≡P2 P1 P0= (I3⊗ T )−Π2 Π1 Π0.
(2.15b)
Then, the eigendecompositions ofC∗C and CC∗ can be summarized as follows.
Theorem 2.4. Assume that the vector k = (k1, k2, k3) in (1.2) is nonzero with
0≤ k1, k2, k3≤12. Then,Q and P are unitary, and
C∗C = Q diag (Λ
q, Λq, 0) Q∗, CC∗=P diag (Λq, Λq, 0) P∗.
(2.16)
Motivated by (2.16), we derive the left and right singular vector matricesP and Q
forC in the following theorem.
Theorem 2.5 (SVD ofC). Let Λq and (Q, P ) be defined in (2.11) and (2.15).
Assume that the vector k = (k1, k2, k3) in (1.2) is nonzero with 0≤ k1, k2, k3≤ 12.
Then, the matrixC in (2.1) has the SVD
C = P diagΛ1/2q , Λ1/2q , 0 Q∗=P rΣrQ∗r, (2.17) where Pr= [P2, P1], Qr= [Q1, Q2], Σr= diag Λ1/2q , Λ1/2q .
Proof. From (2.10) and the definition ofC, it follows that
CT1=− (I3⊗ T ) Λp, C∗T1= (I3⊗ T ) Λp, CQ0= 0, and P0∗C = 0.
(2.18)
Combining (2.18) with (2.13) and (2.14), we have
P∗ 2CQ1=Λ∗pΛp −1/2 T∗ 1C∗C (I − Q0Q∗0)T1Λ∗pΛpΛ−1q −1/2 =Λ∗pΛp−1/2T1∗C∗CT1Λ∗pΛpΛ−1q −1/2 = Λ1/2q , P∗ 1CQ1=Λ∗pΛpΛq−1−1/2T1∗(I − P0P0∗)C (I − Q0Q∗0)T1Λ∗pΛpΛ−1q −1/2 =Λ∗pΛpΛ−1q −1/2T1∗CT1Λ∗pΛpΛ−1q −1/2 =−Λ∗pΛpΛ−1q −1/2T1∗(I3⊗ T ) ΛpΛ∗pΛpΛ−1q −1/2 =−Λ∗pΛpΛ−1q −1/2αI βI IΛp Λ∗pΛpΛ−1q −1/2= 0, P∗ 1CQ2=Λ∗pΛpΛ−1q −1/2T1∗(I − P0P0∗)CC∗T1Λ∗pΛp−1/2 =Λ∗pΛpΛ−1q −1/2T1∗CC∗T1Λ∗pΛp−1/2 = Λ1/2q , P∗ 2CQ2=Λ∗pΛp −1/2T∗ 1C∗CQ2=Λ∗pΛp −1/2T∗ 1Q2Λq =Λ∗pΛp−1/2T1∗C∗T1Λ∗pΛp−1/2Λq = 0.
Therefore, we obtain the SVD ofC described in (2.17).
Now, we consider the case that k = 0 in Lemma 2.6 and Theorem 2.7. Note that,
in this case, the GEP (1.5) hasn + 2 zero eigenvalues [13].
Lemma 2.6. If k = 0, then Λq and Λp in (2.11) have rank n − 1. Furthermore,
Λq(1, 1) = 0 and the first column of Λp is a zero vector.
Proof. If k = 0, it holds that Λa,m(1, 1) = 0 and Λa,m(i, i) = 0 for i = 1 by (2.6).
This fact implies that Λ1(1, 1) = Λ2(1, 1) = Λ3(1, 1) = 0 in (2.10). Therefore, by the
definitions of Λq and Λp in (2.11), Λq and Λp have rankn − 1 and the first columns
of Λq and Λp are zero vectors.
We define two notations that are used in the following theorem. For a given
matrixF ∈ Cn×n, let F
c ∈ Cn×(n−1) (or Frc∈ C(n−1)×(n−1)) be the submatrices of
F in which the first column is deleted (or both the first column and the first row are
deleted). Combining the results in Theorem 2.5 and Lemma 2.6, the SVD ofC with
k = 0 can be obtained in Theorem 2.7.
Theorem 2.7. Let T , (Λ1, Λ2, Λ3), and (Λq, Λp) be defined in (2.8), (2.10), and
(2.11), respectively. Assume k = 0. Then C = PrΣrQ∗r, where Qr = (I3⊗ T ) Π1,c Π2,c, Pr= (I3⊗ T )−Π2,c Π1,c , Σr = diag (Λq)1/2rc , (Λq)1/2rc with Π1,c = ⎡ ⎣((αΛ2− βΛ1)Λ ∗ 2− (Λ1− αΛ3)Λ∗3)c ((βΛ3− Λ2)Λ∗3− (αΛ2− βΛ1)Λ∗1)c ((Λ1− αΛ3)Λ∗1− (βΛ3− Λ2)Λ∗2)c ⎤ ⎦ (Λ∗ pΛp)rc(Λq)−1rc −1/2, Π2,c = ⎡ ⎣(βΛ ∗ 3− Λ∗2)c (Λ∗1− αΛ∗3)c (αΛ∗2− βΛ∗1)c ⎤ ⎦ (Λ∗ pΛp)rc −1/2 .
It is worth mentioning the specific choice of the singular vector matricesP and
Q in (2.15). The choice of these matrices is not unique because the multiplicities of
the nonzero eigenvalues ofC∗C are even and may be large. We have carefully chosen
the P and Q defined in (2.15) to avoid the need to store these two matrices and
the computations involving P and Q can be performed efficiently. We discuss this
computational advantage in section 4.2.
3. The null space free eigenvalue problem. We have derived the SVD ofC
in Theorem 2.5. In this section, we use the SVD to deflate the null space of the GEP obtained by discretizing (1.5) so that we can develop a new solution process for the target eigenvalue problem. The discretization of (1.5) based on Yee’s scheme leads to the following discrete GEP:
C 0 0 C∗ E H =ω i ζd μd −εd −ξd E H . (3.1)
Here,C is the discrete single-curl operator defined in (2.1). The four 3n×3n complex
matrices ζd, ξd, εd, μd are the discrete counterparts of the matrices ζ, ξ, ε, and μ,
respectively.
From Theorem 2.5, we can see that the GEP (3.1) has 2n zero eigenvalues. This null space not only affects the convergence of iterative eigensolvers but also increases the challenge of solving the eigenvalue problem. In this section, we apply the SVD
of C in (2.17) to reduce the GEP (3.1) to the null space free eigenvalue problem
(3.6) equipped with the following two advantages: (i) The dimension of the coefficient
matrix is dramatically reduced from 6n in (3.1) to 4n in (3.6). (ii) The two eigenvalue
problems share the same 4n nonzero eigenvalues.
To derive the NFSEP, we first rewrite (3.1) as
diag (Pr, Qr) diag (Σr, Σr) diag (Q∗r, Pr∗)
E H =ω i ζd μd −εd −ξd E H (3.2)
by applying the SVD of C described in (2.17). We can then explicitly define the
invariant subspace of (3.1) in the following theorem using the matricesPr, Qr, and
Σrdefined in Theorem 2.5.
Theorem 3.1. Assume the matrix
B ≡ i ζd μd −εd −ξd (3.3) is nonsingular. Then span B−1diagP rΣ 1 2 r, QrΣ 1 2 r
is an invariant subspace of (3.1) corresponding to all nonzero eigenvalues. Further-more, it holds that
ω diagΣ12 rQ∗r, Σ 1 2 rPr∗ B−1diagP rΣ 1 2 r, QrΣ 1 2 r y = ωy ={ω | diag (C, C∗)x = ωBx, ω = 0} . (3.4)
Proof. From Theorem 2.5, we have
diag (C, C∗) B−1diagP rΣ 1 2 r, QrΣ 1 2 r (3.5) = {diag (PrΣrQ∗r, QrΣrPr∗)} B−1diagP rΣ 1 2 r, QrΣ 1 2 r =B B−1diagP rΣ 1 2 r, QrΣ 1 2 r diag Σ12 rQ∗r, Σ 1 2 rPr∗ B−1diagP rΣ 1 2 r, QrΣ 1 2 r .
From (3.2), (3.5), and the fact that diag(Σ12
rQ∗r, Σ 1 2 rPr∗)B−1diag(PrΣ 1 2 r, QrΣ 1 2 r) ∈
C4n×4n is nonsingular, we can see that span{B−1diag(P rΣ 1 2 r, QrΣ 1 2 r)} is an
invari-ant subspace of the GEP (3.1) corresponding to all nonzero eigenvalues, and therefore the result in (3.4) holds.
From Theorem 3.1, the null space free eigenvalue problem is derived straightfor-wardly in the following theorem.
Theorem 3.2 (the null space free eigenvalue problem). For any nonsingular
B defined in (3.3), the GEP (3.1) can be reduced to the following null space free eigenvalue problem: diag Σ12 rQ∗r, Σ 1 2 rPr∗ B−1diagP rΣ 1 2 r, QrΣ 1 2 r y = ωy. (3.6)
Furthermore, the GEP (3.1) and the NFSEP (3.6) have the same nonzero eigenvalues.
We have reduced the GEP (3.1) to the NFSEP (3.6), and the NFSEP can be solved by the iterative eigensolvers without being disturbed by the null space. Further computational considerations are discussed in the next section. Finally, we note how theC3n×3nmatricesζd,ξd,εd, andμdare determined. In Yee’s scheme,E and H are
evaluated at the edge centers and the face centers, respectively. However, CE and
C∗H are evaluated at the face centers and the edge centers, respectively. To match
these evaluation points, we can average the corresponding discrete entry values ofζ,
ξ, ε, and μ on the neighbor grid points to form the matrices ζd, ξd, εd, μd in (3.1).
Consequently,ζdE + μdH and εdE + ξdH are evaluated at the face centers and the
edge centers, respectively.
4. Computational and application considerations. The NFSEP defined
in (3.6) can actually be applied to various complex media settings as long as the
corresponding matrix B is nonsingular. Such media include general and Tellegen
bi-isotropic media [30], lossless and reciprocal bianisotropic media [26], and general bianisotropic media [16]. To solve the NFSEP, shift-and-invert type iterative eigen-solvers (e.g., Arnoldi method, Jacobi–Davidson method) can be applied to compute the desired eigenpairs of (3.1) from (3.6) without being affected by zero eigenval-ues. Despite the wide applications on complex media, the process for solving the NFSEP (3.6) can be further accelerated under the mild assumption described in sec-tion 4.1. It is worth mensec-tioning that two important types of media, i.e., chiral media [3, 4, 21, 27, 30, 33] and pseudochiral media [2, 5, 6, 16, 28], satisfy this assumption and can thus be solved by the accelerated eigensolvers. See sections 4.2 and 4.3 for more details and some computational remarks.
4.1. Sufficient conditions for Hermitian and Hermitian positive definite GEPs. The coefficient matrix in the NFSEP (3.6) is in a general form, and the
NFSEP can therefore be solved using, for example, the Arnoldi method. However, under an assumption, we can rewrite the NFSEP (3.6) as a GEP with Hermitian and Hermitian positive definite coefficient matrices. We can then take advantage of the matrix structure to accelerate the solution process by solving this rewritten eigenvalue problem via the invert Lanczos method and the associated conjugate gradient linear system solver.
The acceleration scheme is motivated from the following observations regarding
the matrix B in (3.3). If μd is nonsingular and we let Φ represent the matrixεd−
ξdμ−1d ζd, we have B = i 0 μd −Φ −ξd I3n 0 μ−1 d ζd I3n . Furthermore, if Φ is nonsingular, we have
B−1=−i I3n 0 −μ−1 d ζd I3n −Φ−1ξ dμ−1d −Φ−1 μ−1 d 0 . (4.1)
In other words, the properties of B−1 are closely related to μd, Φ, ξd, and ζd. In
particular, we consider the assumption
(4.2) μd 0, Φ ≡ εd− ξdμd−1ζd 0, and ξd∗=ζd.
The notationμd 0 and Φ 0 suggest that μdand Φ are Hermitian positive definite. Under this assumption, we show that the GEP (3.1) can be transformed to a standard Hermitian eigenvalue problem (so that all the eigenvalues are real) in Theorem 4.1. We then rewrite the NFSEP (3.6) in the new form (4.9) in Theorem 4.2. Consequently, the
corresponding coefficient matrixArto be defined in (4.10) is Hermitian and positive
definite. We can then use the Lanczos method, which consumes less storage and computation than Arnoldi-type methods, to solve (4.9).
We begin the derivation from the following theorem.
Theorem 4.1. Under Assumption (4.2), all eigenvalues ω of the GEP (3.1) are
real. Proof. Let Eζ Hζ = I3n 0 μ−1 d ζd I3n E H . (4.3)
Substituting (4.3) into (3.1) and premultiplying (3.1) by [ξI3n 0
dμ−1d I3n ], it holds that C 0 ξdμ−1d C − C∗μ−1d ζd C∗ Eζ Hζ = iω 0 μd −Φ 0 Eζ Hζ , (4.4)
where Φ is defined in assumption (4.2). By the assumptions thatμd 0 and Φ 0,
we then let
μd=μcμ∗c, Φ = ΦcΦ∗c
(4.5)
be the Cholesky decompositions ofμd and Φ, respectively. Define
E H = Φ∗c 0 0 μ∗c Eζ Hζ . (4.6)
Substituting (4.6) into (4.4) and premultiplying (4.4) by [ 0 −Φ−1c
μ−1 c 0 ], we have (4.7) Ax = ωx, wherex = [ E, H] and A = (−i) −Φ−1 c ξdμ−1d C − C∗μ−1d ζd (Φ∗c)−1 −Φ−1c C∗(μ∗c)−1 μ−1 c C (Φ∗c)−1 0 . (4.8)
BecauseA is Hermitian, all eigenvalues ω of the GEP (3.1) are real.
We have shown that all eigenvalues of (3.1) are real under assumption (4.2). However, the coefficient matrix of the NFSEP in (3.6) is not Hermitian. We reformu-late the NFSEP (3.6) in the following theorem to obtain a Hermitian and Hermitian positive definite GEP (HHPD-GEP).
Theorem 4.2. Under assumption (4.2), the GEP (3.1) can be reduced to a null
space free GEP
i 0 Σ−1r −Σ−1 r 0 yr=ω−1Aryr, (4.9) where Ar≡ diag (Pr∗, Q∗r) μ−1 d ζd −I3n I3n 0 Φ−1 0 0 μ−1d ξdμ−1d I3n −I3n 0 diag (Pr, Qr) (4.10)
is Hermitian and positive definite.
Proof. Let yr= diag Σ1/2r , Σ1/2r y. Rewrite (3.6) as
diag (Q∗r, Pr∗)B−1diag (Pr, Qr)yr=ω diagΣ−1r , Σ−1r yr,
which is equivalent to i diag (Pr∗, Q∗r) 0 I3n −I3n 0 B−1diag (P r, Qr)yr=ω i 0 Σ−1r −Σ−1 r 0 yr. (4.11)
From (4.1), it holds that i 0 I3n −I3n 0 B−1 =μ−1d ζd −I3n I3n 0 Φ−1ξdμ−1d Φ−1 −μ−1 d 0 = μ−1 d ζd −I3n I3n 0 Φ−1 0 0 μ−1d ξdμ−1d I3n −I3n 0 ,
which is Hermitian and positive definite if assumption (4.2) holds. That is, the
coef-ficient matrix on the left-hand side of (4.11) is equal toAr. Therefore, (4.11) can be
rewritten as (4.9).
We have now asserted the sufficient conditions that lead to the Hermitian and Hermitian positive definite null space free GEP (HHPD-NFGEP) (4.9). Next, we discuss some considerations in applying and solving the HHPD-NFGEP.
4.2. The eigenvalue and associated linear system solvers. In the
HHPD-NFGEP (4.9), the coefficient matrix Ar is Hermitian and positive definite. We can
use the generalized Lanczos method to solve (4.9) and obtain the smallest positive eigenvalues that are of interest in complex media. In each step of the generalized Lanczos method, we must solve the linear systems
Aru ≡ P∗ r Q∗ r ζd −I3n I3n 0 Φ−1 0 0 I3n ζ∗ d I3n −I3n 0 Pr Qr u = b (4.12)
for a given vector b. Because Ar is Hermitian positive definite, the linear system
(4.12) can be solved efficiently by using the conjugate gradient method. Furthermore,
the matrix-vector multiplications of the forms (T∗p, T q) for computing (Pr∗pˆ1, Q∗rpˆ2)
and (Prˆq1, Qrqˆ2), which are the most costly parts of solving (4.12), can be computed
efficiently by the 3D FFT because of the periodicity ofT , as shown in (2.8).
4.3. Application remarks. Intensive research has been conducted on chiral
and pseudochiral media. In this section, we assert that the corresponding magneto-electric matrices satisfy the assumptions given in (4.2), and we can therefore solve the HHPD-NFGEP using the Lanczos method. This solution procedure can thus act as a useful numerical tool for simulating 3D chiral and pseudochiral media.
First, we introduce the magnetoelectric matrices in chiral and pseudochiral me-dia. For chiral media (also called Pasteur or reciprocal chiral media), the associated
magnetoelectric matricesζd andξd in (3.1) are of the forms
ξd= iγ ˜I3n andζd=−iγ ˜I3n.
(4.13)
Here,γ is the chirality parameter, and ˜Im∈ Rm×mis a diagonal matrix whose entries
are equal to 0 (outside of the medium) or 1 (within the medium) depending on the corresponding grid point locations. For pseudochiral media, the associated matrices are ξd= ⎡ ⎣ 00 00 iγ ˜I0n iγ ˜In 0 0 ⎤ ⎦ and ζd= ⎡ ⎣ 00 00 −iγ ˜I0 n −iγ ˜In 0 0 ⎤ ⎦. (4.14)
Now, we analyze these magnetoelectric matrices. For chiral and pseudochiral
media, ε are of the form of a 3-by-3 diagonal block matrix where ε = diag( , , ).
Here, is a piecewise constant function that is equal to εiandεoinside and outside the
medium, respectively. Thus, the associated matrixεdin (3.1) is a diagonal matrix with
εoorεion diagonal entries. On the other hand, the permeabilityμ is usually taken as
I3, so the matrixμdin (3.1) is equal to an identity. Combining the diagonal matrices
εdandμdwith (ξd, ζd) in (4.13) and (4.14), we have that Φ =εd−ξdμ−1d ζd=εd−ξdζd
is a diagonal matrix with the entriesεo,εi, orεi− γ2. Becauseεoandεi are positive,
Φ is a positive diagonal matrix, providedγ ∈ (0, √εi).
4.4. A short summary. In Table 1, we summarize all the eigenvalue problems
and eigensolver strategies that have been discussed in the previous sections. From the algorithmic viewpoint, we propose and outline the null space free method in Algorithm 1.
5. Numerical results. In the numerical experiments, we consider 3D
recipro-cal chiral materials, in which ζ = −iγ. The goals of our numerical experiments are
threefold: to verify the correctness of the proposed algorithms and the code imple-mentation, to compare the performance of the proposed algorithm with an existing algorithm, and to study the performance of the proposed method in terms of itera-tion numbers and execuitera-tion time. The details of the numerical experiments are as follows.
Table 1
A summary of the eigenvalue problems and solvers considered in this article. The table lists (a) names of the eigenvalue problems, (b) type of the eigenvalue problems, (c) type of eigenvalues, (d) problem dimensions, (e) number of zero eigenvalues, (f) eigenvalue solver, (g) choice of shiftσ for the smallest eigenvalues, (h) embedded linear systems, (i) linear system solvers, (j) preconditioners for the linear system solvers, and (k) the applicable complex media.
(a) GEP (1.5) NFSEP (3.6) HSEP (4.7) HHPD-NFGEP (4.2) (b) Generalized Standard Standard Generalized
non-Hermitian non-Hermitian Hermitian Hermitian & HPD (c) Complex Complex Real (Thm. 4.1) Real
(d) 6n × 6n 4n × 4n 6n × 6n 4n × 4n
(e) 2n (Thm. 2.5) 0 (Thm. 3.2) 2n 0
(f) S.I. Arnoldi S.I. Arnoldi S.I. Lanczos S.I. Lanczos (g) Hard to choose Zero Hard to choose Zero
(h) (5.1) - - (4.12)
(i) LU or GMRES - - CG with FFT
(not efficient) Mtx-Vec Mult.
(j) Hard to find None (well-cond.) Harder to find None (well-cond.) (k) Media without Media without Media without Media satisfying
Assum. (4.2) Assum. (4.2) Assum. (4.2) Assum. (4.2) (Thm. 4.2)
Algorithm 1. The null space free method for solving the GEP (3.1).
1: Compute Λ1, Λ2, and Λ3 in Theorem 2.1;
2: Compute Λq and Λp in (2.11);
3: Compute Π1and Π2in (2.13a) and (2.14a), respectively;
4: ifμd and Φ are Hermitian positive definite andξd∗=ζd then
5: Solve the HHPD-NFGEP
i 0 Σ−1r −Σ−1 r 0 y = ω−1A ry, whereAr is defined in (4.10); 6: else
7: Solve the NFSEP
diag Σ1/2r Q∗r, Σ1/2r Pr∗ B−1diagP rΣ1/2r , QrΣ1/2r y = ωy. 8: Updatey = diag(Σ1/2r , Σ1/2r )y; 9: end if 10: Compute x = B−1diag (P r, Qr)y.
For the medium structure, we consider a simple cubic lattice consisting of spheres
with radiusr and circular cylinders with radius s, as shown in Figure 1. In particular,
we assume the lattice constanta = 1, r/a = 0.345, and s/a = 0.11. We use the triplet
(εi, εo, γ) to represent the associated permittivity inside the structure, the permittivity
outside the structure, and the chirality parameter, respectively. The perimeter of the irreducible Brillouin zone for the sample cubic lattice is formed by the corners
G = [0, 0, 0],X = 2π
a[12, 0, 0], M =2πa [12,12, 0], andR = 2πa[12,12,12].
For the implementation, the MATLAB functioneigs is used to solve the
HHPD-NFGEP (4.9), andpcg (without preconditioning) is used to solve the associated linear
system (4.12). The stopping criteria foreigs and pcg are 104× /(2
δ−2
x +δy−2+δz−2)
and 10−14, respectively. The constant (≈ 2.2 × 10−16) is the floating-point
rela-tive accuracy in MATLAB. In eigs, the maximal number of Lanczos vectors for the
restart is 3, where = 10 is the number of desired eigenvalues of the GEP (3.1). The
MATLAB functionsfftn and ifftn are applied to compute the matrix-vector
prod-uctsT∗p and T q, respectively. The MATLAB commands tic and toc are used to
measure the elapsed time. All computations are performed in MATLAB 2011b. For the hardware configuration, we use a HP workstation equipped with two Intel Quad-Core Xeon X5687 3.6GHz CPUs, 48 GB of main memory, and RedHat Linux operating system version 5.
5.1. Numerical correctness validation. We validate the correctness of the
proposed algorithm and MATLAB implementation by solving the following three sets of benchmark problems.
First, we consider a special case (εi, εo, γ) = (13, 1, 0), whose corresponding band
structures have been reported in [7, 10, 13, 20]. In this case, because of γ = 0 and
a
Fig. 1. A schema of 3D chiral medium with a simple cubic lattice within a single primitive cell.
G X M R G 0 0.1 0.2 0.3 0.4 0.5 Frequency ( λ a / 2 π ) (a) (εi, εo, γ) = (13, 1, 0) G X M R G 0 0.2 0.4 0.6 0.8 1 1.2 Frequency ( λ a / 2 π ) (b) (εi, εo, γ) = (1, 1, 0)
Fig. 2. The band structure with (εi, εo, γ) = (13, 1, 0) and (1, 1, 0).
(1.4a), we can see that the GEP (3.1) and the eigenvalue problemAE = ω2εdE (for
the photonic crystal as shown in [13]) lead to the same band structure. The computed
band structure due to (3.1) and n1 =n2 = n3 = 50 is shown in Figure 2(a). The
figure is identical (up to numerical precision) to Figure 1(a) in [13].
Second, we consider (εi, εo, γ) = (1, 1, 0), for which some theoretical results are
known. In this case, we know that (3.1) andAE = ω2E have the same band
struc-ture and, from Theorem 2.4, {Λq, Λq} are the nonzero eigenvalues of A. That is,
{Λ1/2q , Λ1/2q } are the eigenvalues of (3.1). Comparing the computed eigenvalues shown
in Figure 2(b) (for n1 = n2 = n3 = 50) with the exact eigenvalues, our numerical
results show that the maximal relative error of all computed eigenvalues in the figure
is 3.65 × 10−14.
104 105 106 107 108 1.64 1.66 1.68 λ1 104 105 106 107 108 1.64 1.66 1.68 λ2 104 105 106 107 108 2.18 2.19 2.2 2.21 λ3 matrix dimension 4n (a) (εi, εo, γ) = (13, 1, 0.5) 104 105 106 107 108 2.71 2.72 2.73 λ1 104 105 106 107 108 2.71 2.715 2.72 λ2 104 105 106 107 108 3.8 3.85 3.9 λ3 matrix dimension 4n (b) (εi, εo, γ) = (1, 1, 0.8)
Fig. 3. The convergent eigenvalues for (εi, εo, γ) = (13, 1, 0.5) and (εi, εo, γ) = (1, 1, 0.8) with
various matrix sizes 4n.
Third, we check the convergence of the eigenvalues in terms of the grid point
num-bers. In particular, we setn1=n2=n3= 2k fork = 3, . . . , 7, and the corresponding
matrix sizes 4n = 4 × 23k of the NFSEP (4.9) range from 2,048 to 8,388,608. The
three smallest positive eigenvaluesλ1,k,λ2,k, andλ3,kfor (εi, εo, γ) = (13, 1, 0.5) and
(εi, εo, γ) = (1, 1, 0.8) are shown in Figure 3 for the wave vector k = [0.5, 0, 0]. The
figure shows that{λ1,k}, {λ2,k}, and {λ3,k} are convergent as k increases.
5.2. Comparison with the shift-and-invert Arnoldi method. The GEP
(3.1) can be solved using the shift-and-invert Arnoldi method. In the shift-and-invert
Arnoldi method, the computational cost is dominated by solving the 6n × 6n linear
system C 0 0 C∗ − σ ζd I3n −εd −ξd y = b (5.1)
for a certain vector b and a shift σ. In contrast, the performance of the null space
free method (Algorithm 1) is dominated by solving the 4n × 4n linear system (4.12).
We thus compare the performance for solving these two linear systems. Here, we
take (εi, εo, γ) = (13, 1, 1) and k = [0.5, 0.5, 0]. To solve (5.1), we use (i) a direct
method based on LU factorization and (ii) the GMRES with SSOR preconditioner.
To solve (4.12), we use the MATLAB pcg without preconditioning. Note that the
chirality parameterγ = 1 implies that the coefficient matrix is Hermitian and positive
definite. The timing results for solving (5.1) and (4.12) are shown in Figure 4. The
results suggest that the performance of thepcg for solving (4.12) outperforms the two
solvers for solving (5.1) remarkably.
In Table 2, we further compare the performance of the two eigenvalue solvers: (i) the shift-and-invert Arnoldi method with direct linear system solver and (ii) the null
space free method (Algorithm 1) withpcg. It is clear that Algorithm 1 is much faster.
Consequently, we do not recommend solving the GEP (3.1) by the shift-and-invert Arnoldi method unless we can develop a good preconditioning scheme for solving the linear system (5.1). However, it is important to note that, even with a good
preconditioner, the effect of a large null space (with rank 2n) can downgrade the
performance significantly [12].
104 105 106 107 108 0 200 400 600 800 1000 1200 1400 1600 1800 2000 6*n CPU times (s) GMRES+SSOR A\b PCG
Fig. 4. The time for solving (5.1) by left matrix divide and thegmres and solving (4.12) by pcg.
Table 2
CPU time in seconds for solving the GEP (3.1). S.I. Arnoldi+LU stands for the shift-and-invert Arnoldi method with LU based linear system solver. We take n1 = n2 = n3 = 32 and (εi, εo, γ) = (13, 1, 0.5).
k (14, 0, 0) (12,12, 0) (12,12,12) (14,14,14) S.I. Arnoldi+LU 18,821 10,533 16,628 20,758
Algorithm 1 155 140 198 155
5.3. Performance of Algorithm 1. We now concentrate on the performance,
in terms of iteration numbers and timing, of the null space free method (Algorithm 1) in finding the 10 smallest positive eigenvalues of the GEP (3.1) with various
combi-nations of the parametersεi,γ, and k. We take n1=n2=n3= 128, and the size of
the coefficient matrix in (4.9) is 4n31= 8,388,608. The Lanczos method is applied to
solve the HHPD-NFGEP (4.9).
In the first test problem set, we vary the wave vector 2πk along the segments
connectingG, X, M, R, and G in the first Brillouin zone to plot the band structure.
In each of the segments, 10 uniformly distributed sampling wave vectors are chosen.
The results are shown in Figure 5 for (εi, εo, γ) = (13, 1, 0.5) or (εi, εo, γ) = (1, 1, 0.8).
In the second test problem set, we change the chirality parameter γ from 0.25 to 3
for εi = 13 (from γ = 0.05 to 0.9 for εi = 1). Note that the condition number of
the linear system (4.12) increases from εi
εo = 13 to ∞ (singular) as γ varies from 0
to √εi ≈ 3.61. We fix k = [0.5, 0, 0]. The results are shown in Figure 6. Based on
Figures 5 and 6, we highlight the following performance results:
• The iteration numbers are very small. Figures 5(c), 5(d) 6(c), and 6(d) show the iteration numbers of the Lanczos method for solving (4.9) with different parameter combinations. In all cases, the iteration numbers are substantially smaller than the matrix size 8,388,608. Note that some higher iteration num-bers in Figures 5(c) and 5(d) are due to the corresponding multiplicity of eigenvalues. The higher iteration numbers in Figures 6(c) and 6(d) are due to the clustering of the eigenvalues.
G X M R G 0 0.1 0.2 0.3 0.4 0.5 Frequency ( λ a / 2 π )
(a) Band structure for (εi, εo, γ) = (13, 1, 0.5)
G X M R G 0 0.2 0.4 0.6 0.8 1 1.2 Frequency ( λ a / 2 π )
(b) Band structure for (εi, εo, γ) = (1, 1, 0.8)
G X M R G 120 130 140 150 160 170 180 190
Iteration numbers of Lanczos method
(c) Iteration numbers ranging from 120 to 190 with average 143 for (εi, εo, γ) = (13, 1, 0.5)
G X M R G 100 110 120 130 140 150 160 170 180 190 200
Iteration numbers of Lanczos method
(d) Iteration numbers ranging from 100 to 200 with average 136 for (εi, εo, γ) = (1, 1, 0.8)
G X M R G
2.5 3 3.5 4
CPU times (hour)
(e) CPU time ranging from 2.75 to 4 hours with average 3.2 hours for (εi, εo, γ) = (13, 1, 0.5)
G X M R G 2 2.5 3 3.5 4 4.5
CPU times (hour)
(f) CPU time ranging from 2.1 to 4.1 hours with average 2.8 hours for (εi, εo, γ) = (1, 1, 0.8)
Fig. 5. The band structures of 3D chiral media, iteration number of the Lanczos method,
and the elapsed times of the NSF (Algorithm 1) for solving (3.1) associated with various wave vectors k.
• The timing is satisfactory. Figures 5(e) and 5(f) show the CPU times for solving (3.1). These results suggest that our approach is efficient for comput-ing 10 (interior) eigenvalue problems as large as 8,388,608. This efficiency is mainly due to the highly efficient linear system solver for (4.12).
0 0.5 1 1.5 2 2.5 3 1 1.5 2 2.5 3 3.5 chirality γ Frequency ( λ / (2 π ))
(a) Band structure for (εi, εo, γ) = (13, 1, γ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 chirality γ Frequency ( λ / (2 π ))
(b) Band structure for (εi, εo, γ) = (1, 1, γ) 0 0.5 1 1.5 2 2.5 3 120 125 130 135 140 145 150 chirality γ Iteration of Lanczos
(c) Iteration numbers ranging from 124 to 147 for (εi, εo, γ) = (13, 1, γ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 120 140 160 180 200 220 240 260 chirality γ Iteration of Lanczos
(d) Iteration numbers ranging from 140 to 260 for (εi, εo, γ) = (1, 1, γ) 0 0.5 1 1.5 2 2.5 3 2.5 3 3.5 4 4.5 5 5.5 6 6.5 chirality γ
CPU times (hour)
(e) CPU time increasing from 2.5 to 6.5 hours for (εi, εo, γ) = (13, 1, γ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.5 2 2.5 3 3.5 4 4.5 5 chirality γ
CPU times (hour)
(f) CPU time increasing from 1.4 to 4.5 hours for (εi, εo, γ) = (1, 1, γ) 0 0.5 1 1.5 2 2.5 3 30 40 50 60 70 80 90 100 chirality γ Average iteration of CG
(g) Iteration numbers of CG for (εi, εo, γ) = (13, 1, γ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 10 20 30 40 50 60 chirality γ Average iteration of CG
(h) Iteration numbers of CG for (εi, εo, γ) = (1, 1, γ)
Fig. 6. The band structures of 3D chiral media, the average of iteration number ofpcg, the
iter-ation numbers of the Lanczos method, and the elapsed times of the NSF for solving (3.1) associated with various chirality parametersγ.
• The linear systems in the form of (4.12) are well-conditioned for the tested γ’s. We take a close look of the behavior of pcg for solving the linear systems
(4.12) with variousγ’s. As shown in Figures 6(g) and 6(h), the (average) pcg
iteration numbers for solving linear system (4.12) in the tested eigenvalue
problems increase from 39 to 90 for εi = 13 and increase from 9 to 59 for
εi = 1. This behavior is parallel to the increase in the condition number
of the linear system (4.12) from εi
εo = 13 to ∞ (singular) as γ varies from
0 to √εi ≈ 3.61. We fix k = [0.5, 0, 0]. A more important observation
is that the timing results are quite satisfactory in the following sense. For
problems as large as 8.4 million and a stopping tolerance as small as 10−14,
these small iteration numbers suggest that the coefficient matrix in (4.12) is quite well-conditioned.
6. Conclusions. In this paper, we focus on the GEPs arising in the source-free
Maxwell equation with magnetoelectric coupling effects in the 3D chiral media. Solv-ing the GEP is a computational challenge. We have proposed a promisSolv-ing theoretical framework for efficiently solving the eigenvalue problem. First, we derive the SVD of the discrete single-curl operator. Using this SVD, we explore an explicit form of the basis for the invariant subspace corresponding to nonzero eigenvalues of the GEP. By applying the basis, the GEP is reduced to an NFSEP, which involves only the eigenspace associated with the nonzero eigenvalues of the GEP and excludes the zero eigenvalues so that they do not degrade the computational efficiency. Next, we show
that all nonzero eigenvalues of the GEP are real ifμdandεd− ξdμ−1d ζdare Hermitian
positive definite andξd∗=ζd. Based on this property, we reformulate the NFSEP to a
null space free GEP whose coefficient matrices are Hermitian and Hermitian positive definite (HHPD-NFGEP). We can then use the invert Lanczos method to solve the HHPD-NFGEP and the CG method to solve the embedded linear systems. The nu-merical results validate the correctness of the proposed algorithms and the computer code implementation. The results also suggest that the proposed methods are efficient in terms of iteration and timing.
Appendix. Proof of Lemma 2.2.
(a) For the sake of contradiction, it is assumed that Λq is singular. By the
definition of Λq, Λq is singular if and only if Λa1,n1, Λa2,n2, and Λa3,n3 are singular if
and only if
eθn1,i+θa1,n1− 1 = 0,
eθn2,j+θa2,n2− 1 = 0,
eθn3,k+θa3,n3− 1 = 0
for somei ∈ {0, . . . , n1− 1}, j ∈ {0, . . . , n2− 1}, and k ∈ {0, . . . , n3− 1}. That is
i + k · a1 n1 = i + k1 n1 , j + k · a2 n2 = j + k2 n2 , k + k · a3 n3 = k + k3 n3
are integers for some i, j, k. By the assumption 0 ≤ k1, k2, k3 ≤ 12, we have i = 0,
j = 0, k = 0, and k1 = k2 = k3 = 0, which contradict k = 0. Therefore, Λq is
nonsingular.
(b) By the definitions of Λ1, Λ2, and Λ3in (2.10), the ((k−1)n1n2+(j−1)n1+i)th
elements of Λ1, Λ2, and Λ3 are
δ−1
x (eθi− 1), δ−1y (eθj− 1), δ−1z (eθk− 1),
respectively, where θi= 2π i + k1 n1 , θj = 2π j + k2 n2 , θk = 2π k + k3 n3 ,
fori = 0, . . . , n1− 1, j = 0, . . . , n2− 1, and k = 0, . . . , n3− 1. Assume that Λp does
not have full column rank. Then there exists somei, j, and k such that
βδ−1 z (eiθk− 1) = δ−1y (eiθj − 1), δ−1 x (eiθi− 1) = αδz−1(eiθk− 1), βδ−1 x (eiθi− 1) = αδy−1(eiθj− 1),
which implies that
β sin θk δz = sinθj δy , sinθi δx =α sin θk δz , β sin θi δx = α sin θj δy and β(cos θk− 1) δz = cosθj− 1 δy , cosθi− 1 δx = α(cos θk− 1) δz , β(cos θi− 1) δx = α(cos θj− 1) δy or equivalent to β sin θi αδx =sinθj δy =β sin θk δz (6.1) and β(cos θi− 1) αδx = cosθj− 1 δy = β(cos θk− 1) δz . (6.2)
From (6.1) and (6.2), it holds that β sin θi αδx 2 + β(cos θi− 1) αδx 2 = sinθj δy 2 + cosθj− 1 δy 2 = β sin θk δz 2 + β(cos θk− 1) δz 2 and then β2(2− 2 cos θ i) α2δ2 x = 2− 2 cos θj δ2 y = β 2(2− 2 cos θ k) δ2 z . Therefore, β(cos θi− 1) αδx =αδx βδy cosθj− 1 δy , β(cos θk− 1) δz = δz βδy cosθj− 1 δy . (6.3)
From (6.2) and (6.3), we can see that ifαδx= βδy andδz= βδy, then
cosθi= cosθj = cosθk = 1.
That is i+k1
n1 ,
j+k2
n2 , and
k+k3
n3 must be integers. This contradicts the assumption for
k. Therefore, Λp has full column rank.
Acknowledgment. The authors appreciate useful comments and suggestions
from the anonymous referees.
REFERENCES
[1] X. Cheng, H. Chen, L. Ran, B.-I. Wu, T. M. Grzegorczyk, and J. A. Kong, Negative
refraction and cross polarization effects in metamaterial realized with bianisotropic s-ring resonator, Phys. Rev. B, 76 (2007), 024402.
[2] R.-L. Chern, Anomalous dispersion in pseudochiral media: Negative refraction and backward
wave, J. Phys. D, 46 (2013), 125307.
[3] R.-L. Chern, Wave propagation in chiral media: Composite Fresnel equations, J. Opt., 15 (2013), 075702.
[4] R.-L. Chern and P.-H. Chang, Negative refraction and backward wave in chiral mediums:
Illustrations of Gaussian beams, J. Appl. Phys., 113 (2013), 153504.
[5] R.-L. Chern and P.-H. Chang, Negative refraction and backward wave in pseudochiral
medi-ums: Illustrations of Gaussian beams, Opt. Express, 21 (2013), pp. 2657–2666.
[6] R.-L. Chern and P.-H. Chang, Wave propagation in pseudochiral media: Generalized Fresnel
equations, J. Opt. Soc. Amer. B, 30 (2013), pp. 552–558.
[7] R. L. Chern, C. C. Chang, C.-C. Chang, and R. R. Hwang, Numerical study of
three-dimensional photonic crystals with large band gaps, J. Phys. Soc. Japan, 73 (2004),
pp. 727–737.
[8] J. Chongjun, Q. Bai, Y. Miao, and Q. Ruhu, Two-dimensional photonic band structure in
the chiral medium-transfer matrix method, Opt. Commun., 142 (1997), pp. 179–183.
[9] C. Engstr¨om and M. Richter, On the spectrum of an operator pencil with applications to
wave propagation in periodic and frequency dependent materials, SIAM J. Appl. Math., 70
(2009), pp. 231–247.
[10] P.-F. Hsieh, T.-T. Wu, and J.-H. Sun, Three-dimensional phononic band gap calculations
using the FDTD method and a pc cluster system, IEEE Trans. Ultrason. Ferroelectr. Freq.
Control, 53 (2006), pp. 148–158.
[11] T.-M. Huang, H.-E. Hsieh, W.-W. Lin, and W. Wang, Eigendecomposition of the discrete
double-curl operator with application to fast eigensolver for three dimensional photonic crystals, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 369–391.
[12] T.-M. Huang, H.-E. Hsieh, W.-W. Lin, and W. Wang, Fast Lanczos Eigenvalue Solvers
for Band Structures of Three Dimensional Photonic Crystals with Face-Centered Cubic Lattice, Technical report 2013-3-003, NCTS Preprints in Mathematics, National Tsing
Hua University, Hsinchu, Taiwan, 2013.
[13] T.-M. Huang, Y.-C. Kuo, and W. Wang, Computing extremal eigenvalues for
three-dimensional photonic crystals with wave vectors near the Brillouin zone center, J. Sci.
Comput., 55 (2013), pp. 529–551.
[14] T.-M. Huang, W.-W. Lin, and W. Wang, Matrix Representation of Discrete Differential
Op-erators and Operations in Electromagnetism, Technical report 2014-8-003, NCTS Preprints
in Mathematics, National Tsing Hua University, Hsinchu, Taiwan, 2014.
[15] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals:
Molding the Flow of Light, Princeton University Press, Princeton, NJ, 2008.
[16] E. O. Kamenetskii, Theory of bianisotropic crystal lattices, Phys. Rev. E, 57 (1998), pp. 3563–3573.
[17] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 2005.
[18] J. A. Kong, Theorems of bianisotropic media, Proc. IEEE, 60 (1972), pp. 1036–1046.
[19] J. Lekner, Optical properties of isotropic chiral media, Pure Appl. Opt., 5 (1996), pp. 417–443. [20] S.-Y. Lin, J. G. Fleming, R. Lin, M. M. Sigalas, R. Biswas, and K. M. Ho, Complete
three-dimensional photonic bandgap in a simple cubic structure, J. Opt. Soc. Amer. B, 18
(2001), pp. 32–35.
[21] Y. Liu and X. Zhang, Metamaterials: A new frontier of science and technology, Chem. Soc. Rev., 40 (2011), pp. 2494–2507.
[22] T. G. Mackay and A. Lakhtakia, Negative refraction, negative phase velocity, and
counter-position in bianisotropic materials and metamaterials, Phys. Rev. B, 79 (2009), 235121.
[23] M. Reed and B. Simon, Methods of modern mathematical physics, in Analysis of Operators IV, Academic Press, San Diego, CA, 1978.
[24] K. Schmidt and R. Kappeler, Efficient computation of photonic crystal waveguide modes
with dispersive material, Opt. Express, 18 (2010), pp. 7307–7322.
[25] K. Schmidt and P. Kauf, Computation of the band structure of two-dimensional photonic
crystals with hp finite elements, Comput. Methods Appl. Mech. Engrg., 198 (2009),
pp. 1249–1259.
[26] A. Serdyukov, I. Semchenko, S. Tretyakov, and A. Sihvola, Electromagnetics of
Bi-anisotropic Materials: Theory and Applications, Gordon and Breach Science, New York,
2001.
[27] A. H. Sihvola, A. J. Viitanen, I. V. Lindell, and S. A. Tretyakov, Electromagnetic Waves
in Chiral and Bi-Isotropic Media, Artech House, Boston, 1994.
[28] S. A. Tretyakov, A. H. Sihvola, A. A. Sochava, and C. R. Simovski,
Magneto-electric interactions in bi-anisotropic media, J. Electromagn. Waves Appl., 12 (1998),
pp. 481–497.
[29] S. A. Tretyakov, C. R. Simovski, and M. Hudliˇcka, Bianisotropic route to the realization
and matching of backward-wave metamaterial slabs, Phys. Rev. B, 75 (2007), 153104.
[30] B. Wang, J. Zhou, T. Koschny, M. Kafesaki, and C. M Soukoulis, Chiral metamaterials:
Simulations and experiments, J. Opt. A, 11 (2009), 114003.
[31] W. S. Weiglhofer and A. Lakhtakia, Introduction to Complex Mediums for Optics and
Electromagnetics, SPIE, Washington, DC, 2003.
[32] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations
in isotropic media, IEEE Trans. Antennas and Propagation, 14 (1966), pp. 302–307.
[33] R. Zhao, T. Koschny, and C. M. Soukoulis, Chiral metamaterials: Retrieval of the effective
parameters with and without substrate, Optics Express, 18 (2010), pp. 14553–14567.
[34] N. Zhou, J. V. Clark, and K. S. J. Pister, Nodal simulation for MEMS design using SUGAR v0.5, in Proceedings of International Conference on Modeling and Simulation of Microsys-tems Semiconductors, Sensors and Actuators, Santa Clara, CA, 1998, pp. 308–313.
