[PPT] - Recent Advances on Mathematical Analysis and Simulations of PowerPoint Presentation

SLIDE 1

Recent Advances on Mathematical Analysis and Simulations of Invisibility Cloaks with Metamaterials

Jichun Li University of Nevada Las Vegas (UNLV)

Jichun Li (UNLV) June 27, 2018 1 / 62

SLIDE 2

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 2 / 62

SLIDE 3

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 3 / 62

SLIDE 4

Invisibility cloak with metamaterials

Science, vol.312 (June 23, 2006): “Controlling Electromagnetic Fields” (by J.B. Pendry, D.

Schurig, D.R. Smith) [Cited 4510 times as 2/25/15; 6115 times as 4/8/17; 6404 times as 8/26/17; 6785 times as 3/15/18]

Science, vol.312 (June 23, 2006): “Optical Conformal Mapping” (by Ulf Leonhardt). [Cited 2369

times as 2/25/15; 3110 times as 4/8/17; 3250 times as 8/26/17; 3401 times as 3/15/18]

Science, vol.314 (Nov. 10, 2006): “Metamaterial Electromagnetic Cloak at Microwave

Frequencies” (by Schurig, Mock, etc.) [Cited 3689 times as 2/25/15; 5037 times as 4/8/17; 5289 times as 8/26/17; 5632 times as 3/15/18]

Science, vol.328 (Apr. 16, 2010): “Three-Dimensional Invisibility Cloak at Optical Wavelengths”

(by Ergin, Stenger, Brenner, Pendry, Wegener) [Cited 918 times as 3/15/18]

Greenleaf, Lassas and Uhlmann (2003) for conductivity equations; Approx/near cloaking:

Milton, Nicorovici (May 3, 2006), Bouchitte, Schweizer (2010), Ammari etc (2013), Kohn, Weinstein etc (2008, 2014), G. Bao, J. Zou, H.Y. Liu, J.Z. Li, ...

Numer. methods: J. Li, Huang, Yang (2011, 2012,...), Z. Xie, C. Chen, etc (CiCP2016), Li-Lian

Wang etc (CiCP2015, CMAME2016), D. Liang etc (JSC2016), Brenner, Gedicke, Sung (JSC2016, M2AS2017), D. Liang etc (JSC2016), J. Liu (OP2014), ...

Jichun Li (UNLV) June 27, 2018 4 / 62

SLIDE 5

Figure 1 : (A) The simulation of the cloak with the exact material properties, (B) the simulation of the cloak with the reduced material properties, (C) the experimental measurement of the bare conducting cylinder, and (D) the experimental measurement of the cloaked conducting cylinder. Source: D. Schurig et al, Science, V.314, Nov. 2006, 977-980. Invisible to an incident plane wave at 8.5 GHz.

Jichun Li (UNLV) June 27, 2018 5 / 62

SLIDE 6

Figure 2 : 2D microwave cloaking structure (background image) with a plot of the material parameters that are implemented. Source: D. Schurig et al, Science, V.314, Nov. 2006, 977-980. Reduced parameters: εz = (

b b−a)2,µr = ( r−a r )2,µθ = 1. Exact parameters:

εz = (

b b−a)2 r−a r ,µr = r−a r ,µθ = r r−a

Jichun Li (UNLV) June 27, 2018 6 / 62

SLIDE 7

Short math history

IMA Hot Topics Workshop on ”Negative Index Materials” during

Oct.2-4, 2006.

CSCAMM at U. of Maryland’s workshop “Electromagnetic

Metamaterials and Their Approximations: Practical and Theoretical Aspects”, Sept. 22-25, 2008. (Organized by E. Tadmor, G. Uhlmann,

M. Vogelius etc)
UCLA’s IPAM workshop on ”Metamaterials: Applications, Analysis

and Modeling”, Jan. 25-29, 2010. (Organized by S. Brenner, M.-C. Calderer, R. Kohn, J. Li, G. Milton, C.-W. Shu, etc)

Tsinghua Sanya International Mathematics Forum (TSIMF)’s

workshop ”Mathematical Analysis of Metamaterials and Applications”, Dec.5-9, 2016. (Organized by J. Li, P . Monk, Y. Huang)

Books on DNG: Caloz, Itoh (2005); Eleftheriades, Balmain (2005); Engheta, Ziolkowski (2006),

Shalaev (2007), Krowne, Zhang (2007), Marques, Martin, Sorolla (2008), Markos, Soukoulis (2008), Y. Hao, R. Mittra (2008), W. Cai, V. Shalaev (2009), Cui, Smith, Liu (2010), U. Leonhardt,

T. Philbin (2010), R.V. Craster, S. Guenneau (2012), D. H. Werner, D.-H. Kwon (2013), ...

Jichun Li (UNLV) June 27, 2018 7 / 62

SLIDE 8

Form invariant property for Maxwell’s equations

Theorem 1

Under a coordinate transformation x′ = x′(x), the Maxwell’s equations ∇×E +jωµH = 0, ∇×H −jωεE = 0, (1) keep the same form in the transformed coordinate system: ∇′ ×E′ +jωµ′H′ = 0, ∇′ ×H′ −jωε′E′ = 0, (2) where all new variables are given by E′(x′) = A−TE(x), H′(x′) = A−TH(x), A = (aij), aij = ∂x′

i

∂xj , (3) and µ′(x′) = Aµ(x)AT/det(A), ε′(x′) = Aε(x)AT/det(A). (4)

Jichun Li (UNLV) June 27, 2018 8 / 62

SLIDE 9

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 9 / 62

SLIDE 10

Cylindrical cloak in time domain

The cloak modeling: Faraday′s Law : ∂B B B ∂t = −∇×E E E, Ampere′s Law : ∂D D D ∂t = ∇×H H H, constitutive relations : D D D = εE E E, B B B = µH H H, E E E and H H H: electric and magnetic fields, D D D and B B B: electric and magnetic flux densities, ε and µ: cloak permittivity and permeability. Cylindrical cloak: Pendry et al (Science 2006): εr = µr = r −R1 r , εφ = µφ = r r −R1 , εz = µz =

R2

R2 −R1 2 r −R1 r , R1 and R2: inner and outer radius of the cloak.

Jichun Li (UNLV) June 27, 2018 10 / 62

SLIDE 11

Cylindrical cloak in time domain

Transforming the polar coordinate system to the Cartesian coordinate system, and using the Drude model for the permittivity: εr(ω) = 1− ω2

p

ω2 −jωγ , we obtain ε0εφ ∂ 2 ∂t2 +γ ∂ ∂t +w2

p

E

E E = ∂ 2 ∂t2 +γ ∂ ∂t +w2

p

MAD

D D +εφ ∂ 2 ∂t2 +γ ∂ ∂t

MBD

D D, where we denote D = (Dx,Dy)′ and MA =

sin2 φ

−sinφ cosφ −sinφ cosφ cos2 φ

, MB =
cos2 φ

sinφ cosφ sinφ cosφ sin2 φ

.

Jichun Li (UNLV) June 27, 2018 11 / 62

SLIDE 12

Cylindrical cloak in time domain

Permeability using the Drude model: µz(ω) = A

1−

ω2

pm

ω2 −jωγm

, A =

R2 R2 −R1 , ωpm > 0 and γm ≥ 0: magnetic plasma and collision frequencies. ∂ 2 ∂t2 +γm ∂ ∂t

Bz = µ0A

∂ 2 ∂t2 +γm ∂ ∂t +ω2

pm

Hz.

In summary,

∂B B B ∂t = −∇×E E E, (5) ∂D D D ∂t = ∇×H H H, (6) ε0εφ ∂ 2 ∂t2 +γ ∂ ∂t +w2

p

E

E E = ∂ 2 ∂t2 +γ ∂ ∂t +w2

p

MAD

D D +εφ ∂ 2 ∂t2 +γ ∂ ∂t

MBD

D D, (7) ∂ 2 ∂t2 +γm ∂ ∂t

Bz = µ0A

∂ 2 ∂t2 +γm ∂ ∂t +ω2

pm

Hz.

(8)

Jichun Li (UNLV) June 27, 2018 12 / 62

SLIDE 13

Li, Huang, Yang: Math Comp, 2015

Assume that γ = γm. µ0Aε0εφ(Et3 +γEt2 +ω2

pEt)

= µ0A(MA +εφMB)∇×(Ht2 +γHt)+ µ0Aω2

pMA∇×H.

(9) To simplify the notation, we denote H = Hz, M = MA +εφMB. Also we have µ0A(Ht2 +γHt +ω2

pmH) = −∇×Et −γ∇×E.

(10) Taking curl of (10) and substituting into (9), we have µ0ε0Aεφ(Et3 +γEt2 +ω2

pEt)+M∇×∇×Et +γM∇×∇×E

= −µ0AM∇×(ω2

pmH)+ µ0Aω2 pMA∇×H.

(11)

Jichun Li (UNLV) June 27, 2018 13 / 62

SLIDE 14

Analysis of the model: cont’d

Lemma 1

Matrix MA is symmetric and non-negative definite, and M is SPD.

Lemma 2

For matrix MC = (MA +εφMB)−1, MC ·MA = MA holds true. Weak formulation: For any φ ∈ H0(curl;Ω),ψ ∈ L2(Ω), (i) ε0µ0A[(εφMCEt3,φ)+γ(εφMCEt2,φ)+(ω2

pεφMCEt,φ)]

+(∇×Et,∇×φ)+γ(∇×E,∇×φ) = −µ0A(ω2

pmH,∇×φ)+ µ0A(ω2 pMCMA∇×H,φ),

(12) (ii) µ0A

(Ht2,ψ)+γ(Ht,ψ)+(ω2

pmH,ψ)

=

−(∇×Et +γ∇×E,ψ). (13)

Jichun Li (UNLV) June 27, 2018 14 / 62

SLIDE 15

Analysis of the model: cont’d

Theorem 1

For the solution of (12)-(13), the following stability holds true: ε0µ0A[(εφMcEt2,Et2)(t)+(ω2

pεφMcEt,Et)(t)]+(∇×Et,∇×Et)(t)

+(∇×E,∇×E)(t)+A(ω2

pεφMcE,E)(t)

+µ0A(||Ht||2

0 +||ωpmH||2 0)(t) ≤ CF(0),

(14) where F(0) depends on initial conditions ∇×E(0),∇×Et(0),E(0), Et(0),Et2(0),H(0),∇×H(0),Ht(0) and D(0).

Jichun Li (UNLV) June 27, 2018 15 / 62

SLIDE 16

Numerical results: Li, Huang, Yang, Math Comp (2015)

Use R1 = 0.1m,R2 = 0.2m, γ = γm = 0 in our Drude model. A plane wave source: specified by Hz = 0.1sin(ωt), where ω = 2πf with operating frequency f = 2.0 GHz. A point source wave (same Hz) at a point (0:078; 0:4). Simulation with 65536 triangles, 28672 rectangles, and time step τ = 0.2ps.

Jichun Li (UNLV) June 27, 2018 16 / 62

SLIDE 17

Numerical results

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 20 / 62

SLIDE 21

Arbitray star-shaped cloak: SINUM 2018

x’=φ(x)

(b)

(x2,y2) (x1,y1) (x2,y2) (x1,y1) k(x2,y2) k(x1,y1)

O

. . . . . . . .

x x’ y’ y

O

For (x,y) = (0,0), arbitrary N-sided polygonal cloaking transformation:

x

′

= kx y1x2 −x1y2 y(x2 −x1)−x(y2 −y1) +(1−k)x, (15) y

′

= ky y1x2 −x1y2 y(x2 −x1)−x(y2 −y1) +(1−k)y. (16)

Jichun Li (UNLV) June 27, 2018 21 / 62

SLIDE 22

Cloaks of arbitray shapes: cont’d

Assume the limit of transformation (16) exists when (x2,y2) → (x1,y1), then the transformation and Jacobian matrix can be written as:

x

′ = kx

K1 y −K2x +(1−k)x, y

′ = ky

K1 y −K2x +(1−k)y, (17) ∂x

′

∂x = (1−k)+ kK1y (y −K2x)2 = (1−k)+ k(1−k)(K1y

′)

(y

′ −K2x ′)2(1−kK3)

= (1−k)+K1K4y

′,

∂x

′

∂y = − kK1x (y −K2x)2 = − k(1−k)K1x

′

(y

′ −K2x ′)2(1−kK3) = −K1K4x ′,

∂y

′

∂x = kK1K2y (y −K2x)2 = k(1−k)K1K2y

′

(y

′ −K2x ′)2(1−kK3) = K1K2K4y ′,

∂y

′

∂y = (1−k)− kK1K2x (y −K2x)2 = (1−k)− k(1−k)K1K2x

′

(y

′ −K2x ′)2(1−kK3)

= (1−k)−K1K2K4x

′,

(18)

where K1 = y1 −K2x1, K2 = P

′(x1), K3 =

K1 y′−K2x′ , K4 = k(1−k) (y′−K2x′)2(1−kK3).

Jichun Li (UNLV) June 27, 2018 22 / 62

SLIDE 23

Cloaks of arbitray shapes: cont’d

The relative permittivity and permeability of the cloak:

ε

′(x ′) =

JJT det(J) =

a0(x

′)

b0(x

′)

b0(x

′)

c0(x

′)

/det(J) =
a(x

′)

b(x

′)

b(x

′)

c(x

′)

,

µ′(x

′) =

1 det(J),

(19)

where

a0(x

′) = (1−k)2 +2(1−k)K1K4y ′ +K 2

1 K 2 4

(x

′)2 +(y ′)2

, b0(x

′) = K 2

1 K 2 4 K2

(x

′)2 −(y ′)2

−(1−k)K1K4(x

′ +K2y ′),

c0(x

′) = (1−k)2 −2(1−k)K1K2K4x ′ +K 2

1 K 2 2 K 2 4

(x

′)2 +(y ′)2

, det(J) = (1−k)2 +(1−k)K1K4(y

′ −K2x ′). Jichun Li (UNLV) June 27, 2018 23 / 62

SLIDE 24

Cloaks of arbitray shapes: cont’d

Theorem 2

Let λ1,λ2 be the eigenvalues of the relative permittivity tensor ε

′, we

have 0 < λ1 ≤ 1, λ2 ≥ 1. Moreover, in order to have the cloaking phenomenon, λ1 = λ2, i.e., the relative permittivity tensor must be anisotropic and one eigenvalue is less than 1. For the real symmetric matrix ε

′, we have decomposition ε ′ = PΛPT,

where Λ = λ1 λ2

,P =
p1

p2 −p2 p1

, and

p1 =

λ2−a(x

′)

λ2−λ1 ,

p2 = sign

b(x

′)

a(x

′)−λ1

λ2−λ1 .

Jichun Li (UNLV) June 27, 2018 24 / 62

SLIDE 25

Cloaks of arbitray shapes: cont’d

Using the Drude dispersion models: λ1(ω) = 1− ω2

e

ω2 −jωγe , (20) and µ′(ω) = µ∞ − ω2

m

ω2 −jωγm , (21) we obtain the governing equations for simulating arbitrary shaped cloaks:

Dt = ∇×H, (22) ε0λ2(Ett +γeEt +ω2

eE) = MADtt +γeMADt +ω2 eMBD,

(23) Bt = −∇×E, (24) µ0(µ∞Htt +γmHt +ω2

mH) = Btt +γmBt,

(25)

subject to PEC BC.

Jichun Li (UNLV) June 27, 2018 25 / 62

SLIDE 26

Cloaks of arbitray shapes: cont’d

Theorem 3

For any t ∈ [0,T], the model problem (22)-(40) exists a unique solution (E(·,t),H(·,t)) ∈ H0(curl;Ω)×H(curl;Ω).

Theorem 4

Denote λ ∗

2 =

sup

0=v∈H(curl;Ω) (Mv,v) (v,v) and (ω∗ e)2 =

inf

0=v∈H(curl;Ω) (ω2

ev,v)

(v,v) . Under

the assumptions that 8γ2

m ≤ (ω∗ e)2/λ ∗ 2, and γm and γe are constants, we

have: For any t ∈ (0,T],

λ ∗

2 µ0(µ∞||Ht||2 0 +||ωmH||2 0)(t)

+ 1 ε0 (||Dt||2

0 +||ωe

MBD||2

0)(t)+ε0(||MEt||2 0 +||ωeME||2 0)(t)

≤ C[λ ∗

2ε0(||

√ MEt||2

0 +||ωe

√ ME||2

0)(0)+λ ∗ 2 µ0(µ∞||Ht||2 0 +||ωmH||2 0)(0)]

+[ 1 2ε0 (||Dt||2

0 +||ωe

MBD||2

0)(0)+ ε0

2 (||MEt||2

0 +||ωeME||2 0)(0)].

(26)

Jichun Li (UNLV) June 27, 2018 26 / 62

SLIDE 27

Mixed FE spaces

Rectangular edge element: U U Uh = {ψh ∈ L2(Ω) : ψh|K ∈ Q0,0, ∀ K ∈ T h}, (27) V V V h = {φh ∈ H(curl;Ω) : φh|K ∈ Q0,1 ×Q1,0, ∀ K ∈ T h}, (28) Triangular edge element: U U Uh = {ψh ∈ L2(Ω) : ψh|K is a piecewise constant, ∀ K ∈ T h}, V V V h = {φh ∈ H(curl;Ω) : φh|K = span{λi∇λj −λj∇λi}, i,j = 1,2,3,∀ K ∈ V V V 0

h = {φh ∈ V

V V h, n ×φh = 0 on ∂Ω}. δτu u un = u u un −u u un−1 τ , δ 2

τ u

u un = u u un −2u u un−1 +u u un−2 τ2 , δ2τu u un = u u un −u u un−2 2τ , u u un−1 = u u un +2u u un−1 +u u un−2 4 , u u un = u u un +u u un−1 2 .

Jichun Li (UNLV) June 27, 2018 27 / 62

SLIDE 28

The FETD scheme

Given initial approximations E

− 3

2

h ,E − 1

2

h ,D − 1

2

h ,H0 h , for any n ≥ 0, find

D

n+ 1

2

h

∈ V 0

h,E n+ 1

2

h

∈ V 0

h,Hn+1 h

∈ Uh such that

(δτD

n+ 1

2

h

,vh) = (Hn

h,∇×vh),

(31) ε0(Mδ 2

τ E n+ 1

2

h

,φh)+ε0γe(Mδ2τE

n+ 1

2

h

,φh)+ε0(ω2

eM ¯

E

n+ 1

2

h

,φh) = (δ 2

τ D n+ 1

2

h

,φh)+γe(δ2τD

n+ 1

2

h

,φh)+(MBω2

e ¯

D

n+ 1

2

h

,φh), (32) µ0µ∞(δ 2

τ Hn+1,ψh)+ µ0γm(δ2τHn+1,ψh)+ µ0(ω2 m ¯

Hn+1

h

,ψh) = −(∇×δτE

n+ 1

2

h

,ψh)−γm

∇×

E

n+ 1

2

h

,ψh

,

(33)

hold true for any vh ∈ V 0

h,φh ∈ V 0 h,ψh ∈ Uh.

Jichun Li (UNLV) June 27, 2018 28 / 62

SLIDE 29

The FETD scheme: discrete stability

Theorem 5

Let (D

n+ 1

2

h

,E

n+ 1

2

h

,Hn+1

h

) be the FE solutions of the FETD scheme, under the assumption 32λ ∗

2γ2 m ≤ (ω∗ e)2 and time step constraint τ < min{ h 12CinvCv , 1 48γe(λ ∗

2)

1 2

, 1 48γm(λ ∗

2)

1 2

, 1 24||ωe||L∞ , 1 12||ωe||L∞(λ ∗

2)

1 2

, h(ω∗

e)2

2γmCinvCv , hµ∞ CinvCvγm }, (34)

then for any m ≥ 1, we have

ε0|| √ MδτE

m+ 1

2

h

||2

0 +ε0||ωe

√ M E

m+ 1

2

h

||2

0 + µ0||δτHm+1 h

||2

0 + µ0||ωm

Hm+1

h

||2 + 1 ε0 ||δτD

m+ 1

2

h

||2

0 + 1

ε0 ||ωe

MB

D

m+ 1

2

h

||2

0 +ε0||δτPhE m+ 1

2

h

||2 ≤ C(ε0|| √ MδτE

1 2

h ||2 0 +ε0||ωe

√ M E

1 2

h ||2 0 +ε0||

E

3 2

h ||2 0 + µ0||δτH1 h||2 0 + µ0||ωm

H1

h||2

+ 1 ε0 ||δτD

1 2

h ||2 0 + 1

ε0 ||ωe

MB

D

1 2

h ||2 0 +ε0||δτPhE

1 2

h ||2 0).

(35)

Jichun Li (UNLV) June 27, 2018 29 / 62

SLIDE 30

A mushroom shaped cloak

The physical domain is Ω = [−2,2]m×[−2,2]m, and the cloaked object is put inside a PEC mushroom enclosed by a contour curve ρ = 0.5+0.1sin(θ)+0.15sin(3θ)+0.15cos(4θ), then the cloaked region is wrapped by a mushroom cloak enclosed by the contour ρ = 1+0.2sin(θ)+0.3sin(3θ)+0.3cos(4θ). In our simulation, the mesh used is the one given in Fig. 8, which has 69930 total number of edges, 36076 total number of triangles, and 7632 total number of rectangles. We choose operating frequency f = 0.4GHz, the time step size τ = 1.5×10−13 s, the final simulation time T = 45ns (i.e, 30000 time steps). A plane wave source is specified by function Hz = 0.1sin(ωt), where ω = 2πf with operating frequency f = 1.0GHz.

Jichun Li (UNLV) June 27, 2018 30 / 62

SLIDE 31

A mushroom shaped cloak: cont’d

Figure 8 : The mesh used for the simulation of mushroom shaped cloak.

Jichun Li (UNLV) June 27, 2018 31 / 62

SLIDE 32

A mushroom shaped cloak: cont’d

Figure 9 : Electric fields Ey at various time steps. Top: t=6 ns, 12 ns, 15 ns; Bottom: t=24 ns, 30 ns, 45 ns.

Jichun Li (UNLV) June 27, 2018 32 / 62

SLIDE 33

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 33 / 62

SLIDE 34

Carpet cloak model: SIAM J Appl Math (2014)

Following Chen, Pendry, et al [Nature Communications, 2 (2011)], a triangular carpet cloak can be achieved with spatially homogeneous anisotropic dielectric materials.

Figure 10 : The physical space of the carpet cloak.

Jichun Li (UNLV) June 27, 2018 34 / 62

SLIDE 35

Carpet cloak modeling equations

Using the transformation optics (x = x′,y = H2−H1

H2

y′ + d−x·sgn(x)

d

H1), the relative permittivity and permeability of the cloak are given by: ε = a b b c

=
H2

H2−H1

−

H1H2 (H2−H1)d sgn(x)

−

H1H2 (H2−H1)d sgn(x) H2−H1 H2

+

H2 H2−H1 ( H1 d )2

,

µ = H2 H2 −H1 , where sgn(x) denotes the standard sign function.

Jichun Li (UNLV) June 27, 2018 35 / 62

SLIDE 36

Carpet cloak modeling equations: cont’d

Diagonalizing the symmetric matrix ε as: ε = PΣPT,, where λ1 = a+c−√

(a−c)2+4b2 2

, λ2 = a+c+√

(a−c)2+4b2 2

, and matrices Σ and P are Σ = λ1 λ2

,

P = p1 p2 p3 p4

,

p2

1 +p2 2 = 1, p1p3 +p2p4 = 0,

and elements pi,1 ≤ i ≤ 4, are given as p1 =

λ2 −a

λ2 −λ1 , p2 =

a−λ1

λ2 −λ1 ·sgn(x), p3 = −

λ2 −c

λ2 −λ1 ·sgn(x), p4 =

c −λ1

λ2 −λ1 .

Jichun Li (UNLV) June 27, 2018 36 / 62

SLIDE 37

Carpet cloak modeling equations: cont’d

It is easy to see that λ2 ≥ a+c+|a−c|

2

≥ a > 1, which leads to λ1 = 1/λ2 < 1. Mapping λ1 by the lossless Drude dispersion model: λ1(ω) = 1− ω2

p

ω2 , where ωp is plasma frequency, and ω is wave frequency. Substituting ε = PΣPT into D = ε0εE, we obtain ε0E = PΣ−1PTD, which can be written in time-domain (assuming eiωt time dependence) as follows: ε0λ2

Et2 +ω2

pE

= MADt2 +MBD,

(36) where matrices MA and MB are MA =

p2

1λ2 +p2 2

p2p4 +p1p3λ2 p2p4 +p1p3λ2 p2

3λ2 +p2 4

,

MB =

p2

2

p2p4 p2p4 p2

4

ω2

p.

Jichun Li (UNLV) June 27, 2018 37 / 62

SLIDE 38

Carpet cloak modeling equations: cont’d

The governing equations for the carpet cloak: Dt = ∇×H, (37) ε0λ2

Et2 +ω2

pE

= MADt2 +MBD,

(38) µ0µHt = −∇×E, (39) supplemented with initial conditions D(x,0) = D0(x), E(x,0) = E0(x), H(x,0) = H0(x), ∀ x ∈ Ω, (40) and the perfect conducting boundary condition (PEC): n ×E = 0

n ∂Ω.

(41) Here H denotes the magnetic field, and 2D vector and scalar curl

perators:

∇×H = (∂H ∂y ,−∂H ∂x )′, ∇×E = ∂Ey ∂x − ∂Ex ∂y , ∀E = (Ex,Ey)′.

Jichun Li (UNLV) June 27, 2018 38 / 62

SLIDE 39

Carpet cloak equations: existence

Lemma 6

The matrix MB is symmetric and non-negative definite, and the matrix MA is symmetric positive definite. Using Laplace transforms, we can prove

Theorem 7

For any t ∈ [0,T], there exists a unique solution (E(·,t),D(·,t),H(·,t)) ∈ (H0(curl;Ω))2 ×H(curl;Ω) of (37)-(41).

Jichun Li (UNLV) June 27, 2018 39 / 62

SLIDE 40

refined stability: Li, Meng, Huang (CMAM2018)

Theorem 8

For any t ∈ [0,T], the following stability holds true:

ε0µ0µλ2(||M

− 1

2

A ∂t2E||2 +ω2 p||M − 1

2

A ∂tE||2)+||∇×∂tE||2

+||M

1 2

A ∂tD||2 +||M

1 2

B D||2 +ε0µ0µλ2ω2 p||H||2

(t)

≤ C

ε0µ0µλ2(||M

− 1

2

A ∂t2E||2 +ω2 p||M − 1

2

A ∂tE||2)+||∇×∂tE||2

+||M

1 2

A ∂tD||2 +||M

1 2

B D||2 +ε0µ0µλ2ω2 p||H||2

(0),

where the constant C > 0 depends on the physical parameters ε0,µ0,d,H1,H2 and ωp.

Jichun Li (UNLV) June 27, 2018 40 / 62

SLIDE 41

The FETD scheme: CMAM2018

Denote difference and average operators: δτun+ 1

2 = un+ 1 2 −un− 1 2

τ , δ 2

τ un+ 1

2 = un+ 1 2 −2un− 1 2 +un− 3 2

τ2 . Construct a leap-frog scheme for the model equations (37)-(39): Given proper initial approximations H0

h,D − 1

2

h ,D − 3

2

h ,E − 1

2

h ,E − 3

2

h , for n ≥ 0 find

D

n+ 1

2

h

,E

n+ 1

2

h

∈ V 0

h, Hn+1 h

∈ Uh such that

δτD

n+ 1

2

h

,φh

= (Hn

h ,∇×φh),

(42) ε0λ2

δ 2

τ E n+ 1

2

h

,ϕh

+ε0λ2ω2

p

E

n+ 1

2

h

,ϕh

=
MAδ 2

τ D n+ 1

2

h

,ϕh

+
MBD

n+ 1

2

h

,ϕh

,

(43) µ0µ

δτHn+1

h

,ψh

= −(∇×E

n+ 1

2

h

,ψh), (44)

hold true for any φ h, ϕh ∈ V 0

h, ψh ∈ Uh.

Jichun Li (UNLV) June 27, 2018 41 / 62

SLIDE 42

The FETD scheme: discrete stability

Theorem 9

For the FETD solution (D

n+ 1

2

h

,E

n+ 1

2

h

), denote the discrete energy at time level n:

ENGn = ε0µ0µλ2(||M

− 1

2

A δ 2 τ E n+ 1

2

h

||2 +ω2

p||M − 1

2

A δτE n+ 1

2

h

||2)+||∇×δτE

n+ 1

2

h

||2 +||M

1 2

A δτD n+ 1

2

h

||2 +||M

1 2

B D n+ 1

2

h

||2 +ε0µ0µλ2ω2

p||Hn h||2.

(45)

Then for any n ≥ 1 and under the constraint τ ≤ min

1,

1 4ε0λ2 , 1 4µ0µ||M−1

A MB||2 , µ0µ

4 , ε0λ2 ||M−1

A MB||2 ,

h2ε0µ0µλ2 8C2

inv||M1/2 A

||2

,

(46) we have ENGn ≤ C ·ENG0, where C > 0 is independent of τ and h.

Jichun Li (UNLV) June 27, 2018 42 / 62

SLIDE 43

Optimal error estimate: a lemma

Lemma 10

(i) ||δτ uk+ 1

2 ||2 := || uk+ 1 2 −uk− 1 2

τ ||2 ≤ 1 τ

tk+ 1 2 tk− 1 2

||∂t u||2ds, ∀ u ∈ H1(0,T;L2(Ω)), (ii) ||δ2

τ uk+ 1 2 ||2 := || uk+ 1 2 −2uk− 1 2 +uk− 3 2

τ2 ||2 ≤ 8 τ

tk+ 1 2 tk− 3 2

||∂t2 u||2ds, ∀ u ∈ H2(0,T;L2(Ω)), (iii) ||δ2

τ uk+ 1 2 −∂t2 u(tk+ 1 2

)||2 ≤ 4τ

tk+ 1 2 tk− 3 2

||∂t3 u||2ds, ∀ u ∈ H3(0,T;L2(Ω)), (iv) ||δ3

τ uk+ 1 2 ||2 ≤ 81

τ

tk+ 1 2 tk− 5 2

||∂t3 u||2ds, ∀ u ∈ H3(0,T;L2(Ω)), (v) ||δ3

τ uk+ 1 2 −δτ ∂t2 u(tk+ 1 2

)||2 ≤ 81τ

tk+ 1 2 tk− 5 2

||∂t4 u||2ds, ∀ u ∈ H4(0,T;L2(Ω)), (vi) ||δτ uk+ 1

2 −∂t u(tk )||2 ≤ τ

12

tk+ 1 2 tk− 1 2

||∂t2 u||2ds, ∀ u ∈ H2(0,T;L2(Ω)), (vii) ||δτ ∂t uk− 1

2 −δ2 τ u(tk )||2 ≤ 2τ

5

tk tk−2

||∂t3 u||2ds, ∀ u ∈ H3(0,T;L2(Ω)), (viii) ||δ2

τ ∂t u(tk )−δ3 τ uk+ 1 2 ||2 ≤ 256τ tk+ 1 2 tk− 5 2

||∂t4 u||2ds, ∀ u ∈ H4(0,T;L2(Ω)). Jichun Li (UNLV) June 27, 2018 43 / 62

SLIDE 44

Optimal error estimate: error equations

The error equations are given as follows:

(δτ D

n+ 1

2

h

,φ h) = ( Hn

h,∇×φ h)+(∂tDn −δτDn+ 1

2 ,φ h),

(47) ε0λ2(δ 2

τ

E

n+ 1

2

h

,ϕh)+ε0λ2ω2

p(

E

n+ 1

2

h

,ϕh)−(MAδ 2

τ

D

n+ 1

2

h

,ϕh)−(MB D

n+ 1

2

h

,ϕh) = ε0λ2(∂t2En+ 1

2 −δ 2

τ En+ 1

2 ,ϕh)−(MA(∂t2Dn+ 1 2 −δ 2

τ Dn+ 1

2 ),ϕh),

(48) µ0µ(δτ Hn+1

h

,ψh) = −(∇× E

n+ 1

2

h

,ψh)+ µ0µ(∂tHn+ 1

2 −δτHn+1,ψh),

(49)

hold true for any φ h,ϕh ∈ V 0

h and ψh ∈ Uh. For simplicity, we assume

that the initial conditions used for the scheme (42)–(44) are as follows:

H0

h = Π2H0, D − 1

2

h

= ΠcD− 1

2 , D

− 3

2

h

= ΠcD− 3

2 , E

− 1

2

h

= ΠcE− 1

2 , E

− 3

2

h

= ΠcE− 3

2 . Jichun Li (UNLV) June 27, 2018 44 / 62

SLIDE 45

Optimal error estimate

Theorem 11

Denote the total errors

Errn =

||M

1 2 A δτ

D

n+ 1 2 h

||2 +||M

1 2 B

D

n+ 1 2 h

||2 +||∇×δτ E

n+ 1 2 h

||2 +||∇× E

n+ 1 2 h

||2 +ε0µ0µλ2(||M

− 1 2 A

δ2

τ

E

n+ 1 2 h

||2 +||ωpM

− 1 2 A

δτ E

n+ 1 2 h

||2 +||ωpM

− 1 2 A

E

n+ 1 2 h

||2) 1/2 . (50)

Then for any m ≥ 1, sufficiently small time step satisfying the constraint (46), and under the following regularity assumptions:

max

t∈[0,T]

||E||Hp(curl;Ω), ||D||Hp(curl;Ω), ||∇×E||Hp(curl;Ω)
< ∞,

T

||∂t2 E||2

Hp(curl;Ω) +||∂t3 E||2 Hp(curl;Ω) +||∂t4 E||2 +||∇×∇×∂t E||2 Hp(curl;Ω) +||∂t3 D||2

+||∂t4 D||2 +||∂t D||2

Hp(curl;Ω) +||∂t2 D||2 Hp(curl;Ω) +||∂t2 H||2 +||∇×∂t3 H||2

ds < ∞,

we have Errm ≤ CT(τ2 +hp), where the constant C > 0 is independent

f h and τ, and p is the degree of the basis functions in Uh and V h.

Jichun Li (UNLV) June 27, 2018 45 / 62

SLIDE 46

Discontinuous Galerkin method: CMAM2018

On triangular meshes, our finite element spaces are given by:

Uh = {ψh ∈ L2(Ω) : ψh|K ∈ Pp, ∀ K ∈ Th}, V V V h = (Uh)2.

Given proper initial approximations H0

h,D − 1

2

h ,D − 3

2

h , E − 1

2

h ,E − 3

2

h , for n ≥ 0

find D

n+ 1

2

h

,E

n+ 1

2

h

∈ V 0

h, Hn+1 h

∈ Uh such that on any element Ki ∈ Th,

Ki

δτ D

n+ 1 2 h

·φh =

Ki

Hn

h ·∇×φh + ∑ K∈νi

aik

φh ×nik ·{{Hn

h }}ik ,

(51) ε0λ2

Ki

δ2

τ E n+ 1 2 h

·ϕh +ε0λ2

Ki

ω2

p E n+ 1 2 h

·ϕh =

Ki

MAδ2

τ D n+ 1 2 h

·ϕh +

Ki

MBD

n+ 1 2 h

·ϕh, (52) µ0µ

Ki

δτ Hn+1

h

ψh = −

Ki

E

n+ 1 2 h

·∇×ψh − ∑

K∈νi

aik

ψh ·nik ×{{E

n+ 1 2 h

}}ik , (53)

hold true for any φ h,ϕh ∈ V0

h,ψh ∈ Uh. {{E n+ 1

2

h

}}ik denotes average of E

n+ 1

2

h

n internal face aik, νi denotes all neighboring elements of K.

Jichun Li (UNLV) June 27, 2018 46 / 62

SLIDE 47

PML

To simulate the cloak phenomenon, we surround the physical domain by a perfectly matched layer (PML), see Fig.10 (Right). In this paper, we use the classical 2D Berenger PML, whose governing equations can be written as ε0 ∂E E E ∂t + σy σx

E

E E = ∇×Hz, (54) µ0 ∂Hzx ∂t +σmxHzx = −∂Ey ∂x , (55) µ0 ∂Hzy ∂t +σmyHzy = ∂Ex ∂y , (56) where Hz = Hzx +Hzy denotes the magnetic field, the parameters σi and σm,i,i = x,y, are the electric and magnetic conductivities in the x- and y- directions, respectively. In our simulation, we use a PML with 12 rectangular cells in thickness around the physical domain.

Jichun Li (UNLV) June 27, 2018 47 / 62

SLIDE 48

Numerical examples: hybrid edge elements

In our simulation, we choose H1 = 0.05m,H2 = 0.2m,d = 0.2m, and the physical domain Ω = [−0.3,0.3] m ×[0,0.3] m, which is partitioned by a uniform triangular mesh with a mesh size h = 0.00625. The PML region surrounding Ω is partitioned by a uniform rectangular mesh. Our final mesh yields 53330 total edges, 26960 total triangular elements, and 6258 total rectangular elements. In the test, we choose the time step size τ = 2∗10−13 s, and the total number of time steps 15000, i.e., the final simulation time T = 3.0 nanosecond (ns). Example 1. The incident wave is generated by a plane wave source Hz = 0.1sin(ωt) imposed at line x = −0.3, where ω = 2πf with frequency f = 3.0 GHz. The numerical magnetic fields Hz at different time steps are shown in Fig.11. Both figures show clearly that the plane wave pattern is recovered very well after passing through the cloaking region, which makes any objects hiden inside the cloaked region invisible to observers at the far end.

Jichun Li (UNLV) June 27, 2018 48 / 62

SLIDE 49

Numerical examples: Ex 1

Figure 11 : Ex1. The Hz fields at 5000, 7000, 10000, 15000 time steps.

Jichun Li (UNLV) June 27, 2018 49 / 62

SLIDE 50

Numerical examples: hybrid edge elements

Example 2. The incident wave is generated by a Gaussian wave Hz(x,y,t) = 0.1e−(y−0.15)2/(60L)2 sin(ωt) imposed along a slanted line y = x +0.45, where L = 0.004 √ 2, and ω = 2πf with frequency f = 6.0 GHz. The numerical magnetic fields Hz at different time steps are presented in Fig.12. To appreciate the cloak phenomenon, in Fig.13 we present the magnetic fields Hz obtained without the cloaking material. It is clear that the cloak phenomenon disappears if the cloaking material is removed.

Jichun Li (UNLV) June 27, 2018 50 / 62

SLIDE 51

Numerical examples: Ex 2

Figure 12 : Ex2. The Hz fields at 5000, 7000, 10000, 15000 time steps.

Jichun Li (UNLV) June 27, 2018 51 / 62

SLIDE 52

Numerical examples: Ex 2

Figure 13 : Ex2. The Hz fields at 5000, 7000, 10000, 15000 time steps

btained with the cloaking material removed.

Jichun Li (UNLV) June 27, 2018 52 / 62

SLIDE 53

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 53 / 62

SLIDE 54

Electromagnetic concentrator: Li etc (CiCP 2018)

Figure 14 : Snapshots of electric field Ey: 2000 time steps (Top Left); 4000 steps (Top Right); 8000 steps (Bottom Left); 10000 steps (Bottom Right).

Jichun Li (UNLV) June 27, 2018 54 / 62

SLIDE 55

Electromagnetic concentrator: cont’d

coordinate transformation: r ′ =

a

br,

0 ≤ r ≤ b,

c−a c−br − b−a c−bc,

b ≤ r ≤ c, (57) φ ′ = φ, 0 ≤ φ ≤ 2π. (58) This mapping compresses the circle r ≤ b into a smaller circle r ′ ≤ a, where a < b, and in the same time expands the annulus b ≤ r ≤ c to a ≤ r ′ ≤ c. The major usage of this transformation is to concentrate the energy in the inner area r ′ ≤ a.

Jichun Li (UNLV) June 27, 2018 55 / 62

SLIDE 56

Electromagnetic concentrator: cont’d

∂tD = ∇×H, (59) ε0εmaxεrMa∂ttE +ε0ω2

eεrMaE = ∂ttD +MbD,

(60) µ0µmax∂ttH + µ0ω2

mH = −∇×∂tE.

(61) εr = r ′+K1

r ′

, where r ′ ∈ [a,c] and K1 = (b−a)c

c−b . The parameters

εmax = c−b

c−a, ωe, µmax = b(c−b) a(c−a) and ωm come from the following Drude

models used to map εφ =

r ′ r ′+K1 and µ′(r ′) =

c−b

c−a

2 r ′+K1

r ′

since they can be less than one: εφ = εmax − ω2

e

ω2 , µ′(r ′) = µmax − ω2

m

ω2 .

Jichun Li (UNLV) June 27, 2018 56 / 62

SLIDE 57

Electromagnetic concentrator: stability

Theorem 12

Denote the energy ENG(t) =

||∂tD||2 +||M1/2

b

D||2 + µ0µmax||∂tH||2 + µ0||ωmH||2 +ε0εmax||ε1/2

r

M1/2

a

∂tE||2 +ε0||ωeε1/2

r

M1/2

a

E||2 +||∇×∂tE||2 +ε0µ0µmax

εmax||ε1/2

r

M1/2

a

∂ttE||2 +||ωeε1/2

r

M1/2

a

∂tE||2 (t). Then for the solution of (59)-(61) satisfying the PEC boundary condition (41) and any t ∈ [0,T], there exists a constant C > 0 such that ENG(t) ≤ C ·ENG(0).

Jichun Li (UNLV) June 27, 2018 57 / 62

SLIDE 58

Electromagnetic concentrator: FETD scheme

2

h ,H0 h, for any

n ∈ [0,N −1], find D

n+ 1

2

h

∈ V 0

h,E n+ 1

2

h

∈ V 0

h, Hn+1 h

∈ Uh such that (δτDn

h,vh) = (Hn h,∇×vh), ∀ vh ∈ V 0 h,

ε0εmax(εrMaδ 2

τ E n− 1

2

Theorem 13

Denote the discrete energy

Edisc(n) = ||δτDn

h||2 +||M1/2 b

¯ D

n h||2 + µ0µmax||δτHn+ 1

2 ||2 + µ0||ωm ¯

H

n+ 1

2

h

||2 +ε0εmax||ε

1 2

r M

1 2

a δτEn h||2 +ε0||ωeε

1 2

r M

1 2

a ¯

E

n h||2 +||∇×δτ ¯

E

n+ 1

2

h

||2 +ε0µ0µmax(εmax||ε

1 2

r M

1 2

a δ 2 τ E n+ 1

2

h

||2 +||ωeε

1 2

r M

1 2

a δτ ¯

E

n+ 1

2

h

||2).

Then for the finite element solution and any n ∈ [1,N −1], and a small enough time step size τ satisfing the CFL condition: τ ≤ Csmallh, (62) there exists a constant C > 0, independent of τ and h, such that Edisc(n) ≤ C ·Edisc(0).

Jichun Li (UNLV) June 27, 2018 59 / 62

SLIDE 60

Electromagnetic concentrator: cont’d

Theorem 14

Denote the fully-discrete error energy

ERRdisc(n) = ||δτ D

n h||2 +||M1/2 b

¯

D

n h||2 + µ0µmax||δτ

Hn+ 1

2 ||2 + µ0||ωm ¯

H

n+ 1

2

h

||2 +ε0εmax||ε

1 2

r M

1 2

a δτ

E

n h||2 +ε0||ωeε

1 2

r M

1 2

a ¯

E

n h||2 +||∇×δτ ¯

E

n+ 1

2

h

||2 +ε0µ0µmax(εmax||ε

1 2

r M

1 2

a δ 2 τ

E

n+ 1

2

h

||2 +||ωeε

1 2

r M

1 2

a δτ ¯

E

n+ 1

2

h

||2).

Then for any n ∈ [1,N −1], under CFL condition (63), there exists a constant C > 0, independent of h and τ, such that ERRdisc(n) ≤ C(τ4 +h2p).

Jichun Li (UNLV) June 27, 2018 60 / 62

SLIDE 61

1

Introduction to electromagenetic cloaking with metamaterials

2

Cylindrical cloak in time-domain

3

Arbitray star-shaped cloak

4

Carpet cloak model

5

Other Applications

6

Summary

Jichun Li (UNLV) June 27, 2018 61 / 62

SLIDE 62

Summary

Li and Huang, Time-Domain Finite Element Methods for Maxwell’s

Equations in Metamaterials, Springer Series in Computational Mathematics, vol.43, Springer, 2013.

Li, Huang, Yang and Wood, SIAM J. Appl. Math. 74(4) (2014)

1136-1151

Li, Huang and Yang, Math. Comp. 84 (2015), no. 292, 543-562
Yang, Li and Huang, SIAM J. Numer. Anal. 56(1) (2018) 136-159
Li, Meng and Huang, Comput. Methods Appl. Math. 2018
well-posedness and regularity; optimal error estimate for DG
a posteriori error estimator; hp-adaptivity (adaptive DG);
potential applications: solar cells, rotator, splitter, Cherenov radiation,
ptical blackhole, ...

Thank you for your attention!

Jichun Li (UNLV) June 27, 2018 62 / 62