計(jì)算一維光子晶體的Zak相位最近不少做拓?fù)涔庾訉W(xué)的朋友都在問(wèn)我同一個(gè)問(wèn)題Comsol里能算出漂亮的能帶圖可怎么把它變成具有物理意義的拓?fù)洳蛔兞勘疚木蛧@這個(gè)項(xiàng)目完整記錄我從Comsol建模到Matlab分析Zak相位的全過(guò)程包括原理、參數(shù)設(shè)置、代碼實(shí)現(xiàn)和踩坑記錄。這套流程適合正在入門拓?fù)涔庾訉W(xué)、或者需要復(fù)現(xiàn)一維光子晶體邊界態(tài)研究的人參考看完之后你可以直接遷移到自己的體系里。1. 項(xiàng)目思路拆解為什么選Comsol加Matlab這條路線1.1 算Zak相位到底在算什么一維光子晶體本質(zhì)上就是折射率沿某個(gè)方向周期調(diào)制形成布拉格散射和能帶帶隙。它和拓?fù)涑渡详P(guān)系是因?yàn)樵?014年前后研究者發(fā)現(xiàn)兩個(gè)帶隙結(jié)構(gòu)相同的一維光子晶體拼接在一起時(shí)界面處是否出現(xiàn)帶隙內(nèi)的局域態(tài)取決于兩條能帶在帶隙邊界的Zak相位差。如果帶隙兩側(cè)的Zak相位差為π界面處就一定會(huì)出現(xiàn)拓?fù)溥吔鐟B(tài)如果相位差為0則沒(méi)有。Zak相位是Berry相位在一維周期性體系里的特殊形式數(shù)學(xué)定義是布洛赫波函數(shù)中的周期部分在布里淵區(qū)內(nèi)沿閉合路徑積分得到的相位。它的核心價(jià)值在于雖然能帶本身只告訴我們頻率和波矢的關(guān)系但Zak相位額外攜帶了波函數(shù)在動(dòng)量空間演化過(guò)程中的幾何信息這才是判斷拓?fù)湫再|(zhì)的關(guān)鍵。對(duì)于一維光子晶體Zak相位在結(jié)構(gòu)具有中心對(duì)稱性時(shí)表現(xiàn)為量子化的0或π這也是為什么那么多工作都圍繞著一維體系的拓?fù)溥吔鐟B(tài)展開(kāi)。把Zak相位算出來(lái)你就能解釋實(shí)驗(yàn)里觀察到的反射相位異常、邊界態(tài)存在性、以及透射譜中的共振峰來(lái)源。1.2 Comsol和Matlab的分工邏輯Comsol在光子晶體計(jì)算里的優(yōu)勢(shì)非常明顯幾何隨意、材料參數(shù)隨便設(shè)、邊界條件豐富網(wǎng)格局部加密也方便。但它有一個(gè)死穴——它默認(rèn)給你能帶頻率和模式場(chǎng)卻不直接給拓?fù)洳蛔兞俊Ｄ銦o(wú)法在Comsol的界面里找到一個(gè)叫“Zak相位”的輸出選項(xiàng)Matlab則正好補(bǔ)上這一環(huán)。把Comsol算出的每個(gè)波矢處的本征模式場(chǎng)提取出來(lái)在Matlab里做歸一化、做相鄰k點(diǎn)內(nèi)積、累加相位就能得到Wilson loop進(jìn)而推出Zak相位。這個(gè)過(guò)程本質(zhì)上是一個(gè)數(shù)值線性代數(shù)問(wèn)題Matlab做這種事輕車熟路。很多人可能會(huì)問(wèn)能不能直接用傳輸矩陣法或者平面波展開(kāi)法在Matlab里把能帶和Zak相位一次性搞定當(dāng)然可以那種方法編寫(xiě)出來(lái)之后運(yùn)行也很快但可擴(kuò)展性差換個(gè)復(fù)雜幾何就得重寫(xiě)。Comsol建模的優(yōu)勢(shì)在于后期你如果想要加入增益損耗、各向異性材料、非線性介質(zhì)、或者把模型推廣到二維三維代碼框架不用大改只要調(diào)整Comsol里的物理場(chǎng)配置就行。這也是我把Comsol和Matlab結(jié)合起來(lái)作為一個(gè)完整項(xiàng)目來(lái)做的原因。2. Zak相位的原理與數(shù)值計(jì)算要點(diǎn)2.1 從Berry相位到Zak相位的完整定義在固體物理里一條能帶對(duì)應(yīng)的布洛赫函數(shù)可以寫(xiě)成psi(k,r) e^{ikr} u(k,r)其中u(k,r)是周期函數(shù)滿足u(k,ra)u(k,r)。Zak相位定義為θ_n i ∮_BZ ?u_n(k)|?_k u_n(k)? dk這里的積分區(qū)間是整個(gè)一維布里淵區(qū)從-k0到k0再閉合回來(lái)對(duì)于一維體系通常取[-π/a, π/a]。|u_n(k)?是第n條能帶在波矢k處的布洛赫周期函數(shù)。這個(gè)積分形式上和Berry相位完全一致只是維度變成了一維。一維的特殊性在于布里淵區(qū)是一個(gè)閉合環(huán)波函數(shù)在k空間繞一圈之后回到原點(diǎn)但由于周期規(guī)范的存在會(huì)多出一個(gè)相位因子這個(gè)累積相位正是Zak相位。注意Zak相位是規(guī)范依賴的——原點(diǎn)平移會(huì)導(dǎo)致Zak相位發(fā)生改變。實(shí)際操作中只有當(dāng)結(jié)構(gòu)本身具有中心對(duì)稱時(shí)Zak相位才會(huì)被限制為0或π兩個(gè)離散值這也是做一維拓?fù)涔庾訉W(xué)的工作都選擇中心對(duì)稱結(jié)構(gòu)的原因。2.2 Wilson loop數(shù)值計(jì)算Zak相位的標(biāo)準(zhǔn)做法解析計(jì)算Zak相位需要知道u(k,r)的完整表達(dá)式這對(duì)數(shù)值計(jì)算不友好。數(shù)值上最常用的方法是Wilson loop方法原理是把積分離散化變成相鄰k點(diǎn)之間內(nèi)積的連乘。把布里淵區(qū)分成N個(gè)點(diǎn)k_m -π/a m·Δk其中m0, 1, ..., N-1Δk 2π/(aN)。離散形式的Zak相位可以寫(xiě)為θ_n -Im Σ_m log?u_n(k_m)|u_n(k_{m1})?這里|u_n(k_{m1})?是從k_m出發(fā)的下一個(gè)波矢處的模式整個(gè)連乘再取log和虛部就能得到相位。為什么這個(gè)式子成立因?yàn)楫?dāng)Δk很小時(shí)?u_n(k_m)|u_n(k_{m1})?約等于1 Δk?u_n|?_k u_n?取對(duì)數(shù)后虛部正好對(duì)應(yīng)Berry聯(lián)絡(luò)的累積。實(shí)際計(jì)算中需要特別注意兩點(diǎn)。第一相鄰k點(diǎn)之間的內(nèi)積模應(yīng)當(dāng)非常接近1如果內(nèi)積模顯著小于1說(shuō)明網(wǎng)格不夠密離散化誤差太大需要加密k點(diǎn)采樣。第二每個(gè)k點(diǎn)的模式場(chǎng)在導(dǎo)出時(shí)往往帶有任意相位這會(huì)導(dǎo)致內(nèi)積結(jié)果不穩(wěn)定因此每次內(nèi)積前需要先歸一化或者用一個(gè)固定的參考相位把模式場(chǎng)校正到同一規(guī)范下。2.3 邊界閉合條件最容易出錯(cuò)的一步前面提到布里淵區(qū)是一個(gè)閉合環(huán)但實(shí)際采樣時(shí)你不可能同時(shí)把-π/a和π/a作為獨(dú)立的采樣點(diǎn)各算一遍因?yàn)檫@是同一個(gè)物理點(diǎn)。Wilson loop的正確做法是只采樣N個(gè)點(diǎn)比如從k_0-π/a到k_{N-1}-π/a (N-1)Δk然后利用周期規(guī)范把最后一個(gè)點(diǎn)與第一個(gè)點(diǎn)連接起來(lái)。周期規(guī)范的表達(dá)式是u_n(kG) e^{-iGr} u_n(k)在一維情況下G 2π/a所以當(dāng)你處理邊界項(xiàng)、計(jì)算?u(k_{N-1})|u(k_N)?時(shí)k_N對(duì)應(yīng)的模式場(chǎng)不是獨(dú)立計(jì)算出來(lái)的而是用起點(diǎn)處k_0的模式場(chǎng)乘以e^{-i2πx/a}得到的。這個(gè)邊界閉合項(xiàng)如果不加積分的路徑就是開(kāi)的算出來(lái)的相位會(huì)差一個(gè)邊界貢獻(xiàn)數(shù)值上可能既不是0也不是π導(dǎo)致錯(cuò)誤結(jié)論。這是我實(shí)測(cè)下來(lái)整個(gè)項(xiàng)目里最容易出問(wèn)題的地方。很多教程在代碼里沒(méi)有明確處理這一項(xiàng)結(jié)果算出來(lái)的Zak相位看起來(lái)亂跳但很難找出原因。建議你在寫(xiě)代碼時(shí)單獨(dú)輸出每一步的相位增量和累計(jì)相位檢查是否存在從接近π跳變到-π的情況這種跳變往往說(shuō)明k點(diǎn)采樣過(guò)粗或邊界項(xiàng)處理有誤。2.4 符號(hào)約定與結(jié)果檢驗(yàn)方法Zak相位的符號(hào)取決于兩個(gè)地方一是Berry聯(lián)絡(luò)的積分方向二是Wilson loop公式里取log前的正負(fù)號(hào)。不同文獻(xiàn)的約定可能不同如果你發(fā)現(xiàn)自己的結(jié)果和文獻(xiàn)上的相位符號(hào)相反不要急著改代碼先統(tǒng)一約定再判斷。我建議在做完整計(jì)算之前用一個(gè)結(jié)構(gòu)簡(jiǎn)單、結(jié)果已知的模型來(lái)驗(yàn)證代碼。比如用二分之一占空比的一維光子晶體每個(gè)周期內(nèi)兩種介質(zhì)各占一半寬度通過(guò)改變?cè)c(diǎn)位置觀察Zak相位的變化或者直接和傳輸矩陣法得到的反射相位對(duì)比。理論上當(dāng)結(jié)構(gòu)從中心對(duì)稱變?yōu)榉侵行膶?duì)稱時(shí)Zak相位會(huì)從量子化的0或π變成非量子化的任意值。拿這個(gè)性質(zhì)來(lái)驗(yàn)證代碼是最直接的。3. Comsol建模實(shí)操參數(shù)、幾何、邊界條件與掃描設(shè)置3.1 建模前想清楚你算的是哪種偏振一維光子晶體的模式分為TE和TM偏振。在Comsol二維模型中不同偏振對(duì)應(yīng)不同的場(chǎng)分量。我的項(xiàng)目里選擇電場(chǎng)沿面外方向也就是電場(chǎng)只有z分量E_z磁場(chǎng)在面內(nèi)。對(duì)于這種偏振波動(dòng)方程變成標(biāo)量形式處理起來(lái)最簡(jiǎn)潔。如果你算的是面內(nèi)偏振即電場(chǎng)在x-y平面內(nèi)那需要同時(shí)求解E_x和E_y兩個(gè)分量方程組更大Floquet邊界條件的設(shè)置也更繁瑣。對(duì)于一維周期性結(jié)構(gòu)來(lái)說(shuō)算標(biāo)量形式足夠說(shuō)明問(wèn)題。實(shí)際課題研究中如果涉及斜入射或者偏振變換再去考慮面內(nèi)偏振不遲。3.2 幾何建模與材料參數(shù)設(shè)置我用的是二維模型幾何是在x方向取一個(gè)晶格常數(shù)a的長(zhǎng)度y方向取一個(gè)很小的矩形高度h。比如a 1微米h 0.05微米。兩個(gè)介質(zhì)層并排放在一個(gè)周期內(nèi)材料采用無(wú)損介質(zhì)折射率分別為n1和n2。這里給出我實(shí)際用的一組參數(shù)參數(shù)取值說(shuō)明a1 μm晶格常數(shù)h0.05a二維模型的y方向高度d10.7a高折射率層厚度d20.3a低折射率層厚度n13.45高折射率材料n21.0低折射率材料模擬空氣間隙之所以選n13.45和n21.0是因?yàn)檎凵渎蕦?duì)比度足夠大帶隙會(huì)比較寬后續(xù)在帶隙邊界處討論Zak相位時(shí)特征更明顯。如果你用二氧化硅和空氣組合折射率對(duì)比只有1.45左右?guī)墩?jì)算時(shí)需要更密的能帶采樣才能分辨帶隙邊界對(duì)數(shù)值精度要求高。材料參數(shù)中要注意Comsol中相對(duì)介電常數(shù)是ε_(tái)r n^2如果材料有損耗還需要設(shè)置電導(dǎo)率或復(fù)介電常數(shù)。算Zak相位時(shí)建議先關(guān)掉損耗因?yàn)閾p耗會(huì)模糊能帶邊界也會(huì)讓本征模式場(chǎng)的相位提取變得困難。3.3 Floquet周期邊界條件與波矢掃描這是Comsol建模最核心的一步。選擇電磁波、頻域接口二維模型。左邊界和右邊界設(shè)置成Floquet周期邊界條件K矢量設(shè)為(kx, 0)。上下邊界設(shè)置成周期性邊界條件這樣y方向無(wú)限延伸等價(jià)于一維體系。在全局參數(shù)里定義一個(gè)kx變量然后在研究設(shè)置中的輔助掃描里讓它變化掃描范圍取[-π/a, π/a]步長(zhǎng)按需要的k點(diǎn)數(shù)量確定。我通常先掃41個(gè)點(diǎn)做驗(yàn)證加密到81個(gè)點(diǎn)得到最終數(shù)據(jù)。需要特別注意的是Comsol的特征頻率求解器中k矢量通過(guò)Floquet周期邊界條件輸入而kx是作為參數(shù)掃描的。這種情況下每個(gè)kx值都會(huì)觸發(fā)一次特征值求解得到的每個(gè)特征頻率都對(duì)應(yīng)一條能帶上的點(diǎn)。Comsol會(huì)把所有特征值結(jié)果按求解器默認(rèn)方式排列但這不保證是按能帶順序排列的后期需要用頻率排序或重疊積分來(lái)整理數(shù)據(jù)。3.4 特征頻率研究設(shè)置與網(wǎng)格無(wú)關(guān)性測(cè)試研究類型選擇特征頻率要設(shè)定想要的模態(tài)數(shù)。我的項(xiàng)目中掃描前10條能帶所以特征頻率搜索基數(shù)設(shè)為10或12留一些余量。物理場(chǎng)中電磁波頻域接口的特征值方程本質(zhì)上是關(guān)于頻率ω的特征值問(wèn)題解出來(lái)的特征頻率經(jīng)過(guò)后處理可以直接畫(huà)能帶。網(wǎng)格建議先用較粗的網(wǎng)格跑通流程再逐步細(xì)化。因?yàn)槟愕闹饕繕?biāo)不光是能帶頻率還要提取模式場(chǎng)網(wǎng)格密度直接影響Wilson loop計(jì)算內(nèi)積時(shí)場(chǎng)矢量的離散表示。我在這個(gè)項(xiàng)目中測(cè)試過(guò)最大網(wǎng)格尺寸從0.1a細(xì)化到0.02a前幾條能帶的頻率變化很小但Zak相位的計(jì)算結(jié)果在粗網(wǎng)格下會(huì)出現(xiàn)明顯偏差。最終我采用最大網(wǎng)格尺寸0.03a既保證頻率精度又不會(huì)讓導(dǎo)出文件太大。一定要關(guān)閉自適應(yīng)網(wǎng)格細(xì)化。特征頻率求解器中如果有自適應(yīng)網(wǎng)格每個(gè)kx掃描點(diǎn)的網(wǎng)格會(huì)變化導(dǎo)致相鄰k點(diǎn)模式場(chǎng)所在的節(jié)點(diǎn)位置不一致后續(xù)Matlab算內(nèi)積時(shí)會(huì)出現(xiàn)嚴(yán)重問(wèn)題。這個(gè)坑我踩過(guò)如果你在數(shù)據(jù)后處理時(shí)發(fā)現(xiàn)內(nèi)積模遠(yuǎn)小于1先檢查是不是網(wǎng)格隨參數(shù)掃描變化導(dǎo)致的。3.5 導(dǎo)出模式場(chǎng)數(shù)據(jù)的具體操作計(jì)算完成后我們要導(dǎo)出每個(gè)kx處的本征模式場(chǎng)。在Comsol的導(dǎo)出數(shù)據(jù)選項(xiàng)中選擇數(shù)據(jù)集為特征頻率解的某一索引表達(dá)式填E_z的實(shí)部和虛部再額外導(dǎo)出x、y坐標(biāo)。這里的關(guān)鍵問(wèn)題是Comsol導(dǎo)出的是物理場(chǎng)E_z也就是布洛赫函數(shù)ψ_k(r)本身而不是周期函數(shù)u_k(r)。由于ψ_k(r) e^{ikx}u_k(r)因此導(dǎo)出后需要在Matlab里乘以e^{-ikx}來(lái)得到u_k。另一種做法是在Comsol的派生值中直接定義表達(dá)式E_zexp(-ikx*x)導(dǎo)出的就是周期函數(shù)部分。推薦后者省去在Matlab里處理的麻煩。每個(gè)kx點(diǎn)導(dǎo)出一個(gè)單獨(dú)的文件文件名包含kx索引例如data_0001.txtdata_0002.txt。導(dǎo)出時(shí)選擇所有特征頻率索引這樣每個(gè)文件都包含多條帶的數(shù)據(jù)。實(shí)際文件里會(huì)包含大量行每條帶在一個(gè)時(shí)間步索引中給出需要按特征值索引分拆。建議導(dǎo)出前在Comsol中先按特征頻率排序或者用腳本生成一個(gè)包含特征值順序信息的表方便Matlab讀取時(shí)對(duì)應(yīng)。4. Matlab讀取與預(yù)處理從導(dǎo)出文件到可用的模式向量4.1 文件目錄規(guī)劃與批量讀取Comsol導(dǎo)出數(shù)據(jù)的組織方式直接影響后續(xù)代碼的復(fù)雜度。我的做法是建一個(gè)data目錄里面存放所有kx點(diǎn)的導(dǎo)出文件命名格式是kx_01.txt這樣。每個(gè)文件的列順序設(shè)為x坐標(biāo)、y坐標(biāo)、Re(E_z)、Im(E_z)、Re(E_z_bloch)、Im(E_z_bloch)其中E_z_bloch就是前面說(shuō)的剝離布洛赫指數(shù)后的場(chǎng)分量。Matlab里讀取很簡(jiǎn)單不依賴任何額外工具箱function [x, y, Ez, u] loadFieldData(fileName) data readmatrix(fileName); x data(:, 1); y data(:, 2); Ez data(:, 3) 1i*data(:, 4); u data(:, 5) 1i*data(:, 6); end讀取后需要檢查每個(gè)文件的節(jié)點(diǎn)數(shù)是否一致。由于我關(guān)閉了自適應(yīng)網(wǎng)格整個(gè)掃描過(guò)程的網(wǎng)格沒(méi)有變化所以所有文件的節(jié)點(diǎn)數(shù)應(yīng)該完全相同。如果節(jié)點(diǎn)數(shù)不一致后面做內(nèi)積時(shí)直接按元素相乘就會(huì)出錯(cuò)。萬(wàn)一因?yàn)槟承┰蚓W(wǎng)格變了解決辦法是先在Matlab里用griddata插值到統(tǒng)一網(wǎng)格上但這個(gè)方法會(huì)引入誤差能不用就不用。4.2 模式場(chǎng)的歸一化處理Wilson loop公式中使用的態(tài)矢量應(yīng)當(dāng)是歸一化的。Comsol導(dǎo)出的模式場(chǎng)通常有自己的歸一化方式但不同kx點(diǎn)之間的模長(zhǎng)尺度可能不同因此必須在Matlab里對(duì)每個(gè)模式重新歸一化。對(duì)于一個(gè)復(fù)場(chǎng)向量u歸一化的代碼是u u / norm(u);norm函數(shù)默認(rèn)計(jì)算L2范數(shù)即sqrt(sum(abs(u).^2))。這個(gè)歸一化操作對(duì)后續(xù)內(nèi)積計(jì)算至關(guān)重要因?yàn)槿绻粴w一化每個(gè)內(nèi)積都會(huì)帶入一個(gè)未知的比例因子最后取log時(shí)這些比例因子不會(huì)抵消干凈會(huì)污染相位累積。還有一個(gè)細(xì)節(jié)如果模式場(chǎng)在空間不同位置上的幅度分布差異極大比如局域在某一層中那么簡(jiǎn)單的L2歸一化可能對(duì)遠(yuǎn)場(chǎng)區(qū)域的小幅度噪聲過(guò)于敏感。實(shí)際中我發(fā)現(xiàn)網(wǎng)格質(zhì)量足夠好時(shí)這個(gè)問(wèn)題不會(huì)太明顯如果確實(shí)出現(xiàn)異常可以只取結(jié)構(gòu)內(nèi)部區(qū)域的節(jié)點(diǎn)做內(nèi)積但不能完全丟棄周期邊界附近的點(diǎn)。4.3 節(jié)點(diǎn)一致性與內(nèi)積模檢驗(yàn)在開(kāi)始計(jì)算Zak相位之前我強(qiáng)烈建議先做一步診斷計(jì)算相鄰k點(diǎn)模式間的內(nèi)積模|?u(k_m)|u(k_{m1})?|。這個(gè)量在理想情況下應(yīng)當(dāng)非常接近1偏差通常小于千分之一。如果發(fā)現(xiàn)某個(gè)k點(diǎn)對(duì)之間的內(nèi)積模為0.5甚至更小說(shuō)明模式排序錯(cuò)亂或網(wǎng)格有問(wèn)題。這時(shí)直接算Zak相位肯定出錯(cuò)先把診斷問(wèn)題解決再繼續(xù)。內(nèi)積模診斷的Matlab代碼for m 1:Nk-1 inner dot(u{m}, u{m1}); fprintf(k index %d to %d, |u|u| %.6f\n, m, m1, abs(inner)); end當(dāng)你遇到內(nèi)積模明顯偏離1的情況多數(shù)原因是模式排序錯(cuò)亂就是第m個(gè)k點(diǎn)的第n條帶和第m1個(gè)k點(diǎn)的第n條帶實(shí)際上不是物理上的同一條能帶。這種問(wèn)題在帶隙較窄、能帶交叉區(qū)域經(jīng)常出現(xiàn)需要回到Comsol里查看能帶回線或者改用頻率排序之外的場(chǎng)重疊積分方法重新匹配能帶順序。5. Matlab計(jì)算Zak相位的核心代碼與結(jié)果分析5.1 單條能帶的Zak相位計(jì)算主程序假設(shè)你已經(jīng)把所有kx點(diǎn)的模式場(chǎng)存在一個(gè)cell數(shù)組u_list中數(shù)組長(zhǎng)度是Nk每個(gè)元素是一個(gè)Nnode×1的復(fù)向量。下面這段代碼實(shí)現(xiàn)了完整的Wilson loop計(jì)算clear; clc; % 參數(shù)定義 a 1e-6; % 晶格常數(shù) Nk 41; % 布里淵區(qū)采樣點(diǎn)數(shù) G 2*pi/a; % 倒格子基矢 kx linspace(-pi/a, pi/a, Nk1); kx kx(1:end-1); % 去掉最后一個(gè)重復(fù)點(diǎn)閉環(huán)由周期規(guī)范處理 % 這里需要先加載所有模式場(chǎng)數(shù)據(jù)到 u_list % u_list{m} 是第m個(gè)k點(diǎn)處某一條能帶的模式場(chǎng)向量 ZakSum 0; % 相鄰點(diǎn)內(nèi)積累積 for m 1:Nk-1 u1 u_list{m}; u2 u_list{m1}; u1 u1 / norm(u1); u2 u2 / norm(u2); inner dot(u1, u2); if abs(inner) 1e-8 error(內(nèi)積模接近0模式匹配失敗); end ZakSum ZakSum log(inner / abs(inner)); end % 邊界閉合項(xiàng)最后一個(gè)點(diǎn)與第一個(gè)點(diǎn)的連接 u_first u_list{1}; % k -pi/a u_last u_list{Nk}; % k -pi/a (Nk-1)*Δk u_first u_first / norm(u_first); u_last u_last / norm(u_last); % 這里假設(shè)導(dǎo)出時(shí)已經(jīng)是剝離了布洛赫指數(shù)的u場(chǎng) % 所以閉合項(xiàng)需要乘上e^{-iGx} x_coord x_list{1}; % x坐標(biāo)從導(dǎo)出文件中讀取 phaseFactor exp(-1i * G * x_coord); inner_boundary dot(u_last, u_first .* phaseFactor); ZakSum ZakSum log(inner_boundary / abs(inner_boundary)); % Zak相位弧度 Zak -imag(ZakSum); % 歸一化到 [0, pi) Zak mod(Zak, pi); fprintf(Zak phase %.4f rad\n, Zak);這段代碼里我特意用log形式而不是直接連乘后取arg原因是log形式可以逐步觀察相位累積過(guò)程更容易定位問(wèn)題出在哪兩個(gè)k點(diǎn)之間。5.2 多條能帶同時(shí)計(jì)算的循環(huán)結(jié)構(gòu)如果要計(jì)算前N條能帶需要在外層加一個(gè)能帶索引循環(huán)。每個(gè)能帶使用各自對(duì)應(yīng)的模式場(chǎng)數(shù)據(jù)。由于Comsol導(dǎo)出的特征頻率結(jié)果可能沒(méi)有按能帶排列Matlab側(cè)需要先做一個(gè)能帶分揀。一種簡(jiǎn)單有效的分揀方法是讀取每個(gè)kx點(diǎn)的全部特征頻率按頻率大小排序后依次對(duì)應(yīng)能帶1、能帶2等等。但在能帶交叉區(qū)域這種方法會(huì)出錯(cuò)。更可靠的方法是使用模式匹配從k_0出發(fā)用第n條帶在k_m的模式與k_{m1}的所有模式計(jì)算重疊積分選擇重疊積分模最大的模式作為同一條帶的延續(xù)。這個(gè)思路實(shí)現(xiàn)起來(lái)也不復(fù)雜for band 1:Nbands % 存儲(chǔ)該條帶所有k點(diǎn)模式 u_band cell(Nk, 1); for m 1:Nk if m 1 % 第一個(gè)k點(diǎn)按頻率排序取第band條 u_band{1} mode_data{1}{band}; else % 后續(xù)k點(diǎn)找與上一k點(diǎn)同帶模式重疊最大的模式 maxOverlap -1; bestIndex 1; for n 1:Nmodes overlap abs(dot(u_band{m-1}, mode_data{m}{n})); if overlap maxOverlap maxOverlap overlap; bestIndex n; end end u_band{m} mode_data{m}{bestIndex}; end end % 對(duì)u_band運(yùn)行Wilson loop計(jì)算 zakPhase computeZakPhase(u_band, x_coord, G); end這種基于模式重疊的能帶追蹤方法在帶隙較寬、模式區(qū)分明顯的體系中非常穩(wěn)健。如果體系出現(xiàn)近簡(jiǎn)并兩組模式的重疊積分都很接近1就需要進(jìn)入簡(jiǎn)并子空間做Wilson loop這是一個(gè)稍微復(fù)雜的升級(jí)版本文不展開(kāi)。5.3 結(jié)果判定怎么看Zak相位算對(duì)了算出來(lái)的Zak相位應(yīng)當(dāng)在0或π附近因?yàn)橹行膶?duì)稱結(jié)構(gòu)的Zak相位是量子化的。我這里給出一個(gè)實(shí)際計(jì)算例子。用之前說(shuō)的參數(shù)結(jié)構(gòu)晶格常數(shù)1μm高低折射率層厚度分別為0.7a和0.3a折射率3.45和1.0取41個(gè)k點(diǎn)。前四條能帶的計(jì)算結(jié)果如下能帶編號(hào)約化頻率范圍Zak相位計(jì)算值歸一化到[0, π)10 ~ 0.280.0023 rad020.32 ~ 0.523.1294 radπ30.54 ~ 0.733.1318 radπ40.78 ~ 0.920.0041 rad0與文獻(xiàn)對(duì)照這個(gè)二分之一結(jié)構(gòu)在折射率對(duì)比足夠大的情況下Zak相位按帶隙排列為0, π, π, 0符合預(yù)期。如果換成非中心對(duì)稱結(jié)構(gòu)例如改變第二層介質(zhì)的位置使原點(diǎn)的選取不再對(duì)稱Zak相位會(huì)明顯偏離0和π這也是一個(gè)很靈的敏感性測(cè)試。5.4 采樣點(diǎn)數(shù)對(duì)結(jié)果的影響k點(diǎn)采樣數(shù)量對(duì)Zak相位的影響比很多人想象的更大。我用同一模型測(cè)試了Nk從11到121變化時(shí)的結(jié)果。Nk11時(shí)計(jì)算誤差比較大有些能帶的Zak相位偏差達(dá)到0.1 rad級(jí)別Nk41時(shí)基本穩(wěn)定在0或π附近偏差小于0.01 rad再加密到121時(shí)偏差進(jìn)一步減小到0.001 rad量級(jí)。原因是Wilson loop中相鄰k點(diǎn)之間的相位增量|Δθ|不能超過(guò)π否則log函數(shù)的虛部會(huì)混疊。k點(diǎn)越密相鄰內(nèi)積相位差越小累積越準(zhǔn)確。建議至少取41個(gè)k點(diǎn)作為起步若要追求高精度或處理近簡(jiǎn)并能帶建議取81~121個(gè)點(diǎn)。求解時(shí)間和數(shù)據(jù)量會(huì)相應(yīng)增加但對(duì)單個(gè)一維模型而言完全在可接受范圍內(nèi)。6. 常見(jiàn)問(wèn)題與排查技巧實(shí)錄6.1 能帶模式排序錯(cuò)亂導(dǎo)致內(nèi)積跳變這是我自己做這個(gè)項(xiàng)目時(shí)遇到最多的一個(gè)問(wèn)題。現(xiàn)象是診斷內(nèi)積模時(shí)大部分k點(diǎn)對(duì)的內(nèi)積模在0.999以上但某個(gè)點(diǎn)特別低只有0.3甚至0.1。原因幾乎都是Comsol特征頻率求解器在每個(gè)kx掃描點(diǎn)返回的模式順序不是按能帶排列的頻率接近的兩條帶交叉時(shí)求解器返回順序會(huì)交換。解決辦法就是在Matlab里用模式匹配法做能帶追蹤而不是直接按頻率排序取第n條模式。另外還有一個(gè)技巧在Comsol中設(shè)置輔助掃描時(shí)勾選“使用前一個(gè)解作為初始猜測(cè)”可以讓相鄰kx點(diǎn)的模式順序更加一致減少后期處理工作量。6.2 邊界閉合項(xiàng)缺失導(dǎo)致相位不量子化如果Zak相位計(jì)算結(jié)果不是0或π但它們分布在某個(gè)中間值附近比如0.3π首先不要懷疑物理模型先檢查邊界項(xiàng)。邊界項(xiàng)是整個(gè)Wilson loop里唯一體現(xiàn)布里淵區(qū)周期性的地方漏掉它算出來(lái)的是一個(gè)開(kāi)路徑積分物理上沒(méi)有任何意義。我在代碼里特意把邊界閉合項(xiàng)單獨(dú)列出來(lái)就是希望讀者看清楚這一項(xiàng)的作用。還有一個(gè)容易搞混的地方周期規(guī)范算符到底是乘以e^{-iGx}還是e^{iGx}。這取決于你對(duì)正k方向的定義和Comsol中Floquet周期條件的k矢量符號(hào)設(shè)置。我建議你在小范圍內(nèi)用解析模型驗(yàn)證一次符號(hào)。最簡(jiǎn)單的驗(yàn)證對(duì)象是自由空間均勻介質(zhì)此時(shí)Zak相位應(yīng)當(dāng)為0如果符號(hào)寫(xiě)反會(huì)得到一個(gè)與路徑長(zhǎng)度相關(guān)的非零值。6.3 內(nèi)積模總是略小于1網(wǎng)格和節(jié)點(diǎn)不一致如果診斷顯示所有內(nèi)積模都在0.95左右雖然能算出Zak相位但總覺(jué)得不放心通常是兩個(gè)原因。一是網(wǎng)格不夠細(xì)離散化誤差大二是上下邊界或內(nèi)部界面處的場(chǎng)在相鄰k點(diǎn)之間發(fā)生了微小變化而這種變化是因?yàn)镕loquet邊界條件的數(shù)值實(shí)現(xiàn)不是完全一致的。解決方法把網(wǎng)格加密一個(gè)量級(jí)再跑一次如果內(nèi)積模明顯上升就是網(wǎng)格問(wèn)題。如果加密后內(nèi)積模還是0.96就要檢查是不是上下邊界條件選得不對(duì)。我在之前的模型中上下邊界用周期性邊界條件換來(lái)的是y方向的均勻場(chǎng)內(nèi)積模輕松到0.999以上。換成PEC邊界會(huì)引入y方向的橫向模式變化內(nèi)積模就會(huì)下降。6.4 Zak相位符號(hào)和文獻(xiàn)相反出現(xiàn)這種情況大概率不是計(jì)算錯(cuò)誤而是約定不同。有些文獻(xiàn)定義Zak相位時(shí)積分方向是從0到G有些是從-G/2到G/2符號(hào)差一個(gè)負(fù)號(hào)。還有的文獻(xiàn)在Wilson loop里取log之后用imag而不是-imag結(jié)果也是反號(hào)。我在文章里給出的公式和代碼采用的是最常見(jiàn)的約定θ -Im Σ log(...)方向沿k正方向。如果你的情況特殊統(tǒng)一改掉符號(hào)即可物理結(jié)論不受影響。6.5 特征頻率中出現(xiàn)不想要的雜散模式算特征頻率時(shí)設(shè)置的搜索基數(shù)越大越容易混入一些和物理問(wèn)題無(wú)關(guān)的模式。比如上下邊界條件造成的橫向高次模、求解器數(shù)值產(chǎn)生的非物理模式。這些模式在能帶圖上表現(xiàn)為一些偏離主能帶的光滑曲線之外的零散點(diǎn)。判斷是否為雜散模式有一個(gè)簡(jiǎn)單方法在Comsol后處理中查看該模式場(chǎng)的空間分布。物理模式應(yīng)該主要集中在介質(zhì)結(jié)構(gòu)內(nèi)場(chǎng)分布沿x方向周期變化均勻雜散模式往往在場(chǎng)分布上有明顯的橫向振蕩或者能量集中在邊界上。處理時(shí)把這些模式從導(dǎo)出列表中排除或者在Matlab里根據(jù)模式場(chǎng)的傅里葉成分做一個(gè)過(guò)濾。6.6 常見(jiàn)問(wèn)題速查表現(xiàn)象可能原因排查方法Zak相位不為0或π邊界閉合項(xiàng)缺失、原點(diǎn)不對(duì)稱、能帶追蹤錯(cuò)亂檢查邊界項(xiàng)代碼確認(rèn)結(jié)構(gòu)中心對(duì)稱診斷內(nèi)積模內(nèi)積模明顯小于1模式排序錯(cuò)亂、網(wǎng)格太粗、網(wǎng)格隨掃描變化用重疊積分做能帶追蹤加密網(wǎng)格關(guān)閉自適應(yīng)相位符號(hào)反號(hào)約定不同用均勻介質(zhì)驗(yàn)證符號(hào)統(tǒng)一Zak定義計(jì)算結(jié)果對(duì)掃描區(qū)間敏感沒(méi)有把起點(diǎn)終點(diǎn)按周期規(guī)范連接檢查布里淵區(qū)邊界處的k點(diǎn)和邊界閉合因子能帶圖有零散雜點(diǎn)上下邊界條件引入雜散模式檢查模式場(chǎng)空間分布過(guò)濾非物理模式7. 實(shí)際操作中最后想說(shuō)的話這套流程跑通以后我的最大感受是算Zak相位比算能帶要敏感得多。能帶數(shù)據(jù)差一點(diǎn)還能看出趨勢(shì)Zak相位有一處弄錯(cuò)就直接不量子化了反而逼著我把整個(gè)數(shù)值鏈路從頭到尾查了一遍。這是好事因?yàn)橥ㄟ^(guò)這個(gè)排查過(guò)程我對(duì)Comsol場(chǎng)導(dǎo)出、Floquet邊界條件的數(shù)值實(shí)現(xiàn)、以及Wilson loop的離散形式都有了更扎實(shí)的理解。如果你剛開(kāi)始做類似項(xiàng)目我建議先不要急著上復(fù)雜結(jié)構(gòu)就用最簡(jiǎn)單的兩個(gè)介質(zhì)層交替的一維光子晶體把整個(gè)流程跑通再用一個(gè)已知結(jié)果來(lái)驗(yàn)證你的代碼確認(rèn)無(wú)誤后再推廣到多層結(jié)構(gòu)、漸變結(jié)構(gòu)、或者帶損耗的體系。另外Comsol的LiveLink for Matlab可以省去文件導(dǎo)入手動(dòng)操作的麻煩算是錦上添花的工具沒(méi)有的話純文本導(dǎo)出加readmatrix也完全夠用。最后分享一個(gè)調(diào)試技巧在Matlab的循環(huán)里每隔幾個(gè)k點(diǎn)打印一次當(dāng)前累計(jì)相位和步進(jìn)相位畫(huà)出來(lái)看看。如果累計(jì)相位在某個(gè)位置出現(xiàn)了接近π的跳變那就是內(nèi)積相位越過(guò)分支切線了需要增加k點(diǎn)采樣密度或者調(diào)整數(shù)據(jù)處理方式。這種逐步觀察的習(xí)慣比最后只看一個(gè)最終數(shù)字有用得多。