計算全流程)
1. 蜂窩晶格為什么值得算狄拉克錐與拓撲光子學的切入點做拓撲光子學研究的人遲早都會碰到蜂窩晶格光子晶體。這個體系幾乎就是為演示“如何從能帶計算走向拓撲不變量”而生的。它的晶格結構和石墨烯一模一樣由兩套三角子格穿插構成所以天然具備狄拉克錐而只要破壞了空間反演對稱性狄拉克點處就會打開一個帶隙進而可以定義陳數(shù)、研究邊緣態(tài)。很多初學者以為能帶圖算出來就萬事大吉了實際上從Comsol模型到Matlab里的陳數(shù)結果中間隔著不少容易踩坑的環(huán)節(jié)。我在這篇文章里會把這套流程完整地拆開講蜂窩晶格單胞怎么建、布洛赫周期邊界怎么設、k路徑怎么掃、本征場怎么導出來最后又怎么用Matlab腳本從實空間場分布出發(fā)算出Berry曲率和陳數(shù)。內(nèi)容適合有有限元基礎、想做拓撲光子學數(shù)值驗證的研究生也適合想快速搭一套“能帶計算拓撲不變量”工作流的工程師。里面所有參數(shù)和腳本邏輯都是我在實際項目中反復驗證過的方案。1.1 從能帶折疊到狄拉克錐蜂窩結構的關鍵是兩套格子蜂窩晶格在實空間里的形象是一張六角網(wǎng)但計算時并不推薦直接畫一個六邊形單胞而是把它看成兩套三角子格。設晶格常數(shù)為a兩個基矢取為a1 (√3/2·a, 1/2·a)a2 (√3/2·a, ?1/2·a)A子格放在原點B子格放在最近鄰位置比如(0, a/√3)。兩套子格在幾何上完全等價時系統(tǒng)的空間反演對稱性沒有被破壞倒易空間第一布里淵區(qū)的K和K′點會出現(xiàn)兩重簡并能帶形成線性交叉也就是狄拉克錐。對應到光子晶體里如果A、B位置放的介質柱半徑相同、介電常數(shù)相同那么TE或TM模式的能帶在K點附近一定會呈錐形。這個結構在Comsol中非常好驗證把A、B柱半徑設成相等掃出一條能帶曲線你會在帶隙關閉的狀態(tài)下看到兩條帶恰好接觸。這里有一個很容易忽略的前提——幾何對稱性必須被數(shù)值網(wǎng)格忠實地保留。如果網(wǎng)格剖分破壞了C3旋轉對稱性K點的簡并度會被人為地打開狄拉克錐就會變成一條假帶隙拓撲結論自然全部作廢。所以后面講網(wǎng)格剖分的時候我專門花了一節(jié)來說這個事。1.2 對稱性破缺決定帶隙能否打開和陳數(shù)直接相關一旦A、B兩套格子的參數(shù)不再相同比如讓A柱半徑rA大于B柱半徑rB空間反演對稱性就被破壞。K點的兩重簡并上升原先線性交叉的兩條能帶被分開形成一段帶隙。這段帶隙的重要性在于它給了我們一個干凈的“頻率窗口”可以定義某條占據(jù)帶通常取帶隙下方的一組帶的拓撲性質。這里需要特別說明光子晶體中“陳數(shù)”的定義方式。光子本身是玻色子不存在電子那樣費米面和占據(jù)態(tài)但我們可以把布洛赫模式當成一個參數(shù)空間中的向量場對某個頻帶的態(tài)進行投影構造出等價于量子Hilbert空間的結構。具體做法是只看帶隙以下那條帶或者那組互相簡并的帶把它當成“占據(jù)帶”在布里淵區(qū)內(nèi)對Berry曲率積分得到的整數(shù)就是陳數(shù)。這個整數(shù)直接和體邊對應關系掛鉤帶隙內(nèi)存在單向傳輸?shù)倪吘墤B(tài)邊緣態(tài)的數(shù)量和方向就由陳數(shù)決定。所以整個項目的邏輯鏈其實很簡潔蜂窩晶格給出狄拉克錐對稱性破缺打開帶隙帶隙的非平庸性質由陳數(shù)描述。下一步Comsol負責給出能帶和本征場Matlab負責把場的相位信息提煉成拓撲判據(jù)。1.3 這個體系的“陳數(shù)”到底在描述什么很多人剛接觸陳數(shù)時總被抽象公式嚇到。我習慣用一個類比來解釋把布里淵區(qū)想象成一個封閉曲面每個k點處都有一個本征向量也就是布洛赫函數(shù)的周期部分。這些向量會隨著k變化而指向不同方向。陳數(shù)描述的就是這個向量場在封閉曲面上纏繞的“圈數(shù)”。就像莫比烏斯帶和普通紙帶的扭轉數(shù)量差一個整數(shù)沒法通過連續(xù)形變改變陳數(shù)也是一個無法連續(xù)變化的整數(shù)拓撲不變量。對于蜂窩晶格光子晶體陳數(shù)通常不止一個值。如果介質柱半徑rA略大于rB某幾條頻帶可能得到1或?1的陳數(shù)而rA小于rB時符號反轉。符號的意義直接映射到邊緣態(tài)的傳播方向——正陳數(shù)的邊界上邊緣態(tài)沿某一方向單向傳播負陳數(shù)則反過來。這個特性就是拓撲保護的根源即使路徑上有彎折、有缺陷只要帶隙沒有被填滿單向傳播就不會被背散射破壞。我在后面會用這個性質和邊緣態(tài)能帶做交叉驗證。2. Comsol建模與MPH工作流單胞怎么畫、邊界怎么設才能不出偽模2.1 先澄清“MPH”這個說法在Comsol生態(tài)里的實際含義很多人以為“MPH算法”是什么高深的學術名詞實際上在Comsol的語境中.mph就是模型文件的擴展名。一份.mph文件包含了整個項目的一切幾何、物理接口、網(wǎng)格、研究步驟、結果和后處理。所謂“含MPH算法”合理的解讀就是依托Comsol的.mph模型文件來組織整套計算流程的方法。我個人喜歡把“MPH”拆成三個動作Modeling建模與參數(shù)化、Processing求解與數(shù)據(jù)提取、Harvesting后處理與拓撲判據(jù)計算。這樣工程文件的結構會很清晰模型目錄下放幾何和物理研究目錄下放參數(shù)掃描結果目錄下放導出數(shù)據(jù)Matlab腳本單獨放一個文件夾。如果你后續(xù)要換體系、換晶格類型只需要改參數(shù)和幾何研究步驟幾乎不用動。2.2 幾何搭建和參數(shù)化用晶格矢量而不是普通直角坐標在Comsol中創(chuàng)建一個二維模型物理接口選擇“電磁波頻域”Electromagnetic Waves, Frequency Domain。幾何上我們建一個平行四邊形單胞而不是正六邊形單胞。原因是平行四邊形單胞的兩條邊正好對應a1、a2兩個平移矢量Floquet周期邊界條件的兩個相對邊就沿這兩個方向施加周期關系一目了然。具體參數(shù)可以這樣設置全局參數(shù)a 800 nm作為晶格常數(shù)A柱半徑rA 0.2aB柱半徑rB 0.1a對稱時令兩者相等背景介電常數(shù)eps_bg 1柱介電常數(shù)eps_rod 12這是硅在近紅外的典型值平行四邊形單胞的兩個相鄰邊向量就是a1、a2繪制時將左下角定在原點。介質柱在單胞內(nèi)的位置需要嚴格對應子格坐標。我推薦用“幾何”里的“圓”來添加兩個圓域圓心分別放在(0,0)和(0, a/√3)。注意不要用普通坐標要用參數(shù)表達這樣以后改a時整個模型自動更新。材料域分配上背景域設為空氣圓形域設為高介電材料。為了后面網(wǎng)格收斂性驗證我還會單獨加一個參數(shù)mesh_size來控制網(wǎng)格最大單元尺寸。還有一個細節(jié)是單位問題。Comsol的幾何尺寸一般帶單位你用nm畫圖時特征頻率結果會是Hz級別的高頻數(shù)看起來很不直觀。我習慣把幾何單位設成μm微米頻率出來后換算到常用的歸一化頻率f a / c也就是在Matlab后處理里統(tǒng)一處理這樣能帶圖的縱軸才具有普適性。2.3 布洛赫邊界條件的關鍵周期性結對與網(wǎng)格一致性這一步是整個Comsol建模中最容易出問題的地方。在物理接口中添加“周期條件”類型選擇“Floquet周期”然后為兩對相對邊界分別指定波矢分量。Comsol中通常需要輸入兩個方向的波數(shù)kx、ky我建議先在全局參數(shù)中定義kx k1 * (2pi/(sqrt(3)a)) k2 * (2pi/(sqrt(3)a))ky k1 * (2pi/a) - k2 * (2pi/a)這里的k1、k2是倒格矢約化坐標后續(xù)k路徑掃描會以它們?yōu)閰?shù)。這里有一個到k空間直角坐標的換算必須和你的基矢定義嚴格對應。最穩(wěn)妥的做法是用約化坐標(k1, k2)作為掃描變量kx、ky都通過參數(shù)表達式實時計算這樣最后導出的數(shù)據(jù)和Matlab腳本里的k網(wǎng)格也是一一對應的。周期性邊界條件還有一個坑相對邊界上的網(wǎng)格必須一一配對否則Comsol會在邊界插值上引入誤差產(chǎn)生大量偽模。你需要進入網(wǎng)格步驟打開“周期性結對”Periodic Pair功能讓左右邊界、上下邊界的網(wǎng)格節(jié)點精確對應。我自己的經(jīng)驗是如果忘了配對特征頻率會多出很多能帶曲線看起來就像一團毛線完全沒有規(guī)律而一旦配對正確哪怕網(wǎng)格密度一般低頻段的能帶也會非常清爽。另外在Floquet邊界條件下Comsol內(nèi)部的場變量自動帶有布洛赫相位exp(ik·r)的約定。這一點對后面的Matlab相位修正特別重要——我建議在建模完成后先做一個自檢算一個均勻背景沒有介質柱的平行平板模式把解析色散關系和Comsol結果對比確認k空間相位約定到底差一個正負號。這個自檢兩小時就能做完能省掉后面調試陳數(shù)的三天時間。2.4 求解器設置特征頻率研究、搜索基準、所需模態(tài)數(shù)研究類型選擇“特征頻率”Eigenfrequency。在研究設置里需要指定“所需特征頻率數(shù)”我一般取8到12個覆蓋帶隙附近主要頻段即可。搜索基準值也要設置如果目標帶隙在歸一化頻率0.3附近可以根據(jù)公式f (0.3*c/a)換算成Hz填入。這個基準值不需要非常精確但能極大提升特征值求解器的收斂速度也能避免求解器抓住一堆高頻模式。參數(shù)掃描要在“研究擴展”里使用“輔助掃描”Auxiliary sweep而不是簡單地在研究里嵌套掃描。輔助掃描可以保證每個參數(shù)組合都執(zhí)行完整的特征值求解步驟而且結果數(shù)據(jù)會用一個額外的維度標記每個掃描點方便后面一維繪圖。我們把s設為掃描變量范圍0到3步長0.02左右然后用參數(shù)表達式把s映射到k1、k2上。求解器設置里還有一個容易忽略的點特征值求解器默認使用SPOOLES或MUMPS對于二維小模型都夠用。如果模式數(shù)多、網(wǎng)格細建議切換成MUMPS并打開“稀疏直接求解器”的“行預排序”選項速度通常能改善不少。對蜂窩晶格這個體系模型規(guī)模其實不大Intel八代以上處理器一般一兩分鐘就能掃完一個k點整條路徑掃描在半小時到一小時之間。3. 能帶路徑掃描與數(shù)據(jù)導出把k空間“走”成一條連續(xù)的曲線3.1 波矢路徑和參數(shù)映射表Γ-K-M-Γ的換算蜂窩晶格的倒空間高對稱路徑通常取?！鶮→M→Γ路徑在倒格矢坐標下可以用一個統(tǒng)一的路徑參數(shù)s來定義。s從0到3每段對應一個高對稱段。我建議把k1、k2寫成s的分段表達式s范圍k1k20到1s/3s/31到21/3 (s?1)/61/3 ? (s?1)/32到31/2 ? (s?2)/20比如s0代表Γ點s1代表K點s2代表M點s3回到Γ點。把這些表達式直接寫進Comsol的全局參數(shù)定義里k1、k2就變成s的函數(shù)了。此時也可以定義kx、ky的表達式如前文所述。需要注意的是這里的s的物理含義是無量綱路徑比例而不是k的大小所以能帶圖的橫軸直接用s或者累計路徑長度都行。我個人習慣在Matlab里把s映射成實際倒空間距離這樣橫軸單位是1/a看起來更專業(yè)。3.2 在Comsol中設置輔助掃描的完整步驟我在Comsol界面的實際操作流程如下在“全局參數(shù)”中定義a、rA、rB、eps等幾何和材料參數(shù)再定義s、k1、k2、kx、ky。在“研究1”中選擇特征頻率研究在“研究設置”里填所需特征頻率數(shù)。展開“研究擴展”勾選“輔助掃描”在掃描參數(shù)里選擇s填入范圍0到3、步長0.02。確?！霸谒袇?shù)組合上求解”被勾選。運行研究Comsol會自動對每個s做一次特征值求解。這里有個經(jīng)驗數(shù)值步長0.02對應的路徑段上有大約150個k點已經(jīng)足夠畫出一條光滑的能帶曲線。如果你后面要做精確的Wilson loop建議路徑掃描的s步長再細化到0.005不過這會讓總求解時間拉長一倍以上前期調模型時保持0.02就好。3.3 導出頻率和場數(shù)據(jù)的推薦方式數(shù)據(jù)導出有兩種主流路線取決于你手頭有沒有LiveLink for MATLAB。如果你有LiveLink那最理想的做法是直接在Matlab里通過COMSOL Java API控制模型掃描s后把頻率和本征場一次性讀入工作區(qū)所有數(shù)據(jù)都不落盤。這個方案最適合反復迭代缺點是前期需要寫一段比較長的調用腳本。我自己的項目用的是這個方案后面給Matlab腳本時也會給出對應的數(shù)據(jù)接口設計。如果你沒有LiveLink也不用慌Comsol GUI導出的數(shù)據(jù)足夠完成陳數(shù)計算。需要導出兩個東西頻率數(shù)據(jù)通過“派生值 → 全局計算”得到每個s下的所有特征頻率把表格導出為文本文件每行是“s, f1, f2, …”。場數(shù)據(jù)在“數(shù)據(jù)集”中選擇某個參數(shù)點的主特征模式然后“導出 → 解”中選擇你關心的域勾選坐標列以及Ez的實部、虛部導出為CSV文件。文件格式大致是x, y, Re(Ez), Im(Ez)0.000, 0.000, 2.34e-1, -1.22e-3...要特別提醒的是Comsol導出本征場時實部和虛部通常是分離的而且歸一化到能量或某種絕對尺度。這意味著不同k點導出的場幅度不一定可比但對Wilson loop重疊積分來說每個k點內(nèi)部的歸一化已經(jīng)足夠——重疊積分的值會在適當?shù)貧w一化后變?yōu)橛暇仃囋?。我建議在Matlab腳本中顯式地對每條能帶的場向量做一次單位化避免幅度不匹配導致的重疊矩陣退化。3.4 網(wǎng)格收斂性校驗一個系數(shù)就能判斷結果可不可信能帶計算和有限元的所有計算一樣必須先確認結果不受網(wǎng)格影響。蜂窩晶格光子晶體有一個天然的探針對稱狀態(tài)下K點的狄拉克頻率。它必須和網(wǎng)格密度無關地收斂到某個值。我的做法是分別用最大單元尺寸為a/10、a/20、a/30三套網(wǎng)格各算一次對稱狀態(tài)下的K點頻率和帶隙關閉程度。如果兩次加密之間狄拉克頻率變化小于0.5%就認為網(wǎng)格足夠。另一個更苛刻的驗證是rA≠rB時帶隙寬度的收斂。網(wǎng)格太粗通常會高估帶隙因為數(shù)值色散會人為加大模式分裂網(wǎng)格加細之后帶隙會單調下降到一個平臺。另外圓柱邊界處的網(wǎng)格需要單獨加密。圓邊界如果只用粗網(wǎng)格相當于用一個多邊形去近似圓弧這等于在幾何層面額外引入了對稱性破缺。我一般會在“網(wǎng)格”里添加一個“邊界層”沿圓柱邊界設置5到8層邊界層網(wǎng)格厚度因子設為0.2左右。這一步對狄拉克錐頻率的精度影響非常大值得每次建模都加上。4. 陳數(shù)計算的Matlab腳本詳解從本征場到Berry曲率積分4.1 核心公式和數(shù)值算法為什么用plaquette Wilson loop陳數(shù)的定義是Berry曲率在布里淵區(qū)上的積分C (1/2π) ∫ F(kx, ky) d2kBerry曲率F由Berry聯(lián)絡A i?u|?_k u?的旋度給出。數(shù)值實現(xiàn)時我們不會去顯式求導而是使用Wilson loop方法把k空間劃成很多小方格plaquette每個小格子的四個角點各自對應一個本征場u(k)。然后構造相鄰k點之間的重疊矩陣M(k, k′) ?u(k)|u(k′)?由四個角點圍成閉環(huán)的Wilson loop相位就是該小格上的Berry曲率U_total M(k1→k2) M(k2→k3) M(k3→k4) M(k4→k1)F ≈ arg(det(U_total))把所有小格子的相位加起來除以2π就得到了陳數(shù)。這個算法的好處是回避了本征態(tài)相位規(guī)范選擇的問題因為每條能帶的全局相位在共軛內(nèi)積中抵消了。4.2 Matlab讀取Comsol數(shù)據(jù)的格式約定我建議建立一個統(tǒng)一的數(shù)據(jù)目錄結構comsol_output/ klist.txt % 每一行: kx, ky, s freq_list.txt % 每一行: s, f1, f2, f3... fields/ k001.txt % 每個k點的場文件 k002.txt ...場文件的格式為x, y, Re(Ez_band1), Im(Ez_band1), Re(Ez_band2), Im(Ez_band2), ... 0.0, 0.0, 2.3e-2, -1.1e-3, ...Matlab讀取腳本里我建議固定按行讀入并reshape成網(wǎng)格矩陣然后用實際網(wǎng)格坐標做數(shù)值積分。為了讓Wilson loop計算省事我推薦在Comsol導出場時使用統(tǒng)一的矩形網(wǎng)格坐標這樣不同k點之間的場可以直接配對。如果拿到的場分布是三角形網(wǎng)格上的散點就要先在Matlab里用griddata插值到統(tǒng)一網(wǎng)格。這一步我強烈建議在Comsol導出時就做好在導出節(jié)點中多花幾分鐘指定網(wǎng)格坐標比之后在Matlab里插值省事得多也避免插值誤差積累進Wilson loop。4.3 Wilson loop矩陣構造和相位追蹤下面是核心的Matlab代碼片段。我先給出主循環(huán)部分然后單獨說明兩個重要細節(jié)。% 載入k列表和場數(shù)據(jù) [kx, ky] load_klist(comsol_output/klist.txt); fields load_fields(comsol_output/fields/, length(kx)); bands 1:2; % 占據(jù)帶索引例如帶1~2 Nx 30; Ny 30; % k網(wǎng)格數(shù)實際使用需根據(jù)數(shù)據(jù)覆蓋調整 C 0; for i 1:Nx for j 1:Ny % 四個角點的索引 idxA sub2ind([Nx, Ny], i, j); idxB sub2ind([Nx, Ny], i1, j); idxC sub2ind([Nx, Ny], i1, j1); idxD sub2ind([Nx, Ny], i, j1); % 構造各邊重疊矩陣 Uab wloop_matrix(fields(idxA), fields(idxB), bands); Ubc wloop_matrix(fields(idxB), fields(idxC), bands); Ucd wloop_matrix(fields(idxC), fields(idxD), bands); Uda wloop_matrix(fields(idxD), fields(idxA), bands); % 閉合環(huán)路相位 Uloop Uab * Ubc * Ucd * Uda; F angle(det(Uloop)); C C F; end end C C / (2*pi);wloop_matrix函數(shù)的實現(xiàn)如下function U wloop_matrix(f1, f2, bands) nb length(bands); M zeros(nb, nb); for n 1:nb for m 1:nb u1 f1.u{bands(n)}; u2 f2.u{bands(m)}; M(n, m) sum(conj(u1(:)) .* u2(:)); end end % M是方陣直接作為Wilson線算符 % 必要時可以做逆矩陣修正M / (M*M)^0.5 U M; end如果計算中curly的M矩陣不等于酉矩陣需要做極分解修正。我自己的經(jīng)驗是當場數(shù)據(jù)歸一化且平面波近似時M矩陣都會很接近酉矩陣但嚴格起見可以加上奇異值分解修正[Uq, ~, Vq] svd(M); U Uq * Vq;這一步能有效消除數(shù)值噪聲特別是網(wǎng)格較粗時。4.4 最終陳數(shù)計算與結果驗證直接算出來的C未必是漂亮整數(shù)。比如你可能得到0.934或者1.06這是k網(wǎng)格離散誤差的正常表現(xiàn)加密網(wǎng)格會逼近整數(shù)。我自己通常畫一條“網(wǎng)格加密曲線”把k網(wǎng)格從10×10逐步增加到50×50觀察C怎么收斂。如果C穩(wěn)定收斂到1或者?1就可以放心如果C在某個網(wǎng)格密度下反復跳變大概率是能帶交叉或者相位修正出了問題。另外一個重要的驗證手段是和已知對稱性對照在rArB的極限帶隙關閉陳數(shù)失去定義在rA和rB的關系相反時陳數(shù)應當變號。用Matlab跑一遍參數(shù)掃描畫出“rA/rB ? C”曲線如果看到清晰的平臺區(qū)比如比值大于1的地帶一直保持C1小于1的地帶C?1這個結果就和物理圖像完全吻合。4.5 腳本中兩個容易寫錯的細節(jié)第一個細節(jié)是布洛赫相位修正。Comsol的Floquet邊界條件解出來的物理場是完整的布洛赫函數(shù)E(r) exp(ik·r)u(r)但Wilson loop中需要的是周期部分u(r)。所以從Comsol導出Ez后必須做一次換算u E·exp(?ik·r)。這里的正負號以上面提到的自檢結果為準。如果在Matlab里把這個符號搞反重疊矩陣會多出一個隨位移指數(shù)變化的相位因子直接導致Wilson loop相位不閉合陳數(shù)算出個非整數(shù)。第二個細節(jié)是能帶排序。參數(shù)掃描時Comsol輸出的特征頻率是按大小排列的但相鄰k點之間的模式身份可能互換尤其在高對稱點附近。做Wilson loop時必須確?！暗趎條帶”在相鄰k點上是同一條物理能帶否則重疊矩陣的行列會錯位。解決的辦法就是前面提到的用上一個k點的模式與當前k點所有模式做重疊選擇重疊最大的索引作為當前能帶順序。我通常在load_fields之后就完成按順序重排再送進Wilson loop避免在循環(huán)里反復處理。5. 踩坑實錄假帶隙、邊界模式混疊和陳數(shù)符號問題5.1 狄拉克錐沒有打開先查網(wǎng)格對稱和邊界條件如果你設置了rArB理論預期K點應該是狄拉克錐可是算出來的能帶在K點附近有明顯分裂這就是一個典型警告信號幾何或者網(wǎng)格對稱性被破壞了。第一個可查的項目是單胞的幾何坐標確認A柱圓心在原點、B柱圓心在準確的(0, a/√3)不要有小數(shù)點誤差。第二個可查的項目是網(wǎng)格剖分。我遇到過最隱蔽的一次問題出在周期性結對邊界網(wǎng)格雖然節(jié)點數(shù)一致但兩對邊界的“源/目標”配對方向反了導致周期性條件施加的是反對稱關系狄拉克錐被強行打開成一條假帶隙。排查方法很簡單——把網(wǎng)格隱藏掉單獨看能帶如果帶隙依然存在就關掉周期條件改做鏡像模式或者查看邊界兩側的場確認相位關系。這種問題純靠看頻率數(shù)值很難發(fā)現(xiàn)關鍵時刻還是要看場分布。5.2 布里淵區(qū)路徑端點不連續(xù)模式互換與能帶穿越能帶圖上最讓人頭疼的現(xiàn)象是某條能帶在跨過高對稱點后“跳”到了另一條顏色對應的曲線上。這通常不是物理問題而是特征頻率排序造成的。在K點附近群速度為零模式簡并度高微小數(shù)值誤差會讓求解器把模式順序打亂。處理思路有兩條。第一條是把s步長加密比如從0.02改成0.005讓每個模式在相鄰兩點間的頻率差足夠小排序更穩(wěn)定。第二條是在后處理時按重疊積分重排方法在4.5節(jié)里已經(jīng)說了。如果項目只關心帶隙位置和陳數(shù)其實能帶圖橫軸上的這種跳躍不影響最終結果但如果你要畫漂亮的能帶圖給論文用就必須把重排邏輯加進去。5.3 陳數(shù)算出0.3或?0.7這不是非整數(shù)是網(wǎng)格和相位問題當你第一次跑通整個腳本算出來的陳數(shù)往往是0.31、?0.72這種不尷不尬的值。有時還會出現(xiàn)正負號漂移。出現(xiàn)這個問題的原因幾乎總是以下三個之一布洛赫相位修正的正負號錯了k網(wǎng)格太粗Berry曲率集中在K/K′尖點附近沒有被充分采樣場數(shù)據(jù)沒有統(tǒng)一網(wǎng)格導致重疊積分誤差大。我的建議是先花半天時間用緊束縛模型驗證Matlab腳本本身沒有Bug。寫一個2×2或者4×4的緊束縛哈密頓量把本征向量按解析公式算好直接喂給同一個Wilson loop函數(shù)如果算出來結果與解析陳數(shù)一致那腳本就沒問題剩下的鍋全甩給Comsol導出和相位修正。這個方法幫我節(jié)省了大量排查時間。5.4 與其他結果交叉驗證邊緣態(tài)和群速度方向即使陳數(shù)算出了漂亮的整數(shù)我也建議做一步額外的驗證計算超胞的邊緣態(tài)能帶。把單胞改成條帶結構沿一個方向有限另一個方向仍用周期邊界在同一個Comsol模型里再算一次特征頻率。如果體態(tài)帶隙內(nèi)確實出現(xiàn)了跨過整個帶隙的邊緣態(tài)色散曲線而且邊緣態(tài)在同一邊界的兩個方向的群速度方向相反對應不同手性那陳數(shù)結果基本可以蓋章確認。這一步還有個額外好處可以直接觀察邊緣態(tài)場分布來判斷拓撲保護是否成立。比如在帶隙頻率范圍內(nèi)激勵一個源看波是否能繞過缺陷轉角繼續(xù)單向傳播。我在實際項目中就是靠這個最終驗證了設計的拓撲波導整個過程比單純相信一個整數(shù)踏實得多。最后分享一點個人的體會能帶計算和陳數(shù)計算雖然是兩套東西但它們的錯誤模式高度耦合——能帶有一點點不干凈陳數(shù)就會顯現(xiàn)成莫名其妙的數(shù)值。所以做這個項目時不要急著沖到最后一步先把“對稱狀態(tài)下狄拉克錐正確關閉”這個最基本的結果反復算穩(wěn)、看穩(wěn)再談拓撲不變量。把這套流程走通之后以后再換三角晶格、Kagome晶格或者其他復雜體系框架完全可以直接復用只是換參數(shù)、換波矢路徑而已。 蜂窩晶格光子晶體：從Comsol能帶到Matlab陳數(shù)計算全流程
1. 蜂窩晶格為什么值得算狄拉克錐與拓撲光子學的切入點做拓撲光子學研究的人遲早都會碰到蜂窩晶格光子晶體。這個體系幾乎就是為演示“如何從能帶計算走向拓撲不變量”而生的。它的晶格結構和石墨烯一模一樣由兩套三角子格穿插構成所以天然具備狄拉克錐而只要破壞了空間反演對稱性狄拉克點處就會打開一個帶隙進而可以定義陳數(shù)、研究邊緣態(tài)。很多初學者以為能帶圖算出來就萬事大吉了實際上從Comsol模型到Matlab里的陳數(shù)結果中間隔著不少容易踩坑的環(huán)節(jié)。我在這篇文章里會把這套流程完整地拆開講蜂窩晶格單胞怎么建、布洛赫周期邊界怎么設、k路徑怎么掃、本征場怎么導出來最后又怎么用Matlab腳本從實空間場分布出發(fā)算出Berry曲率和陳數(shù)。內(nèi)容適合有有限元基礎、想做拓撲光子學數(shù)值驗證的研究生也適合想快速搭一套“能帶計算拓撲不變量”工作流的工程師。里面所有參數(shù)和腳本邏輯都是我在實際項目中反復驗證過的方案。1.1 從能帶折疊到狄拉克錐蜂窩結構的關鍵是兩套格子蜂窩晶格在實空間里的形象是一張六角網(wǎng)但計算時并不推薦直接畫一個六邊形單胞而是把它看成兩套三角子格。設晶格常數(shù)為a兩個基矢取為a1 (√3/2·a, 1/2·a)a2 (√3/2·a, ?1/2·a)A子格放在原點B子格放在最近鄰位置比如(0, a/√3)。兩套子格在幾何上完全等價時系統(tǒng)的空間反演對稱性沒有被破壞倒易空間第一布里淵區(qū)的K和K′點會出現(xiàn)兩重簡并能帶形成線性交叉也就是狄拉克錐。對應到光子晶體里如果A、B位置放的介質柱半徑相同、介電常數(shù)相同那么TE或TM模式的能帶在K點附近一定會呈錐形。這個結構在Comsol中非常好驗證把A、B柱半徑設成相等掃出一條能帶曲線你會在帶隙關閉的狀態(tài)下看到兩條帶恰好接觸。這里有一個很容易忽略的前提——幾何對稱性必須被數(shù)值網(wǎng)格忠實地保留。如果網(wǎng)格剖分破壞了C3旋轉對稱性K點的簡并度會被人為地打開狄拉克錐就會變成一條假帶隙拓撲結論自然全部作廢。所以后面講網(wǎng)格剖分的時候我專門花了一節(jié)來說這個事。1.2 對稱性破缺決定帶隙能否打開和陳數(shù)直接相關一旦A、B兩套格子的參數(shù)不再相同比如讓A柱半徑rA大于B柱半徑rB空間反演對稱性就被破壞。K點的兩重簡并上升原先線性交叉的兩條能帶被分開形成一段帶隙。這段帶隙的重要性在于它給了我們一個干凈的“頻率窗口”可以定義某條占據(jù)帶通常取帶隙下方的一組帶的拓撲性質。這里需要特別說明光子晶體中“陳數(shù)”的定義方式。光子本身是玻色子不存在電子那樣費米面和占據(jù)態(tài)但我們可以把布洛赫模式當成一個參數(shù)空間中的向量場對某個頻帶的態(tài)進行投影構造出等價于量子Hilbert空間的結構。具體做法是只看帶隙以下那條帶或者那組互相簡并的帶把它當成“占據(jù)帶”在布里淵區(qū)內(nèi)對Berry曲率積分得到的整數(shù)就是陳數(shù)。這個整數(shù)直接和體邊對應關系掛鉤帶隙內(nèi)存在單向傳輸?shù)倪吘墤B(tài)邊緣態(tài)的數(shù)量和方向就由陳數(shù)決定。所以整個項目的邏輯鏈其實很簡潔蜂窩晶格給出狄拉克錐對稱性破缺打開帶隙帶隙的非平庸性質由陳數(shù)描述。下一步Comsol負責給出能帶和本征場Matlab負責把場的相位信息提煉成拓撲判據(jù)。1.3 這個體系的“陳數(shù)”到底在描述什么很多人剛接觸陳數(shù)時總被抽象公式嚇到。我習慣用一個類比來解釋把布里淵區(qū)想象成一個封閉曲面每個k點處都有一個本征向量也就是布洛赫函數(shù)的周期部分。這些向量會隨著k變化而指向不同方向。陳數(shù)描述的就是這個向量場在封閉曲面上纏繞的“圈數(shù)”。就像莫比烏斯帶和普通紙帶的扭轉數(shù)量差一個整數(shù)沒法通過連續(xù)形變改變陳數(shù)也是一個無法連續(xù)變化的整數(shù)拓撲不變量。對于蜂窩晶格光子晶體陳數(shù)通常不止一個值。如果介質柱半徑rA略大于rB某幾條頻帶可能得到1或?1的陳數(shù)而rA小于rB時符號反轉。符號的意義直接映射到邊緣態(tài)的傳播方向——正陳數(shù)的邊界上邊緣態(tài)沿某一方向單向傳播負陳數(shù)則反過來。這個特性就是拓撲保護的根源即使路徑上有彎折、有缺陷只要帶隙沒有被填滿單向傳播就不會被背散射破壞。我在后面會用這個性質和邊緣態(tài)能帶做交叉驗證。2. Comsol建模與MPH工作流單胞怎么畫、邊界怎么設才能不出偽模2.1 先澄清“MPH”這個說法在Comsol生態(tài)里的實際含義很多人以為“MPH算法”是什么高深的學術名詞實際上在Comsol的語境中.mph就是模型文件的擴展名。一份.mph文件包含了整個項目的一切幾何、物理接口、網(wǎng)格、研究步驟、結果和后處理。所謂“含MPH算法”合理的解讀就是依托Comsol的.mph模型文件來組織整套計算流程的方法。我個人喜歡把“MPH”拆成三個動作Modeling建模與參數(shù)化、Processing求解與數(shù)據(jù)提取、Harvesting后處理與拓撲判據(jù)計算。這樣工程文件的結構會很清晰模型目錄下放幾何和物理研究目錄下放參數(shù)掃描結果目錄下放導出數(shù)據(jù)Matlab腳本單獨放一個文件夾。如果你后續(xù)要換體系、換晶格類型只需要改參數(shù)和幾何研究步驟幾乎不用動。2.2 幾何搭建和參數(shù)化用晶格矢量而不是普通直角坐標在Comsol中創(chuàng)建一個二維模型物理接口選擇“電磁波頻域”Electromagnetic Waves, Frequency Domain。幾何上我們建一個平行四邊形單胞而不是正六邊形單胞。原因是平行四邊形單胞的兩條邊正好對應a1、a2兩個平移矢量Floquet周期邊界條件的兩個相對邊就沿這兩個方向施加周期關系一目了然。具體參數(shù)可以這樣設置全局參數(shù)a 800 nm作為晶格常數(shù)A柱半徑rA 0.2aB柱半徑rB 0.1a對稱時令兩者相等背景介電常數(shù)eps_bg 1柱介電常數(shù)eps_rod 12這是硅在近紅外的典型值平行四邊形單胞的兩個相鄰邊向量就是a1、a2繪制時將左下角定在原點。介質柱在單胞內(nèi)的位置需要嚴格對應子格坐標。我推薦用“幾何”里的“圓”來添加兩個圓域圓心分別放在(0,0)和(0, a/√3)。注意不要用普通坐標要用參數(shù)表達這樣以后改a時整個模型自動更新。材料域分配上背景域設為空氣圓形域設為高介電材料。為了后面網(wǎng)格收斂性驗證我還會單獨加一個參數(shù)mesh_size來控制網(wǎng)格最大單元尺寸。還有一個細節(jié)是單位問題。Comsol的幾何尺寸一般帶單位你用nm畫圖時特征頻率結果會是Hz級別的高頻數(shù)看起來很不直觀。我習慣把幾何單位設成μm微米頻率出來后換算到常用的歸一化頻率f a / c也就是在Matlab后處理里統(tǒng)一處理這樣能帶圖的縱軸才具有普適性。2.3 布洛赫邊界條件的關鍵周期性結對與網(wǎng)格一致性這一步是整個Comsol建模中最容易出問題的地方。在物理接口中添加“周期條件”類型選擇“Floquet周期”然后為兩對相對邊界分別指定波矢分量。Comsol中通常需要輸入兩個方向的波數(shù)kx、ky我建議先在全局參數(shù)中定義kx k1 * (2pi/(sqrt(3)a)) k2 * (2pi/(sqrt(3)a))ky k1 * (2pi/a) - k2 * (2pi/a)這里的k1、k2是倒格矢約化坐標后續(xù)k路徑掃描會以它們?yōu)閰?shù)。這里有一個到k空間直角坐標的換算必須和你的基矢定義嚴格對應。最穩(wěn)妥的做法是用約化坐標(k1, k2)作為掃描變量kx、ky都通過參數(shù)表達式實時計算這樣最后導出的數(shù)據(jù)和Matlab腳本里的k網(wǎng)格也是一一對應的。周期性邊界條件還有一個坑相對邊界上的網(wǎng)格必須一一配對否則Comsol會在邊界插值上引入誤差產(chǎn)生大量偽模。你需要進入網(wǎng)格步驟打開“周期性結對”Periodic Pair功能讓左右邊界、上下邊界的網(wǎng)格節(jié)點精確對應。我自己的經(jīng)驗是如果忘了配對特征頻率會多出很多能帶曲線看起來就像一團毛線完全沒有規(guī)律而一旦配對正確哪怕網(wǎng)格密度一般低頻段的能帶也會非常清爽。另外在Floquet邊界條件下Comsol內(nèi)部的場變量自動帶有布洛赫相位exp(ik·r)的約定。這一點對后面的Matlab相位修正特別重要——我建議在建模完成后先做一個自檢算一個均勻背景沒有介質柱的平行平板模式把解析色散關系和Comsol結果對比確認k空間相位約定到底差一個正負號。這個自檢兩小時就能做完能省掉后面調試陳數(shù)的三天時間。2.4 求解器設置特征頻率研究、搜索基準、所需模態(tài)數(shù)研究類型選擇“特征頻率”Eigenfrequency。在研究設置里需要指定“所需特征頻率數(shù)”我一般取8到12個覆蓋帶隙附近主要頻段即可。搜索基準值也要設置如果目標帶隙在歸一化頻率0.3附近可以根據(jù)公式f (0.3*c/a)換算成Hz填入。這個基準值不需要非常精確但能極大提升特征值求解器的收斂速度也能避免求解器抓住一堆高頻模式。參數(shù)掃描要在“研究擴展”里使用“輔助掃描”Auxiliary sweep而不是簡單地在研究里嵌套掃描。輔助掃描可以保證每個參數(shù)組合都執(zhí)行完整的特征值求解步驟而且結果數(shù)據(jù)會用一個額外的維度標記每個掃描點方便后面一維繪圖。我們把s設為掃描變量范圍0到3步長0.02左右然后用參數(shù)表達式把s映射到k1、k2上。求解器設置里還有一個容易忽略的點特征值求解器默認使用SPOOLES或MUMPS對于二維小模型都夠用。如果模式數(shù)多、網(wǎng)格細建議切換成MUMPS并打開“稀疏直接求解器”的“行預排序”選項速度通常能改善不少。對蜂窩晶格這個體系模型規(guī)模其實不大Intel八代以上處理器一般一兩分鐘就能掃完一個k點整條路徑掃描在半小時到一小時之間。3. 能帶路徑掃描與數(shù)據(jù)導出把k空間“走”成一條連續(xù)的曲線3.1 波矢路徑和參數(shù)映射表Γ-K-M-Γ的換算蜂窩晶格的倒空間高對稱路徑通常取?！鶮→M→Γ路徑在倒格矢坐標下可以用一個統(tǒng)一的路徑參數(shù)s來定義。s從0到3每段對應一個高對稱段。我建議把k1、k2寫成s的分段表達式s范圍k1k20到1s/3s/31到21/3 (s?1)/61/3 ? (s?1)/32到31/2 ? (s?2)/20比如s0代表Γ點s1代表K點s2代表M點s3回到Γ點。把這些表達式直接寫進Comsol的全局參數(shù)定義里k1、k2就變成s的函數(shù)了。此時也可以定義kx、ky的表達式如前文所述。需要注意的是這里的s的物理含義是無量綱路徑比例而不是k的大小所以能帶圖的橫軸直接用s或者累計路徑長度都行。我個人習慣在Matlab里把s映射成實際倒空間距離這樣橫軸單位是1/a看起來更專業(yè)。3.2 在Comsol中設置輔助掃描的完整步驟我在Comsol界面的實際操作流程如下在“全局參數(shù)”中定義a、rA、rB、eps等幾何和材料參數(shù)再定義s、k1、k2、kx、ky。在“研究1”中選擇特征頻率研究在“研究設置”里填所需特征頻率數(shù)。展開“研究擴展”勾選“輔助掃描”在掃描參數(shù)里選擇s填入范圍0到3、步長0.02。確?！霸谒袇?shù)組合上求解”被勾選。運行研究Comsol會自動對每個s做一次特征值求解。這里有個經(jīng)驗數(shù)值步長0.02對應的路徑段上有大約150個k點已經(jīng)足夠畫出一條光滑的能帶曲線。如果你后面要做精確的Wilson loop建議路徑掃描的s步長再細化到0.005不過這會讓總求解時間拉長一倍以上前期調模型時保持0.02就好。3.3 導出頻率和場數(shù)據(jù)的推薦方式數(shù)據(jù)導出有兩種主流路線取決于你手頭有沒有LiveLink for MATLAB。如果你有LiveLink那最理想的做法是直接在Matlab里通過COMSOL Java API控制模型掃描s后把頻率和本征場一次性讀入工作區(qū)所有數(shù)據(jù)都不落盤。這個方案最適合反復迭代缺點是前期需要寫一段比較長的調用腳本。我自己的項目用的是這個方案后面給Matlab腳本時也會給出對應的數(shù)據(jù)接口設計。如果你沒有LiveLink也不用慌Comsol GUI導出的數(shù)據(jù)足夠完成陳數(shù)計算。需要導出兩個東西頻率數(shù)據(jù)通過“派生值 → 全局計算”得到每個s下的所有特征頻率把表格導出為文本文件每行是“s, f1, f2, …”。場數(shù)據(jù)在“數(shù)據(jù)集”中選擇某個參數(shù)點的主特征模式然后“導出 → 解”中選擇你關心的域勾選坐標列以及Ez的實部、虛部導出為CSV文件。文件格式大致是x, y, Re(Ez), Im(Ez)0.000, 0.000, 2.34e-1, -1.22e-3...要特別提醒的是Comsol導出本征場時實部和虛部通常是分離的而且歸一化到能量或某種絕對尺度。這意味著不同k點導出的場幅度不一定可比但對Wilson loop重疊積分來說每個k點內(nèi)部的歸一化已經(jīng)足夠——重疊積分的值會在適當?shù)貧w一化后變?yōu)橛暇仃囋?。我建議在Matlab腳本中顯式地對每條能帶的場向量做一次單位化避免幅度不匹配導致的重疊矩陣退化。3.4 網(wǎng)格收斂性校驗一個系數(shù)就能判斷結果可不可信能帶計算和有限元的所有計算一樣必須先確認結果不受網(wǎng)格影響。蜂窩晶格光子晶體有一個天然的探針對稱狀態(tài)下K點的狄拉克頻率。它必須和網(wǎng)格密度無關地收斂到某個值。我的做法是分別用最大單元尺寸為a/10、a/20、a/30三套網(wǎng)格各算一次對稱狀態(tài)下的K點頻率和帶隙關閉程度。如果兩次加密之間狄拉克頻率變化小于0.5%就認為網(wǎng)格足夠。另一個更苛刻的驗證是rA≠rB時帶隙寬度的收斂。網(wǎng)格太粗通常會高估帶隙因為數(shù)值色散會人為加大模式分裂網(wǎng)格加細之后帶隙會單調下降到一個平臺。另外圓柱邊界處的網(wǎng)格需要單獨加密。圓邊界如果只用粗網(wǎng)格相當于用一個多邊形去近似圓弧這等于在幾何層面額外引入了對稱性破缺。我一般會在“網(wǎng)格”里添加一個“邊界層”沿圓柱邊界設置5到8層邊界層網(wǎng)格厚度因子設為0.2左右。這一步對狄拉克錐頻率的精度影響非常大值得每次建模都加上。4. 陳數(shù)計算的Matlab腳本詳解從本征場到Berry曲率積分4.1 核心公式和數(shù)值算法為什么用plaquette Wilson loop陳數(shù)的定義是Berry曲率在布里淵區(qū)上的積分C (1/2π) ∫ F(kx, ky) d2kBerry曲率F由Berry聯(lián)絡A i?u|?_k u?的旋度給出。數(shù)值實現(xiàn)時我們不會去顯式求導而是使用Wilson loop方法把k空間劃成很多小方格plaquette每個小格子的四個角點各自對應一個本征場u(k)。然后構造相鄰k點之間的重疊矩陣M(k, k′) ?u(k)|u(k′)?由四個角點圍成閉環(huán)的Wilson loop相位就是該小格上的Berry曲率U_total M(k1→k2) M(k2→k3) M(k3→k4) M(k4→k1)F ≈ arg(det(U_total))把所有小格子的相位加起來除以2π就得到了陳數(shù)。這個算法的好處是回避了本征態(tài)相位規(guī)范選擇的問題因為每條能帶的全局相位在共軛內(nèi)積中抵消了。4.2 Matlab讀取Comsol數(shù)據(jù)的格式約定我建議建立一個統(tǒng)一的數(shù)據(jù)目錄結構comsol_output/ klist.txt % 每一行: kx, ky, s freq_list.txt % 每一行: s, f1, f2, f3... fields/ k001.txt % 每個k點的場文件 k002.txt ...場文件的格式為x, y, Re(Ez_band1), Im(Ez_band1), Re(Ez_band2), Im(Ez_band2), ... 0.0, 0.0, 2.3e-2, -1.1e-3, ...Matlab讀取腳本里我建議固定按行讀入并reshape成網(wǎng)格矩陣然后用實際網(wǎng)格坐標做數(shù)值積分。為了讓Wilson loop計算省事我推薦在Comsol導出場時使用統(tǒng)一的矩形網(wǎng)格坐標這樣不同k點之間的場可以直接配對。如果拿到的場分布是三角形網(wǎng)格上的散點就要先在Matlab里用griddata插值到統(tǒng)一網(wǎng)格。這一步我強烈建議在Comsol導出時就做好在導出節(jié)點中多花幾分鐘指定網(wǎng)格坐標比之后在Matlab里插值省事得多也避免插值誤差積累進Wilson loop。4.3 Wilson loop矩陣構造和相位追蹤下面是核心的Matlab代碼片段。我先給出主循環(huán)部分然后單獨說明兩個重要細節(jié)。% 載入k列表和場數(shù)據(jù) [kx, ky] load_klist(comsol_output/klist.txt); fields load_fields(comsol_output/fields/, length(kx)); bands 1:2; % 占據(jù)帶索引例如帶1~2 Nx 30; Ny 30; % k網(wǎng)格數(shù)實際使用需根據(jù)數(shù)據(jù)覆蓋調整 C 0; for i 1:Nx for j 1:Ny % 四個角點的索引 idxA sub2ind([Nx, Ny], i, j); idxB sub2ind([Nx, Ny], i1, j); idxC sub2ind([Nx, Ny], i1, j1); idxD sub2ind([Nx, Ny], i, j1); % 構造各邊重疊矩陣 Uab wloop_matrix(fields(idxA), fields(idxB), bands); Ubc wloop_matrix(fields(idxB), fields(idxC), bands); Ucd wloop_matrix(fields(idxC), fields(idxD), bands); Uda wloop_matrix(fields(idxD), fields(idxA), bands); % 閉合環(huán)路相位 Uloop Uab * Ubc * Ucd * Uda; F angle(det(Uloop)); C C F; end end C C / (2*pi);wloop_matrix函數(shù)的實現(xiàn)如下function U wloop_matrix(f1, f2, bands) nb length(bands); M zeros(nb, nb); for n 1:nb for m 1:nb u1 f1.u{bands(n)}; u2 f2.u{bands(m)}; M(n, m) sum(conj(u1(:)) .* u2(:)); end end % M是方陣直接作為Wilson線算符 % 必要時可以做逆矩陣修正M / (M*M)^0.5 U M; end如果計算中curly的M矩陣不等于酉矩陣需要做極分解修正。我自己的經(jīng)驗是當場數(shù)據(jù)歸一化且平面波近似時M矩陣都會很接近酉矩陣但嚴格起見可以加上奇異值分解修正[Uq, ~, Vq] svd(M); U Uq * Vq;這一步能有效消除數(shù)值噪聲特別是網(wǎng)格較粗時。4.4 最終陳數(shù)計算與結果驗證直接算出來的C未必是漂亮整數(shù)。比如你可能得到0.934或者1.06這是k網(wǎng)格離散誤差的正常表現(xiàn)加密網(wǎng)格會逼近整數(shù)。我自己通常畫一條“網(wǎng)格加密曲線”把k網(wǎng)格從10×10逐步增加到50×50觀察C怎么收斂。如果C穩(wěn)定收斂到1或者?1就可以放心如果C在某個網(wǎng)格密度下反復跳變大概率是能帶交叉或者相位修正出了問題。另外一個重要的驗證手段是和已知對稱性對照在rArB的極限帶隙關閉陳數(shù)失去定義在rA和rB的關系相反時陳數(shù)應當變號。用Matlab跑一遍參數(shù)掃描畫出“rA/rB ? C”曲線如果看到清晰的平臺區(qū)比如比值大于1的地帶一直保持C1小于1的地帶C?1這個結果就和物理圖像完全吻合。4.5 腳本中兩個容易寫錯的細節(jié)第一個細節(jié)是布洛赫相位修正。Comsol的Floquet邊界條件解出來的物理場是完整的布洛赫函數(shù)E(r) exp(ik·r)u(r)但Wilson loop中需要的是周期部分u(r)。所以從Comsol導出Ez后必須做一次換算u E·exp(?ik·r)。這里的正負號以上面提到的自檢結果為準。如果在Matlab里把這個符號搞反重疊矩陣會多出一個隨位移指數(shù)變化的相位因子直接導致Wilson loop相位不閉合陳數(shù)算出個非整數(shù)。第二個細節(jié)是能帶排序。參數(shù)掃描時Comsol輸出的特征頻率是按大小排列的但相鄰k點之間的模式身份可能互換尤其在高對稱點附近。做Wilson loop時必須確?！暗趎條帶”在相鄰k點上是同一條物理能帶否則重疊矩陣的行列會錯位。解決的辦法就是前面提到的用上一個k點的模式與當前k點所有模式做重疊選擇重疊最大的索引作為當前能帶順序。我通常在load_fields之后就完成按順序重排再送進Wilson loop避免在循環(huán)里反復處理。5. 踩坑實錄假帶隙、邊界模式混疊和陳數(shù)符號問題5.1 狄拉克錐沒有打開先查網(wǎng)格對稱和邊界條件如果你設置了rArB理論預期K點應該是狄拉克錐可是算出來的能帶在K點附近有明顯分裂這就是一個典型警告信號幾何或者網(wǎng)格對稱性被破壞了。第一個可查的項目是單胞的幾何坐標確認A柱圓心在原點、B柱圓心在準確的(0, a/√3)不要有小數(shù)點誤差。第二個可查的項目是網(wǎng)格剖分。我遇到過最隱蔽的一次問題出在周期性結對邊界網(wǎng)格雖然節(jié)點數(shù)一致但兩對邊界的“源/目標”配對方向反了導致周期性條件施加的是反對稱關系狄拉克錐被強行打開成一條假帶隙。排查方法很簡單——把網(wǎng)格隱藏掉單獨看能帶如果帶隙依然存在就關掉周期條件改做鏡像模式或者查看邊界兩側的場確認相位關系。這種問題純靠看頻率數(shù)值很難發(fā)現(xiàn)關鍵時刻還是要看場分布。5.2 布里淵區(qū)路徑端點不連續(xù)模式互換與能帶穿越能帶圖上最讓人頭疼的現(xiàn)象是某條能帶在跨過高對稱點后“跳”到了另一條顏色對應的曲線上。這通常不是物理問題而是特征頻率排序造成的。在K點附近群速度為零模式簡并度高微小數(shù)值誤差會讓求解器把模式順序打亂。處理思路有兩條。第一條是把s步長加密比如從0.02改成0.005讓每個模式在相鄰兩點間的頻率差足夠小排序更穩(wěn)定。第二條是在后處理時按重疊積分重排方法在4.5節(jié)里已經(jīng)說了。如果項目只關心帶隙位置和陳數(shù)其實能帶圖橫軸上的這種跳躍不影響最終結果但如果你要畫漂亮的能帶圖給論文用就必須把重排邏輯加進去。5.3 陳數(shù)算出0.3或?0.7這不是非整數(shù)是網(wǎng)格和相位問題當你第一次跑通整個腳本算出來的陳數(shù)往往是0.31、?0.72這種不尷不尬的值。有時還會出現(xiàn)正負號漂移。出現(xiàn)這個問題的原因幾乎總是以下三個之一布洛赫相位修正的正負號錯了k網(wǎng)格太粗Berry曲率集中在K/K′尖點附近沒有被充分采樣場數(shù)據(jù)沒有統(tǒng)一網(wǎng)格導致重疊積分誤差大。我的建議是先花半天時間用緊束縛模型驗證Matlab腳本本身沒有Bug。寫一個2×2或者4×4的緊束縛哈密頓量把本征向量按解析公式算好直接喂給同一個Wilson loop函數(shù)如果算出來結果與解析陳數(shù)一致那腳本就沒問題剩下的鍋全甩給Comsol導出和相位修正。這個方法幫我節(jié)省了大量排查時間。5.4 與其他結果交叉驗證邊緣態(tài)和群速度方向即使陳數(shù)算出了漂亮的整數(shù)我也建議做一步額外的驗證計算超胞的邊緣態(tài)能帶。把單胞改成條帶結構沿一個方向有限另一個方向仍用周期邊界在同一個Comsol模型里再算一次特征頻率。如果體態(tài)帶隙內(nèi)確實出現(xiàn)了跨過整個帶隙的邊緣態(tài)色散曲線而且邊緣態(tài)在同一邊界的兩個方向的群速度方向相反對應不同手性那陳數(shù)結果基本可以蓋章確認。這一步還有個額外好處可以直接觀察邊緣態(tài)場分布來判斷拓撲保護是否成立。比如在帶隙頻率范圍內(nèi)激勵一個源看波是否能繞過缺陷轉角繼續(xù)單向傳播。我在實際項目中就是靠這個最終驗證了設計的拓撲波導整個過程比單純相信一個整數(shù)踏實得多。最后分享一點個人的體會能帶計算和陳數(shù)計算雖然是兩套東西但它們的錯誤模式高度耦合——能帶有一點點不干凈陳數(shù)就會顯現(xiàn)成莫名其妙的數(shù)值。所以做這個項目時不要急著沖到最后一步先把“對稱狀態(tài)下狄拉克錐正確關閉”這個最基本的結果反復算穩(wěn)、看穩(wěn)再談拓撲不變量。把這套流程走通之后以后再換三角晶格、Kagome晶格或者其他復雜體系框架完全可以直接復用只是換參數(shù)、換波矢路徑而已。