簡介圍繞通用自旋動量鎖定光學(xué)力整理的這份項目資料面向光學(xué)、光子學(xué)及微納操控方向的學(xué)習(xí)者與研究人員聚焦光與物質(zhì)相互作用中自旋動量鎖定與光學(xué)力建模這一較抽象主題。壓縮包共13個文件以4個Jupyter Notebook交互腳本為核心輔以6張結(jié)果示意圖、1份說明文檔、1個HTML頁面及1個YAML配置文件整體約465KB便于邊運(yùn)行代碼邊對照圖示理解。內(nèi)容涉及自旋與軌道動量關(guān)系、輻射壓力與梯度力并延伸到表面等離激元場、倏逝場及HE11模式受力隨角動量變化等案例可通過Python數(shù)值模擬與可視化復(fù)現(xiàn)光場分布和受力分析。已有190人學(xué)習(xí)適合希望從理論概念推進(jìn)到可執(zhí)行案例并借此掌握建模思路、繪圖方法與結(jié)果解讀的讀者也可作為相關(guān)課程設(shè)計或課題入門的參考。1. 通用自旋動量鎖定光學(xué)力一份 Jupyter Notebook 能算到哪一步“通用自旋動量鎖定光學(xué)力”這幾個字?jǐn)[在一起很容易被當(dāng)成純理論名詞一眼掃過去。我第一次打開這份 Jupyter Notebook 時也是這個反應(yīng)直到把一顆亞波長粒子放進(jìn)全內(nèi)反射產(chǎn)生的倏逝場里看著橫向力的符號跟著傳播方向一起翻面才意識到它更像一套能逐格執(zhí)行、能改參數(shù)、能出圖的力場計算流程而不是公式陳列。它做的事很具體給定一束場倏逝場、聚焦光束、渦旋光都行和一顆小粒子把梯度力、輻射壓、以及和自旋鎖定相關(guān)的橫向力分別算出來讓“自旋—動量鎖定”這件事在數(shù)值上看得見。做光鑷、手性分選、近場操控的人可以拿它當(dāng)自己的場型試驗臺。2. 自旋動量鎖定的物理骨架力為什么必須拆成三份2.1 倏逝場里那點(diǎn)看不見的橫向自旋先把場擺出來。界面在 z0上方介質(zhì)折射率 n1下方 n2 n1入射角 θ 超過臨界角 θc arcsin(n2/n1)界面下方就只剩一層貼著界面的倏逝場。它的波矢在傳播方向是實的在垂直方向是虛的k_x n1·k0·sinθ 是實數(shù)衰減常數(shù) κ sqrt(k_x2 ? n22k02) 決定穿透深度 δ 1/κ。場沿 x 振蕩、沿 z 指數(shù)衰減這是任何光學(xué)教材都會畫的圖。真正有意思的是偏振。TE 偏振下電場只有 E_y 一個分量E 和 E* 平行叉乘為零一點(diǎn)自旋都沒有。換成 TM 偏振情況立刻變了由 ?·(εE)0 可以解出 E_z 與 E_x 差一個 90° 相位寫成復(fù)振幅就是 E_z i(k_x/κ)·E_x。這個純虛系數(shù)是整個自旋故事的源頭——兩個正交分量相位差四分之一周期場在局部就是橢圓的。把復(fù)振幅代進(jìn)自旋角動量密度 s (ε0/4ω)·Im(E* × E)取 E (1, 0, i·k_x/κ)·E0·exp(i·k_x·x ? κz)算出來只有 y 分量不為零s_y ?ε0·k_x / (2ω·κ)負(fù)號來自我這里選的時間因子 e^(?iωt) 和 E_z 的相位寫法換個約定符號就翻過來不用糾結(jié)。關(guān)鍵在結(jié)構(gòu)s_y 的大小只由 k_x 和 κ 決定符號完全跟著 k_x 走。傳播方向反了橫向自旋立刻反號傳播方向定了自旋方向就定了沒有第二個自由度。這就是“自旋—動量鎖定”四個字的全部含義——它是幾何約束不是某種需要特殊激發(fā)才能出現(xiàn)的效應(yīng)。這也解釋了為什么很多人第一次算橫向力會算出零拿線偏振平面波去測當(dāng)然什么都測不到只有倏逝場這種“波矢有虛部”的結(jié)構(gòu)才會強(qiáng)制出來一個垂直于傳播平面、又鎖定于傳播方向的自旋分量。2.2 力分解梯度項、殘余項與自旋相關(guān)項場有了接下來是粒子。在瑞利近似下粒子半徑遠(yuǎn)小于波長粒子被當(dāng)成一個點(diǎn)偶極整個粒子的響應(yīng)壓縮成一個復(fù)極化率 α。時間平均光學(xué)力寫成?F? (1/2)·Re[ α · Σ_i E_i · ?_i E* ]這里的 E 是復(fù)振幅時間因子取 e^(?iωt)下標(biāo) i 遍歷 x、y、z 三個方向。這個式子是整個 Notebook 的計算核心一行 einsum 就能實現(xiàn)但它的符號約定極其敏感——后面避坑那一章會專門講如果整體符號反了無論怎么調(diào)參數(shù)都救不回來。把實部項單獨(dú)拎出來得到的是大家最熟悉的梯度力F_grad (1/4)·Re(α)·?|E|2注意這里是 1/4 而不是 1/2因為 |E|2 用的是峰值平方時間平均之后要除以 2兩個因子疊在一起就成了 1/4。這個系數(shù)是新手最容易翻車的地方一旦你習(xí)慣性地把場換成 RMS 值1/4 就得改成 1/2改錯一處力的量級直接差一倍而且圖上看起來“形狀很對”非常難發(fā)現(xiàn)。剩下那部分我習(xí)慣叫殘余項F_res ?F? ? F_grad。它包含輻射壓也包含自旋相關(guān)的貢獻(xiàn)。在均勻平面波里殘余項就沿傳播方向大小是 (1/2)·k·Im(α)·|E|2可以拿來做符號自檢。在倏逝場里殘余項多出一個沿 x 的分量這個分量才是這套 Notebook 真正想讓人看到的東西。橫向力可以再分兩支。一支是推著粒子沿傳播方向走的輻射壓所有粒子都有正比于 Im(α)方向只由能流方向決定。另一支是自旋鎖定力只有當(dāng)粒子對左旋場和右旋場的響應(yīng)不一樣時才會顯著——手性粒子、帶磁響應(yīng)的粒子、或者被人為賦予圓極化率的模型都屬于這一類。它的方向由 s_y 的符號決定大小和場自旋的空間旋度掛鉤。因為 s_y 鎖定于 k_x這一支力的方向也就鎖定于傳播方向翻轉(zhuǎn)傳播方向粒子被推向相反的一側(cè)。這就是“自旋動量鎖定光學(xué)力”在實驗上真正能被觀測到的形式不是力的大小變了是力的方向被場的手性單方面指定了。Notebook 里這兩支共用同一個 α 入口手性那一支通過一個額外的圓極化率參數(shù)打開。這樣安排的好處是只研究普通介電小球的人把那個參數(shù)置零就行研究手性分選的人把它打開兩套物理在同一份代碼里跑不用維護(hù)兩個分支。2.3 “通用”體現(xiàn)在哪場函數(shù)與極化率解耦很多光學(xué)力腳本寫死在某一個場型上換個光束就得重寫一遍梯度代碼。這份 Notebook 的結(jié)構(gòu)上做了一個很樸素的切分把“通用”落在了三個互相不知道對方存在的模塊上模塊輸入輸出換掉它意味著什么場模塊坐標(biāo)網(wǎng)格、波長、角度等E(x,z) 復(fù)振幅形狀 (Nz, Nx, 3)換光束類型力計算代碼一行不動粒子模塊半徑、介電常數(shù)、背景折射率、波長復(fù)極化率 α標(biāo)量或 3×3 張量換材料、換手性場代碼一行不動力模塊E、數(shù)值梯度 dE、α三分量力場只有物理公式變了才動它這個切分看著平平無奇實際用起來差別很大。想看貝塞爾光束的橫向力只要寫一個新的場函數(shù)返回同樣形狀的數(shù)組想比較金納米球和聚苯乙烯球的區(qū)別只要換 α 的輸入?yún)?shù)想驗證某個解析解直接把解析梯度喂進(jìn)力模塊跳過數(shù)值差分。真正的復(fù)用發(fā)生在接口層而不是在公式層。我自己改這類腳本時有個習(xí)慣先把場模塊單獨(dú)跑一遍畫 |E|2 和 s_y 兩張圖確認(rèn)自旋的符號和傳播方向?qū)Φ蒙显偃ソ恿δK。如果跳過這一步直接看力一旦符號反了你根本分不清是場寫錯了還是力的約定錯了——兩個錯誤疊在一起有時候結(jié)果反而“看起來合理”這是最難查的一類。3. Notebook 拆解環(huán)境、單元格順序與第一次跑通3.1 依賴與環(huán)境比你想象的少但版本要對這類計算 Notebook 的依賴很少核心就三件numpy 做數(shù)組和復(fù)數(shù)運(yùn)算scipy 用來做可選的插值和特殊函數(shù)matplotlib 負(fù)責(zé)出圖。真正會影響結(jié)果的是 numpy 的大版本——1.20 之后 np.gradient 的 axis 參數(shù)行為和更早版本有細(xì)微差別如果你拿到的是老 Notebook在新環(huán)境里跑出來邊界有異常條紋先懷疑這個。我一般這么建環(huán)境# 建議單獨(dú)建環(huán)境避免和已有科學(xué)計算棧打架 conda create -n spinforce python3.11 -y conda activate spinforce pip install numpy scipy matplotlib jupyterlab這里不裝 numba 也能跑網(wǎng)格在 512×512 這個量級時純 numpy 的 einsum 單次計算在百毫秒級交互式調(diào)參完全夠用。等到你要掃幾十組參數(shù)做收斂性測試再考慮用 numba 加速梯度那一段但別一上來就加——加速會掩蓋掉公式層面的錯誤你連中間變量都看不見。啟動就用 jupyter lab在瀏覽器里逐格執(zhí)行。不建議用 jupyter notebook 老界面這份資源里有不少需要反復(fù)回看中間數(shù)組形狀的格子lab 的多標(biāo)簽和變量檢查器省事很多。3.2 三個入口函數(shù)與單元格的組織順序一份組織得好的計算 Notebook單元格順序應(yīng)該和物理推導(dǎo)順序一致先定義網(wǎng)格再構(gòu)造函數(shù)再給參數(shù)最后才算和畫。亂序執(zhí)行是這類 Notebook 最常見的“玄學(xué) bug”來源——上一格改了 α 但沒重跑力計算圖上顯示的還是舊結(jié)果你還以為是物理效應(yīng)。三個入口函數(shù)的骨架大致是這樣import numpy as np def build_field(x, z, kx, kappa, E01.0): 構(gòu)造 TM 偏振倏逝場的復(fù)振幅。 x, z: 一維坐標(biāo)數(shù)組由 meshgrid 生成 kx: 傳播方向波矢分量實數(shù) kappa: 垂直方向衰減常數(shù)實數(shù)單位 1/m 返回: E, 形狀 (Nz, Nx, 3)最后一維是 x/y/z 分量 phase np.exp(1j * kx * x - kappa * z) Ex E0 * phase Ez 1j * (kx / kappa) * E0 * phase # 由 div(eps*E)0 解出 Ey np.zeros_like(Ex) return np.stack([Ex, Ey, Ez], axis-1) def numeric_gradient(E, dx, dz): 對每個分量做中心差分返回 dE形狀 (Nz, Nx, 3, 3)。 索引約定: dE[..., i, j] d(E_j)/d(x_i)i0 對應(yīng) xi2 對應(yīng) z。 dE np.zeros(E.shape (3,), dtypecomplex) for j in range(3): gz, gx np.gradient(E[..., j], dz, dx, axis(0, 1)) dE[..., 0, j] gx dE[..., 2, j] gz return dE def optical_force(E, dE, alpha): 時間平均光學(xué)力約定時間因子 e^(-i*omega*t)。 F_j 0.5 * Re[ alpha * sum_i E_i * conj(dE_i_j) ] return 0.5 * np.real(alpha * np.einsum(...i,...ij-...j, E, np.conj(dE)))第一個函數(shù)對應(yīng) 2.1 節(jié)那套推導(dǎo)Ez 前面的 i 系數(shù)就是自旋的來源改掉它自旋立刻消失。第二個函數(shù)把梯度寫成了 (3,3) 張量看著有點(diǎn)繞但這樣做的好處是力模塊里不用關(guān)心坐標(biāo)系順序einsum 一行就能把 Σ_i E_i·?_i 收掉。第三個函數(shù)就是 2.2 節(jié)那個總力公式注意 conj 加在 dE 上而不是 E 上——這是 e^(?iωt) 約定下的正確寫法加反了散射力的方向會整體倒過來。提示如果你的場是解析可導(dǎo)的完全可以把 numeric_gradient 換成解析表達(dá)式。數(shù)值差分在邊界只有單邊精度倏逝場在 z 方向衰減很快邊界處的梯度誤差會被放大這一點(diǎn)在第 4 章的收斂性測試?yán)飼炕?.3 第一次跑通單粒子在倏逝場中的三分量力參數(shù)給一組最常見的波長 633 nm棱鏡折射率 1.5水 1.33入射角 70°粒子半徑 25 nm介電常數(shù)取玻璃的 2.25。先算臨界角確認(rèn)在超臨界區(qū)再生成網(wǎng)格再算力。lam 633e-9 n1, n2 1.5, 1.33 theta np.deg2rad(70.0) k0 2 * np.pi / lam theta_c np.arcsin(n2 / n1) # 約 62.5 度確認(rèn) 70 度超臨界 kx n1 * k0 * np.sin(theta) kappa np.sqrt(kx**2 - (n2 * k0)**2) print(f臨界角 {np.rad2deg(theta_c):.1f} deg, 穿透深度 {1/kappa*1e9:.1f} nm) # 網(wǎng)格: x 方向多給一點(diǎn)z 方向貼著界面切 Nx, Nz 400, 200 dx lam / (20 * n2) x (np.arange(Nx) - Nx // 2) * dx z np.arange(Nz) * dx dx # 從界面上方一格開始避開 z0 X, Z np.meshgrid(x, z, indexingxy) E build_field(X, Z, kx, kappa) dE numeric_gradient(E, dx, dx) a, eps_p, eps_m 25e-9, 2.25, n2**2 eps0 8.854e-12 alpha0 4 * np.pi * eps0 * a**3 * (eps_p - eps_m) / (eps_p 2 * eps_m) alpha alpha0 / (1 - 1j * (n2*k0)**3 * alpha0 / (6 * np.pi * eps0)) F optical_force(E, dE, alpha) iz, ix Nz // 2, Nx // 2 print(Fx, Fy, Fz , F[iz, ix, :])跑完這一格你會看到 F_y 嚴(yán)格為零——場在 y 方向均勻沒有任何梯度這是對的。F_z 是負(fù)的粒子被往界面方向拉這是倏逝場梯度力的正常表現(xiàn)衰減常數(shù)越大拉得越狠。F_x 不為零方向沿傳播方向由輻射壓貢獻(xiàn)普通介電小球沒有手性自旋鎖定那一支在這里不顯著。想看自旋鎖定得打開手性。把粒子模塊換成帶圓極化率的版本或者干脆用雙各向異性極化率再跑一次F_x 的大小會隨手性參數(shù)的符號變化甚至可能反向。這時候你再把 kx 取負(fù)等價于反向傳播會看到 F_x 跟著翻面——這一個動作就能把“鎖定”兩個字坐實。把 F 的三個分量畫成三張圖疊加在 |E|2 的灰度底圖上是這套 Notebook 最有說服力的一頁輸出。梯度力的圖沿 z 方向分層輻射壓的圖沿 x 方向單調(diào)自旋鎖定力的圖則在界面附近出現(xiàn)一條明顯的條帶越靠近界面越強(qiáng)因為自旋密度正比于 1/κ而 κ 越大入射角越靠近 90°條帶越薄。4. 參數(shù)怎么定波長、半徑、極化率與網(wǎng)格4.1 復(fù)極化率從 Clausius-Mossotti 到輻射修正極化率是整份代碼里唯一描述“粒子是什么”的量取錯了后面所有圖都沒有意義。起點(diǎn)是 Clausius-Mossotti 形式α0 4πε0·a3·(εp ? εm)/(εp 2εm)這里 a 是半徑εp 是粒子介電常數(shù)εm 是背景介質(zhì)介電常數(shù)。這個式子在靜場下是準(zhǔn)的但光頻下必須加輻射修正否則實部和虛部都不對尤其是虛部——它決定吸收直接決定輻射壓的大小α α0 / (1 ? i·k3·α0/(6πε0))k n2·k0 是背景介質(zhì)里的波數(shù)。修正項里的 k3 意味著波長越短、修正越強(qiáng)在 400 nm 附近這項能改變百分之十幾不做修正力的大小會和文獻(xiàn)值對不上而你又找不到哪里錯了。參數(shù)取值上有個容易忽略的點(diǎn)εp 和 εm 都要用相對介電常數(shù)即折射率的平方不要混用絕對介電常數(shù)。水的 n1.33 對應(yīng) εm≈1.77玻璃 n1.5 對應(yīng) εp2.25這兩個數(shù)代進(jìn)去分子是 0.48分母是 5.79α0 是正的粒子被拉向高場強(qiáng)區(qū)。如果換成金屬粒子εp 在光頻下是負(fù)數(shù)且?guī)Т筇摬喀?的實部可能變號梯度力的方向就反了——這不是 bug是物理。這個式子只在瑞利區(qū)成立判據(jù)是 ka ? 1工程上一般要求 ka 0.3。633 nm 波長、n21.33 時k≈1.32×10? m?1a25 nm 對應(yīng) ka≈0.33已經(jīng)踩在邊界上了。把半徑調(diào)到 50 nmka≈0.66嚴(yán)格說不該用這個公式但很多人照用不誤結(jié)果就是力和半徑的標(biāo)度關(guān)系從 a3 變成了 a? 甚至更亂。真要算大粒子得換米氏級數(shù)那是另一個量級的工程量。4.2 場參數(shù)臨界角、衰減常數(shù)與穿透深度場這一側(cè)的自由度其實只有三個波長、入射角、折射率對比。波長決定 k0折射率對比決定臨界角入射角決定 κ 和穿透深度。三者不是獨(dú)立的——給定材料和波長入射角是唯一能連續(xù)調(diào)節(jié) κ 的旋鈕。參數(shù)符號本組取值依賴關(guān)系越界后果真空波長λ0633 nm獨(dú)立改變 k0全盤縮放棱鏡折射率n11.5材料決定影響臨界角背景折射率n21.33材料決定影響臨界角與 α入射角θ70°必須 θc低于臨界角無倏逝場衰減常數(shù)κ由 θ 算出κ sqrt(kx2?n22k02)越大穿透越淺邊界誤差越敏感穿透深度δ1/κ反比于 κ小于網(wǎng)格步長時數(shù)值發(fā)散70° 時算出來的穿透深度在 100 nm 量級比波長小得多這意味著場在 z 方向不到兩個波長就衰減到可以忽略。好處是 z 方向網(wǎng)格不用給太寬壞處是邊界效應(yīng)被放大如果 z 方向只切了 200 格、每格 20 nm 出頭總深度才 4 μm看起來夠但最上面幾格的場已經(jīng)小到接近浮點(diǎn)誤差那里的梯度基本是噪聲。入射角往 90° 推κ 迅速變大穿透深度掉到幾十納米這時候網(wǎng)格步長必須跟著縮小否則連一個衰減長度都采樣不到幾個點(diǎn)。反過來入射角剛過臨界角時 κ 趨近于零穿透深度可以到微米量級場幾乎不衰減自旋密度 s_y 正比于 1/κ 會變得很大但梯度力很小兩者的相對權(quán)重會發(fā)生劇烈變化。這兩個極端都值得單獨(dú)掃一遍比在中間取一個“好看”的角度信息量大得多。4.3 網(wǎng)格與收斂步長、padding 與自檢數(shù)值差分最怕兩件事步長太大梯度的截斷誤差壓不住邊界沒有 padding單邊差分把假的梯度引進(jìn)結(jié)果。這兩件事在倏逝場里都會被放大因為場的二階導(dǎo)數(shù)本身就很大。判斷步長夠不夠不要靠肉眼做一個循環(huán)# 收斂性自檢逐步加密網(wǎng)格看中心點(diǎn)力是否穩(wěn)定 base None for npl in [8, 12, 16, 24, 32, 48]: dx lam / (npl * n2) Nx, Nz int(6 * lam / dx), int(1.5 * lam / dx) x (np.arange(Nx) - Nx // 2) * dx z np.arange(Nz) * dx dx X, Z np.meshgrid(x, z, indexingxy) E build p a hrefhttps://download.csdn.net/download/weixin_42164702/15355077 stylecolor:#ec7500;font-size:14px; 本文還有配套的精品資源點(diǎn)擊獲取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p