简介围绕通用自旋动量锁定光学力整理的这份项目资料面向光学、光子学及微纳操控方向的学习者与研究人员聚焦光与物质相互作用中自旋动量锁定与光学力建模这一较抽象主题。压缩包共13个文件以4个Jupyter Notebook交互脚本为核心辅以6张结果示意图、1份说明文档、1个HTML页面及1个YAML配置文件整体约465KB便于边运行代码边对照图示理解。内容涉及自旋与轨道动量关系、辐射压力与梯度力并延伸到表面等离激元场、倏逝场及HE11模式受力随角动量变化等案例可通过Python数值模拟与可视化复现光场分布和受力分析。已有190人学习适合希望从理论概念推进到可执行案例并借此掌握建模思路、绘图方法与结果解读的读者也可作为相关课程设计或课题入门的参考。1. 通用自旋动量锁定光学力一份 Jupyter Notebook 能算到哪一步“通用自旋动量锁定光学力”这几个字摆在一起很容易被当成纯理论名词一眼扫过去。我第一次打开这份 Jupyter Notebook 时也是这个反应直到把一颗亚波长粒子放进全内反射产生的倏逝场里看着横向力的符号跟着传播方向一起翻面才意识到它更像一套能逐格执行、能改参数、能出图的力场计算流程而不是公式陈列。它做的事很具体给定一束场倏逝场、聚焦光束、涡旋光都行和一颗小粒子把梯度力、辐射压、以及和自旋锁定相关的横向力分别算出来让“自旋—动量锁定”这件事在数值上看得见。做光镊、手性分选、近场操控的人可以拿它当自己的场型试验台。2. 自旋动量锁定的物理骨架力为什么必须拆成三份2.1 倏逝场里那点看不见的横向自旋先把场摆出来。界面在 z0上方介质折射率 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(x,z) 复振幅形状 (Nz, Nx, 3)换光束类型力计算代码一行不动粒子模块半径、介电常数、背景折射率、波长复极化率 α标量或 3×3 张量换材料、换手性场代码一行不动力模块E、数值梯度 dE、α三分量力场只有物理公式变了才动它这个切分看着平平无奇实际用起来差别很大。想看贝塞尔光束的横向力只要写一个新的场函数返回同样形状的数组想比较金纳米球和聚苯乙烯球的区别只要换 α 的输入参数想验证某个解析解直接把解析梯度喂进力模块跳过数值差分。真正的复用发生在接口层而不是在公式层。我自己改这类脚本时有个习惯先把场模块单独跑一遍画 |E|² 和 s_y 两张图确认自旋的符号和传播方向对得上再去接力模块。如果跳过这一步直接看力一旦符号反了你根本分不清是场写错了还是力的约定错了——两个错误叠在一起有时候结果反而“看起来合理”这是最难查的一类。3. Notebook 拆解环境、单元格顺序与第一次跑通3.1 依赖与环境比你想象的少但版本要对这类计算 Notebook 的依赖很少核心就三件numpy 做数组和复数运算scipy 用来做可选的插值和特殊函数matplotlib 负责出图。真正会影响结果的是 numpy 的大版本——1.20 之后 np.gradient 的 axis 参数行为和更早版本有细微差别如果你拿到的是老 Notebook在新环境里跑出来边界有异常条纹先怀疑这个。我一般这么建环境# 建议单独建环境避免和已有科学计算栈打架 conda create -n spinforce python3.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, 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 节那套推导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 # 从界面上方一格开始避开 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|² 的灰度底图上是这套 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 都要用相对介电常数即折射率的平方不要混用绝对介电常数。水的 n1.33 对应 εm≈1.77玻璃 n1.5 对应 εp2.25这两个数代进去分子是 0.48分母是 5.79α0 是正的粒子被拉向高场强区。如果换成金属粒子εp 在光频下是负数且带大虚部α 的实部可能变号梯度力的方向就反了——这不是 bug是物理。这个式子只在瑞利区成立判据是 ka ≪ 1工程上一般要求 ka 0.3。633 nm 波长、n21.33 时k≈1.32×10⁷ m⁻¹a25 nm 对应 ka≈0.33已经踩在边界上了。把半径调到 50 nmka≈0.66严格说不该用这个公式但很多人照用不误结果就是力和半径的标度关系从 a³ 变成了 a⁴ 甚至更乱。真要算大粒子得换米氏级数那是另一个量级的工程量。4.2 场参数临界角、衰减常数与穿透深度场这一侧的自由度其实只有三个波长、入射角、折射率对比。波长决定 k0折射率对比决定临界角入射角决定 κ 和穿透深度。三者不是独立的——给定材料和波长入射角是唯一能连续调节 κ 的旋钮。参数符号本组取值依赖关系越界后果真空波长λ0633 nm独立改变 k0全盘缩放棱镜折射率n11.5材料决定影响临界角背景折射率n21.33材料决定影响临界角与 α入射角θ70°必须 θc低于临界角无倏逝场衰减常数κ由 θ 算出κ sqrt(kx²−n2²k0²)越大穿透越浅边界误差越敏感穿透深度δ1/κ反比于 κ小于网格步长时数值发散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, 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
