今天看到有些同學在求chern number的計算方法，正好前段時間我有花時間做過相關的工作，本著科研共享的精神，來分享一下我的理解。

由於我不是專業做這個方向的研究人員，所以我就不能做更多理論方面的解釋了，直接給大家一步步說明我的演算法就好了。0.0

準備工具：COMSOL +MATLAB

參考文獻：PRL 100,013905 (2008) + J.phys.Soc.Jpn. vol,74 no,6 (2005)

1、利用COMSOL求解光子晶體結構的本徵頻率：把第一布里淵區等分成n1*n2個點，點命名為N（ki,kj）（i=1:n1,j=1:n2）好了，於是可以得到n1*n2個本徵頻率解，保存他們的電場或磁場分布(按介質的不同保存，比如說是Si和Air就分別保存兩個區域的場分布)。特別注意的是由於求解器的問題有可能有些結構在算本徵頻率的時候，同一條能帶的解的位置會不一樣，例如我想求4條能帶的chern number，但是N14和N20（只是舉例）的解並不是正好跟能帶位置對應的地方（解的位置~=能帶位置的話會影響後面的計算）。所以建議每次求解完之後簡單的在COMSOL中把這些解放在圖中畫一下，仔細比對一下（不用把n1*n2個點都畫出來，只用畫簡約布里淵區邊界的那些點就足夠精確了（從效率上來說，當然時間夠有興趣可以導出來在其他工具中畫出三維的（k1,k2,freq）圖來確認））。

2、得到了n1*n2個本徵頻率下的電場/磁場強度分布---->導出保存為EH1.txt+EH2.txt(兩種介質)文件（注意每個介質的實部虛部的導出，我是比較麻煩每次分實部虛部導出兩個文件（總共4個文件）再相加的）。這裡感覺我可以多說一句，對於那些COMSOL不是很熟悉的人。首先COMSOL導出的電場數據我假設是4條能帶11*11=121個點即會有484個不同的場分布，每個場分布中假設有10000個點（這個跟網格劃分有關），則我導出的數據為10000*484的矩陣，加上COMSOL的位置（x,y）信息則導出數據為10000*486的矩陣。

3、在得到EH1.txt+EH2.txt文件之後下一步就是對這兩個10000*486的矩陣進行插值處理了。首先在進行插值之前要把各個能帶之間的數據分割開來，在這個案例中能帶數為4，則對3：end列的數據每隔4個為一個能帶的場分布（即3,7,11，...為第一能帶，4,8,12，...為第二能帶）命名為Eh1,Eh2,Eh3,Eh4（Air）+Eh5,Eh6,Eh7,Eh8(Si)。（Eh為10000*121的矩陣對應了121個不同的點N（ki,kj））。

以求第一個能帶的chern number為例。首先對Eh1,Eh5進行插值，griddata命令即可。插值以後的場分布分別稱為eH1和eH5。

4、理論上的東西大家可以參考上面推薦的兩個文獻資料。簡單來說是， Int1:conj(eH1_{(i,j)})*eH1_{(i+1,j)}+conj(eH5_{(i,j)})*eH5_{(i+1,j)} 做面積分，

同理

Int2:conj(eH1_{(i+1,j)})*eH1_{(i+1,j+1)}+conj(eH5_{(i+1,j)})*eH5_{(i+1,j+1)} 做面積分， Int3:conj(eH1_{(i,j+1)})*eH1_{(i+1,j+1)}+conj(eH5_{(i,j+1)})*eH5_{(i+1,j+1)} 做面積分，

Int4:conj(eH1_{(i,j)})*eH1_{(i,j+1)}+conj(eH5_{(i,j)})*eH5_{(i,j+1)} 做面積分，

然後 C_{(i,j)}=ln(frac{Int1}{abs(Int1)}*frac{Int2}{abs(Int2)}/frac{Int3}{abs(Int3)}/frac{Int4}{abs(Int4)}) ，

最後 Chern number=sum_{1}^{10}sum_{1}^{10}{C_{(a,b)}}/2/pi/1i

很多不夠準確的地方大家見諒，如果有什麼錯誤後續會改，希望大家科研順利。合作共享是科研的出路！

推薦閱讀：

相关文章