简介:压缩包内整合了“通用自旋动量锁定光学力”主题的交互式学习与仿真资源,面向光学工程、光子学等相关专业学习者。包内以 Jupyter Notebook 为核心,共13个文件、总大小465KB,包含4个可运行的 ipynb 脚本(如 SPP_Field、Evanescent Field、HE11_General 等)、6张电场/受力分布图、1份 Markdown 说明文档、1份 HTML 页面及1份 YAML 配置,便于对照图文理解建模过程。已有189人学习下载,适合希望通过数值模拟掌握自旋动量锁定光学力原理的读者。资源通过电磁场仿真展示自旋动量锁定现象,覆盖表面等离激元、倏逝场及 HE11 模式等典型场景,可观察光学力分布并理解自旋-轨道耦合机制。结合源码与结果图,读者能完成从理论模型到数值实验的完整复现,为光镊、微纳操控、量子光学等应用的进一步探索提供基础。
1. 通用自旋动量锁定光学力:从“横向力”到可交互的Jupyter仿真
第一次接触到“通用自旋动量锁定光学力”这个概念时报错,是因为它在光镊、微流控和量子弱测量方向都被人反复引用,但真正给出可运行计算的代码却很少。这份资源提供了一整套基于JupyterNotebook的仿真流程,把光子自旋动量与光力场之间的锁定关系拆成可调参数的数值实验:你可以随意切换光束类型、粒子半径、介质折射率,实时看到法向力、散射力和横向力分量的变化。它解决的是“为什么自旋动量方向变化会导致粒子出现侧向位移”这个问题,也顺带把角度谱法、瑞利散射模型、麦克斯韦应力张量的数值实现全串起来了。适合正在做光力计算、微颗粒操控或者矢量光束传输研究的从业者。
2. 把物理模型落成代码:瑞利粒子受力分解与角谱法选型
2.1 瑞利粒子的受力分解:梯度力、散射力与自旋力
我最早踩的坑是把光学力当成一个整体去算,结果发现不管怎么调网格,横向分量的方向总是和文献对不上。后来在某个平台开发光力模块时,某开发者提醒我:瑞利粒子受力必须拆成三部分才能和“自旋动量锁定”挂上钩。
对于半径远小于波长的球形粒子,通常采用的分解方式是:
- 梯度力:正比于电场强度梯度的实部,方向指向光强极大处。
- 散射力:正比于坡印廷矢量时间平均值,方向沿光传播方向。
- 自旋力:来自电磁场自旋动量密度的旋度分量,会产生垂直于传播方向的横向力。
三者的相对大小随光束偏振态变化。尤其要注意的是,圆偏振涡旋光束下,自旋力方向和涡旋拓扑荷的符号强相关。这个“方向锁定”关系就是标题里“自旋动量锁定”的物理核心。我在代码里把三种力分开输出,就是为了能直接验证这条锁定关系。
2.2 用角谱法计算聚焦光场:为什么选角谱而不是FDTD
在对比了几种场计算方法之后,我最终选择角谱法作为这份Notebook的主算法。原因有三条:
第一,光学力仿真中关心的是焦平面附近的场分布,角谱法把光场分解成平面波叠加,可以直接通过二维傅里叶变换计算,不用搭建三维立体网格,内存开销比FDTD低一个量级。
第二,角谱法天然满足矢量光束的偏振需求。径向偏振光、角向偏振光、涡旋光都可以通过在频谱面上定义不同的偏振初始相位来生成,不需要额外引入边界条件。
第三,后续计算自旋动量密度时需要用到场的空间导数,角谱法在频域里乘以对应的波矢分量就能得到导数,数值上比有限差分平滑许多。
当然角谱法也有一个使用边界:它适用于折射率变化不剧烈的均匀介质场景。如果粒子周围有金属纳米结构或者强梯度折射率层,我一般会改用FDTD单独验证关键节点。
2.3 Jupyter中的力场计算流程
把受力模型和场算法合到一起,这份Notebook里最核心的代码块是这样的:
def compute_optical_force(E_field, H_field, alpha_polar, k0, n_m, grid_dx): # 计算梯度力 E_sq = np.abs(E_field) ** 2 grad_E_sq = np.gradient(E_sq, grid_dx, axis=(0, 1, 2)) F_grad = (alpha_polar.real / 4.0) * grad_E_sq # 计算散射力(沿坡印廷方向) S = 0.5 * np.real(np.cross(E_field, np.conj(H_field), axis=-1)) F_scatt = (k0 ** 4 * np.abs(alpha_polar) ** 2) / (12.0 * np.pi * n_m ** 2) * S # 自旋力通过自旋动量密度旋度计算 p_spin = compute_spin_momentum(E_field, H_field, k0, grid_dx) F_spin = k0 * np.imag(np.gradient(p_spin, grid_dx, axis=(0, 1, 2))) return F_grad, F_scatt, F_spin这一段的关键点在于三个力都用了不同的物理来源:梯度力取电场模平方的实部梯度,散射力取坡印廷矢量,自旋力取自旋动量密度的旋度。很多人一次性把三个力加起来当总数用,我建议还是分开跑一遍看分量图,否则一旦横向力方向和理论对不上,根本不知道是哪一项算错了。
alpha_polar是粒子的极化率复数,对介质球用 Clausius-Mossotti 公式计算;grid_dx是网格步长,直接影响所有空间梯度的精度,网格设置需要和聚焦数值孔径匹配。这里最容易被忽视的是H_field相位基准,和电场必须使用同一个时间相位约定,否则坡印廷矢量的实部方向会翻转。
3. 自旋动量密度的数值实现:旋度算子与离散化误差控制
3.1 自旋动量密度的标准公式与符号约定
自旋动量密度是理解横向光学力的关键中间量。在均匀介质中,电磁场自旋密度矢量通常写成:
p_spin = (1 / (4 * omega)) * Im[ E* x (nabla x E) + H* x (nabla x H) ]其中 omega 是角频率,E* 和 H* 表示复共轭,x 是叉积。“锁定”这个词在这里的意思是:涡旋光束中自旋密度的符号,决定了横向力的方向;一旦涡旋拓扑荷或者圆偏振手性反转,横向力方向也跟着翻转。这就是为什么在数值模拟中自旋动量的离散化必须精确到符号层面,不能只关心幅值。
实现时要注意,这个公式里用的是“共轭后的场和场的旋度做叉积再取虚部”,不是直接对坡印廷矢量取旋度。二者的数值结果有系统性差异,方向极容易差到相反。
3.2 频域法计算旋度:比有限差分更稳的离散策略
在Jupyter中逐点做有限差分旋度,高频噪声会明显放大。我在这份资源里用的是频域算子方案:
import numpy as np def curl_fourier(field, dx): # field: (Ny, Nx, Nz, 3),最后一个维度是E或H的三个分量 Nx, Ny, Nz = field.shape[:3] kx = np.fft.fftfreq(Nx, dx) * 2.0 * np.pi ky = np.fft.fftfreq(Ny, dx) * 2.0 * np.pi kz = np.fft.fftfreq(Nz, dx) * 2.0 * np.pi KX, KY, KZ = np.meshgrid(kx, ky, kz, indexing='ij') def curl_component(comp_vec): Fx_fft = np.fft.fftn(comp_vec[..., 0]) Fy_fft = np.fft.fftn(comp_vec[..., 1]) Fz_fft = np.fft.fftn(comp_vec[..., 2]) cx = np.fft.ifftn(1j * KY * Fz_fft - 1j * KZ * Fy_fft) cy = np.fft.ifftn(1j * KZ * Fx_fft - 1j * KX * Fz_fft) cz = np.fft.ifftn(1j * KX * Fy_fft - 1j * KY * Fx_fft) return np.stack([cx, cy, cz], axis=-1) return curl_component这里的核心逻辑在频域完成:傅里叶变换后把空间偏导变成乘上对应的波矢分量,再逆变换回来。这样得到的旋度在整个计算区域内保持全局一致的精度,不会因为个别点噪声产生毛刺。
用这个函数时要特别注意网格是和波矢分辨率绑定的:dx越小,fftfreq返回的波矢范围越大,高频成分越容易被引入。所以网格密度要配合低通滤波使用,尤其是在粒子内部折射率突变的地方。
我在实际跑偏的时候还发现一个细节:np.fft.fftn默认把变换放在多维数组的前三维,field的第4维是三个分量,这个维度不能参与变换。不少刚上手的人会把field[..., 0]和整个field混在一起做变换,结果得到的是完全错误的复数场,力方向也随之漂移。
4. 拿起这份Notebook直接改参数:光束、粒子与介质的三方配合
4.1 文件结构与入口函数速览
解压之后,核心文件有四个,全在同一个主目录下:
| 文件 | 作用 | 主要输出 |
|---|---|---|
spin_force_core.py | 角谱生成聚焦场、三种力计算 | 力分量数组 |
spin_momentum.py | 自旋动量密度与旋度计算 | 自旋矢量场 |
demo_vortex_beam.ipynb | 涡旋光束下的横向力演示 | 交互式图表 |
params.yaml | 所有物理参数与网格参数入口 | 配置字典 |
入口函数写在demo_vortex_beam.ipynb的第一个cell里,函数名是run_simulation(params_dict)。只需要保证params_dict里的键名和你改动的yaml字段一致,调用后会把三种力分量写进results/目录下的npz文件。
4.2 参数设置表与推荐起点
在params.yaml里,我一般只动这几个字段:
| 参数名 | 推荐初值 | 作用 | 调整边界 |
|---|---|---|---|
wavelength_nm | 800 | 入射光波长 | 400-1000 |
numerical_aperture | 1.2 | 物镜数值孔径 | 0.6-1.4 |
n_medium | 1.33 | 介质折射率 | 1.0-1.6 |
n_particle | 1.59 | 粒子折射率 | 1.4-2.0 |
particle_radius_nm | 120 | 粒子半径 | 50-250 |
grid_points | 128 | 单轴网格点数 | 64-256 |
grid_spacing_nm | 12 | 网格步长 | 8-20 |
beam_type | lg01 | 光束模式 | gauss/lg01/radial |
特别提醒一点:particle_radius_nm不要超过波长的十分之一除以介质折射率。一旦粒子半径进入米氏散射区间,瑞利分解模型会整体失效,这时候无论如何调网格,横向力曲线都会在零点附近乱抖。
4.3 跑通第一个示例:涡旋光束下的横向光学力
打开demo_vortex_beam.ipynb,逐格执行到标记“Force calculation”的代码块,看到的是下面这段:
import yaml from spin_force_core import run_simulation with open('params.yaml', 'r') as f: params = yaml.safe_load(f) params['beam_type'] = 'lg01' params['topological_charge'] = 1 params['polarization'] = 'circular_right' F_grad, F_scatt, F_spin = run_simulation(params) print("梯度力均值:", np.abs(F_grad).mean()) print("散射力均值:", np.abs(F_scatt).mean()) print("自旋力均值:", np.abs(F_spin).mean())跑通后你会在输出图里看到三个分量的空间分布。如果topological_charge改成1再改成-1,F_spin的横向分量方向会跟着翻转,而梯度力方向保持不变。这个对称性破缺现象,就是自旋动量锁定最直观的体现。
这里参数设置的逻辑是:lg01表示拉盖尔高斯光束径向指数为0、拓扑荷为1,在焦平面上光场中心是暗的;circular_right指定右旋圆偏振入射,让自旋角动量方向和轨道角动量叠加。两者共同作用,产生稳定的横向自旋力。
我建议第一次跑时把三个力度量都打开,只看总力的大小,会掩盖掉自旋力的符号变化。自旋力通常比梯度力小一个量级,只有分开看才能确认它确实存在。
5. 避坑指南:六个高频翻车点与排查经验
5.1 横向力方向翻车:相位约定不统一
现象:把拓扑荷从+1改成-1,横向力方向不变,和理论预期相反。
原因:电场分量里混用了时间相位约定,有的分使用了exp(-i*omega*t),有的用了exp(+i*omega*t),复振幅的虚部符号被翻转。
解决:写死一个全局约定,我强制在compute_optical_force里加入相位标志位,所有场分量统一使用exp(-i*omega*t)约定,再在结果里自检:拓扑荷反号时自旋力必须跟着反号。
5.2 力曲线震荡:网格太粗或波矢混叠
现象:自旋力分布出现棋盘状花纹,数值比相邻格点交替跳变。
原因:网格间距dx超过焦场最高空间频率对应采样的一半,高空间频率成分混叠到低频,梯度算子产生虚假震荡。
解决:用dx < wavelength / (4 * n_medium * numerical_aperture)作为强约束。我的经验是 NA=1.2、波长800nm 时,dx至少要取 12nm,再小到 8nm 会更稳,但内存占用会指数上升。
5.3 归一化不正确:力值整体偏大一截
现象:梯度力和散射力比文献预期大了几十倍,分布形状却正确。
原因:入射功率算成了峰值强度乘以无穷大积分面积,没有限定光束功率为 1mW。
解决:在角谱生成阶段对总功率做归一化,即把电场模平方在焦平面的积分设为1μW。参数表里加上input_power_mw字段,每次改光束模式后都重新跑一遍归一化。
5.4 粒子折射率和介质折射率太接近时极化率虚部被吞掉
现象:散射力分量趋近于零,自旋力被梯度力淹没。
原因:alpha_polar的虚部在折射率差趋近零时也趋近零,散射力在数值上变成微小量减微小量,正负号噪声主导。
解决:把折射率差设置到 0.01 以上,或者改用n_particle = 1.59配合n_medium = 1.0这种空气环境抓典型的散射力曲线。边界上算力时,我习惯把折射率差小于 0.005 的配置直接标记为无效输入。
5.5 计算窗口边界截断导致横向力整体偏离
现象:粒子正好位于计算窗口边缘附近时,横向力方向有时会倒转。
原因:角谱法默认周期延拓,窗口边界处场被周期镜像污染,等效于旁边多了一个假粒子。
解决:在粒子周围至少保留 8 个波长的空白区域,或者计算前先对场做空间窗函数截断。我习惯在run_simulation里打开一个window_padding_ratio参数,默认设为 1.5,也就是把粒子周围的空白区拉伸到粒子区域的1.5倍。
5.6 Jupyter内核重启后结果和缓存对不上
现象:改了某个参数重跑,画出来的图还是旧数据。
原因:Notebook所有cell都依赖内存变量,重启内核后变量还在旧值,重新执行也会把旧结果覆盖写进新路径,但图里缓存没更新。
解决:在保存结果时强制带上params_hash,每次参数变动都会生成新的文件名,避免覆盖旧结果。改完参数先执行一次params_yaml_verify(),确认变更被真正加载。
6. 进阶用法:把单点力扫描改成参数扫描和并行加速
6.1 用ipywidgets做实时交互调参
跑通基础Notebook之后,我第一个加的增强是做实时调参面板。核心是用ipywidgets.interact包装仿真函数,滑块参数直接绑定到params.yaml的字段上:
import ipywidgets as widgets from IPython.display import display def interactive_force(radius_nm, n_particle, charge): params['particle_radius_nm'] = radius_nm params['n_particle'] = n_particle params['topological_charge'] = charge F_results = run_simulation(params, save_result=False) plot_force_components(F_results) widgets.interact( interactive_force, radius_nm=widgets.FloatSlider(min=50, max=250, step=10, value=120), n_particle=widgets.FloatSlider(min=1.4, max=2.0, step=0.01, value=1.59), charge=widgets.IntSlider(min=-2, max=2, step=1, value=1) )这段代码的价值在于它能实时显示横向力方向随拓扑荷翻转的动态过程。跑过一次后你就能直观理解“自旋动量锁定”不是抽象概念,而是参数图上清晰可见的方向反转。
6.2 用numba给自旋动量计算提速
角谱计算和自旋动量旋度都在频域做傅里叶变换,这部分提速空间有限。真正耗时的是三维网格上的逐点操作,我一般会用numba加速:
from numba import jit @jit(nopython=True) def compute_spin_force_numba(p_spin_x, p_spin_y, p_spin_z, dx): Nx, Ny, Nz = p_spin_x.shape F_spin = np.zeros_like(p_spin_x) for i in range(1, Nx-1): for j in range(1, Ny-1): for k in range(1, Nz-1): F_spin[i, j, k] = (p_spin_z[i+1, j, k] - p_spin_z[i-1, j, k] - p_spin_y[i, j+1, k] + p_spin_y[i, j-1, k]) / (2 * dx) return F_spin用@jit装饰后,网格 128³ 的力计算可以从几秒压到几百毫秒。注意numba不支持复数域的全部numpy功能,所以自旋动量密度的计算仍然保留在numpy里,只有最后的实部旋度交给numba。
6.3 批量扫描拓扑荷和粒子半径
最后做成参数扫描的脚本,输出二维热力图:
charges = [-2, -1, 1, 2] radii = np.arange(80, 201, 20) force_heatmap = np.zeros((len(charges), len(radii))) for i, ch in enumerate(charges): for j, r in enumerate(radii): params['topological_charge'] = ch params['particle_radius_nm'] = r _, _, F_spin = run_simulation(params, save_result=False) # 提取横向力在粒子中心的值 center = tuple(np.array(F_spin.shape[:3]) // 2) force_heatmap[i, j] = F_spin[center]这段批量扫描代码帮助我确认了一个已知规律:自旋力横向分量的最大值出现在粒子半径 120-160nm 区间,且拓扑荷为正时横向力偏向正方向。参数扫描比单点调试更有价值的地方在于,它是验证“通用”性的关键手段——不是某单一参数组合下偶然出现的现象,而是在多个参数组合下呈现统一规律,才能叫“锁定”。
我现在的使用习惯是:每次拿到一份新光学力资源,先跑通它的默认和复现结果,再改参数扫描一次,而且会把扫描结果保存成带参数哈希的npz文件。从那以后我每次跑新的仿真案例都强制走一遍这套流程:核对相位约定、检查网格分辨率、确认归一化、最后做参数扫描验证方向律。希望这种排查顺序也能帮你在自旋动量锁定光学力仿真上少走些弯路。
本文还有配套的精品资源,点击获取