色散:交界面(2DTM)
前面範例有金屬材料交界面的反射率範例(交界面interface)，在這裡我們試著解說如何自行輸入電性金屬材料參數、及如何利用CompactFDTD求得色散曲線。
例如若有一組材料是eV單位，要在風行裡使用必須經過單位轉換，參考
S.-H. Chang, S. Gray, and G. Schatz, “Surface plasmon generation and light transmission by isolated nanoholes and arrays of nanoholes in thin metal films,” Opt. Express 13, 3150-3165 (2005).
http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-13-8-3150
裡有提到模擬金材料的Drude參數 500nm-1000nm 是
epsilon_infinite = 11.4577;
omegad =9.40274 eV;
gamad =0.0831375 ev;
轉成風行可用單位
ehbar = 1.51926751447914e+015; % e/hbar where hbar=h/(2*pi) and e=1.6e-19
omegad= 9.40274/27.2214/(2.41888*10^-17)=9.40274*1.51926751447914e+015= 1.428527742909359e+016 rad/s
gamad =0.0831375/27.2214/(2.41888*10^-17)=0.0831375*1.51926751447914e+015=1.263081029850095e+014 rad/s
那麼將這些值輸入進電性材料模型如下圖
這個範例、是僅只使用一個Drude模型。而內建的金材料如下圖，是使用了一個Drude模型加五個Lorentz模型。在有限差分時域法裡，當模型使用的越多，需更新方程將越多(請參考FDTD相關書籍)，做二維計算可能影響尚不明顯，但在三維計算中，需更新的方程式越多、記憶體與計算時間將大大拉長。因此，對您所感興趣的波段做模擬(可能需重新擬合色散曲線)，使用較少的模擬的模型，將能大大的減少計算時間。
設定為softsource 做為空間一點的激發源
從設定圖可知我們所輸入的波數k(wave number)從k=1.571e6 rad/m 至 k=20.944e6 rad/m 之間切20個等分
因為是要計算單一圓柱體散射吸收頻譜，因此四周圍皆為使用吸收層PML
勾選 ☑ 設計結構 (Design Structures)
勾選 ☑ E,H,V,I (FFT) : 量測電場、磁場、電壓、電流
【指令輸入方塊】輸入以下程式碼，然後按【執行】，確定好結構無誤的話後再按【輸出結構】
%=================矩形 Brick===========================
xstart=ib/2*gdx;
xend=ib*gdx;
ystart=1*dy;
yend=1*dy;
zstart=1*gdz;
zend=1*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===========================
%------------探測磁場 1 Probe Magnetic 1------------
% direction 1(Hx),2(Hy),3(Hz),4(|H|)
direction=2;
x1=(ib/2+1)*gdx;
x2=(ib/2+1)*gdx;
y1=(jcenter)*gdy;
y2=(jcenter)*gdy;
z1=(kcenter)*gdz;
z2=(kcenter)*gdz;
probe_magnetic(1,:)=[x1,x2,y1,y2,z1,z2,direction];
%----------------------------------
%------------探測磁場 2 Probe Magnetic 1------------
% direction 1(Hx),2(Hy),3(Hz),4(|H|)
direction=2;
x1=(ib/3+1)*gdx;
x2=(ib/3+1)*gdx;
y1=(jcenter)*gdy;
y2=(jcenter)*gdy;
z1=(kcenter)*gdz;
z2=(kcenter)*gdz;
probe_magnetic(2,:)=[x1,x2,y1,y2,z1,z2,direction];
%----------------------------------
%------------探測磁場 3 Probe Magnetic 1------------
% direction 1(Hx),2(Hy),3(Hz),4(|H|)
direction=2;
x1=(ib/4+1)*gdx;
x2=(ib/4+1)*gdx;
y1=(jcenter)*gdy;
y2=(jcenter)*gdy;
z1=(kcenter)*gdz;
z2=(kcenter)*gdz;
probe_magnetic(3,:)=[x1,x2,y1,y2,z1,z2,direction];
%----------------------------------
%------------探測電場 1 Probe Electric 1------------
% direction 1(Ex),2(Ey),3(Ez),4(|E|)
direction=1;
x1=(icenter-1)*gdx;
x2=(icenter-1)*gdx;
y1=(jcenter-1)*gdy;
y2=(jcenter)*gdy;
z1=(kcenter-1)*gdz;
z2=(kcenter+1)*gdz;
probe_electric(1,:)=[x1,x2,y1,y2,z1,z2,direction];
%----------------------------------
前面範例有金屬材料交界面的反射率範例(交界面interface)，在這裡我們試著解說如何自行輸入電性金屬材料參數、及如何利用CompactFDTD求得色散曲線。
例如若有一組材料是eV單位，要在風行裡使用必須經過單位轉換，參考
S.-H. Chang, S. Gray, and G. Schatz, “Surface plasmon generation and light transmission by isolated nanoholes and arrays of nanoholes in thin metal films,” Opt. Express 13, 3150-3165 (2005).
http://www.opticsinfobase.org/oe/abstract.cfm?uri=oe-13-8-3150
裡有提到模擬金材料的Drude參數 500nm-1000nm 是
epsilon_infinite = 11.4577;
omegad =9.40274 eV;
gamad =0.0831375 ev;
轉成風行可用單位
ehbar = 1.51926751447914e+015; % e/hbar where hbar=h/(2*pi) and e=1.6e-19
omegad= 9.40274/27.2214/(2.41888*10^-17)=9.40274*1.51926751447914e+015= 1.428527742909359e+016 rad/s
gamad =0.0831375/27.2214/(2.41888*10^-17)=0.0831375*1.51926751447914e+015=1.263081029850095e+014 rad/s
那麼將這些值輸入進電性材料模型如下圖
要輸出檔案則勾選 ☑穩態場輸出 (Output Steady State - Fields) -(Poynting DFT)。
運算結束後
【分析結果】
參考網頁：可用變數
N_Probed_E_Td(觀測位置編號,計算步數,分割數) %量測電場(時域) Probed Electric (Time-domain)
N_Probed_H_Td(觀測位置編號,計算步數,分割數) %量測磁場(時域) Probed Magnetic (Time-domain)
N_Probed_E_Fd(觀測位置編號,計算步數,分割數) %量測電場(頻域) Probed Electric (Frequency-domain)
N_Probed_H_Fd(觀測位置編號,計算步數,分割數) %量測磁場(頻域) Probed Magnetic (Frequency-domain)
★這裡您只需要調整兩個參數即可，第一個是觀測場的輸入，我們前面設定四個觀測場，其中三個磁場位置，一個電場位置
★要做模態判定，一次只能輸入電場或磁場其中一種。
★Probes=N_Probed_H_Td; 這個方式，是將三個觀測磁場都輸入，Mode_Analysis會將三個觀測的磁場做平均化來分析尋找模態
★Probes=N_Probed_H_Td(3,:,:); 這個方式，是將第三個觀測磁場都輸入，Mode_Analysis只會用第三個觀測的磁場來做分析尋找模態
★FieldMax 設定臨界值，高於此值的將視為找到的模態，小於此值的將過濾掉。若找到太多雜訊模態，調整前面的係數0.01變大來過濾
★max_band_number 在一輸入波數k(wave number)位置上，最多找幾個模態
★lambda0 中心入射波長
Probes=N_Probed_H_Td; %這裡將我們觀測的電場或磁場輸入
FieldMax=0.025*max(Probes(:)); %這裡我們設定臨界值，
max_band_number=6;
Mode_Analysis(Probes,FieldMax,max_band_number);
figure(1);
plot(betaz,Eigen_omega(:,:),'r*'); ylabel('ω');xlabel('k')
figure(2)
pcolor(Eigen_omega_record');shading flat;colorbar; %colormap('hot');
%===理論解析解: 金屬介電質色散關係 metal-dielectri interface dispersion relation==============
RunNumber=100;
epsInf =11.4577;
WDp = 9.40274/27.2214/(2.41888*10^-17);
gama = 0.0831375/27.2214/(2.41888*10^-17);
freq0=3e8/lambda0;
omegas=linspace(0.1,1.5,RunNumber)*2*pi*freq0;
eps2= 1;
eps1= epsInf-WDp^2./(omegas.^2+omegas*gama*j);
kxx=(omegas).*((eps2*eps1)./(eps1+eps2) ).^0.5;
figure(1);hold on;title('金屬介電質色散關係')
plot(kxx/cc,kxx,'k','LineWidth',1,'MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',1);
plot(kxx/cc,omegas,'b','LineWidth',1,'MarkerEdgeColor','k','MarkerFaceColor','k','MarkerSize',1);
hold off;
%===================================
我們可比較一下上圖左與右，分別調整臨界值如圖所示，可看出調低臨界值會找到更多的模態。有些單個波數k找到對應的兩個ω，但不見得都是對的。這個範例右邊顯示有些對應的模態出現在light-line色散曲線上(黑色實線)，但金屬與介電質交界面的色散並不是在這條線上。我們再調低一點如下圖，會將許多雜訊判定為模態。因此所過濾出的模態要很小心，不見得都能用。
下面這我們看一下
cc = 299792458;
N_spare=1; %k=2*pi/300e-9=1.571x10^6; lamda_start=300nm,
t1=strcat('k=',num2str(betaz(N_spare)));
figure(2)
subplot(2,1,1);
plot(1*dt:dt:nmax*dt,N_Probed_H_Td(1,:,N_spare));
subplot(2,1,2);
plot(freqFFT,N_Probed_H_Fd(1,:,N_spare));
N_spare=20; %k=2*pi/4e-6=20.944x10^6; lamda_end=4um,
t20=strcat('k=',num2str(betaz(N_spare)));
figure(2);
subplot(2,1,1);
hold on;
plot(1*dt:dt:nmax*dt,N_Probed_H_Td(1,:,N_spare),'r');
title('Time-Domain');
xlabel('s');
legend(t1,t20);
hold off
subplot(2,1,2);
hold on;
plot(freqFFT,N_Probed_H_Fd(1,:,N_spare),'r');
title('Frequency-Domain');
xlabel('Hz');
legend(t1,t20);
hold off
