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