Finite-Difference Time-Domain Method for Complex Media Jinjie Liu - - PowerPoint PPT Presentation

Finite-Difference Time-Domain Method for Complex Media

Jinjie Liu

Delaware state university Collaborators:

– Arizona: Profs. J. Moloney, M. Brio, C. Dineen, and Ph.D. students J. Nehls and A. Ha – Delaware State: Dr. J. Zhang, Dr. P. Xu, J. Cornelius, A. Strong

ICERM Workshop “Computational Aspects of Time Dependent Electromagnetic Wave Problems in Complex Materials” June 26, 2018

1

Outline

• Maxwell’s equations in complex media

– Anisotropic, Dispersive, and Nonlinear

• Anisotropic FDTD

– Grid Transformation (Transformation Optics)

• Dual grid FDTD method

– Magneto-electric materials (spacetime cloak) – Pulse propagation in Dispersive and Nonlinear media

• Nonlinear Drude model for Second Harmonic Generation

(SHG) from metallic structures

2

The Maxwell’s Equations of electrodynamics

Gauss' laws: 𝛼 ⋅ 𝐶 = 0 𝛼 ⋅ 𝐸 = 0 Faraday’s Law: 𝜖𝐶 𝜖𝑢 = −𝛼 × 𝐹 Ampère's Law: 𝜖𝐸 𝜖𝑢 = 𝛼 × 𝐼 Constitutive Relations: 𝐸 = 𝜗 𝐹 + 𝑸 𝐶 = 𝜈 𝐼

3

polarization 𝑄 = 𝑄𝑀 + 𝑄𝑂𝑀

𝜈 𝜖𝐼 𝜖𝑢 = −𝛼 × 𝐹 𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 × 𝐼 − 𝑲

current density: 𝐾 = 𝑒𝑸𝑴

𝑒𝑢

Complex Media

𝐸 = 𝝑 𝐹 + 𝑄𝑀 + 𝑄𝑂𝑀

• Dielectric: 𝝑 is constant, 𝐸 = 𝝑 𝐹
• Anisotropic: 𝝑 is a tensor, 𝐸 = ത

𝝑 𝐹 ത 𝝑 = 𝜗11 𝜗12 𝜗13 𝜗21 𝜗22 𝜗23 𝜗31 𝜗32 𝜗33

• Dispersive: 𝝑(𝝏) complex and frequency dependent, 𝑄𝑀 = 𝜗 𝜕 𝐹

– Example: Drude 𝜗 𝜕 = 1 − 𝜕𝑞

2

𝜕2 + 𝑗𝛿𝜕

• Nonlinear: 𝐸 = 𝜗 𝐹 + 𝑸𝑶𝑴

– Example: Kerr 𝑄𝑂𝑀 = 𝜓 𝐹 2𝐹

• Magneto-electric material: 𝐸 = 𝜗𝐹 + 𝛾𝐼

– Examples: Space-time cloak, photonic topological insulator

4

Finite-Difference Time-Domain (FDTD)

• Finite-Difference Time-Domain (FDTD) method

𝜈

𝜖𝐼 𝜖𝑢 = −𝛼 × 𝐹

𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 × 𝐼

• Yee ’66, Taflove ’75
• Staggered Cartesian grid in space & time
• Centered Finite Difference and Leapfrog
• Non-dissipative
• Divergence free
• Robust and easy to implement
• Variety of materials: Dielectric, Anisotropic, Dispersive, Nonlinear, etc

5

Ey Ez Ex Hy Hx Hz

Electrically and Magnetically Anisotropic media

• Use Finite-Difference Time-Domain (FDTD) method to solve Ampère and

Faraday Laws to update 𝐸 & 𝐶 𝜖𝐸 𝜖𝑢 = 𝛼 × 𝐼, 𝜖𝐶 𝜖𝑢 = −𝛼 × 𝐹

• Use Constitutive Relations to update 𝐹 & 𝐼:

𝐹 = 𝜗−1𝐸, 𝐼 = 𝜈−1𝐶

𝐹𝑦 𝐹𝑧 𝐹𝑨 = 𝜊𝑦𝑦 𝜊𝑦𝑧 𝜊𝑦𝑨 𝜊𝑧𝑦 𝜊𝑧𝑧 𝜊𝑧𝑨 𝜊𝑨𝑦 𝜊𝑨𝑧 𝜊𝑨𝑨 𝐸𝑦 𝐸𝑧 𝐸𝑨 𝐼𝑦 𝐼𝑧 𝐼𝑨 = 𝜃𝑦𝑦 𝜃𝑦𝑧 𝜃𝑦𝑨 𝜃𝑧𝑦 𝜃𝑧𝑧 𝜃𝑧𝑨 𝜃𝑨𝑦 𝜃𝑨𝑧 𝜃𝑨𝑨 𝐶𝑦 𝐶𝑧 𝐶𝑨 𝜗−1 = (𝜊), and 𝜈−1 = 𝜃 .

6

Anisotropic material

Constitutive Equations: 𝑭 = 𝝑−𝟐𝑬, 𝐼 = 𝜈−1𝐶

𝐹𝑦 = 𝜊𝑦𝑦𝐸𝑦 + 𝜊𝑦𝑧𝐸𝑧 + 𝜊𝑦𝑨𝐸𝑨 𝐹𝑦,𝑗+1

2,𝑘,𝑙 = 𝜊𝑦𝑦𝐸𝑦,𝑗+1 2,𝑘,𝑙 + 𝜊𝑦𝑧𝑬𝒛,𝒋+𝟐 𝟑,𝒌,𝒍 + 𝜊𝑦𝑨𝑬𝒜,𝒋+𝟐 𝟑,𝒌,𝒍

7

𝑭/𝑬𝒛,𝒋+𝟐,𝒌+𝟐

𝟑,𝒍

𝑭/𝑬𝒜,𝒋,𝒌,𝒍+𝟐

𝟑

𝑭/𝑬𝒚,𝒋+𝟐

𝟑,𝒌,𝒍

Staggered Yee Lattice

Update equation for anisotropic material

Constitutive Equations: 𝐹 = 𝜗−1𝐸, 𝐼 = 𝜈−1𝐶

𝐹𝑦 = 𝜊𝑦𝑦𝐸𝑦 + 𝜊𝑦𝑧𝐸𝑧 + 𝜊𝑦𝑨𝐸𝑨 𝐹𝑦,𝑗+1

2,𝑘,𝑙 = 𝜊𝑦𝑦𝐸𝑦,𝑗+1 2,𝑘,𝑙 + 𝝄𝒚𝒛𝑬𝒛,𝒋+𝟐 𝟑,𝒌,𝒍 + 𝜊𝑦𝑨𝐸𝑨,𝑗+1 2,𝑘,𝑙

• Non-averaging (1st order)

𝐸𝑧,𝑗+1

2,𝑘,𝑙 = 𝐸𝑧,𝑗,𝑘+1 2,𝑙

• Averaging method 1 (unstable)

𝐸𝑧,𝑗+1

2,𝑘,𝑙 =

𝐸𝑧,𝑗,𝑘+1

2,𝑙 + 𝐸𝑧,𝑗,𝑘−1 2,𝑙 + 𝐸𝑧,𝑗+1,𝑘+1 2,𝑙 + 𝐸𝑧+1,𝑗,𝑘−1 2,𝑙

4

8

Update equation for anisotropic material (cont.)

Constitutive Equations: 𝐹 = 𝜗−1𝐸, 𝐼 = 𝜈−1𝐶

𝐹𝑦 = 𝜊𝑦𝑦𝐸𝑦 + 𝜊𝑦𝑧𝐸𝑧 + 𝜊𝑦𝑨𝐸𝑨 𝐹𝑦,𝑗+1

2,𝑘,𝑙 = 𝜊𝑦𝑦𝐸𝑦,𝑗+1 2,𝑘,𝑙 + 𝝄𝒚𝒛𝑬𝒛,𝒋+𝟐 𝟑,𝒌,𝒍 + 𝜊𝑦𝑨𝐸𝑨,𝑗+1 2,𝑘,𝑙

• Averaging method 2 (Werner & Cary ’07,’13)

𝜊𝑦𝑧,𝑗+1

2,𝑘,𝑙𝐸𝑧,𝑗+1 2,𝑘,𝑙 = 𝜊𝑦𝑧,𝑗,𝑘,𝑙𝐸𝑀 + 𝜊𝑦𝑧,𝑗+1,𝑘,𝑙𝐸𝑆

2 𝐸𝑀 = 𝐸𝑧,𝑗,𝑘+1

2,𝑙 + 𝐸𝑧,𝑗,𝑘−1 2,𝑙

2 , 𝐸𝑆 = 𝐸𝑧,𝑗+1,𝑘+1

2,𝑙 + 𝐸𝑧,𝑗+1,𝑘−1 2,𝑙

2

9

Invariant Coordinate Transformation (Transformation Optics)

𝜈 𝜖𝐼 𝜖𝑢 = −𝛼 × 𝐹 𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 × 𝐼

𝜈′ 𝜖𝐼′ 𝜖𝑢 = −𝛼′ × 𝐹′ 𝜗′ 𝜖𝐹′ 𝜖𝑢 = 𝛼′ × 𝐼′ Coordinate Transformation from (x,y,z) to (x’,y’,z’)

𝐹′ = 𝛭𝑈 𝐹, 𝜗′ = 𝛭 𝛭−1𝜗 𝛭−𝑈 𝐼′ = 𝛭𝑈𝐼, 𝜈′ = 𝛭 𝛭−1𝜈𝛭−𝑈 𝛭 is the Jacobian matrix

10

isotropic anisotropic

Teixeira ‘99, Pendry ’06, Leonhardt ‘06

Coordinate Stretching

11

𝑠′ = 𝑔 𝑠

Computational domain Physical domain

Mapping function

12

𝑠′ = 𝑔(𝑠) 𝑆1′ 𝑆2 𝑆1 𝑆2 𝑝

𝑔 maps 0 < 𝑠 ≤ 𝑆1 to 0 < 𝑠′ ≤ 𝑆′1 and 𝑆1 ≤ 𝑠 ≤ 𝑆2 to 𝑆1

′ ≤ 𝑠 ≤ 𝑆2

𝑔 continuous and 𝑔’ continuous, for example, by using cubic spline

𝑠 𝑠′

Transformed mesh and material

13

𝜗11 not constant Physical domain Computational domain

𝜗 = 1

mapping

Transformation based Maxwell Solver

14

Output Post-Processing Maxwell Solver Pre-Processing Input

Grid Mapping: from physical domain to a computational domain inverse mapping back to physical domain

EM Anisotropic material

Liu, et. al. J. Comput. Phys. 2014

TO-FDTD for dispersive material: the 𝐾 version

𝜈 𝜖𝐼 𝜖𝑢 = −𝛼 × 𝐹 𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 × 𝐼 − 𝐾 𝜖𝐾 𝜖𝑢 + 𝛿𝐾 = 𝜗0𝜕𝑞

2𝐹

𝝂′ 𝜖𝐼′ 𝜖𝑢 = −𝛼′ × 𝐹′ 𝝑′ 𝜖𝐹′ 𝜖𝑢 = 𝛼′ × 𝐼′ − 𝐾′ 𝜖𝐾′ 𝜖𝑢 + 𝛿𝐾′ = Λ −1Λ𝜗0 Λ𝑈𝜕𝑞

2𝐹′

𝐹′ = 𝛭𝑈𝐹, 𝝑′ = 𝛭 𝛭−1𝜗 𝛭−𝑈 𝐼′ = 𝛭𝑈𝐼, 𝝂′ = 𝛭 𝛭−1𝜈𝛭−𝑈 𝐾′ = Λ 𝛭−1𝐾, 𝛭 is the Jacobian matrix from (𝑦’, 𝑧’, 𝑨’) 𝑢𝑝 (𝑦, 𝑧, 𝑨)

15

(𝑦, 𝑧, 𝑨) (𝑦’, 𝑧’, 𝑨’) mapping

Coordinate Transformation from (x,y,z) to (x’,y’,z’)

TO-FDTD for dispersive material: the 𝑄 version

𝜖𝐶 𝜖𝑢 = −𝛼 × 𝐹 𝜖𝐸 𝜖𝑢 = 𝛼 × 𝐼 𝐶 = 𝜈𝐼 𝐸 = 𝜗0𝜗∞ 𝐹 + 𝑄 𝜖2𝑄 𝜖𝑢2 + 𝛿 𝜖𝑄 𝜖𝑢 = 𝜕𝑞

2

𝜗∞ 𝐹 𝜖𝐶′ 𝜖𝑢 = −𝛼′ × 𝐹′ 𝜖𝐸′ 𝜖𝑢 = 𝛼′ × 𝐼′ 𝐶′ = 𝝂′𝐼′ 𝐸′ = 𝝑′ 𝐹′ + 𝑄′ 𝜖2𝑄′ 𝜖𝑢2 + 𝛿 𝜖𝑄′ 𝜖𝑢 = 𝜕𝑞

2

𝜗∞ 𝐹′ 𝐹′ = 𝛭𝑈𝐹, 𝝑′ = 𝛭 𝛭−1𝜗0𝜗∞ 𝛭−𝑈 𝐼′ = 𝛭𝑈𝐼, 𝝂′ = 𝛭 𝛭−1𝜈 𝛭−𝑈 𝑄′ = 𝛭 𝛭−1𝑄, 𝛭 is the Jacobian matrix from (𝑦’, 𝑧’, 𝑨’) 𝑢𝑝 (𝑦, 𝑧, 𝑨)

16

(𝑦, 𝑧, 𝑨) (𝑦’, 𝑧’, 𝑨’) mapping

Coordinate Transformation from (x,y,z) to (x’,y’,z’)

Ring Simulation using grid transformation

FDTD mesh (Cartesian, isotropic) TO Physical Domain (Non-Rectangular mesh) TO Comput. domain (Cartesian, anisotropic)

17

Grid Rotation

18

ϵ > 1 ϵ’ 𝜄′ = 𝑔(𝜄) Computational domain Physical domain

Metallic Bowtie simulation

1600 nm = 2 𝜇

FDTD: Δx = 4 nm FDTD: Δx = 2 nm TO-FDTD: Δx = 20 nm

19 20 nm gap

Mesh 𝑶𝒛 CPU Time FDTD 400

(Δ = 4𝑜𝑛)

40 FDTD 800

(Δ = 2𝑜𝑛)

320 TO-FDTD 80 1

Subgridding/Adaptive Mesh Refinement

20

Beam Propagation

21

Spatial and Temporal Subgridding

22

𝚬𝒚, 𝚬𝒛, 𝚬𝒖 𝚬𝒚 𝟑 , 𝚬𝒛 𝟑 , 𝚬𝒖 𝟑

Temporal Subgridding + irregular mesh

23

nonrectangular mesh,

𝚬𝒖 𝟑

𝚬𝒚, 𝚬𝒛, 𝚬𝒖

Temporal subcycling with iterations

Update coarse mesh for Δ𝑢 Update fine mesh for Δ𝑢/2 Update fine mesh for Δ𝑢/2

24

several times (near interface)

𝐼/𝐹𝑔𝑗𝑜𝑓 𝐼/𝐹𝑑𝑝𝑏𝑠𝑡𝑓

Update coarse mesh by 𝚬𝐮

• Update 𝐹𝑑𝑝𝑏𝑠𝑡𝑓 from 𝑜 −

1 2 to 𝑜 + 1 2

• Update 𝐼𝑑𝑝𝑏𝑠𝑡𝑓 from 𝑜 to 𝑜 + 1

Update fine mesh by 𝚬𝐮/𝟑

• Update 𝐹

𝑔𝑗𝑜𝑓 from 𝑜 − 1 2 to 𝑜

• Interpolate 𝐹

𝑔𝑗𝑜𝑓 at mesh boundary using 𝐹𝑑𝑝𝑏𝑠𝑡𝑓

• Update 𝐼

𝑔𝑗𝑜𝑓 from 𝑜 − 1 4 to 𝑜 + 1 4

• Interpolate 𝑰𝒅𝒑𝒃𝒔𝒕𝒇 at 𝒐 using 𝑰𝒈𝒋𝒐𝒇 at 𝒐 −

𝟐 𝟓 and 𝒐 + 𝟐 𝟓

Update fine mesh by 𝚬𝐮/𝟑

• Update 𝐹

𝑔𝑗𝑜𝑓 from 𝑜 to 𝑜 + 1 2

• Interpolate 𝐹

𝑔𝑗𝑜𝑓 at mesh boundary using 𝐹𝑑𝑝𝑏𝑠𝑡𝑓

• Update 𝐼

𝑔𝑗𝑜𝑓 from 𝑜 + 1 4 to 𝑜 + 3 4

Temporal subcycling FDTD with iteration

25

iteration

Stability test of temporal subcycling

3D, coarse mesh size: 6x6x6, fine mesh size: 2x2x2

26

• P. Xu, PhD dissertation, 2017

Numerical Example: Cylinder Scattering

27

Using TO to enlarge the cylinder twice

TMz

Convergence test (TEz mode)

28

FDTD TO-FDTD 𝜇/Δ 40 20 𝑀1 Error 0.091 0.087 CPU time 3 1

Dual Grid FDTD

• The FDTD Yee lattice provides: 𝐹𝑜, 𝐼𝑜+1

2

• The Second Yee lattice provides: 𝐹𝑜+1

2, 𝐼𝑜

• Two grids coupled together in constitutive relations

– Magneto-electric material – Nonliearity

29

Spacetime Cloak

McCall & Kinsler, 2011

30

Free space Magneto-electric material 𝐹 = 𝛽𝐸 + 𝛾𝐶 𝐼 = 𝛾𝐸 + 𝛿𝐶 𝐹 = 1 𝜗0 𝐸, 𝐼 = 1 𝜈0 𝐶

Cloaked region

Magneto-electric material

𝐸𝑜 = 𝜗𝐹𝑜 + 𝛾𝐼𝑜 𝐶𝑜+1

2 = 𝛾𝐹𝑜+1 2 + 𝜈𝐼𝑜+1 2

𝐸𝑜 = 𝜗𝐹𝑜 + 𝛾 (3𝐼𝑜−1

2−𝐼𝑜−3 2)

2 𝐶𝑜+1

2 = 𝛾 (3𝐹𝑜 − 𝐹𝑜−1)

2 + 𝜈𝐼𝑜+1

2

Yee (with extrapolation) Dual Yee grids

magneto-electric material: Unstable

𝐸𝑜 = 𝜗𝐹𝑜 + 𝛾 ෡ 𝐼𝑜 𝐶𝑜+1

2 = 𝛾 ෠

𝐹𝑜+1

2 + 𝜈𝐼𝑜+1 2

Stable ෡ 𝐼 and ෠ 𝐹 are provided by the 2nd Yee grid (offset by

Δ𝑢 2 )

31

Spacetime cloak simulation

t x

Cornelius, Biro, & Liu, Optics Express, 2014

t x

Yee (with extrapolation) Dual Grid FDTD

32

Dual Grid FDTD for Nonlinear Case

𝐸 = 𝜁0(𝜁𝐹 + 𝜓 𝐹 2𝐹) = 𝜁0(𝜁 + 𝜓 𝐹 2)𝐹 (1)

• On the FDTD Yee lattice, we have 𝐹𝑜, 𝐸𝑜, 𝐼𝑜+1

2

• Use a second Yee lattice to get 𝐹𝑜+1

2, 𝐸𝑜+1 2, 𝐼𝑜

• Discretize constitutive equation at time level 𝑜 +

1 2:

𝑬𝒐+𝟐

𝟑 = 𝜁0 𝜁 + 𝜓 𝐹𝑜+1 2

2

𝑭𝒐+𝟐

𝟑

(2)

𝐸𝑜+1+𝐸𝑜 2

= 𝜁0 𝜁 + 𝜓 𝐹𝑜+1

2

2 𝐹𝑜+1+𝐹𝑜 2

(3)

• Then solve (3) for 𝐹𝑜+1:

𝐹𝑜+1 =

𝐸𝑜+1+𝐸𝑜 𝜁0 𝜁+𝜓 𝐹𝑜+1

2 2 − 𝐹𝑜

(4)

33

Soliton propagation in dispersive and nonlinear media

34

Lorentz Dispersion + Kerr Nonlinearity 𝐸 = 𝜁0𝜁𝐹 + 𝑄𝑀 + 𝑄𝑂𝑀 𝑄𝑂𝑀 = 𝜓 𝐹 2𝐹 𝑄𝑀 = 𝜕𝑞

2

𝜕𝑞

2 − 𝜕2 − 𝑗𝛿𝜕 𝐹

PML

Simulation results of Soliton in Lorentz + Kerr media

• FDTD + Newton iteration
• FDTD using Dual Grid (Explicit)
• Kerr effect is weak

35

Mesh size 1600x800

The Nonlinear Drude Model

Numerical Method:

(1) FDTD for E and B (the Maxwell equations); (2) a time-splitting method for J equation together with semi-implicit finite difference method and Lax-Wendroff.

Hydrodynamic Nonlinear Drude Model

• Liu, et. al, JCP, 2010

36

Electron Fluid Maxwell System

Fluid Equations for Electron Gas

Electron density

𝑣 velocity

Electron mass (constant) Electron charge (constant) Current density J = 𝑟𝑓𝑛𝑓𝑣

Maxwell’s Equations for Electromagnetic Wave

𝛿 electron collision rate

𝑞 electron pressure

37

Zhao, et. al., Photonics, 2015

Second Harmonic generation from Metamaterials

Log-log plot of the Fourier transform of the time history of Ex and Ey as functions of the wavelength.

SHG Conversion Efficiency ( × 10−11 ) Experimental Numerical ~ 2.00 1.02

38

Incident pulse Ex (Red) propagates along z Second Harmonic Wave Ey (Blue) Boundaries: x, y: Periodic; z: PML “U” shape metal: Nonlinear Drude

Summary

• Anisotropic media FDTD

– Fully Anisotropic, Coordinate Mapping, Temporal subcycling

• Dual grid FDTD

– Magneto-electric material, Dispersive and Nonlinear Kerr media

• Hydrodynamic Nonlinear Drude model

– Future work: metal-air interface, higher order solver

39

Thank You!