FAME-matlab Package: Fast Algorithm for Maxwell Equations

(1)

FAME-matlab Package: Fast Algorithm for Maxwell Equations

Tsung-Ming Huang 

Modelling, Simulation and Analysis of

Nonlinear Optics, NUK, September, 4-8, 2015

1

(2)

2

(3)

3 FAME group

2 Wen-Wei Lin

Department of Applied Mathematics National Chiao-Tung University

Weichung Wang

Department of Mathematics National Taiwan University

Han-En Hsieh(謝函恩)

Department of Mathematics National Taiwan University

Chien-Chih Huang(黃建智)

Department of Mathematics

National Taiwan Normal University

(4)

F ast A lgorithm for M axwell E q

4 FAME

FAME. ^m

FAME. ^mpi

FAME. ^GPU

(5)

FAME

5 Eigen-solvers (JD, SIRA)

Dispersive Metallic

materials Complex materials Photonic Crystals

Simple Cubic Face-Centered Cubic

∇ × ∇ × E(r) = λε ^(r)E(r) ∇ × ∇ × E(r) = ω ² ε ^(r, ω ^)E(r) ^⎡ _⎣ ^⎢ ^{∇ ×} ₀ _{∇ ×} ⁰ ^⎤ _⎦ ^⎥ ^⎡ _⎣ ^⎢ _H ^E ^⎤ _⎦ ^{⎥ = i} ^ω ^⎡ ₋ ^ζ _ε ₋ ^µ _ξ

⎣ ⎢

⎢

⎤

⎦ ⎥

⎥ E H

⎡

⎣ ⎢ ⎤

⎦ ⎥

employing projections onto the range space of the discrete matrix. (However, it is noted that the problem of many zeros does not occur in the PWE as in each plane wave component the number of polarizations is chosen to be the independent two rather three.) (iii) As the third issue of importance, it is of great interest to develop a fast algorithm for the eigenvalue problem. An inverse method, accelerated by multigrid technique with use of projection is proposed for this purpose. The method exploits the sparsity of the matrix for the eigenvalue problem in the finite diﬀerence formulation. Because of the above mentioned diﬃculties, the present method is a nontrivial extension of a similar method, recently developed by the authors²⁴⁾for computing photonic band structures in two dimensions.

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.²⁵⁾Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.²⁶⁾The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,^27,28)as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.²⁹⁾The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp³-like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.^12,28)In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp³structure of the electrons of diamond atoms.

Recently, submicron diamond-lattice photonic crystals have been successfully produced by two-photon laser nanofabri-

cation (photopolymerization).³²⁾

The order of presentation of the paper is organized as follows. In §2, we show how to correctly formulate the finite diﬀerence method for the double curl operator of the photonic eigenvalue problem. In §3, we develop the numerical method (inverse iteration with the full multigrid acceleration) and present the fast algorithm, in which two alternative methods of projection are proposed to avoid the necessity of deflating zeros). In §4, we first present numerical results that illustrate the eﬃciency of the presently developed method. Then, the band structures are computed for the modified simple cubic lattice, the tetragonal square spiral structure (direct and inverse structure) and the diamond structure with sp³-like configuration. Finally, concluding remarks with a summary of results are drawn in

§5.

2. Basic Equations and Finite Diﬀerence Formulation The electromagnetic waves propagating in the photonic crystals are well described by Maxwell’s equations. For linear isotropic and frequency-independent dielectric materials with permeability close to one, the time-harmonic Fig. 1. Modified simple cubic lattice comprising dielectric spheres and

connecting thin circular cylinders.²⁵⁾

Fig. 2. Tetragonal square spiral structure comprising circular cylinders.^27,28)

Fig. 3. diamond structure with sp³-like configuration comprising dielectric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.²⁵⁾Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.²⁶⁾ The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,^27,28) as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.²⁹⁾The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp³-like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.^12,28)In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp³structure of the electrons of diamond atoms.

§5.

connecting thin circular cylinders.²⁵⁾

Fig. 2. Tetragonal square spiral structure comprising circular cylinders.^27,28)

(6)

Generalized eigenvalue problems for 3D

photonic crystal

6

(7)

Curl operator

Central edge points 

Central face points  

where 

Resulting generalized eigenvalue problem 

with diagonal B 

7 ∇ × ∇ × E(r) = ω ² ε ^(r)E(r)

∇ × E(r) = H(r) ⇒ Ce = h

∇ × H(r) = ω ² ε (r)E(r) ⇒ C ^∗ h = ω ² B e

∇ × E =

0 − ∂

∂z

∂

∂y

∂

∂z 0 − ∂

∂x

− ∂ ∂y

∂

∂x 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ E

₁

E

₂

E

₃

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ ∈! ^3n×3n

C ₁ = I _n ₂ _n ₃ ⊗ K ₁ ∈! ⁿ ^×n , C ₂ = I _n ₃ ⊗ K ₂ ∈! ⁿ ^×n , C ₃ = K ₃ ∈! ⁿ ^×n

C ^∗ C − ω ² B

( ) ^x ^{≡ A −} ⁽ ^λ ^B ⁾ ^x ^{= 0}

(8)

Finite Diff. Assoc. with Quasi-Periodic Cond.

8 K ₁ = 1 δ _x

−1 1

! ! −1 1 e ^ı2 ^π ^k⋅a ¹ −1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁿ ¹ ^×n ¹ ,

K ₂ = 1 δ _y

−I _n ₁ I _n ₁

! !

−I _n ₁ I _n

1 e ^ı2 ^π ^k⋅a ² J ₂ −I _n ₁

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾ ,

K ₃ = 1 δ _z

−I _n ₁ _n ₂ I _n ₁ _n ₂

! !

−I _n ₁ _n ₂ I _n

1 n ₂

e ^ı2 ^π ^k⋅a ³ J ₃ −I _n ₁ _n ₂

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ^n×n

(9)

Finite Diff. Assoc. with Quasi-Periodic Cond.

8 K ₁ = 1 δ _x

−1 1

! ! −1 1 e ^ı2 ^π ^k⋅a ¹ −1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁿ ¹ ^×n ¹ ,

K ₂ = 1 δ _y

−I _n ₁ I _n ₁

! !

−I _n ₁ I _n

1 e ^ı2 ^π ^k⋅a ² J ₂ −I _n ₁

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾ ,

K ₃ = 1 δ _z

−I _n ₁ _n ₂ I _n ₁ _n ₂

! !

−I _n ₁ _n ₂ I _n

1 n ₂

e ^ı2 ^π ^k⋅a ³ J ₃ −I _n ₁ _n ₂

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ^n×n

E(r + a _ℓ ) = e ⁱ² ^π ^k⋅a ^ℓ E(r)

(10)

Finite Diff. Assoc. with Quasi-Periodic Cond.

8 K ₁ = 1 δ _x

−1 1

! ! −1 1 e ^ı2 ^π ^k⋅a ¹ −1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁿ ¹ ^×n ¹ ,

K ₂ = 1 δ _y

−I _n ₁ I _n ₁

! !

−I _n ₁ I _n

1 e ^ı2 ^π ^k⋅a ² J ₂ −I _n ₁

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾ ,

K ₃ = 1 δ _z

−I _n ₁ _n ₂ I _n ₁ _n ₂

! !

−I _n ₁ _n ₂ I _n

1 n ₂

e ^ı2 ^π ^k⋅a ³ J ₃ −I _n ₁ _n ₂

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ^n×n

(11)

Finite Diff. Assoc. with Quasi-Periodic Cond.

8 K ₁ = 1 δ _x

−1 1

! ! −1 1 e ^ı2 ^π ^k⋅a ¹ −1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁿ ¹ ^×n ¹ ,

K ₂ = 1 δ _y

−I _n ₁ I _n ₁

! !

−I _n ₁ I _n

1 e ^ı2 ^π ^k⋅a ² J ₂ −I _n ₁

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾ ,

K ₃ = 1 δ _z

−I _n ₁ _n ₂ I _n ₁ _n ₂

! !

−I _n ₁ _n ₂ I _n

1 n ₂

e ^ı2 ^π ^k⋅a ³ J ₃ −I _n ₁ _n ₂

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ^n×n

J ₂ = I _n ₁ , J ₃ = I _n ₁ _n ₂

For SC lattice

(12)

Finite Diff. Assoc. with Quasi-Periodic Cond.

8 K ₁ = 1 δ _x

−1 1

! ! −1 1 e ^ı2 ^π ^k⋅a ¹ −1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁿ ¹ ^×n ¹ ,

K ₂ = 1 δ _y

−I _n ₁ I _n ₁

! !

−I _n ₁ I _n

1 e ^ı2 ^π ^k⋅a ² J ₂ −I _n ₁

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾ ,

K ₃ = 1 δ _z

−I _n ₁ _n ₂ I _n ₁ _n ₂

! !

−I _n ₁ _n ₂ I _n

1 n ₂

e ^ı2 ^π ^k⋅a ³ J ₃ −I _n ₁ _n ₂

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

∈! ^n×n

J ₂ = I _n ₁ , J ₃ = I _n ₁ _n ₂

For SC lattice

J ₂ = 0 e ^−ı2 ^πk⋅a ¹ I _n

1 /2

I _n

1 /2 0

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥ ∈! ⁿ ¹ ^×n ¹ ,

J ₃ =

0 e ^{−ı2πk⋅a} ² I ₁

3 n ₂ ⊗ I _n

1 I ₂

3 n ₂ ⊗ J ₂ 0

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ⁽ⁿ ¹ ⁿ ² ^)×(n ¹ ⁿ ² ⁾

For FCC lattice

(13)

Power method

Let be the eigenpairs of A where  is linearly independent

For any nonzero vector u 

Since , we have 

If for i >1 and , then 

Given shift value

9 ( λ _i ^{, x} _i ^{) for i} = 1,…,n

x ₁ , …, x _n

u = α ₁ ^x ₁ +!+ α _n ^x _n

A ^k x _i = λ _i ^k ^x _i

A ^k u = α ₁ λ ₁ ^k ^x ₁ +!+ α _n λ _n ^k ^x _n

λ ₁ > λ _i α ₁ ≠ 0

1 λ ₁ ^k ^A

k u = α ₁ ^x ₁ + ( λ ₂ λ ₁ ⁾

k α ₂ ^x ₂ +!+ α _n ⁽ λ _n λ ₁ ⁾

k x _n → α ₁ ^x ₁ ^{as k} → ∞

(A − σ ^{I )} ⁻¹

{ } ^k ^u ⁼ ^α ¹ { ⁽ ^λ ¹ ⁻ ^σ ⁾ ⁻¹ } ^k ^x ¹ ^+!+ ^α ⁿ { ⁽ ^λ ⁿ ⁻ ^σ ⁾ ⁻¹ } ^k ^x ⁿ

(14)

Solving

Use shift-and-invert Lanczos method

In each iteration of shift-and-invert Lanczos method, we need to solve 

How to efficiently solve this linear system? 

10 A − λ ^B

( ) ^x ^{= 0}

(A − σ ^B)y = b

(15)

Solving linear system

11 (A − σ ^B)y = b

(16)

Solve

Direct method (Gaussian elimination) 

Iterative method

Matrix vector multiplication with Preconditioner M 

sol = bicgstabl(coef_mtx, rhs, tol, maxit,

@(x)SSOR_prec(x, diag_coef_mtx, lower_L));

12 (A − σ ^B)y = b

y = (A − σ ^{B) \ b}

A − σ ^B

(17)

Solve

Direct method (Gaussian elimination) 

Iterative method

Matrix vector multiplication with Preconditioner M 

sol = bicgstabl(coef_mtx, rhs, tol, maxit,

@(x)SSOR_prec(x, diag_coef_mtx, lower_L));

12 (A − σ ^B)y = b

y = (A − σ ^{B) \ b}

A − σ ^B

Demo performance

(18)

Eigen-decomp. of C ₁ , C ₂ , C ₃ for SC lattice

Define 

Define unitary matrix T as 

Then it holds that 

13 D _a,m = diag 1,e ( ^θ ^{a ,m} , !,e ^(m ⁻¹⁾ ^θ ^{a ,m} ) ^, ^Λ ^a,m ^{= diag e} ( ^θ ^{m ,1} ⁺ ^θ ^{a ,m} ^{−1 ! e} ^θ ^{m ,m} ⁺ ^θ ^{a ,m} ⁻¹ ) ^,

U _m =

1 1 ! 1

e ^θ ^{m ,1} e ^θ ^{m ,2} ! 1

! ! !

e ^(m−1) ^θ ^{m ,1} e ^(m−1) ^θ ^{m ,2} ! 1

⎡

⎣

⎢ ⎢

⎤

⎦

⎥ ⎥

∈! ^m×m , θ _a,m = ı2 π ^k ⋅a

m , θ _m,i = ı2 π ⁱ m

T = 1

n D _a

3 ,n ₃ ⊗ D _a ₂ _,n ₂ ⊗ D _a ₁ _,n ₁

( ) ( ^U ⁿ ³ ^⊗U ⁿ ² ^⊗U ⁿ ¹ )

C ₁ T = δ _x ⁻¹ ^{T I} ( _n ₃ ⊗ I _n ₂ ⊗ Λ _a ₁ _,n ₁ ) ^{≡ T Λ} ¹ ^,

C ₂ T = δ _y ⁻¹ ^{T I} ( _n ₃ ⊗ Λ _a ₂ _,n ₂ ⊗ I _n ₁ ) ^{≡ T Λ} ² ^,

C ₃ T = δ _z ⁻¹ ^T ( Λ _a ₃ _,n ₃ ⊗ I _n ₂ ⊗ I _n ₁ ) ^{≡ T Λ} ³

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.²⁵⁾Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.²⁶⁾The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,^27,28)as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.²⁹⁾The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp³-like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.^12,28)In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp³structure of the electrons of diamond atoms.

§5.

connecting thin circular cylinders.²⁵⁾ Fig. 2. Tetragonal square spiral structure comprising circular cylinders.^27,28)

(19)

Eigen-decomp. of C ₁ , C ₂ , C ₃ for FCC lattice

Define  

Define unitary matrix T as 

Then it holds that 

14 x _i = D _x U _n ₁ (:,i), y _{i, j} = D _y,i U _n ₂ (:, j)

T = 1

n ⎡ T ₁ T ₂ ! T _n ₁

⎣⎢ ⎤

⎦⎥ ∈! ⁿ ^×n , T _i = T ⎡ _i,1 T _i,2 ! T _i,n ₂

⎣⎢ ⎤

⎦⎥ ∈! ⁿ ^×(n ² ⁿ ³ ⁾ , T _{i, j} = D _z,i _{+ j} U _n

( 3 ) ^{⊗ y} ⁽ ^{i, j} ^{⊗ x} ⁱ ⁾

C ₁ T = T Λ ( _n ₁ ⊗ I _n ₂ _n ₃ ) ^{≡ T Λ} ¹ ^,

C ₂ T = T ⊕ ( ( _i=1 ⁿ ¹ Λ _i,n ₂ ) ^{⊗ I} ⁿ ³ ) ^{≡ T Λ} ² ^,

C ₃ T = T ⊕ ( _i=1 ⁿ ¹ ⊕ ⁿ _j ₌₁ ² Λ _{i, j,n} ₃ ) ^{≡ T Λ} ³

ψ _x = ı2 π ^k ⋅a ₁

n ₁ , D _x = diag 1,e ( ^ψ

^x

, !,e ⁽ⁿ

¹

⁻¹⁾ ^ψ

^x

) ^,

ψ _y,i = ı2 π

n ₂ k ⋅ a ₂ − a ₁ 2

⎛ ⎝⎜ ⎞

⎠⎟ − i 2

⎧ ⎨

⎩

⎫ ⎬

⎭ , D _y,i = diag 1,e ( ^ψ

^{y ,i}

, !,e ⁽ⁿ

²

⁻¹⁾ ^ψ

^{y ,i}

) ^,

ψ _{z,i+ j} = ı2 π

n ₃ k ⋅ a ₃ − a ₁ + a ₂ 3

⎛ ⎝⎜ ⎞

⎠⎟ −

i + j 3

⎧ ⎨

⎩

⎫ ⎬

⎭ , D _{z,i+ j} = diag 1,e ( ^ψ

^{y ,i+ j}

, !,e ⁽ⁿ

³

⁻¹⁾ ^ψ

^{y ,i+ j}

)

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.²⁵⁾Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.²⁶⁾ The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,^27,28)as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.²⁹⁾The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp³-like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.^12,28)In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp³structure of the electrons of diamond atoms.

§5.

(20)

Eigen-decomp. of C ₁ , C ₂ , C ₃ for FCC lattice

Define  

Define unitary matrix T as 

Then it holds that 

14 x _i = D _x U _n ₁ (:,i), y _{i, j} = D _y,i U _n ₂ (:, j)

T = 1

n ⎡ T ₁ T ₂ ! T _n ₁

⎣⎢ ⎤

⎦⎥ ∈! ⁿ ^×n , T _i = T ⎡ _i,1 T _i,2 ! T _i,n ₂

⎣⎢ ⎤

⎦⎥ ∈! ⁿ ^×(n ² ⁿ ³ ⁾ , T _{i, j} = D _z,i _{+ j} U _n

( 3 ) ^{⊗ y} ⁽ ^{i, j} ^{⊗ x} ⁱ ⁾

C ₁ T = T Λ ( _n ₁ ⊗ I _n ₂ _n ₃ ) ^{≡ T Λ} ¹ ^,

C ₂ T = T ⊕ ( ( _i=1 ⁿ ¹ Λ _i,n ₂ ) ^{⊗ I} ⁿ ³ ) ^{≡ T Λ} ² ^,

C ₃ T = T ⊕ ( _i=1 ⁿ ¹ ⊕ ⁿ _j ₌₁ ² Λ _{i, j,n} ₃ ) ^{≡ T Λ} ³

ψ _x = ı2 π ^k ⋅a ₁

n ₁ , D _x = diag 1,e ( ^ψ

^x

, !,e ⁽ⁿ

¹

⁻¹⁾ ^ψ

^x

) ^,

ψ _y,i = ı2 π

n ₂ k ⋅ a ₂ − a ₁ 2

⎛ ⎝⎜ ⎞

⎠⎟ − i 2

⎧ ⎨

⎩

⎫ ⎬

⎭ , D _y,i = diag 1,e ( ^ψ

^{y ,i}

, !,e ⁽ⁿ

²

⁻¹⁾ ^ψ

^{y ,i}

) ^,

ψ _{z,i+ j} = ı2 π

n ₃ k ⋅ a ₃ − a ₁ + a ₂ 3

⎛ ⎝⎜ ⎞

⎠⎟ −

i + j 3

⎧ ⎨

⎩

⎫ ⎬

⎭ , D _{z,i+ j} = diag 1,e ( ^ψ

^{y ,i+ j}

, !,e ⁽ⁿ

³

⁻¹⁾ ^ψ

^{y ,i+ j}

)

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the lattice sites, each connected to its six nearest neighbors by thin circular cylinders, which was proposed by Biswas et al.²⁵⁾Figure 1 shows the modified simple cubic lattice. It is noted that the original simple cubic structure comprising a lattice of rods has been fabricated recently with advanced silicon processing techniques.²⁶⁾ The second one is the tetragonal square spiral structure comprising a lattice of circular or square cylinders, which was proposed by Toader and John,^27,28)as shown in Fig. 2. Spiral structure was discussed previously by Chutinan and Noda.²⁹⁾The square spiral structure is arranged to connect the lattice points of diamond structure with specific order, and is amenable to the current technique of fabrication GLAD (GLancing Angle Deposition) as discussed in refs. 30 and 31. As a third example, we propose a diamond structure that has sp³-like configuration, composed of dielectric spheres with connecting spheroids, as shown in Fig. 3. Diamond structures are known to have large band gaps between relatively lower branches either in diamond network or inverse diamond structure.^12,28)In the present study, the spheroids, instead of circular cylinders, take the positions of ‘‘valence bonds’’ to imitate the sp³structure of the electrons of diamond atoms.

§5.

Demo performance

(21)

**CPU Times for T*p and Tq with FCC**

15 Eigendecomposition of Double-Curl Operator 19

0 2 4 6 8 10

x 10 ⁷ 0

1 2 3 4 5 6 7 8 9 10

CPU times (sec.)

n Tq

T ^* p

(a) n

0 1 2 3 4 5 6 7 8

x 10 ⁸ 0

1 2 3 4 5 6 7 8 9 10

CPU times (sec.)

n log(n) Tq

T ^* p

(b) n log(n) Fig. 6.1. CPU time for computing T ^⇤ p and T q with various n.

X U L G X W K

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

frequency

Fig. 6.2. Band structure of the 3D photonic crystals with FCC lattice. The vectors k’s along the boundary of first Brillouin zone. The frequency ! = a p

/(2⇡) is shown on the y-axis. The radius of the sphere is r = 0.12a and the connecting spheroid has minor axis length s = 0.11a.

in MATLAB are used for the IPL and CG methods, respectively. The stopping tolerances for eigs and pcg are set to be 10 ⁴ ⇥ ✏/(2

q 2

x + _y ² + _z ² ) and ^✏ 2

x ⇥ 10 ³ , respectively, where ✏ ⇡ 2 ⁵² is the floating-point relative accuracy in MATLAB. The maximal number of Lanczos vectors for the restart in eigs is set to be 20. All computations are carried out on a workstation with two Intel Quad-Core Xeon X5687 3.6 GHz CPUs, 48 GB main memory, the RedHat Linux operation system, and IEEE double-precision floating-point arithmetic operations.

Figure 6.1(a) shows the timing results for computing T ^⇤ p and T q by Algorithm 2 and 3, respectively. The matrix size of T ranges from 884, 736 to 94, 818, 816. In particular, the dimension of T is ¯ n ³ _j , where ¯ n _j = 96 + 24j = n ₁ = n ₂ = n ₃ for j = 0, 1, . . . , 15. The average CPU time out of ten trials for each j is then plotted in the figure. We can see Algorithms 2 and 3 are extraordinarily efficient. They take less than 10 seconds to finish a T ^⇤ p or T q matrix-vector multiplication even for the matrix T whose dimension is as large as 95 million. Figure 6.1(b) shows that the complexity of T ^⇤ p and T q is O(n log(n)).

Being equipped with these fast T ^⇤ p or T q computational kernels, we evaluate how the IPL method (Algorithm 1) performs, in terms of CPU time and iteration numbers, to solve the eigenvalue problems for the band structure of the target photonic crystals.

MATLAB **T*p : ﬀt**

Tq : iﬀt

(22)

Solving preconditioning linear system

16 (C ^∗ C − τ ^{I )y} ^{= d}

(23)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

(C ^∗ C − τ ^{I )y} ^{= d}

(24)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

(25)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ CG = 0 _GG ^∗ _y _{= −} _τ ⁻¹ _GG ^∗ _d

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

(26)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ CG = 0 _GG ^∗ _y _{= −} _τ ⁻¹ _GG ^∗ _d

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d −} ^τ ⁻¹ ^GG ^∗ ^d

(27)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ CG = 0 _GG ^∗ _y _{= −} _τ ⁻¹ _GG ^∗ _d

C ₁ T = T Λ ₁ , C ₂ T = T Λ ₂ , C ₃ T = T Λ ₃ Λ _q = Λ ₁ ^∗ Λ ₁ + Λ ^∗ ₂ Λ ₂ + Λ ^∗ ₃ Λ ₃

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d −} ^τ ⁻¹ ^GG ^∗ ^d

(28)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ CG = 0 _GG ^∗ _y _{= −} _τ ⁻¹ _GG ^∗ _d

C ₁ T = T Λ ₁ , C ₂ T = T Λ ₂ , C ₃ T = T Λ ₃ Λ _q = Λ ₁ ^∗ Λ ₁ + Λ ^∗ ₂ Λ ₂ + Λ ^∗ ₃ Λ ₃

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d −} ^τ ⁻¹ ^GG ^∗ ^d

I ₃ ⊗ Λ _q − τ ^I

( ) ^y! ^{= I −} ^τ ⁻¹ ^Λ ^Λ ¹ ²

Λ ₃

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ ⎡ Λ ₁ ^* Λ ₂ ^* Λ ^* ₃

⎣⎢ ⎤

⎦⎥

⎛

⎝

⎜ ⎜

⎜

⎞

⎠

⎟ ⎟

⎟ ( I ₃ ⊗T ) ^* ^{d, y} ^{= I} ( ³ ^⊗T ) ^y!

(29)

Solving preconditioning linear system

16 C ^∗ C = I ₃ ⊗ G ( ) ^∗ G ^{− GG} ^∗

G = [C ₁ ^⊤ ,C ₂ ^⊤ ,C ₃ ^⊤ ] ^⊤

C =

0 −C ₃ C ₂ C ₃ 0 −C ₁

−C ₂ C ₁ 0

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ CG = 0 _GG ^∗ _y _{= −} _τ ⁻¹ _GG ^∗ _d

C ₁ T = T Λ ₁ , C ₂ T = T Λ ₂ , C ₃ T = T Λ ₃ Λ _q = Λ ₁ ^∗ Λ ₁ + Λ ^∗ ₂ Λ ₂ + Λ ^∗ ₃ Λ ₃

(C ^∗ C − τ ^{I )y} ^{= d}

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d + GG} ^∗ ^y

I ₃ ⊗ (G ^∗ G) − τ ^I

{ } ^y ^{= d −} ^τ ⁻¹ ^GG ^∗ ^d

I ₃ ⊗ Λ _q − τ ^I

( ) ^y! ^{= I −} ^τ ⁻¹ ^Λ ^Λ ¹ ²

Λ ₃

⎡

⎣

⎢ ⎢

⎢

⎤

⎦

⎥ ⎥

⎥ ⎡ Λ ₁ ^* Λ ₂ ^* Λ ^* ₃

⎣⎢ ⎤

⎦⎥

⎛

⎝

⎜ ⎜

⎜

⎞

⎠

⎟ ⎟

⎟ ( I ₃ ⊗T ) ^* ^{d, y} ^{= I} ( ³ ^⊗T ) ^y!

Demo performance

(30)

Preconditioner

Iterative solver with preconditioner M: 

sol = bicgstabl(coef_mtx, rhs, tol, maxit, @(vec)FFT_based_precond(vec, Lambda, tau, EigDecompDoubCurl_cell, fun_mtx_TH_prod_vec,

fun_mtx_T_prod_vec)); 

17 M = C ^∗ C − τ ^I

(31)

Preconditioner

Iterative solver with preconditioner M: 

sol = bicgstabl(coef_mtx, rhs, tol, maxit, @(vec)FFT_based_precond(vec, Lambda, tau, EigDecompDoubCurl_cell, fun_mtx_TH_prod_vec,

fun_mtx_T_prod_vec)); 

17 M = C ^∗ C − τ ^I

Since 

we have 

No need to compute the matrix-vector multiplication involving A: 

sol = bicgstabl(@(vec)mtx_prod_vec_shift_invert_LS(vec, tau, Lambda_new, EigDecompDoubCurl_cell, mtx_B_sigma, fun_mtx_TH_prod_vec,

fun_mtx_T_prod_vec), rhs, tol, maxit); 

M ⁻¹ (A − σ ^B) = M ⁻¹ (A- τ Î + τ Î − σ ^B) = I + M ⁻¹ ( τ Î − σ ^B)

I + M ⁻¹ ( τ I − σ B)

{ } ^y ^{= M} ⁻¹ ^b

(32)

Challenge in Solving Linear System

SC lattice (dim = 46875)

FCC lattice

18

(33)

Null-space free

eigenvalue problem

19

(34)

Huge zero eigenvalues

Eigen-decomposition 

where 

is unitary and

20 Q ₀ Q

⎡ ⎣ ⎤

⎦ := I ( ³ ⊗T ) ^⎡ _⎣ ^Π ⁰ ^Π ¹ ^⎤ _{⎦ ≡ I} ( ³ ^⊗T ) ^Π ^Π ^0,2 ^0,1 ^Π ^Π ^1,3 ^1,1 ^Π ^Π ^1,4 ^1,2

FAME-matlab Package: Fast Algorithm for Maxwell Equations