简介：围绕通用自旋动量锁定光学力整理的这份项目资料，面向光学、光子学及微纳操控方向的学习者与研究人员，聚焦光与物质相互作用中自旋动量锁定与光学力建模这一较抽象主题。压缩包共13个文件，以4个Jupyter Notebook交互脚本为核心，辅以6张结果示意图、1份说明文档、1个HTML页面及1个YAML配置文件，整体约465KB，便于边运行代码边对照图示理解。内容涉及自旋与轨道动量关系、辐射压力与梯度力，并延伸到表面等离激元场、倏逝场及HE11模式受力随角动量变化等案例，可通过Python数值模拟与可视化复现光场分布和受力分析。已有190人学习，适合希望从理论概念推进到可执行案例，并借此掌握建模思路、绘图方法与结果解读的读者，也可作为相关课程设计或课题入门的参考。

1. 通用自旋动量锁定光学力：一份 Jupyter Notebook 能算到哪一步

“通用自旋动量锁定光学力”这几个字摆在一起，很容易被当成纯理论名词一眼扫过去。我第一次打开这份 Jupyter Notebook 时也是这个反应，直到把一颗亚波长粒子放进全内反射产生的倏逝场里，看着横向力的符号跟着传播方向一起翻面，才意识到它更像一套能逐格执行、能改参数、能出图的力场计算流程，而不是公式陈列。它做的事很具体：给定一束场（倏逝场、聚焦光束、涡旋光都行）和一颗小粒子，把梯度力、辐射压、以及和自旋锁定相关的横向力分别算出来，让“自旋—动量锁定”这件事在数值上看得见。做光镊、手性分选、近场操控的人，可以拿它当自己的场型试验台。

2. 自旋动量锁定的物理骨架：力为什么必须拆成三份

2.1 倏逝场里那点看不见的横向自旋

先把场摆出来。界面在 z=0，上方介质折射率 n1，下方 n2 < n1，入射角 θ 超过临界角 θc = arcsin(n2/n1)，界面下方就只剩一层贴着界面的倏逝场。它的波矢在传播方向是实的，在垂直方向是虚的：k_x = n1·k0·sinθ 是实数，衰减常数 κ = sqrt(k_x² − n2²k0²) 决定穿透深度 δ = 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 走。传播方向反了，横向自旋立刻反号；传播方向定了，自旋方向就定了，没有第二个自由度。这就是“自旋—动量锁定”四个字的全部含义——它是几何约束，不是某种需要特殊激发才能出现的效应。

这也解释了为什么很多人第一次算横向力会算出零：拿线偏振平面波去测，当然什么都测不到；只有倏逝场这种“波矢有虚部”的结构，才会强制出来一个垂直于传播平面、又锁定于传播方向的自旋分量。

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|²

注意这里是 1/4 而不是 1/2，因为 |E|² 用的是峰值平方，时间平均之后要除以 2，两个因子叠在一起就成了 1/4。这个系数是新手最容易翻车的地方：一旦你习惯性地把场换成 RMS 值，1/4 就得改成 1/2，改错一处，力的量级直接差一倍，而且图上看起来“形状很对”，非常难发现。

剩下那部分我习惯叫残余项：F_res = ⟨F⟩ − F_grad。它包含辐射压，也包含自旋相关的贡献。在均匀平面波里，残余项就沿传播方向，大小是 (1/2)·k·Im(α)·|E|²，可以拿来做符号自检。在倏逝场里，残余项多出一个沿 x 的分量，这个分量才是这套 Notebook 真正想让人看到的东西。

横向力可以再分两支。一支是推着粒子沿传播方向走的辐射压，所有粒子都有，正比于 Im(α)，方向只由能流方向决定。另一支是自旋锁定力，只有当粒子对左旋场和右旋场的响应不一样时才会显著——手性粒子、带磁响应的粒子、或者被人为赋予圆极化率的模型都属于这一类。它的方向由 s_y 的符号决定，大小和场自旋的空间旋度挂钩。因为 s_y 锁定于 k_x，这一支力的方向也就锁定于传播方向，翻转传播方向，粒子被推向相反的一侧。这就是“自旋动量锁定光学力”在实验上真正能被观测到的形式：不是力的大小变了，是力的方向被场的手性单方面指定了。

Notebook 里这两支共用同一个 α 入口，手性那一支通过一个额外的圆极化率参数打开。这样安排的好处是：只研究普通介电小球的人把那个参数置零就行，研究手性分选的人把它打开，两套物理在同一份代码里跑，不用维护两个分支。

2.3 “通用”体现在哪：场函数与极化率解耦

很多光学力脚本写死在某一个场型上，换个光束就得重写一遍梯度代码。这份 Notebook 的结构上做了一个很朴素的切分，把“通用”落在了三个互相不知道对方存在的模块上：

这个切分看着平平无奇，实际用起来差别很大。想看贝塞尔光束的横向力，只要写一个新的场函数返回同样形状的数组；想比较金纳米球和聚苯乙烯球的区别，只要换 α 的输入参数；想验证某个解析解，直接把解析梯度喂进力模块，跳过数值差分。真正的复用发生在接口层，而不是在公式层。

我自己改这类脚本时有个习惯：先把场模块单独跑一遍，画 |E|² 和 s_y 两张图，确认自旋的符号和传播方向对得上，再去接力模块。如果跳过这一步直接看力，一旦符号反了，你根本分不清是场写错了还是力的约定错了——两个错误叠在一起，有时候结果反而“看起来合理”，这是最难查的一类。

3. Notebook 拆解：环境、单元格顺序与第一次跑通

3.1 依赖与环境：比你想象的少，但版本要对

这类计算 Notebook 的依赖很少，核心就三件：numpy 做数组和复数运算，scipy 用来做可选的插值和特殊函数，matplotlib 负责出图。真正会影响结果的是 numpy 的大版本——1.20 之后 np.gradient 的 axis 参数行为和更早版本有细微差别，如果你拿到的是老 Notebook，在新环境里跑出来边界有异常条纹，先怀疑这个。

我一般这么建环境：

# 建议单独建环境，避免和已有科学计算栈打架 conda create -n spinforce python=3.11 -y conda activate spinforce pip install numpy scipy matplotlib jupyterlab

这里不装 numba 也能跑，网格在 512×512 这个量级时，纯 numpy 的 einsum 单次计算在百毫秒级，交互式调参完全够用。等到你要扫几十组参数做收敛性测试，再考虑用 numba 加速梯度那一段，但别一上来就加——加速会掩盖掉公式层面的错误，你连中间变量都看不见。

启动就用 jupyter lab，在浏览器里逐格执行。不建议用 jupyter notebook 老界面，这份资源里有不少需要反复回看中间数组形状的格子，lab 的多标签和变量检查器省事很多。

3.2 三个入口函数与单元格的组织顺序

一份组织得好的计算 Notebook，单元格顺序应该和物理推导顺序一致：先定义网格，再构造函数，再给参数，最后才算和画。乱序执行是这类 Notebook 最常见的“玄学 bug”来源——上一格改了 α 但没重跑力计算，图上显示的还是旧结果，你还以为是物理效应。

三个入口函数的骨架大致是这样：

import numpy as np def build_field(x, z, kx, kappa, E0=1.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)，i=0 对应 x，i=2 对应 z。 """ dE = np.zeros(E.shape + (3,), dtype=complex) 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 节那套推导，Ez 前面的 i 系数就是自旋的来源，改掉它自旋立刻消失。第二个函数把梯度写成了 (3,3) 张量，看着有点绕，但这样做的好处是力模块里不用关心坐标系顺序，einsum 一行就能把 Σ_i E_i·∂_i 收掉。第三个函数就是 2.2 节那个总力公式，注意 conj 加在 dE 上而不是 E 上——这是 e^(−iωt) 约定下的正确写法，加反了散射力的方向会整体倒过来。

提示：如果你的场是解析可导的，完全可以把 numeric_gradient 换成解析表达式。数值差分在边界只有单边精度，倏逝场在 z 方向衰减很快，边界处的梯度误差会被放大，这一点在第 4 章的收敛性测试里会量化。

3.3 第一次跑通：单粒子在倏逝场中的三分量力

参数给一组最常见的：波长 633 nm，棱镜折射率 1.5，水 1.33，入射角 70°，粒子半径 25 nm，介电常数取玻璃的 2.25。先算临界角确认在超临界区，再生成网格，再算力。

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 # 从界面上方一格开始，避开 z=0 X, Z = np.meshgrid(x, z, indexing='xy') 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|² 的灰度底图上，是这套 Notebook 最有说服力的一页输出。梯度力的图沿 z 方向分层，辐射压的图沿 x 方向单调，自旋锁定力的图则在界面附近出现一条明显的条带，越靠近界面越强，因为自旋密度正比于 1/κ，而 κ 越大（入射角越靠近 90°）条带越薄。

4. 参数怎么定：波长、半径、极化率与网格

4.1 复极化率：从 Clausius-Mossotti 到辐射修正

极化率是整份代码里唯一描述“粒子是什么”的量，取错了，后面所有图都没有意义。起点是 Clausius-Mossotti 形式：

α0 = 4πε0·a³·(εp − εm)/(εp + 2εm)

这里 a 是半径，εp 是粒子介电常数，εm 是背景介质介电常数。这个式子在静场下是准的，但光频下必须加辐射修正，否则实部和虚部都不对，尤其是虚部——它决定吸收，直接决定辐射压的大小：

α = α0 / (1 − i·k³·α0/(6πε0))

k = n2·k0 是背景介质里的波数。修正项里的 k³ 意味着波长越短、修正越强，在 400 nm 附近这项能改变百分之十几，不做修正，力的大小会和文献值对不上，而你又找不到哪里错了。

参数取值上有个容易忽略的点：εp 和 εm 都要用相对介电常数（即折射率的平方），不要混用绝对介电常数。水的 n=1.33 对应 εm≈1.77，玻璃 n=1.5 对应 εp=2.25，这两个数代进去，分子是 0.48，分母是 5.79，α0 是正的，粒子被拉向高场强区。如果换成金属粒子，εp 在光频下是负数且带大虚部，α 的实部可能变号，梯度力的方向就反了——这不是 bug，是物理。

这个式子只在瑞利区成立，判据是 ka ≪ 1，工程上一般要求 ka < 0.3。633 nm 波长、n2=1.33 时，k≈1.32×10⁷ m⁻¹，a=25 nm 对应 ka≈0.33，已经踩在边界上了。把半径调到 50 nm，ka≈0.66，严格说不该用这个公式，但很多人照用不误，结果就是力和半径的标度关系从 a³ 变成了 a⁴ 甚至更乱。真要算大粒子，得换米氏级数，那是另一个量级的工程量。

4.2 场参数：临界角、衰减常数与穿透深度

场这一侧的自由度其实只有三个：波长、入射角、折射率对比。波长决定 k0，折射率对比决定临界角，入射角决定 κ 和穿透深度。三者不是独立的——给定材料和波长，入射角是唯一能连续调节 κ 的旋钮。

70° 时算出来的穿透深度在 100 nm 量级，比波长小得多，这意味着场在 z 方向不到两个波长就衰减到可以忽略。好处是 z 方向网格不用给太宽，坏处是边界效应被放大：如果 z 方向只切了 200 格、每格 20 nm 出头，总深度才 4 μm，看起来够，但最上面几格的场已经小到接近浮点误差，那里的梯度基本是噪声。

入射角往 90° 推，κ 迅速变大，穿透深度掉到几十纳米，这时候网格步长必须跟着缩小，否则连一个衰减长度都采样不到几个点。反过来，入射角刚过临界角时 κ 趋近于零，穿透深度可以到微米量级，场几乎不衰减，自旋密度 s_y 正比于 1/κ 会变得很大，但梯度力很小，两者的相对权重会发生剧烈变化。这两个极端都值得单独扫一遍，比在中间取一个“好看”的角度信息量大得多。

4.3 网格与收敛：步长、padding 与自检

数值差分最怕两件事：步长太大，梯度的截断误差压不住；边界没有 padding，单边差分把假的梯度引进结果。这两件事在倏逝场里都会被放大，因为场的二阶导数本身就很大。

判断步长够不够，不要靠肉眼，做一个循环：

# 收敛性自检：逐步加密网格，看中心点力是否稳定 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, indexing='xy') E = build <p> <a href="https://download.csdn.net/download/weixin_42164702/15355077" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>