The computation of photonic band gaps Christian Wieners Institut fr - - PowerPoint PPT Presentation

the computation of photonic band gaps
SMART_READER_LITE
LIVE PREVIEW

The computation of photonic band gaps Christian Wieners Institut fr - - PowerPoint PPT Presentation

The computation of photonic band gaps Christian Wieners Institut fr Angewandte und Numerische Mathematik, Karlsruhe www.kit.edu Light propagation in periodic structures We consider light propagation in periodic structures. We assume for


slide-1
SLIDE 1

www.kit.edu

The computation of photonic band gaps

Christian Wieners

Institut für Angewandte und Numerische Mathematik, Karlsruhe

slide-2
SLIDE 2

Light propagation in periodic structures

We consider light propagation in periodic structures. We assume for simplicity that the electro-magnetic fields are described by the linear Maxwell system; the permitivity µ = µ0 is constant; the permeability ε is periodic; the structure is infinite without boundary conditions; we consider only the time-harmonic case without sources. Optics Express, 16, 14812 (2008)

1

slide-3
SLIDE 3

Outline

  • 1. The Mathematics of Photonic Band Gaps (revisited)
  • 2. The Approximation of Photonic Band Gaps (revisited)
  • 3. The Computation of Photonic Band Gaps (state-of-the-art)
  • 4. The Verification of Photonic Band Gap (summarized)
  • 5. The Computation of Photonic Band Gaps (alternative approach)

The results are contained in the PhD thesis of A. Bulovyatov and

  • V. Huang, M. Plum, and C. Wieners: A computer-assisted proof for photonic band gaps.
  • Z. angew. Math. Phys. (ZAMP) 60 (2009) 1035-1052
  • C. Wieners: Calculation of the photonic band structure In: Photonic Crystals. Mathematical

Analysis and Numerical Approximation (ed. W. Dörfler) Birkhäuser, New York (2011) 40-62

  • C. Wieners, J. Xin: Boundary Element Approximation for Maxwell’s Eigenvalue Problem.

Mathematical Methods in the Applied Sciences 36 (2013) 2524-2539

2

slide-4
SLIDE 4

The Maxwell Eigenvalue Problem

We consider the Maxwell eigenvalue problem for the magnetic field H ∇ × ε−1∇ × H = λH (1a) ∇ · H = (1b) in the case that ε is a periodic function with ε(x) = ε(x + z) for all z ∈ Z3 and µ ≡ 1. Let Ω = (0, 1)3 be a fundamental cell, and K = [−π, π]3 be the Brillouin zone. The ansatz H(x) = eik·x H(x) with x ∈ Ω, k ∈ K and H periodic yields ∇k × ε−1∇k × H = λ H (2a) ∇k · H = (2b) in R3/Z3, where ∇k = ∇ + ik. The reduced problem for k ∈ K has a discrete spectrum with eigenvalues 0 ≤ λk,1 ≤ λk,1 ≤ · · · ≤ λk,N ≤ · · · Theorem Problem (1) has the spectrum σ =

  • n∈N

[ inf

k∈K

λk,n, sup

k∈K

λk,n].

(application of the Floquet-Bloch theory, for photonic crystals see Kuchment and Figotin)

If (supk∈K λk,n, inf k∈K λk,n+1) is non-empty for some n ∈ N, this is a band gap.

3

slide-5
SLIDE 5

Approximation of the Maxwell Eigenvalue Problem

Let X = Hper(curl; Ω) ⊂ H(curl; Ω) be the subspace of periodic vector fields. For u, v ∈ X we define the Hermitian bilinear forms ak(u, v) =

ε−1∇k × u · ∇k × v dx , m(u, v) =

u · v dx Let Q = H1

per(Ω) ⊂ H1(Ω) be the subspace of periodic scalar functions. Let

bk(v, q) =

v · ∇kq dx , ck(p, q) =

∇k ∇k ∇kp · ∇k ∇k ∇kq dx , v ∈ X , p, q ∈ Q , and Vk = {v ∈ X: bk(v, q) = 0 for all q ∈ Q}. Weak Form Find (u, λ) ∈ Vk × R such that ak(u, v) = λ m(u, v) for all v ∈ Vk . Discrete Approximation Let Xh,k ⊂ X and Qh,k ⊂ Q be discrete subspaces. Set Vh,k = {vh,k ∈ Xh,k : bk(vh,k, qh,k) = 0 for all qh,k ∈ Qh,k}. Find (uh,k, λh,k) ∈ Vh,k × R such that ak(uh,k, vh,k) = λh,k m(uh,k, vh,k) for all vh,k ∈ V V V h,k .

4

slide-6
SLIDE 6

Lowest Order Conforming Finite Elements

Let Ω =

  • T∈T h

T be a decomposition into cells T ∈ Th, and let ϕT : ˆ T = [0, 1]3 − → T be the transformation to the reference cell. In case of hexahedral cells we define Xh,0 = {u ∈ X: Dϕ−1

T u ◦ ϕT ∈ P0,1,1ex + P1,0,1ey + P1,1,0ez for all T ∈ Th},

Qh,0 = {q ∈ Q : q ◦ ϕT ∈ P1,1,1 for all T ∈ Th} . Let Vh ⊂ ¯ Ω be the set of all vertices z, and let Eh be the set of all edges e = (xe, ye) with midpoint ze = 0.5(xe + ye) and tangent te = ye − xe. A nodal basis {ψe,0 : e ∈ Eh} ⊂ Xh,0 exists such that for all uh,0 ∈ Xh,0 uh,0 =

  • e∈Eh

ψ′

e,0, uh,0ψe,0 ,

ψ′

e,0, uh,0 =

ye

xe

uh,0 · teds , and a nodal basis {φz,0 : z ∈ Vh} ⊂ Qh,0 exists such that for all qh,0 ∈ Qh,0 qh,0 =

  • z∈Vh

φ′

z,0, qh,0φz,0 ,

φ′

z,0, qh,0 = qh,0(z) . (curl conforming elements were introduced by Nédélec)

5

slide-7
SLIDE 7

Conforming Elements with Shifted Basis

For k ∈ K we define modified elements with a phase shift Xh,k = span{ψe,k : e ∈ Eh}, ψe,k(x) = e−ik·(x−ze)ψe,0(x) Qh,k = span{φz,k : z ∈ Vh}, φz,k(x) = e−ik·(x−z)φz,0(x) with uh,k =

  • e∈Eh

ψ′

e,k, uh,kψe,k ,

ψ′

e,k, uh,k =

ye

xe

eik·(x−ze)uh,0 · teds , qh,k =

  • z∈Vh

φ′

z,k, qh,kφz,k ,

φ′

z,k, qh,0 = qh,k(z) .

We have ∇k φz,k(x) = e−ik·(x−z)∇φz,0(x) ∇k × ψe,k(x) = e−ik·(x−ze)∇ × ψe,0(x) ∇k · ψe,k(x) = e−ik·(x−ze)∇ · ψe,0(x) ak(ψe1,k, ψe2,k) = eik·(ze1 −ze2 )a0(ψe1,0, ψe2,0) m(ψe1,k, ψe2,k) = eik·(ze1 −ze2 )m(ψe1,0, ψe2,0) bk(ψe,k, φz,k) = eik·(ze−z)b0(ψe,0, φz,0) ck(φz1,k, φz2,k) = eik·(z1−z2)c0(φz1,0, φz2,0) .

6

slide-8
SLIDE 8

Finite Element Convergence

The finite element spaces Xh,k, Qh,k, and Vh,k satisfy Ellipticity on Vh,k There exists C > 0 s.t. ak(uh,k, uh,k) ≥ C uh,k2

L2

for all uh,k ∈ Vh,k . Weak approximability of Q There exists ρ1(h) > 0, tending to zero as h goes to zero such that sup

vh,k∈Vh,k

bk(vh,k, q) vh,kcurl ≤ ρ1(h) qH1 for all q ∈ Q . Strong approximability of V For some r > 0 there exists ρ2(h) > 0, tending to zero as h goes to zero such that for any u ∈ V ∩ H1+r(Ω, C3) there exists uh,k ∈ Vh,k satisfying u − uh,kcurl ≤ ρ2(h) uH1+r .

(for coercivity and regularity see Dauge-Norton-Scheichl 2015)

Theorem Let (u, λ) a solution of the continuous eigenvalue problem, and let (uh,k, λh,k) be the corresponding discrete solution. Then, we have (uh,k, λh,k) → (u, λ) for h → 0 .

(Boffi-Conforti-Gastaldi 2006)7

slide-9
SLIDE 9

The Projection

For uh,k ∈ Xh,k we construct ph,k ∈ Qh,k such that uh,k − ∇kph,k ∈ Vh,k, i.e., = bk(uh,k − ∇kph,k, qh,k) = bk(uh,k, qh,k) − ck(ph,k, qh,k) for all qh,k ∈ Qh,k Thus, we have Bh,kuh,k = Ch,kph,k and uh,k = Sh,kph,k with operators "div": Bh,k : Xh,k − → Q′

h,k is defined by

Bh,kvh,k, φz,k = bk(vh,k, φz,k), z ∈ Vh "Laplace": Ch,k : Qh,k − → Q′

h,k is defined by

Ch,kqh,k, φz,k = ck(qh,k, φz,k), z ∈ Vh "grad": Sh,k : Qh,k − → Xh,k is given by nodal evaluation Sh,kqh,k =

  • e=(xe,ye)∈Eh
  • qh,k(ye)eik·(ye−ze) − qh,k(xe)eik·(xe−ze)

ψe,k. This defines a projection Ph,k = id −Sh,k ◦ C−1

h,k ◦ Bh,k : Xh,k −

→ Vh,k .

(special care is required for k = 0)

8

slide-10
SLIDE 10

Modified LOBPCG Method (including projection)

Let Th,k : X′

h,k −

→ Xh,k be a preconditioner for Aδ

h,k = Ah,k + δMh,k : X

X X h,k − → X′

h,k.

S0) Choose randomly u1

h,k, ..., uN h,k ∈ Xh,k. Compute vn h,k = Ph,kun h,k ∈ Vh,k.

S1) Ritz-step: Set up Hermitian matrices ˆ A =

  • ak(vm

h,k, vn h,k)

  • m,n=1,...,N ,

ˆ M =

  • m(vm

h,k, vn h,k)

  • m,n=1,...,N ∈ CN×N

and solve the matrix eigenvalue problem ˆ Aˆ zn = λn ˆ Mˆ zn. S2) Compute yn

h,k = N

  • n=1

ˆ zn

mvn h,k ∈ Vh,k.

S3) Compute rn

h,k = Ah,kyn h,k − λnMh,kyn h,k ∈ X′ h,k, check for convergence.

S4) Compute un

h,k := Th,krn h,k ∈ Xh,k and wn h,k = Ph,kuh,k ∈ Vh,k.

S5) Perform Ritz-step for {v1

h,k, ..., vN h,k, w1 h,k, ..., wN h,k} ⊂ Vh,k of size 2N.

S6) Go to step S2). The full algorithm uses orthogonalization, new random vectors, and a Ritz-step of size 3N.

(the LOBPCG method was introduced by Knyazev 2001)

9

slide-11
SLIDE 11

Multigrid Preconditioner for the Maxwell Operator

Let X0,k ⊂ X1,k ⊂ · · · ⊂ XJ,k be finite element space of mesh size hj = 2−jh0. The multigrid preconditioner Tj,k : X′

j,k −

→ Xj,k is defined recursively: For j = 0, define T0,k =

0,k

−1 . For j > 0, the definition of Tj,k requires: 1) a prolongation operator Ij,k : Xj−1,k − → Xj,k; 2) the adjoint operator I′

j,k : X′ j,k −

→ X′

j−1,k;

3) a smoother Rj,k : X′

j,k −

→ Xj,k for Aj,k; 4) a smoother Dj,k : Q′

j,k −

→ Qj,k for Cj,k; 5) a transfer operator Sj,k : Qj,k − → Xj,k; 6) the adjoint transfer operator S′

j,k : X′ j,k −

→ Q′

j,k.

Now, define Tj,k by id −Aδ

j,kTj,k =

  • id −Aδ

j,kIj,kTj−1,kI′ j,k

  • id −δ−1Aδ

j,kSj,kDj−1,kS′ j,k

  • id −Aδ

j,kRj,k

  • .

(this multigrid variant was introduced and analyzed by Hiptmair 1998)

10

slide-12
SLIDE 12

Eigenvalue Computation for a Photonic Crystal

We compute Block-Floquet modes for k = (3, 1, −1) in the unit cube Ω = (0, 1)3 with periodic boundary conditions ∇k × ε−1∇k × H = λ H , ∇k · H = 0 .

Band structure of a photonic crystal Eigenfunction – (Bloch-Floquet mode)

(configuration from Dobson/Pasciak, Comp. Meth. Appl. Math. 1 (2001) 138–153)

11

slide-13
SLIDE 13

Photonic Crystals: Parallel Eigenvalue Computation

Calculation of the first 7 eigenvalues (≈ 1 million d.o.f.) N processor kernels Total time Speed up factor 32 5:11 min. 1.79 64 2:53 min. 1.85 128 1:33 min. 1.42 256 1:05 min. 1.39 512 0:47 min. Calculation of the first 7 eigenvalues on 512 processor kernels Refinement level d.o.f. Total time Scale up factor j = 4 13 872 00:18 min. 1.32 j = 5 104 544 00:24 min. 1.91 j = 6 811 200 00:47 min. 3.98 j = 7 6 390 144 03:08 min. 6.85 j = 8 50 725 632 21:31 min. (every refinement level increases the problem size by a factor of 8)

12

slide-14
SLIDE 14

Photonic Crystals: Eigenvalue Convergence

Calculation of the first 7 eigenfunctions (up to ≈ 50 million d.o.f.) j λj

1

λj

2

λj

3

λj

4

λj

5

λj

6

λj

7

3 3.4818 3.7189 8.4285 9.0216 13.4387 13.7301 14.2479 4 3.4449 3.6828 7.9698 8.5039 12.4049 12.8074 13.3287 5 3.4310 3.6689 7.8410 8.3580 12.0515 12.5959 13.1179 6 3.4258 3.6637 7.8026 8.3145 11.9466 12.5339 13.0613 7 3.4239 3.6618 7.7906 8.3008 11.9136 12.5156 13.0451 8 3.4232 3.6612 7.7868 8.2964 11.9027 12.5100 13.0402 Convergence of the first 7 eigenfunctions meassured by |λj+1

n

− λj

n|

|λj

n − λj−1 n

|

|λ4

n − λ3 n|

0.0369 0.0361 0.4587 0.5177 1.0338 0.9227 0.9192 2.65 2.60 3.56 3.55 2.93 4.36 4.36 |λ5

n − λ4 n|

0.0139 0.0139 0.1288 0.1459 0.3534 0.2115 0.2108 2.67 2.67 3.35 3.35 3.37 3.41 3.72 |λ6

n − λ5 n|

0.0052 0.0052 0.0384 0.0435 0.1049 0.0620 0.0566 2.74 2.74 3.20 3.18 3.18 3.39 3.49 |λ7

n − λ6 n|

0.0019 0.0019 0.0120 0.0137 0.0330 0.0183 0.0162 2.71 3.17 3.16 3.11 3.03 3.27 3.31 |λ8

n − λ7 n|

0.0007 0.0006 0.0038 0.0044 0.0109 0.0056 0.0049

13

slide-15
SLIDE 15

Photonic Crystals: Eigenmodes

Illustration of an isosurface of the field intensity and the field directions.

14

slide-16
SLIDE 16

Band Gap Verification

Analytically, it can be shown that photonic crystals with band gaps exist (by asymptotic investigations). Numerically,

  • ne can check a given photonic crystals

for band gap candidates. Our aim is—combining analysis and numerics—to provide a rigorous proof for the existence of a band gap for a specific test configuration. This requires the following techniques: Close finite element approximations of eigenmodes; Upper spectral bounds (Rayleigh-Ritz); Lower spectral bounds (Weinstein, Lehmann/Goerisch); Spectral exclosure for perturbed operators (Hoang, Plum); Spectral homotopy (Goerisch, Plum).

15

slide-17
SLIDE 17

Polarized Waves in 2D

Now, we consider the Maxwell eigenvalue problem for the electic field E ∇ × ∇ × E = λεE , ∇ · εE = 0 in special case that ε is periodic in x1 and x2, constant in x3 and that the wave E = (0, 0, u) is polarized, which leads to the Helmholtz problem −∆u = λεu . We compute in the periodicity cell Ω for finitely many k in the Brillouin zone (uk,n, λk,n) ∈ H1

per(Ω) × R:

∇kuk,n · ∇kv dx = λk,n

εuk,nv dx , v ∈ H1

per(Ω)

(n = 1, ..., N) with guaranteed accuracy. We define ak(u, v) =

∇ku · ∇kv dx , (u, v)ε =

εuv dx .

16

slide-18
SLIDE 18

A Candidate

We set ε(x) = 1 for x ∈ [1/16, 15/16]2 and ε(x) = 5 else. By symmetry we have the same spectrum for k = (k1, k2), (π − k1, k2), (k1, π − k2), ... (−π, −π) (π, −π) (π, π) (−π, π) Brillouin zone K periodic ε distribution

17

slide-19
SLIDE 19

A Candidate

(0, 0) (π/2, 0) (π, 0) (π, π/2) (π, π) (π/2, π/2) (0, 0) λk,n 35 30 25 20 band gap 15 10 5 eigenvalues λk,1, ..., λk,6 along the test path k ∈ P ⊂ K

18

slide-20
SLIDE 20

A Candidate

  • 4
  • 3
  • 2
  • 1

1 2 3 4

  • 4
  • 3
  • 2
  • 1

1 2 3 4 5 10 15 20 25 30

eigenvalue λk,3 for all k ∈ K

19

slide-21
SLIDE 21

A Candidate

  • 4
  • 3
  • 2
  • 1

1 2 3 4

  • 4
  • 3
  • 2
  • 1

1 2 3 4 5 10 15 20 25 30

eigenvalue λk,4 for all k ∈ K

20

slide-22
SLIDE 22

A Candidate

5 10 15 20 25 30

  • 4
  • 3
  • 2
  • 1

1 2 3 4

  • 4
  • 3
  • 2
  • 1

1 2 3 4

eigenvalues λk,1, ..., λk,5 for all k ∈ K

21

slide-23
SLIDE 23

A Candidate

Bloch modes uk,1, ..., uk,6 for k = (π, π)

22

slide-24
SLIDE 24

Upper Spectral Bounds (Ritz-Galerkin)

Theorem For eigenfunction approximations ˜ uk,1, ..., ˜ uk,N ∈ H1(Ω/Λ, C) define A =

  • ak(˜

uk,m, ˜ uk,n)

  • m,n=1,...,N ,

B =

uk,m, ˜ uk,n)ε

  • m,n=1,...,N ∈ CN,N

and let ψk,1 ≤ ψk,2 ≤ · · · ≤ ψk,N be the eigenvalues of the problem A x = ψ B x. Then we have λk,n ≤ ψk,n , n = 1, ..., N .

23

slide-25
SLIDE 25

Lower Spectral Bounds (Lehmann-Goerisch)

Theorem Let γ > 0. For eigenfunction approximations ˜ uk,n ∈ V and scaled dual approximations ˆ σk,n ≈

1 γ+˜ λk,n ∇k ˜

uk,n ∈ W define S =

σk,m, ˆ σk,n)

  • m,n=1,...,N ,

T = 1 γ

uk,m + 1 ε ∇k · ˆ σk,m, ˜ uk,n + 1 ε ∇k · ˆ σk,n)ε

  • m,n=1,...,N .

If β > γ satisfies β − γ ≤ λk,N+1, if N = A + (γ − 2β)B + β2(S + T T T) is positive definite, and if the eigenvalues θ1 ≥ θ2 ≥ · · · ≥ θN

  • f the eigenvalue problem
  • A + (γ − β)B
  • x = θ N x

are negative, we have the lower eigenvalue bound β − γ − β 1 − θn ≤ λn,k, n = 1, ..., N. Remark: A suitable β is obtained by a spectral homotopy.

24

slide-26
SLIDE 26

Eigenvalue Exclosure by Perturbation Analysis

Theorem For some λgap, n > 1, and k ∈ K assume that δ > 0 exits such that λk,n−1 + δ < λgap < λk,n − δ . Then we have λk+h,n−1 < λgap < λk+h,n for all h ∈ R2 with |h| <

  • min

x∈Ω ε(x)

  • λk,n + δ −
  • λk,n
  • .

For the application we compute a guaranteed lower bound λk,n,inf ≤ λk,n and upper bounds λk,n−1,sup ≥ λk,n−1, λk,n,sup ≥ λk,n, and rk =

  • min

x∈Ω ε(x)

  • λk,n,sup + δ −
  • λk,n,sup
  • .

This ensures that the interval (λk,n−1,sup + δ, λk

k k,n,inf − δ) is contained in the resolvent

set of the elliptic operator associated to the eigenvalue problem for |k′ − k| < rk.

25

slide-27
SLIDE 27

A Verified Band Gap

This figure illustrates Brillouin vectors k′ ∈

k∈K

Ball(k, rk) where λgap ∈ (18.2, 18.25) is not a spectral value. By symmetry this proves that for all k ∈ K the existence of a band gap (18.2, 18.25) ⊂ (λmax,3, λmin,4) for the eigenvalue problem −∆u = λεu in R2.

The proof requires the close approximation of more than 5 000 eigenvalues and eigenfunctions (for 100 vectors k ∈ K and for a homotopy to bound the rest of the spectrum) and takes about 90 h computing time.

26

slide-28
SLIDE 28

The Representation Formula for Maxwell’s Equations

The fundamental solution for the Maxwell the cavity problem ∇ × ∇ × u − κ2u = 0 and ∇ · u = 0 is given by Gκ(x, y) = Eκ(x, y)I + κ−2∇x

  • ∇yEκ(x, y)

T , Eκ(x, y) = exp(iκ|x − y|) 4π|x − y| . This defines the Maxwell single-layer and double-layer potentials Ψκ

SL(v)(x)

= −

  • Γ

Eκ(x, y)v(y) dsy − 1 κ2 ∇x

  • Γ

Eκ(x, y) divΓ(v(y)) dsy Ψκ

DL(w)(x)

= −∇x × Vκ(w)(x) = −∇x ×

  • Γ

Eκ(x, y)w(y) dsy , x ∈ Ω for tangential boundary data v, w where v · n = w · n = 0. The potentials are Maxwell solutions, and the representation formula u(x) = Ψκ

SL(n × (∇ × u))(x) + Ψκ DL(n × u)(x)

x ∈ Ω , determines the Maxwell solution u in terms of its boundary traces.

27

slide-29
SLIDE 29

Trace Operators and Function Spaces

Let Γ be composed of smooth manifolds Γj with exterior normals nj such that Γ =

j Γj and Γj ∩ Γm = ∅, j = m. For smooth vector fields v in Ω and x ∈ Γj define

γt(v)(x) = lim

y∈Ω→x∈Γj nj(x) × v(y) ,

γκ

N(v)(x) = 1

κγt(∇ × v)(x) . Let ejm = ∂Γj ∩ ∂Γm be a boundary edge with tangent tjm, and let n∂Γj ∈ affine(Γj) be the tangential normal along ∂Γj, i.e., n∂Γj = nj × tjm on the edge ejm. From γt(v)|¯

Γj · n∂Γj = (nj × v|¯ Γj ) · (nj × tjm) = tjm · v|¯ Γj

  • n ejm

(3) we obtain

  • j
  • ∂Γj

γt(v)|¯

Γj · n∂Γj w dl =

  • j=m
  • ejm

tjm · v|¯

Γj w dl = 0. This yields

∇ × v · ∇w dx = −

  • Γ

divΓ γt(v)w ds . (4) For v ∈ H(curl, Ω) the weak Dirichlet trace operator γt(v) ∈ H−1/2(Γ) is defined by

  • γt(v), w|Γ
  • =

∇ × v · w dx −

v · ∇ × w dx , w ∈ C∞(Ω) . Extending (4) to v ∈ H(curl, Ω) shows divΓ γt(v) ∈ H−1/2(Γ) defined by

  • divΓ γt(v), w|Γ
  • = −

∇ × v · ∇w dx , w ∈ C∞(Ω) .

28

slide-30
SLIDE 30

Boundary Integral Operators

On the trace space W−1/2(Γ) = γt

  • H(curl, Ω)
  • ⊂ H−1/2(Γ) we define

Sκ =γtΨκ

SL : W−1/2(Γ) → W−1/2(Γ) ,

1 2I + Cκ =γtΨκ

DL : W−1/2(Γ) → W−1/2(Γ) .

Note that we have the identities γtΨκ

SL = γκ NΨκ DL and γtΨκ DL = γκ NΨκ SL.

For smooth tangential vector fields we define the anti-linear pairing v, wτ,Γ =

  • Γ

(v × n) · w ds which extends to v ∈ H(curl2, Ω) and v ∈ H(curl, Ω) by

  • γκ

N(v), γt(w)

  • τ,Γ = 1

κ

  • ∇ × ∇ × v · w − ∇ × v · ∇ × w
  • dx .

For tangential vector fields v, w ∈ L∞(Γ) with divΓ v, divΓ w ∈ L∞(Γ)

  • Sκ(v), w
  • τ,Γ =
  • Γ
  • Γ

1 κ divΓ v(y) divΓ w(x) − κv(y) · w(x)

  • Eκ(x, y) dsydsx ,
  • Cκ(w), v
  • τ,Γ = −
  • Γ
  • Γ

∇xEκ(x, y) ·

  • w(y) × v(x)
  • dsydsx .

... more details by Buffa, Costabel, Hiptmair, Schwab, ...

29

slide-31
SLIDE 31

The Maxwell Eigenvalue Problem

Find a non-trivial solution pair (u, κ) with u × n = 0 on Γ and ∇ × ∇ × u + κ2u = 0 and ∇ · u = 0 in Ω . A suitable normalization of the eigenfunction u is required. (BIE) Find (σ σ σ, κ) ∈ W−1/2(Γ) × R+ with σ σ σ = 0 satisfying Sκ(σ σ σ),χ χ χτ,Γ = 0 , χ χ χ ∈ W−1/2(Γ) . (BEM) Find (σ σ σh, κh) ∈ W−1/2

h

(Γ) × R+ with σ σ σh = 0 such that Sκh(σ σ σh),χ χ χhτ,Γ = 0 , χ χ χh ∈ W−1/2

h

(Γ) . ... more details by Steinbach and Unger ... Note that S(κ) = Sκ : Λ ⊂ C → L(W −1/2(Γ), W −1/2(Γ)) is holomorphic and satisfies a generalized Garding inequality.

30

slide-32
SLIDE 32

The Contour Integral Method

Eigenvalues are singular points λ ∈ C of the holomorphic matrix function A: C − → CN×N. Thus, the matrix function A(·)−1 is meromorphic. Fix D ⊂ C such that A(z)−1 is regular for z ∈ ∂D, and assume that λ1, ..., λn are the eigenvalues in D with eigenvectors ξj ∈ CN, i.e., A(λj)ξj = 0 and ξH

j ξj = 1.

If all eigenvalues λj are simple, we have A(z)−1 =

n

  • j=1

1 z − λj ξjξH

j + B(z) ,

where B(z) is a holomorphic matrix function in D. The residual theorem yields for the contour integral matrices A0 := 1 2πi

  • ∂D

A(z)−1R dz =

n

  • j=1

ξjξH

j R ,

A1 := 1 2πi

  • ∂D

zA(z)−1R dz =

n

  • j=1

λjξjξH

j R ,

where R ∈ CN×n are suitable test matrices such that A0 and A1 have full rank. By construction we have A1A+

0 = n j=1 λjξjξH j ∈ CN×N which allows to recover the

eigenvalues from A+

0 A1 ∈ Cn×n.

... more details by W.-D. Beyn ...

31

slide-33
SLIDE 33

A Domain Decomposition Method

Let Ω =

P

  • p=1

Ωp with permittivity εp and permeability µp in Ωp. In Ωp solve ∇ × ∇ × up − (ω/cp)2up = 0 and ∇ · up = 0 . Set βp =

  • εp/µp and Γp = ∂Ωp. On Γpq = Γp ∩ Γq holds

γp

t (up) + γq t (uq) = 0 ,

βpγ

κp,p N

(up) + βqγ

κq,q N

(uq) = 0 . Dirichlet traces on interfaces and Neumann traces on subdomain boundaries Find ω, ϕ = (ϕpq)p<q ∈

p<q W−1/2(Γpq) and σ

σ σ = (σ σ σp) ∈

p W−1/2(Γp) with

  • q

1

2I − CΓp ω/cp

  • (ϕpq), τ p

τ,Γp +

  • SΓp

ω/cp(σ

σ σp), τ p

τ,Γp

= for all τ p ,

  • p<q
  • βpSΓp

ω/cp + βqSΓq ω/cq

  • (ϕpq),χ

χ χpq

τ,Γpq

+

  • βp

1

2 + CΓp ω/cp

σ σp) + βq 1

2 + CΓq ω/cq

σ σq),χ χ χpq

τ,Γpq

= for all χ χ χpq . ... more details by Langer, Of, Steinbach, Zulehner ...

32

slide-34
SLIDE 34

A Test Example with two subdomains

σ σ σ2

h

ϕ12

h

σ σ σ1

h

33

slide-35
SLIDE 35

Band Structure Computation for Photonic Crystals

We consider a periodic medium with ε(x + z) = ε(x) for all z ∈ Z3 and µ ≡ µ0. As a consequence of the Floquet theory, the spectrum of the Maxwell operator is the union of all eigenvalues in the periodicity cell Ω = (0, 1)3 with quasi-periodic boundary conditions γt(u)(x + ej) + exp(ik · ej)γt(u)(x) = 0 , γκ

N(u)(x + ej) + exp(ik · ej)γκ N(u)(x) = 0

  • n the faces G1 = [0, 1]2 × {0}, G2 = [0, 1] × {0} × [0, 1] and G3 = {0} × [0, 1]2,

where ej are the unit vectors and k ∈ (−π, π]3 is a parameter in the Brillouin zone. For two subdomains (with ε1 = 1 and ε2 = 13) this leads to the problem to find (ϕ12

h ,σ

σ σ2

h|Γ,σ

σ σ2|Γ12, ϕ2

h|Γ,σ

σ σ1) ∈ W

− 1

2

h

(Γ12) × W

− 1

2

h

(Γ) × W

− 1

2

h

(Γ12) × W

− 1

2

h

(Γ) × W

− 1

2

h

(Γ1) and ωh such that       A11(ωh) A12(ωh)Bk A13(ωh) A14(ωh)Bk A15(ωh) BH

k A21(ωh)

BH

k A22(ωh)Bk

BH

k A23(ωh)

BH

k A24(ωh)Bk

A31(ωh) A32(ωh)Bk A33(ωh) A34(ωh)Bk BH

k A41(ωh)

BH

k A42(ωh)Bk

BH

k A43(ωh)

BH

k A44(ωh)Bk

A51(ωh) A55(ωh)       is singular, where the quasi-periodic boundary conditions are realized by the matrix Bk extending the degrees of freedom on G1 ∪ G2 ∪ G3 to Γ = ∂Ω.

34

slide-36
SLIDE 36

Band Structure Computation for Photonic Crystals

We consider a periodic medium with ε(x + z) = ε(x) for all z ∈ Z3 and µ ≡ µ0. As a consequence of the Floquet theory, the spectrum of the Maxwell operator is the union of all eigenvalues in the periodicity cell Ω = (0, 1)3 with quasi-periodic boundary conditions γt(u)(x + ej) + exp(ik · ej)γt(u)(x) = 0 , γκ

N(u)(x + ej) + exp(ik · ej)γκ N(u)(x) = 0

  • n the faces G1 = [0, 1]2 × {0}, G2 = [0, 1] × {0} × [0, 1] and G3 = {0} × [0, 1]2,

where ej are the unit vectors and k ∈ (−π, π]3 is a parameter in the Brillouin zone. boundary elements finite elements band structure along a path in the Brillouin zone (−π, π]3 for a photonic crystal

35