Metal-Dielectric-Interface

中文/English

Problem:


To simulate a plane wave incident into the metal layer, and observe the reflection spectrum (wavelength = 100 ~ 1000 nm)

And compared with analytical solution

Results :



Define Parameters

A.Main-Setup

1. Define total x, y, z length that equal to CAD, and then set a appropriate resolution.

2. Press【建立網格 Cad to Grids (Create)】button to create simulation grids.


Here demonstrate how to simulate by using the 2D TM (Ex Ez Hy),

define total size x= 50 nm, z= 500 nm; delta= 5 nm




B.Boundary Conditions


A. Define the boundary conditions

B. Press the【創建 (Create)】button => Create the total grid size (Include the boundary condition & add space)



C.Wave & Observation

The incident wave is linear polarization in y axis (pis=0) from bottom side(theta =0).

Reflection Spectrum observation : - Z direction

Transmission Spectrum observation : + Z direction


instruction webpage: DFT (Observation Range)



D. Geometries. Design


Checked out the simulation geometries, and then press 【輸出 Output】button to output the geometries files

Data_Materials.csv // (simulation index of structures)


In this example the "矩形 Brick" structure is using the built-in function.

instruction webpage: Built-in Function 


%=================矩形 Brick===========================
xstart=1*gdx;
xend=10*gdx;
ystart=1*dy;
yend=1*dy;
zstart=50*gdz;
zend=100*gdz;
nindex=1^2;
sigma=0;
choice='E_Model1'; %E_Iso,PEC,M_Iso,PMC,E_Model1,M_Model1,EM_Model1
gridtype=-1;    %
Iso_Brick(choice,gridtype,nindex,sigma,xstart,xend,ystart,yend,zstart,zend)

%=================矩形 Brick===========================




E. Electric & Magnetic Models

The dielectric constant of gold in this simulation can be obtain by press "預覽 Preview" button in above figure




F.Spectrum & Field Pattern

1. Set the Spectrum Analysis & Wavelength range

2. Check out the source




G.Calculation

1. Set the simulation time setp

2. Save parameters

3. Calculation



H.Result Analysis

★(Result Analysis)：

頻譜 (Spectrum) : Analysis the spectrum


This is the FDTD 2D TM simulation result, 3D, SplitField Method can get the same results.



Theoretical (Analytical) Solution code :


clear ; clc;
%讀取電性色散參數檔案
a=textread('Data.Mats','%s','delimiter', 's');
Debyemodel1=str2num(char(a((1))));
Debyemodel2=str2num(char(a((2))));
Debyemodel3=str2num(char(a((3))));
Debyemodel4=str2num(char(a((4))));
Debyemodel5=str2num(char(a((5))));
Drudemodel1=str2num(char(a((6))));
Drudemodel2=str2num(char(a((7))));
Drudemodel3=str2num(char(a((8))));
Drudemodel4=str2num(char(a((9))));
Drudemodel5=str2num(char(a((10))));
Lorentzmodel1=str2num(char(a((11))));
Lorentzmodel2=str2num(char(a((12))));
Lorentzmodel3=str2num(char(a((13))));
Lorentzmodel4=str2num(char(a((14))));
Lorentzmodel5=str2num(char(a((15))));
Pademodel1=str2num(char(a((16))));
Pademodel2=str2num(char(a((17))));
Pademodel3=str2num(char(a((18))));
Pademodel4=str2num(char(a((19))));
Pademodel5=str2num(char(a((20))));
 
 
 
hDebyemodel1=str2num(char(a((21))));
hDebyemodel2=str2num(char(a((22))));
hDebyemodel3=str2num(char(a((23))));
hDebyemodel4=str2num(char(a((24))));
hDebyemodel5=str2num(char(a((25))));
hDrudemodel1=str2num(char(a((26))));
hDrudemodel2=str2num(char(a((27))));
hDrudemodel3=str2num(char(a((28))));
hDrudemodel4=str2num(char(a((29))));
hDrudemodel5=str2num(char(a((30))));
hLorentzmodel1=str2num(char(a((31))));
hLorentzmodel2=str2num(char(a((32))));
hLorentzmodel3=str2num(char(a((33))));
hLorentzmodel4=str2num(char(a((34))));
hLorentzmodel5=str2num(char(a((35))));
Pademodel1=str2num(char(a((36))));
Pademodel2=str2num(char(a((37))));
Pademodel3=str2num(char(a((38))));
Pademodel4=str2num(char(a((39))));
Pademodel5=str2num(char(a((40))));
 

epsinf1=str2num(char(a((71))));
sigmametal1=str2num(char(a((72))));
epsinf2=str2num(char(a((73))));
sigmametal2=str2num(char(a((74))));
epsinf3=str2num(char(a((75))));
sigmametal3=str2num(char(a((76))));
epsinf4=str2num(char(a((77))));
sigmametal4=str2num(char(a((78))));
epsinf5=str2num(char(a((79))));
sigmametal5=str2num(char(a((80))));
 
hepsinf1=str2num(char(a((81))));
hsigmametal1=str2num(char(a((82))));
hepsinf2=str2num(char(a((83))));
hsigmametal2=str2num(char(a((84))));
hepsinf3=str2num(char(a((85))));
hsigmametal3=str2num(char(a((86))));
hepsinf4=str2num(char(a((87))));
hsigmametal4=str2num(char(a((88))));
hepsinf5=str2num(char(a((89))));
hsigmametal5=str2num(char(a((90))));
 
gg=111;
%將各項係數分配入變數
 
for i=1:Debyemodel1
Detp1(i)=str2num(char(a((gg))));
gamabye1(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Debyemodel2
Detp2(i)=str2num(char(a((gg))));
gamabye2(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Debyemodel3
Detp3(i)=str2num(char(a((gg))));
gamabye3(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Debyemodel4
Detp4(i)=str2num(char(a((gg))));
gamabye4(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Debyemodel5
Detp5(i)=str2num(char(a((gg))));
gamabye5(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
 
for i=1:Drudemodel1
Wd1(i)=str2num(char(a((gg))));
gamad1(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Drudemodel2
Wd2(i)=str2num(char(a((gg))));
gamad2(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Drudemodel3
Wd3(i)=str2num(char(a((gg))));
gamad3(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Drudemodel4
Wd4(i)=str2num(char(a((gg))));
gamad4(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
for i=1:Drudemodel5
Wd5(i)=str2num(char(a((gg))));
gamad5(i)=str2num(char(a((gg+1))));
gg=gg+2;
end
 

for i=1:Lorentzmodel1
eps_p1(i)=str2num(char(a((gg))));
Wp1(i)=str2num(char(a((gg+1))));
gamap1(i)=str2num(char(a((gg+2))));
gg=gg+3;
end
for i=1:Lorentzmodel2
eps_p2(i)=str2num(char(a((gg))));
Wp2(i)=str2num(char(a((gg+1))));
gamap2(i)=str2num(char(a((gg+2))));
gg=gg+3;
end
for i=1:Lorentzmodel3
eps_p3(i)=str2num(char(a((gg))));
Wp3(i)=str2num(char(a((gg+1))));
gamap3(i)=str2num(char(a((gg+2))));
gg=gg+3;
end
for i=1:Lorentzmodel4
eps_p4(i)=str2num(char(a((gg))));
Wp4(i)=str2num(char(a((gg+1))));
gamap4(i)=str2num(char(a((gg+2))));
gg=gg+3;
end
for i=1:Lorentzmodel5
eps_p5(i)=str2num(char(a((gg))));
Wp5(i)=str2num(char(a((gg+1))));
gamap5(i)=str2num(char(a((gg+2))));
gg=gg+3;
end
 
cc=2.99792458e8;
muz=4.0*pi*1.0D-7;
epsz=1.D0/(muz*cc*cc);
 
lambda=linspace(100,1000,500)*1e-9;
twopic = 1.883651567308853e+09; % twopic=2*pi*c where c is speed of light
omega = twopic.*(lambda.^(-1)); % angular frequency of light (rad/s)
 

%=================================================
epsilon=epsinf1-sigmametal1./(j*omega*epsz);
%debye
for i=1:Debyemodel1
epsilon=epsilon-Detp1(i)./(1-j*omega*gamabye1(i));
end
%drude
for i=1:Drudemodel1
epsilon=epsilon-Wd1(i)^2./(omega.^2+j*omega*gamad1(i));
end
%lorentz
for i=1:Lorentzmodel1
    epsilon=epsilon-eps_p1(i)^2./(omega.^2+j.*omega*gamap1(i)-Wp1(i)^2);
end
 
% epsilon_r 對 lambda
figure(1);hold on;
plot(lambda,real(epsilon),'r','linewidth',3);
plot(lambda,imag(epsilon),'b','linewidth',3);
hold off;
title('Au');
xlabel('\lambda (nm)');
ylabel('\epsilon r');
 
% 反射率頻譜
n=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)+real(epsilon))/2 );
k=sqrt( (sqrt(real(epsilon).^2+imag(epsilon).^2)-real(epsilon))/2 );
r=(n0-n1)./(n0+n1);
figure(1);
plot(lambda,abs(r).^2,'^b','linewidth',1);
title('R');
xlabel('\lambda (nm)');
ylabel('R');