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: 𝐸 = 𝜗 𝐹 + 𝑸 current density: 𝐾 = 𝑒𝑸 𝑴 𝑒𝑢 𝐶 = 𝜈 𝐼 polarization 𝑄 = 𝑄 𝑀 + 𝑄 𝑂𝑀 3

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 2 𝜕 𝑞 𝜗 𝜕 = 1 − 𝜕 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 𝜖𝐼 Hz 𝜈 𝜖𝑢 = −𝛼 × 𝐹 𝜗 𝜖𝐹 Hx 𝜖𝑢 = 𝛼 × 𝐼 Ez Hy • Y ee ’66, Taflove ’75 • Staggered Cartesian grid in space & time Ey • Centered Finite Difference and Leapfrog Ex • Non-dissipative • Divergence free • Robust and easy to implement • Variety of materials: Dielectric, Anisotropic, Dispersive, Nonlinear, etc 5

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,𝑘,𝑙 + 𝜊 𝑦𝑧 𝑬 𝒛,𝒋+𝟐 𝟑,𝒌,𝒍 + 𝜊 𝑦𝑨 𝑬 𝒜,𝒋+𝟐 𝟑,𝒌,𝒍 𝑭/𝑬 𝒜,𝒋,𝒌,𝒍+ 𝟐 𝟑 𝑭/𝑬 𝒛,𝒋+𝟐,𝒌+ 𝟐 Staggered Yee Lattice 𝟑 ,𝒍 𝑭/𝑬 𝒚,𝒋+𝟐 7 𝟑,𝒌,𝒍

Update equation for anisotropic material Constitutive Equations: 𝐹 = 𝜗 −1 𝐸, 𝐼 = 𝜈 −1 𝐶 𝐹 𝑦 = 𝜊 𝑦𝑦 𝐸 𝑦 + 𝜊 𝑦𝑧 𝐸 𝑧 + 𝜊 𝑦𝑨 𝐸 𝑨 𝐹 𝑦,𝑗+1 2,𝑘,𝑙 = 𝜊 𝑦𝑦 𝐸 𝑦,𝑗+1 2,𝑘,𝑙 + 𝝄 𝒚𝒛 𝑬 𝒛,𝒋+𝟐 𝟑,𝒌,𝒍 + 𝜊 𝑦𝑨 𝐸 𝑨,𝑗+1 2,𝑘,𝑙 • Non-averaging (1 st order) 𝐸 𝑧,𝑗+1 2,𝑘,𝑙 = 𝐸 𝑧,𝑗,𝑘+1 2,𝑙 • Averaging method 1 (unstable) 𝐸 𝑧,𝑗,𝑘+1 2,𝑙 + 𝐸 𝑧,𝑗,𝑘−1 2,𝑙 + 𝐸 𝑧,𝑗+1,𝑘+1 2,𝑙 + 𝐸 𝑧+1,𝑗,𝑘−1 2,𝑙 𝐸 𝑧,𝑗+ 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) 2,𝑘,𝑙 = 𝜊 𝑦𝑧,𝑗,𝑘,𝑙 𝐸 𝑀 + 𝜊 𝑦𝑧,𝑗+1,𝑘,𝑙 𝐸 𝑆 𝜊 𝑦𝑧,𝑗+1 2,𝑘,𝑙 𝐸 𝑧,𝑗+1 2 𝐸 𝑧,𝑗,𝑘+1 2,𝑙 + 𝐸 𝑧,𝑗,𝑘−1 𝐸 𝑧,𝑗+1,𝑘+1 2,𝑙 + 𝐸 𝑧,𝑗+1,𝑘−1 2,𝑙 2,𝑙 𝐸 𝑀 = , 𝐸 𝑆 = 2 2 9

Invariant Coordinate Transformation (Transformation Optics) Coordinate 𝜈 ′ 𝜖𝐼 ′ 𝜈 𝜖𝐼 Transformation 𝜖𝑢 = −𝛼 ′ × 𝐹 ′ 𝜖𝑢 = −𝛼 × 𝐹 from (x,y,z) to ( x’,y’,z’) 𝜗 ′ 𝜖𝐹 ′ 𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 ′ × 𝐼 ′ 𝜖𝑢 = 𝛼 × 𝐼 anisotropic isotropic 𝐹 ′ = 𝛭 𝑈 𝐹 , 𝜗 ′ = 𝛭 𝛭 −1 𝜗 𝛭 −𝑈 𝐼 ′ = 𝛭 𝑈 𝐼 , 𝜈 ′ = 𝛭 𝛭 −1 𝜈𝛭 −𝑈 𝛭 is the Jacobian matrix Teixeira ‘99, Pendry ’06, Leonhardt ‘06 10

Coordinate Stretching 𝑠′ = 𝑔 𝑠 Physical domain Computational domain 11

Mapping function 𝑠′ 𝑠′ = 𝑔(𝑠) 𝑆 2 𝑆 1 ′ 𝑠 𝑝 𝑆 1 𝑆 2 ′ ≤ 𝑠 ≤ 𝑆 2 𝑔 maps 0 < 𝑠 ≤ 𝑆 1 to 0 < 𝑠 ′ ≤ 𝑆′ 1 and 𝑆 1 ≤ 𝑠 ≤ 𝑆 2 to 𝑆 1 𝑔 continuous and 𝑔’ continuous, for example, by using cubic spline 12

Transformed mesh and material mapping Physical domain Computational domain 𝜗 11 not 𝜗 = 1 constant 13

Transformation based Maxwell Solver Input Grid Mapping: from physical Pre-Processing domain to a computational domain Maxwell Solver EM Anisotropic material inverse mapping back to physical Post-Processing domain Output Liu, et. al. J. Comput. Phys. 2014 14

TO-FDTD for dispersive material: the 𝐾 version (𝑦, 𝑧, 𝑨) (𝑦’, 𝑧’, 𝑨’) Coordinate 𝝂 ′ 𝜖𝐼 ′ 𝜈 𝜖𝐼 𝜖𝑢 = −𝛼 ′ × 𝐹 ′ Transformation 𝜖𝑢 = −𝛼 × 𝐹 from (x,y,z) to ( x’,y’,z’) mapping 𝝑 ′ 𝜖𝐹 ′ 𝜗 𝜖𝐹 𝜖𝑢 = 𝛼 ′ × 𝐼 ′ − 𝐾 ′ 𝜖𝑢 = 𝛼 × 𝐼 − 𝐾 𝜖𝐾 𝜖𝐾′ 2 𝐹 𝜖𝑢 + 𝛿𝐾 = 𝜗 0 𝜕 𝑞 𝜖𝑢 + 𝛿𝐾′ = Λ −1 Λ𝜗 0 Λ 𝑈 𝜕 𝑞 2 𝐹′ 𝐹 ′ = 𝛭 𝑈 𝐹 , 𝝑 ′ = 𝛭 𝛭 −1 𝜗 𝛭 −𝑈 𝐼 ′ = 𝛭 𝑈 𝐼 , 𝝂 ′ = 𝛭 𝛭 −1 𝜈𝛭 −𝑈 𝐾 ′ = Λ 𝛭 −1 𝐾 , 𝛭 is the Jacobian matrix from (𝑦’, 𝑧’, 𝑨’) 𝑢𝑝 (𝑦, 𝑧, 𝑨) 15

TO-FDTD for dispersive material: the 𝑄 version (𝑦, 𝑧, 𝑨) (𝑦’, 𝑧’, 𝑨’) 𝜖𝐶 ′ 𝜖𝐶 𝜖𝑢 = −𝛼 ′ × 𝐹 ′ 𝜖𝑢 = −𝛼 × 𝐹 Coordinate Transformation 𝜖𝐸 𝜖𝐸′ from (x,y,z) to ( x’,y’,z’) 𝜖𝑢 = 𝛼 ′ × 𝐼 ′ 𝜖𝑢 = 𝛼 × 𝐼 mapping 𝐶 = 𝜈𝐼 𝐶′ = 𝝂′𝐼′ 𝐸 = 𝜗 0 𝜗 ∞ 𝐹 + 𝑄 𝐸′ = 𝝑′ 𝐹′ + 𝑄′ 𝜖 2 𝑄 2 𝜖𝑢 2 + 𝛿 𝜖𝑄 𝜖𝑢 = 𝜕 𝑞 𝜖 2 𝑄′ 2 𝜖𝑢 2 + 𝛿 𝜖𝑄′ 𝜖𝑢 = 𝜕 𝑞 𝐹 𝐹′ 𝜗 ∞ 𝜗 ∞ 𝐹 ′ = 𝛭 𝑈 𝐹 , 𝝑 ′ = 𝛭 𝛭 −1 𝜗 0 𝜗 ∞ 𝛭 −𝑈 𝐼 ′ = 𝛭 𝑈 𝐼 , 𝝂 ′ = 𝛭 𝛭 −1 𝜈 𝛭 −𝑈 𝑄 ′ = 𝛭 𝛭 −1 𝑄 , 16 𝛭 is the Jacobian matrix from (𝑦’, 𝑧’, 𝑨’) 𝑢𝑝 (𝑦, 𝑧, 𝑨)

Ring Simulation using grid transformation FDTD mesh TO Physical Domain TO Comput. domain (Cartesian, isotropic) (Non-Rectangular mesh) (Cartesian, anisotropic) 17

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

Metallic Bowtie simulation Mesh 𝑶 𝒛 CPU Time FDTD 400 40 ( Δ = 4𝑜𝑛) FDTD 800 320 ( Δ = 2𝑜𝑛) 20 nm gap TO-FDTD 80 1 1600 nm = 2 𝜇 19 FDTD: Δ x = 4 nm FDTD: Δ x = 2 nm TO-FDTD: Δ x = 20 nm

Subgridding/Adaptive Mesh Refinement 20

Beam Propagation 21

Spatial and Temporal Subgridding 𝚬𝒚, 𝚬𝒛, 𝚬𝒖 𝚬𝒚 𝟑 , 𝚬𝒛 𝟑 , 𝚬𝒖 𝟑 22

Temporal Subgridding + irregular mesh 𝚬𝒚, 𝚬𝒛, 𝚬𝒖 𝚬𝒖 nonrectangular mesh, 𝟑 23

Temporal subcycling with iterations 𝐼/𝐹 𝑑𝑝𝑏𝑠𝑡𝑓 𝐼/𝐹 𝑔𝑗𝑜𝑓 Update coarse mesh for Δ𝑢 several times (near interface) Update fine mesh for Δ𝑢/2 Update fine mesh for Δ𝑢/2 24