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