域束縛態(tài)遠場偏振與空間能帶值計算指南)
簡介圍繞光子晶體連續(xù)域束縛態(tài)BIC的遠場偏振計算這份資料包面向從事COMSOL光學(xué)仿真、拓撲光子學(xué)及微納光子器件研究的高年級本科生、研究生和科研人員內(nèi)容涵蓋BIC理論基礎(chǔ)、k空間能帶與Q值分析、遠場偏振特性計算等整合了仿真模型思路與Matlab腳本可直接作為課題入門和復(fù)現(xiàn)計算的參考。包內(nèi)共15個文件包括6篇docx研究文檔理論梳理與計算說明、4個html過程記錄仿真分析步驟與結(jié)果討論、4張jpg結(jié)果圖像能帶/偏振圖等及1份txt說明文件壓縮包整體大小2.83MB結(jié)構(gòu)清晰便于按需查閱。目前已有42人學(xué)習(xí)瀏覽借助其中的文檔可掌握BIC參數(shù)掃描、Q值提取與遠場偏振投影方法Matlab腳本提供能帶繪制和偏振分析的具體代碼圖像則直觀展示計算結(jié)果適合正在搭建BIC仿真框架、需要對比驗證遠場偏振計算路徑的讀者。1. 光子晶體連續(xù)域束縛態(tài)這個“不輻射”的模式為什么值得算遠場偏振“光子晶體連續(xù)域束縛態(tài)遠場偏振計算含空間能帶值”這串關(guān)鍵詞拆開看就是三件事設(shè)計一個周期性結(jié)構(gòu)讓某個模式泡在連續(xù)輻射譜里卻不向外泄漏在離這個模式很近的動量點上算出它的遠場偏振方向和旋轉(zhuǎn)規(guī)律再把結(jié)果放回能帶圖上確認它落在空間能帶值的哪個位置。做光芯片激光選模、超表面窄帶濾波和二次諧波增強的人會被這串問題卡住BIC 的遠場輻射嚴格為零直接在奇點上算偏振只會得到一堆數(shù)值噪聲能帶圖也會失真。這里直接給出一條可復(fù)現(xiàn)的路線先用平面波展開法定位候選能帶再用 Stokes 參數(shù)做遠場偏振提取最后繞奇點走一圈驗證拓撲荷。2. 能帶計算路線從結(jié)構(gòu)建模到 BIC 候選模式的完整流程2.1 為什么選 PWE 找能帶、選 FDTD 算遠場兩條路線各管一段光子晶體能帶計算最常見的兩套工具是平面波展開法PWEMPB 就是代表實現(xiàn)和時域有限差分法FDTD。PWE 在周期性結(jié)構(gòu)里求解本征頻率非?？煲淮涡园颜麄€布里淵區(qū)路徑上的能帶算完特別適合先摸清模式分布但它默認結(jié)構(gòu)無損、周期無限拿不到模式泄漏率自然也算不出 Q 值。FDTD 能直接激勵結(jié)構(gòu)再測衰減Q 值給得很實但掃一個 k 點就要一個完整仿真把整條能帶掃出來成本太高。所以我的常規(guī)做法是兩段式用 PWE 圈出疑似 BIC 的能帶區(qū)間再用 FDTD 在那個 k 點附近做衰減測量驗證 Q 值是否隨分辨率增加而發(fā)散。物理模型也要先定住。這里采用二維方形晶格介質(zhì)柱結(jié)構(gòu)晶格常數(shù)設(shè)為 a圓柱半徑 r折射率 n。取二維模型而不是三維板是為了先把 xy 平面內(nèi)的輻射通道算清楚避開 z 方向 slab 模式帶來的額外泄漏項。對這個模型BIC 的判據(jù)很直接某個布洛赫模式的頻率落在背景模式的連續(xù)譜區(qū)間內(nèi)但它的對稱性與所有可輻射通道不匹配于是能量只被局域在結(jié)構(gòu)里無法耦合到遠場。2.2 最小建模腳本方形晶格介質(zhì)柱算到 BIC 候選點用一個可跑的 MPB Python 腳本把結(jié)構(gòu)建起來。這里以 meep.mpb 的 Python 接口為例不同小版本的參數(shù)名略有差異但主體流程一致。import meep as mp from meep import mpb resolution 32 # 網(wǎng)格分辨率單位是每個晶格常數(shù)內(nèi)的采樣點數(shù) r 0.2 # 介質(zhì)柱半徑單位是晶格常數(shù) a n_si 3.5 # 柱體折射率硅在近紅外波段約 3.5 geometry_lattice mp.Lattice( basis1mp.Vector3(1.0, 0.0), basis2mp.Vector3(0.0, 1.0) ) geometry [ mp.Cylinder( radiusr, heightmp.inf, # 二維模型認為柱體沿 z 向無限長 materialmp.Medium(indexn_si), axismp.Vector3(0, 0, 1) ) ] ms mpb.ModeSolver( geometrygeometry, geometry_latticegeometry_lattice, resolutionresolution, num_bands6 ) ms.run_tm() # 計算 TM 極化電場沿 z 方向的能帶這里有兩個參數(shù)直接影響后續(xù)找 BIC 的難度。第一個是resolution建議從 16 開始一旦看到某個能帶交叉點疑似出現(xiàn)“嵌入連續(xù)譜”的特征再提到 64 復(fù)核 Q 值趨勢網(wǎng)格不夠時BIC 會被數(shù)值泄漏偽裝成普通泄漏模。第二個是r半徑 0.2 只是一個起點不同半徑會移動能帶位置和對稱性匹配條件后面做參數(shù)掃描時它是最重要的變量。run_tm()只算 TM 模是因為介質(zhì)柱結(jié)構(gòu)里 BIC 通常更容易出現(xiàn)在 TM 能帶的交叉點附近TE 模也可以算但建議分開跑避免簡并模式混在一起難判斷。2.3 從能帶圖看 Q 值連續(xù)域里“不泄漏”的模式長什么樣MPB 輸出的是實頻率不能直接讀出 Q 值。判斷一個點是 BIC 候選先看它的頻率是否落在連續(xù)譜區(qū)間內(nèi)再看它能否輻射。能帶圖上最直觀的信號是某條能帶的頻率穿過其他模式形成的連續(xù)譜帶但它所在的分支沒有表現(xiàn)出明顯的帶隙或反交叉顏色/對稱性標記與相鄰分支不同。這時候才值得對它做 FDTD 復(fù)核。復(fù)核 Q 值的做法是在結(jié)構(gòu)中心放一個窄帶偶極子源頻率對準候選模式的頻點然后在結(jié)構(gòu)內(nèi)放一個探針記錄電場衰減曲線。衰減滿足 |E(t)| ∝ exp(-ωt/2Q)對包絡(luò)做對數(shù)線性擬合就能得到 Q。真正的 BIC 在數(shù)值上表現(xiàn)為 Q 隨網(wǎng)格分辨率持續(xù)上升沒有收斂上限的趨勢如果 Q 穩(wěn)定在某個有限值說明它只是一個高 Q 準 BIC而不是嚴格 BIC。以下表格是我用來區(qū)分兩者的快速對照。特征準 BIC嚴格 BICQ 隨分辨率變化先上升后收斂持續(xù)上升無上限遠場偏振幅值有限值但很小嚴格趨近于零繞奇點拓撲荷整數(shù)但可能混入雜散相位清晰的整數(shù)對稱性要求弱對稱破缺即可需要嚴格匹配這一步容易忽略的是對稱性條件。BIC 之所以不輻射本質(zhì)上是模式的遠場輻射通道因?qū)ΨQ性被禁止網(wǎng)格一旦把結(jié)構(gòu)的旋轉(zhuǎn)對稱性破壞這條禁令就失效了。所以建模時要把介質(zhì)柱的中心放在網(wǎng)格元胞的對稱中心上還要保證圓柱的圓形輪廓不被矩形網(wǎng)格過度“切”出棱角后面第五章會專門展開這個坑。3. 遠場偏振計算用 Stokes 參數(shù)從復(fù)電場里提取偏振態(tài)3.1 遠場偏振為什么不能直接看近場橫波條件先滿足近場里的電場分量和遠場傳播方向并不垂直這是計算偏振時最容易犯的錯誤。結(jié)構(gòu)表面的近場電場里有很強的 Ez 分量而遠場是橫波電場必須垂直于傳播方向。直接拿近場復(fù)矢量去算偏振得到的結(jié)果和遠場實際輻射偏振完全對不上。更麻煩的是 BIC 附近的近場振幅非常大因為它把能量緊緊鎖在結(jié)構(gòu)里而同一模式的遠場振幅趨近于零。用幅度大的近場去推幅度小的遠場數(shù)值噪聲會徹底淹沒真實偏振信息。正確做法是先選定遠場觀察方向單位波矢記為 k?。把 k? 與全局坐標的 z? 做叉積得到投影平面的一組正交基 e1、e2再把復(fù)雜電場矢量投影到這兩個基上。這一步完成后剩下的就是標準偏振分析。實際項目里BIC 的遠場偏振通常是在泄漏輻射顯微鏡里測的測量方向正好對應(yīng)布里淵區(qū)里的某個動量方向數(shù)值計算時我們就在該動量方向構(gòu)造 k?。3.2 用 Stokes 參數(shù)做偏振提取一個可以直接抄的函數(shù)Stokes 參數(shù)的優(yōu)點是它完整描述偏振態(tài)包括強度、線偏振方位角和橢圓率而且不需要假設(shè)電場是線偏振的。這里給出一個已經(jīng)封裝好的函數(shù)輸入某個 k 方向的復(fù)電場矢量返回 Stokes 分量和偏振方位角。import numpy as np def far_field_stokes(E_cplx, kvec): # E_cplx: 長度為 3 的復(fù)電場矢量 [Ex, Ey, Ez] # kvec: 長度為 3 的實數(shù)波矢無需歸一化 E np.asarray(E_cplx, dtypecomplex) k np.asarray(kvec, dtypefloat) k_hat k / np.linalg.norm(k) # 構(gòu)造與傳播方向正交的兩個基矢量 e1, e2 z_hat np.array([0.0, 0.0, 1.0]) e1 np.cross(k_hat, z_hat) if np.linalg.norm(e1) 1e-12: # 當 k 與 z 軸平行時退化為固定取 x 方向 e1 np.array([1.0, 0.0, 0.0]) e1 e1 / np.linalg.norm(e1) e2 np.cross(k_hat, e1) # 投影到橫向平面 E1 np.dot(E, e1) E2 np.dot(E, e2) # Stokes 參數(shù) S0 np.abs(E1)**2 np.abs(E2)**2 S1 np.abs(E1)**2 - np.abs(E2)**2 S2 2.0 * np.real(E1 * np.conj(E2)) S3 -2.0 * np.imag(E1 * np.conj(E2)) # 偏振方位角范圍在 0 到 pi 之間 psi 0.5 * np.arctan2(S2, S1) return S0, S1, S2, S3, psi調(diào)用這個函數(shù)之前要確認 E_cplx 來自同一個時間切片的復(fù)場相位關(guān)系不能亂。S0 表示總強度S1 和 S2 決定線偏振方向S3 表示圓偏振分量。BIC 點附近 S0 趨近于零此時 psi 會因為數(shù)值噪聲到處亂跳這不是程序 bug而是物理上該點沒有遠場輻射、偏振本來就沒有定義。所以偏振計算要刻意避開嚴格 BIC 點取它周圍一圈動量點來計算。3.3 從近場到遠場傅里葉變換與相位歸一化的兩個參數(shù)遠場復(fù)電場可以從 FDTD 的近場監(jiān)測面做空間傅里葉變換得到。近場平面上每個點的復(fù)電場 E(x, y) 做一個二維 FFT得到的頻域分量就對應(yīng)不同遠場方向的幅值和相位。這一步有兩個參數(shù)直接決定結(jié)果質(zhì)量近場監(jiān)測面的寬度必須是整數(shù)個晶格周期否則 FFT 會引入 sinc 振鈴遠場角度譜上出現(xiàn)一堆偽峰時間窗口長度要覆蓋模式衰減到足夠小否則頻率分辨率不夠BIC 附近極窄的偏振變化會被抹掉。def near_to_far_spectrum(E_near): # E_near: 二維復(fù)數(shù)數(shù)組形狀為 (ny, nx) # 先做空間 FFT再平移到低頻居中 E_far np.fft.fftshift(np.fft.fft2(E_near)) # 返回的每個元素對應(yīng)一個 (kx, ky) 采樣點 return E_farFFT 輸出的索引坐標與動量 (kx, ky) 的對應(yīng)關(guān)系取決于采樣網(wǎng)格間距。實空間采樣點間距 dx 決定動量空間總范圍采樣點數(shù)量決定動量分辨率。用這個結(jié)果作為far_field_stokes的輸入時要注意 E_near 必須包含完整的相位信息不能用強度值替代。取相位時所有監(jiān)測點要從同一個時間參考點導(dǎo)出通常先讓結(jié)構(gòu)響應(yīng)進入穩(wěn)態(tài)再記錄一個時間快照。提示離 BIC 太近時S0 會被數(shù)值噪聲主導(dǎo)。實際計算中我一般取布里淵區(qū)半徑的 1%~5% 作為偏移量并在圓上取 8~12 個點這樣既能看到完整的偏振旋轉(zhuǎn)又不會踩進零輻射區(qū)。4. 空間能帶值提取把能帶曲面、BIC 位置與偏振分布對應(yīng)起來4.1 把能帶鋪成空間曲面空間能帶值指的是哪一組數(shù)據(jù)“空間能帶值”這個說法在光子晶體語境里容易產(chǎn)生歧義。它不是指實空間坐標上每個點的能帶頻率那是聲子晶體局域態(tài)密度的處理方式在光子晶體 BIC 計算里它通常指動量空間里的能帶取值把第一布里淵區(qū)當作坐標平面每個點 (kx, ky) 對應(yīng)某個模式的一個本征頻率把這一組頻率鋪成曲面就是空間能帶值。BIC 就藏在這張曲面上位置往往在能帶交叉形成的錐形點附近。構(gòu)建這張曲面時不要只沿 Γ-X-M 路徑算一維能帶那樣會漏掉 BIC 所在的方向。正確做法是在二維動量網(wǎng)格上做掃描kx 和 ky 都取第一布里淵區(qū)一半的范圍網(wǎng)格數(shù)從 21×21 起步。BIC 本身是一個零維奇點網(wǎng)格太粗會直接把它跳過去所以第一輪粗掃只是定位可疑區(qū)域第二輪要在可疑區(qū)域附近做局部加密掃描。4.2 從這里判斷 BIC 位置能帶轉(zhuǎn)折點與 Q 值不發(fā)散的模式把能帶曲面畫出來后BIC 候選點有幾個可辨識的特征頻率等高線圍繞某個點出現(xiàn)渦旋狀分布能帶曲面在該點附近呈現(xiàn)尖銳的漏斗形狀沿任意方向穿過該點頻率連續(xù)但模式的遠場輻射突然消失。下面這段代碼把從求解器導(dǎo)出的頻率表填充成能帶曲面再做可視化。import numpy as np import matplotlib.pyplot as plt def build_band_surface(kx_grid, ky_grid, band_index, freq_table): # kx_grid, ky_grid: 由 np.meshgrid 生成的動量坐標網(wǎng)格 # freq_table: 求解器導(dǎo)出的數(shù)組形狀為 (n_k, n_band) # band_index: 要繪制的能帶序號從 0 開始計數(shù) surface np.zeros_like(kx_grid) nky, nkx kx_grid.shape for i in range(nky): for j in range(nkx): k_idx i * nkx j surface[i, j] freq_table[k_idx, band_index] return surface # 示例在 21x21 網(wǎng)格上計算第三條能帶 kx np.linspace(-0.5, 0.5, 21) ky np.linspace(-0.5, 0.5, 21) KX, KY np.meshgrid(kx, ky) # freq_table 由能帶求解器輸出每個 k 點取若干條能帶頻率 # surface build_band_surface(KX, KY, band_index2, freq_tablefreq_table) plt.contourf(KX, KY, surface, levels40, cmapviridis) plt.colorbar(labelnormalized frequency) plt.xlabel(r$k_x$ (units of $2\pi/a$)) plt.ylabel(r$k_y$ (units of $2\pi/a$))注意band_index的取值要和第二、三章的能帶順序保持一致。TM 模式下第三條能帶通常是最容易出 BIC 的位置因為前兩條能帶往往落在空氣光錐之下輻射通道本身就不存在。判斷曲面上的漏斗點是不是真 BIC還要回到 2.3 節(jié)的 Q 值復(fù)核單靠頻率曲面只能確定“能帶在這里發(fā)生了交叉”不能確定“交叉點不輻射”。4.3 熱圖輔助把偏振方位角疊到能帶曲面上空間能帶值除了頻率還可以把遠場偏振方位角作為第二層信息疊加上去。做法是在能帶曲面的等高線圖上把每個 k 點的偏振方位角畫成短線短線方向代表該動量的線偏振方向長度代表 S0 的相對強度。BIC 附近會看到清晰的旋渦狀短線分布繞奇點走一圈短線方向旋轉(zhuǎn) 2π 的整數(shù)倍。這個旋渦結(jié)構(gòu)直接對應(yīng) BIC 的拓撲荷比單純看能帶交叉更可靠。疊加圖還能幫你判斷 BIC 的對稱性類型。如果旋渦旋轉(zhuǎn)方向相反說明拓撲荷符號不同對應(yīng)的偏振構(gòu)型也不同。這類信息在做超表面偏振調(diào)控時非常關(guān)鍵目標偏振轉(zhuǎn)換器件的功能本質(zhì)上就是讓遠場偏振沿特定路徑繞 BIC 演化。繪圖時要注意角度顯示的范圍方位角是 0 到 π 的周期量直接用原始值畫會看到劇烈跳變建議先對角度做 unwrap 處理再映射到顏色或短線方向。5. 常見問題排查BIC 偏振與能帶計算最容易翻車的四個坑5.1 BIC 點上算偏振結(jié)果全是噪聲現(xiàn)象把 k 點正好設(shè)在 BIC 位置調(diào)用far_field_stokes后 S1、S2 在每次運行之間完全不一致方位角隨機跳動。原因BIC 的遠場輻射嚴格為零求解器算出來的只是浮點噪聲。Stokes 參數(shù)的 S1、S2 都由遠場復(fù)電場乘積構(gòu)成真實信號趨近零時噪聲就成了主導(dǎo)項。解決永遠不要在嚴格奇點上算偏振。繞 BIC 取一個半徑為 Δk 的小圓Δk 取布里淵區(qū)半徑的 1% 到 5%圓上相鄰點之間夾角均勻。設(shè)定幾個不同半徑各算一圈如果方位角旋轉(zhuǎn)角一致說明結(jié)果收斂如果旋轉(zhuǎn)角隨半徑變化就要加密動量網(wǎng)格可能是附近還有其他奇點干擾。5.2 能帶曲面出現(xiàn)多余的折疊分支現(xiàn)象空間能帶曲面圖上有大量重復(fù)的等高線條紋看起來像是多了一整條能帶。原因掃描動量網(wǎng)格時沒有使用不可約布里淵區(qū)而是把整個方形布里淵區(qū)掃了一遍。晶格對稱性導(dǎo)致某些能帶在多個區(qū)域內(nèi)重復(fù)出現(xiàn)疊加在一起就成了偽造的多分支結(jié)構(gòu)。解決先確認結(jié)構(gòu)屬于哪個點群方形晶格通常取不可約布里淵區(qū)的三角區(qū)域掃描即可。如果做全布里淵區(qū)掃描后處理時按點群對稱性把等效點合并或者只在等高線圖里保留頻率最低的那條分支。這個坑在第四章代碼里特別明顯我當時第一次畫出“八爪魚”狀的曲面排查了半天才發(fā)現(xiàn)是掃描范圍選錯了。5.3 Q 值上不去總是封頂在某個有限值現(xiàn)象FDTD 算出的 Q 值隨分辨率從 16 升到 64漲到幾千后就不動了再怎么加密網(wǎng)格也上不去。原因結(jié)構(gòu)本身是準 BIC不是嚴格 BIC或者介質(zhì)柱的圓形邊界被直角網(wǎng)格離散后引入了非物理的對稱性破缺。網(wǎng)格一粗圓柱輪廓變成鋸齒狀旋轉(zhuǎn)對稱性被破壞本應(yīng)禁止的輻射通道被重新打開。解決先檢查結(jié)構(gòu)是否滿足嚴格對稱條件。圓柱中心要落在網(wǎng)格元胞中心半徑在網(wǎng)格上的離散誤差用亞像素平滑處理MEEP 里可以開啟 subpixel 默認平均。如果已經(jīng)做了這些Q 值還封頂說明這個幾何參數(shù)點本身就不是 BIC 條件沿能帶方向掃幾個相鄰參數(shù)讓 r 或折射率微調(diào)再重新看 Q 值趨勢。5.4 近場到遠場頻譜出現(xiàn)非物理的旁瓣現(xiàn)象near_to_far_spectrum得到的遠場角度譜里目標方向附近出現(xiàn)一串規(guī)則的偽峰間距相等。原因近場監(jiān)測面寬度不是晶格周期的整數(shù)倍。監(jiān)測面被邊界截斷后FFT 把邊界的跳變卷進頻譜產(chǎn)生 sinc 振鈴。偽峰間距恰好對應(yīng)截斷窗口的尺寸。解決把監(jiān)測面寬度精確設(shè)置為整數(shù)個周期。如果結(jié)構(gòu)本身有缺陷或需要局部激勵先在監(jiān)測面加一個窗函數(shù)做平滑再用 FFT。加窗會犧牲一點角分辨率但能換來干凈的頻譜。防止這類問題的習(xí)慣做法是在設(shè)計仿真時就記錄監(jiān)測面的起點坐標保證它正好落在周期邊界上。提示以上四個坑都會在“先算能帶、再算遠場”的流程里暴露。每改一個參數(shù)先記錄網(wǎng)格分辨率和 k 點偏移量不然翻車后根本不知道是幾何參數(shù)變了還是數(shù)值設(shè)置變了。6. 兩個驗證技巧對稱性自檢與繞 BIC 追蹤偏振旋轉(zhuǎn)6.1 用對稱性做第一次自檢拿到能帶曲面和遠場偏振結(jié)果后先別急著往下做器件設(shè)計。跑一個對稱性檢查把整個結(jié)構(gòu)旋轉(zhuǎn) 90 度重新算一遍能帶和 Stoke驗證 KX、KY 坐標相應(yīng)旋轉(zhuǎn)后頻率沒有變化。如果旋轉(zhuǎn)后的偏振方位角沒有跟著旋轉(zhuǎn)說明投影坐標系定義有問題要回 3.2 節(jié)檢查 e1、e2 的構(gòu)造——這類問題不是物理錯了是坐標系沒對齊。這個檢查只需要多跑一次能帶掃描成本很低卻能把一半的坐標系 bug 擋在門外。6.2 繞 BIC 走一圈追蹤偏振方向驗證拓撲荷繞 BIC 繞一圈看偏振方位角的總旋轉(zhuǎn)量這是 BIC 所有數(shù)值結(jié)果里最有說服力的一致性驗證。拓撲荷由總旋轉(zhuǎn)角除以 2π 得到理論上必須是整數(shù)。如果算出 0.5 或 1.7 這種值說明網(wǎng)格、動量偏移或相位歸一化至少有一項出了問題。下面這個函數(shù)直接用小圓采樣計算繞數(shù)。def winding_number(k_center, dk, n_points12, get_field): # k_center: BIC 所在 (kx, ky) # dk: 繞圈半徑 # get_field: 輸入 (kx, ky)返回該方向遠場復(fù)電場 psi_list [] for idx in range(n_points): angle 2.0 * np.pi * idx / n_points kx k_center[0] dk * np.cos(angle) ky k_center[1] dk * np.sin(angle) E get_field(kx, ky) S0, S1, S2, S3, psi far_field_stokes(E, [kx, ky, 0.0]) psi_list.append(psi) # unwrap 防止方位角在 ±pi/2 處跳變 unwrapped np.unwrap(2.0 * np.array(psi_list)) total (unwrapped[-1] - unwrapped[0]) / 2.0 return int(round(total / (2.0 * np.pi)))這個小函數(shù)現(xiàn)在是我所有 BIC 計算的標配。每次換結(jié)構(gòu)參數(shù)、換網(wǎng)格密度我都會先用它驗證一遍繞數(shù)再去做投影器件設(shè)計。如果繞數(shù)不是整數(shù)我會回到第一章的建模流程重新檢查不會繼續(xù)往下算。代價只是多花十分鐘但能省下一個晚上重算的時間——這是我在這類計算里最實在的血淚經(jīng)驗希望幫到你。本文還有配套的精品資源點擊獲取 光子晶體連續(xù)域束縛態(tài)遠場偏振與空間能帶值計算指南
簡介圍繞光子晶體連續(xù)域束縛態(tài)BIC的遠場偏振計算這份資料包面向從事COMSOL光學(xué)仿真、拓撲光子學(xué)及微納光子器件研究的高年級本科生、研究生和科研人員內(nèi)容涵蓋BIC理論基礎(chǔ)、k空間能帶與Q值分析、遠場偏振特性計算等整合了仿真模型思路與Matlab腳本可直接作為課題入門和復(fù)現(xiàn)計算的參考。包內(nèi)共15個文件包括6篇docx研究文檔理論梳理與計算說明、4個html過程記錄仿真分析步驟與結(jié)果討論、4張jpg結(jié)果圖像能帶/偏振圖等及1份txt說明文件壓縮包整體大小2.83MB結(jié)構(gòu)清晰便于按需查閱。目前已有42人學(xué)習(xí)瀏覽借助其中的文檔可掌握BIC參數(shù)掃描、Q值提取與遠場偏振投影方法Matlab腳本提供能帶繪制和偏振分析的具體代碼圖像則直觀展示計算結(jié)果適合正在搭建BIC仿真框架、需要對比驗證遠場偏振計算路徑的讀者。1. 光子晶體連續(xù)域束縛態(tài)這個“不輻射”的模式為什么值得算遠場偏振“光子晶體連續(xù)域束縛態(tài)遠場偏振計算含空間能帶值”這串關(guān)鍵詞拆開看就是三件事設(shè)計一個周期性結(jié)構(gòu)讓某個模式泡在連續(xù)輻射譜里卻不向外泄漏在離這個模式很近的動量點上算出它的遠場偏振方向和旋轉(zhuǎn)規(guī)律再把結(jié)果放回能帶圖上確認它落在空間能帶值的哪個位置。做光芯片激光選模、超表面窄帶濾波和二次諧波增強的人會被這串問題卡住BIC 的遠場輻射嚴格為零直接在奇點上算偏振只會得到一堆數(shù)值噪聲能帶圖也會失真。這里直接給出一條可復(fù)現(xiàn)的路線先用平面波展開法定位候選能帶再用 Stokes 參數(shù)做遠場偏振提取最后繞奇點走一圈驗證拓撲荷。2. 能帶計算路線從結(jié)構(gòu)建模到 BIC 候選模式的完整流程2.1 為什么選 PWE 找能帶、選 FDTD 算遠場兩條路線各管一段光子晶體能帶計算最常見的兩套工具是平面波展開法PWEMPB 就是代表實現(xiàn)和時域有限差分法FDTD。PWE 在周期性結(jié)構(gòu)里求解本征頻率非?？煲淮涡园颜麄€布里淵區(qū)路徑上的能帶算完特別適合先摸清模式分布但它默認結(jié)構(gòu)無損、周期無限拿不到模式泄漏率自然也算不出 Q 值。FDTD 能直接激勵結(jié)構(gòu)再測衰減Q 值給得很實但掃一個 k 點就要一個完整仿真把整條能帶掃出來成本太高。所以我的常規(guī)做法是兩段式用 PWE 圈出疑似 BIC 的能帶區(qū)間再用 FDTD 在那個 k 點附近做衰減測量驗證 Q 值是否隨分辨率增加而發(fā)散。物理模型也要先定住。這里采用二維方形晶格介質(zhì)柱結(jié)構(gòu)晶格常數(shù)設(shè)為 a圓柱半徑 r折射率 n。取二維模型而不是三維板是為了先把 xy 平面內(nèi)的輻射通道算清楚避開 z 方向 slab 模式帶來的額外泄漏項。對這個模型BIC 的判據(jù)很直接某個布洛赫模式的頻率落在背景模式的連續(xù)譜區(qū)間內(nèi)但它的對稱性與所有可輻射通道不匹配于是能量只被局域在結(jié)構(gòu)里無法耦合到遠場。2.2 最小建模腳本方形晶格介質(zhì)柱算到 BIC 候選點用一個可跑的 MPB Python 腳本把結(jié)構(gòu)建起來。這里以 meep.mpb 的 Python 接口為例不同小版本的參數(shù)名略有差異但主體流程一致。import meep as mp from meep import mpb resolution 32 # 網(wǎng)格分辨率單位是每個晶格常數(shù)內(nèi)的采樣點數(shù) r 0.2 # 介質(zhì)柱半徑單位是晶格常數(shù) a n_si 3.5 # 柱體折射率硅在近紅外波段約 3.5 geometry_lattice mp.Lattice( basis1mp.Vector3(1.0, 0.0), basis2mp.Vector3(0.0, 1.0) ) geometry [ mp.Cylinder( radiusr, heightmp.inf, # 二維模型認為柱體沿 z 向無限長 materialmp.Medium(indexn_si), axismp.Vector3(0, 0, 1) ) ] ms mpb.ModeSolver( geometrygeometry, geometry_latticegeometry_lattice, resolutionresolution, num_bands6 ) ms.run_tm() # 計算 TM 極化電場沿 z 方向的能帶這里有兩個參數(shù)直接影響后續(xù)找 BIC 的難度。第一個是resolution建議從 16 開始一旦看到某個能帶交叉點疑似出現(xiàn)“嵌入連續(xù)譜”的特征再提到 64 復(fù)核 Q 值趨勢網(wǎng)格不夠時BIC 會被數(shù)值泄漏偽裝成普通泄漏模。第二個是r半徑 0.2 只是一個起點不同半徑會移動能帶位置和對稱性匹配條件后面做參數(shù)掃描時它是最重要的變量。run_tm()只算 TM 模是因為介質(zhì)柱結(jié)構(gòu)里 BIC 通常更容易出現(xiàn)在 TM 能帶的交叉點附近TE 模也可以算但建議分開跑避免簡并模式混在一起難判斷。2.3 從能帶圖看 Q 值連續(xù)域里“不泄漏”的模式長什么樣MPB 輸出的是實頻率不能直接讀出 Q 值。判斷一個點是 BIC 候選先看它的頻率是否落在連續(xù)譜區(qū)間內(nèi)再看它能否輻射。能帶圖上最直觀的信號是某條能帶的頻率穿過其他模式形成的連續(xù)譜帶但它所在的分支沒有表現(xiàn)出明顯的帶隙或反交叉顏色/對稱性標記與相鄰分支不同。這時候才值得對它做 FDTD 復(fù)核。復(fù)核 Q 值的做法是在結(jié)構(gòu)中心放一個窄帶偶極子源頻率對準候選模式的頻點然后在結(jié)構(gòu)內(nèi)放一個探針記錄電場衰減曲線。衰減滿足 |E(t)| ∝ exp(-ωt/2Q)對包絡(luò)做對數(shù)線性擬合就能得到 Q。真正的 BIC 在數(shù)值上表現(xiàn)為 Q 隨網(wǎng)格分辨率持續(xù)上升沒有收斂上限的趨勢如果 Q 穩(wěn)定在某個有限值說明它只是一個高 Q 準 BIC而不是嚴格 BIC。以下表格是我用來區(qū)分兩者的快速對照。特征準 BIC嚴格 BICQ 隨分辨率變化先上升后收斂持續(xù)上升無上限遠場偏振幅值有限值但很小嚴格趨近于零繞奇點拓撲荷整數(shù)但可能混入雜散相位清晰的整數(shù)對稱性要求弱對稱破缺即可需要嚴格匹配這一步容易忽略的是對稱性條件。BIC 之所以不輻射本質(zhì)上是模式的遠場輻射通道因?qū)ΨQ性被禁止網(wǎng)格一旦把結(jié)構(gòu)的旋轉(zhuǎn)對稱性破壞這條禁令就失效了。所以建模時要把介質(zhì)柱的中心放在網(wǎng)格元胞的對稱中心上還要保證圓柱的圓形輪廓不被矩形網(wǎng)格過度“切”出棱角后面第五章會專門展開這個坑。3. 遠場偏振計算用 Stokes 參數(shù)從復(fù)電場里提取偏振態(tài)3.1 遠場偏振為什么不能直接看近場橫波條件先滿足近場里的電場分量和遠場傳播方向并不垂直這是計算偏振時最容易犯的錯誤。結(jié)構(gòu)表面的近場電場里有很強的 Ez 分量而遠場是橫波電場必須垂直于傳播方向。直接拿近場復(fù)矢量去算偏振得到的結(jié)果和遠場實際輻射偏振完全對不上。更麻煩的是 BIC 附近的近場振幅非常大因為它把能量緊緊鎖在結(jié)構(gòu)里而同一模式的遠場振幅趨近于零。用幅度大的近場去推幅度小的遠場數(shù)值噪聲會徹底淹沒真實偏振信息。正確做法是先選定遠場觀察方向單位波矢記為 k?。把 k? 與全局坐標的 z? 做叉積得到投影平面的一組正交基 e1、e2再把復(fù)雜電場矢量投影到這兩個基上。這一步完成后剩下的就是標準偏振分析。實際項目里BIC 的遠場偏振通常是在泄漏輻射顯微鏡里測的測量方向正好對應(yīng)布里淵區(qū)里的某個動量方向數(shù)值計算時我們就在該動量方向構(gòu)造 k?。3.2 用 Stokes 參數(shù)做偏振提取一個可以直接抄的函數(shù)Stokes 參數(shù)的優(yōu)點是它完整描述偏振態(tài)包括強度、線偏振方位角和橢圓率而且不需要假設(shè)電場是線偏振的。這里給出一個已經(jīng)封裝好的函數(shù)輸入某個 k 方向的復(fù)電場矢量返回 Stokes 分量和偏振方位角。import numpy as np def far_field_stokes(E_cplx, kvec): # E_cplx: 長度為 3 的復(fù)電場矢量 [Ex, Ey, Ez] # kvec: 長度為 3 的實數(shù)波矢無需歸一化 E np.asarray(E_cplx, dtypecomplex) k np.asarray(kvec, dtypefloat) k_hat k / np.linalg.norm(k) # 構(gòu)造與傳播方向正交的兩個基矢量 e1, e2 z_hat np.array([0.0, 0.0, 1.0]) e1 np.cross(k_hat, z_hat) if np.linalg.norm(e1) 1e-12: # 當 k 與 z 軸平行時退化為固定取 x 方向 e1 np.array([1.0, 0.0, 0.0]) e1 e1 / np.linalg.norm(e1) e2 np.cross(k_hat, e1) # 投影到橫向平面 E1 np.dot(E, e1) E2 np.dot(E, e2) # Stokes 參數(shù) S0 np.abs(E1)**2 np.abs(E2)**2 S1 np.abs(E1)**2 - np.abs(E2)**2 S2 2.0 * np.real(E1 * np.conj(E2)) S3 -2.0 * np.imag(E1 * np.conj(E2)) # 偏振方位角范圍在 0 到 pi 之間 psi 0.5 * np.arctan2(S2, S1) return S0, S1, S2, S3, psi調(diào)用這個函數(shù)之前要確認 E_cplx 來自同一個時間切片的復(fù)場相位關(guān)系不能亂。S0 表示總強度S1 和 S2 決定線偏振方向S3 表示圓偏振分量。BIC 點附近 S0 趨近于零此時 psi 會因為數(shù)值噪聲到處亂跳這不是程序 bug而是物理上該點沒有遠場輻射、偏振本來就沒有定義。所以偏振計算要刻意避開嚴格 BIC 點取它周圍一圈動量點來計算。3.3 從近場到遠場傅里葉變換與相位歸一化的兩個參數(shù)遠場復(fù)電場可以從 FDTD 的近場監(jiān)測面做空間傅里葉變換得到。近場平面上每個點的復(fù)電場 E(x, y) 做一個二維 FFT得到的頻域分量就對應(yīng)不同遠場方向的幅值和相位。這一步有兩個參數(shù)直接決定結(jié)果質(zhì)量近場監(jiān)測面的寬度必須是整數(shù)個晶格周期否則 FFT 會引入 sinc 振鈴遠場角度譜上出現(xiàn)一堆偽峰時間窗口長度要覆蓋模式衰減到足夠小否則頻率分辨率不夠BIC 附近極窄的偏振變化會被抹掉。def near_to_far_spectrum(E_near): # E_near: 二維復(fù)數(shù)數(shù)組形狀為 (ny, nx) # 先做空間 FFT再平移到低頻居中 E_far np.fft.fftshift(np.fft.fft2(E_near)) # 返回的每個元素對應(yīng)一個 (kx, ky) 采樣點 return E_farFFT 輸出的索引坐標與動量 (kx, ky) 的對應(yīng)關(guān)系取決于采樣網(wǎng)格間距。實空間采樣點間距 dx 決定動量空間總范圍采樣點數(shù)量決定動量分辨率。用這個結(jié)果作為far_field_stokes的輸入時要注意 E_near 必須包含完整的相位信息不能用強度值替代。取相位時所有監(jiān)測點要從同一個時間參考點導(dǎo)出通常先讓結(jié)構(gòu)響應(yīng)進入穩(wěn)態(tài)再記錄一個時間快照。提示離 BIC 太近時S0 會被數(shù)值噪聲主導(dǎo)。實際計算中我一般取布里淵區(qū)半徑的 1%~5% 作為偏移量并在圓上取 8~12 個點這樣既能看到完整的偏振旋轉(zhuǎn)又不會踩進零輻射區(qū)。4. 空間能帶值提取把能帶曲面、BIC 位置與偏振分布對應(yīng)起來4.1 把能帶鋪成空間曲面空間能帶值指的是哪一組數(shù)據(jù)“空間能帶值”這個說法在光子晶體語境里容易產(chǎn)生歧義。它不是指實空間坐標上每個點的能帶頻率那是聲子晶體局域態(tài)密度的處理方式在光子晶體 BIC 計算里它通常指動量空間里的能帶取值把第一布里淵區(qū)當作坐標平面每個點 (kx, ky) 對應(yīng)某個模式的一個本征頻率把這一組頻率鋪成曲面就是空間能帶值。BIC 就藏在這張曲面上位置往往在能帶交叉形成的錐形點附近。構(gòu)建這張曲面時不要只沿 Γ-X-M 路徑算一維能帶那樣會漏掉 BIC 所在的方向。正確做法是在二維動量網(wǎng)格上做掃描kx 和 ky 都取第一布里淵區(qū)一半的范圍網(wǎng)格數(shù)從 21×21 起步。BIC 本身是一個零維奇點網(wǎng)格太粗會直接把它跳過去所以第一輪粗掃只是定位可疑區(qū)域第二輪要在可疑區(qū)域附近做局部加密掃描。4.2 從這里判斷 BIC 位置能帶轉(zhuǎn)折點與 Q 值不發(fā)散的模式把能帶曲面畫出來后BIC 候選點有幾個可辨識的特征頻率等高線圍繞某個點出現(xiàn)渦旋狀分布能帶曲面在該點附近呈現(xiàn)尖銳的漏斗形狀沿任意方向穿過該點頻率連續(xù)但模式的遠場輻射突然消失。下面這段代碼把從求解器導(dǎo)出的頻率表填充成能帶曲面再做可視化。import numpy as np import matplotlib.pyplot as plt def build_band_surface(kx_grid, ky_grid, band_index, freq_table): # kx_grid, ky_grid: 由 np.meshgrid 生成的動量坐標網(wǎng)格 # freq_table: 求解器導(dǎo)出的數(shù)組形狀為 (n_k, n_band) # band_index: 要繪制的能帶序號從 0 開始計數(shù) surface np.zeros_like(kx_grid) nky, nkx kx_grid.shape for i in range(nky): for j in range(nkx): k_idx i * nkx j surface[i, j] freq_table[k_idx, band_index] return surface # 示例在 21x21 網(wǎng)格上計算第三條能帶 kx np.linspace(-0.5, 0.5, 21) ky np.linspace(-0.5, 0.5, 21) KX, KY np.meshgrid(kx, ky) # freq_table 由能帶求解器輸出每個 k 點取若干條能帶頻率 # surface build_band_surface(KX, KY, band_index2, freq_tablefreq_table) plt.contourf(KX, KY, surface, levels40, cmapviridis) plt.colorbar(labelnormalized frequency) plt.xlabel(r$k_x$ (units of $2\pi/a$)) plt.ylabel(r$k_y$ (units of $2\pi/a$))注意band_index的取值要和第二、三章的能帶順序保持一致。TM 模式下第三條能帶通常是最容易出 BIC 的位置因為前兩條能帶往往落在空氣光錐之下輻射通道本身就不存在。判斷曲面上的漏斗點是不是真 BIC還要回到 2.3 節(jié)的 Q 值復(fù)核單靠頻率曲面只能確定“能帶在這里發(fā)生了交叉”不能確定“交叉點不輻射”。4.3 熱圖輔助把偏振方位角疊到能帶曲面上空間能帶值除了頻率還可以把遠場偏振方位角作為第二層信息疊加上去。做法是在能帶曲面的等高線圖上把每個 k 點的偏振方位角畫成短線短線方向代表該動量的線偏振方向長度代表 S0 的相對強度。BIC 附近會看到清晰的旋渦狀短線分布繞奇點走一圈短線方向旋轉(zhuǎn) 2π 的整數(shù)倍。這個旋渦結(jié)構(gòu)直接對應(yīng) BIC 的拓撲荷比單純看能帶交叉更可靠。疊加圖還能幫你判斷 BIC 的對稱性類型。如果旋渦旋轉(zhuǎn)方向相反說明拓撲荷符號不同對應(yīng)的偏振構(gòu)型也不同。這類信息在做超表面偏振調(diào)控時非常關(guān)鍵目標偏振轉(zhuǎn)換器件的功能本質(zhì)上就是讓遠場偏振沿特定路徑繞 BIC 演化。繪圖時要注意角度顯示的范圍方位角是 0 到 π 的周期量直接用原始值畫會看到劇烈跳變建議先對角度做 unwrap 處理再映射到顏色或短線方向。5. 常見問題排查BIC 偏振與能帶計算最容易翻車的四個坑5.1 BIC 點上算偏振結(jié)果全是噪聲現(xiàn)象把 k 點正好設(shè)在 BIC 位置調(diào)用far_field_stokes后 S1、S2 在每次運行之間完全不一致方位角隨機跳動。原因BIC 的遠場輻射嚴格為零求解器算出來的只是浮點噪聲。Stokes 參數(shù)的 S1、S2 都由遠場復(fù)電場乘積構(gòu)成真實信號趨近零時噪聲就成了主導(dǎo)項。解決永遠不要在嚴格奇點上算偏振。繞 BIC 取一個半徑為 Δk 的小圓Δk 取布里淵區(qū)半徑的 1% 到 5%圓上相鄰點之間夾角均勻。設(shè)定幾個不同半徑各算一圈如果方位角旋轉(zhuǎn)角一致說明結(jié)果收斂如果旋轉(zhuǎn)角隨半徑變化就要加密動量網(wǎng)格可能是附近還有其他奇點干擾。5.2 能帶曲面出現(xiàn)多余的折疊分支現(xiàn)象空間能帶曲面圖上有大量重復(fù)的等高線條紋看起來像是多了一整條能帶。原因掃描動量網(wǎng)格時沒有使用不可約布里淵區(qū)而是把整個方形布里淵區(qū)掃了一遍。晶格對稱性導(dǎo)致某些能帶在多個區(qū)域內(nèi)重復(fù)出現(xiàn)疊加在一起就成了偽造的多分支結(jié)構(gòu)。解決先確認結(jié)構(gòu)屬于哪個點群方形晶格通常取不可約布里淵區(qū)的三角區(qū)域掃描即可。如果做全布里淵區(qū)掃描后處理時按點群對稱性把等效點合并或者只在等高線圖里保留頻率最低的那條分支。這個坑在第四章代碼里特別明顯我當時第一次畫出“八爪魚”狀的曲面排查了半天才發(fā)現(xiàn)是掃描范圍選錯了。5.3 Q 值上不去總是封頂在某個有限值現(xiàn)象FDTD 算出的 Q 值隨分辨率從 16 升到 64漲到幾千后就不動了再怎么加密網(wǎng)格也上不去。原因結(jié)構(gòu)本身是準 BIC不是嚴格 BIC或者介質(zhì)柱的圓形邊界被直角網(wǎng)格離散后引入了非物理的對稱性破缺。網(wǎng)格一粗圓柱輪廓變成鋸齒狀旋轉(zhuǎn)對稱性被破壞本應(yīng)禁止的輻射通道被重新打開。解決先檢查結(jié)構(gòu)是否滿足嚴格對稱條件。圓柱中心要落在網(wǎng)格元胞中心半徑在網(wǎng)格上的離散誤差用亞像素平滑處理MEEP 里可以開啟 subpixel 默認平均。如果已經(jīng)做了這些Q 值還封頂說明這個幾何參數(shù)點本身就不是 BIC 條件沿能帶方向掃幾個相鄰參數(shù)讓 r 或折射率微調(diào)再重新看 Q 值趨勢。5.4 近場到遠場頻譜出現(xiàn)非物理的旁瓣現(xiàn)象near_to_far_spectrum得到的遠場角度譜里目標方向附近出現(xiàn)一串規(guī)則的偽峰間距相等。原因近場監(jiān)測面寬度不是晶格周期的整數(shù)倍。監(jiān)測面被邊界截斷后FFT 把邊界的跳變卷進頻譜產(chǎn)生 sinc 振鈴。偽峰間距恰好對應(yīng)截斷窗口的尺寸。解決把監(jiān)測面寬度精確設(shè)置為整數(shù)個周期。如果結(jié)構(gòu)本身有缺陷或需要局部激勵先在監(jiān)測面加一個窗函數(shù)做平滑再用 FFT。加窗會犧牲一點角分辨率但能換來干凈的頻譜。防止這類問題的習(xí)慣做法是在設(shè)計仿真時就記錄監(jiān)測面的起點坐標保證它正好落在周期邊界上。提示以上四個坑都會在“先算能帶、再算遠場”的流程里暴露。每改一個參數(shù)先記錄網(wǎng)格分辨率和 k 點偏移量不然翻車后根本不知道是幾何參數(shù)變了還是數(shù)值設(shè)置變了。6. 兩個驗證技巧對稱性自檢與繞 BIC 追蹤偏振旋轉(zhuǎn)6.1 用對稱性做第一次自檢拿到能帶曲面和遠場偏振結(jié)果后先別急著往下做器件設(shè)計。跑一個對稱性檢查把整個結(jié)構(gòu)旋轉(zhuǎn) 90 度重新算一遍能帶和 Stoke驗證 KX、KY 坐標相應(yīng)旋轉(zhuǎn)后頻率沒有變化。如果旋轉(zhuǎn)后的偏振方位角沒有跟著旋轉(zhuǎn)說明投影坐標系定義有問題要回 3.2 節(jié)檢查 e1、e2 的構(gòu)造——這類問題不是物理錯了是坐標系沒對齊。這個檢查只需要多跑一次能帶掃描成本很低卻能把一半的坐標系 bug 擋在門外。6.2 繞 BIC 走一圈追蹤偏振方向驗證拓撲荷繞 BIC 繞一圈看偏振方位角的總旋轉(zhuǎn)量這是 BIC 所有數(shù)值結(jié)果里最有說服力的一致性驗證。拓撲荷由總旋轉(zhuǎn)角除以 2π 得到理論上必須是整數(shù)。如果算出 0.5 或 1.7 這種值說明網(wǎng)格、動量偏移或相位歸一化至少有一項出了問題。下面這個函數(shù)直接用小圓采樣計算繞數(shù)。def winding_number(k_center, dk, n_points12, get_field): # k_center: BIC 所在 (kx, ky) # dk: 繞圈半徑 # get_field: 輸入 (kx, ky)返回該方向遠場復(fù)電場 psi_list [] for idx in range(n_points): angle 2.0 * np.pi * idx / n_points kx k_center[0] dk * np.cos(angle) ky k_center[1] dk * np.sin(angle) E get_field(kx, ky) S0, S1, S2, S3, psi far_field_stokes(E, [kx, ky, 0.0]) psi_list.append(psi) # unwrap 防止方位角在 ±pi/2 處跳變 unwrapped np.unwrap(2.0 * np.array(psi_list)) total (unwrapped[-1] - unwrapped[0]) / 2.0 return int(round(total / (2.0 * np.pi)))這個小函數(shù)現(xiàn)在是我所有 BIC 計算的標配。每次換結(jié)構(gòu)參數(shù)、換網(wǎng)格密度我都會先用它驗證一遍繞數(shù)再去做投影器件設(shè)計。如果繞數(shù)不是整數(shù)我會回到第一章的建模流程重新檢查不會繼續(xù)往下算。代價只是多花十分鐘但能省下一個晚上重算的時間——這是我在這類計算里最實在的血淚經(jīng)驗希望幫到你。本文還有配套的精品資源點擊獲取