光力场中粒子运动轨迹的计算

谢海磊 杜骏杰

摘要： 在采用米氏散射理论严格计算微纳粒子所受光力的基础上 ， 研究了基于欧拉-理查森算法计算粒子运动轨迹的问题. 相比欧拉算法和欧拉-克罗默算法 ， 欧拉-理查森算法精度更高且收敛速度更快 ， 是非常适合描绘粒子运动轨迹的方法.纳米粒子在周期性保守光力场中的运动轨迹与物理分析完全吻合 ， 验证了该方法的有效性和稳定性.给出的计算方法 ， 可用于更高效地研究光学微操控中胶体粒子和生物大分子的囚禁、输运、分类 ， 以及宏观粒子的冷却等.

关键词：欧拉-理查森算法;  光力场;  粒子运动轨迹

中图分类号： O436.2    文献标志码： A    DOI： 10.3969/j.issn.1000-5641.2022.02.012

Calculation of particle motion trajectories in optical force fields

XIE Hailei，  DU Junjie

（School of Physics and Electronic Science， East China Normal University， Shanghai  200241 China）

Abstract： In this paper， the motion trajectory of micro-nanoparticles is calculated based on the Euler- Richardson algorithm after the optical force exerted on the particles is determined using Mie scattering theory. The Euler-Richardson algorithm has better calculation accuracy and faster convergence speed than the Euler algorithm and the Euler-Kromer algorithm， and thus is an appropriate approach to describe the trajectory of particles. Hence， the motion trajectory of a nanoparticle in a periodic conservative optical force field is calculated based on the Euler-Kromer algorithm; the results confirm consistency with the physical analysis， further verifying the effectiveness and stability of the approach. The calculation method shown in this paper provides a high-efficiency approach to study optical trapping， transport， sorting of colloidal particles， and biological macromolecules as well as the cooling of macroscopic particles in optical micro-manipulation.

Keywords： Euler-Richardson algorithm;  optical force field;  particle motion trajectory

0  引言

400年前开普勒将彗星尾的形成解释为光对物质的机械效应. 随着麦克斯韦电磁理论的发展 ， 20世纪初 Lebedew[1]、Nichols 和 Hull[2]建立了光压的概念. 激光产生后 ， 光学微操控即借用光力俘获和控制原子、分子、生物与胶体粒子等中性粒子成为可能 ， 其也成为光学研究中的重要方向. 最早的光学微操控实验由 Ashkin 于1970年完成[3-4] ， 该项工作实现了利用单光束或双光束牵引微粒子 ， 同时还得到了两个反向传播的光场俘获粒子的结果 ， 这就是最早的用于俘获粒子的光阱. 随后 ， 1986年出现了单光束梯度阱 ， 它利用了强聚焦光束提供的梯度力[5-10] ， 捕获和移动尺寸在几十纳米到几十微米范围内的粒子[11-12] ， 该梯度阱后来被称为“光镊”. 光镊目前已成为生物学[13-15]、物理、化学和软凝聚态物理领域的重要研究工具 ， 它能够在没有物理接触的情况下俘获和移动微纳粒子以及活体细胞.除此之外 ， 光力还可以实现对微粒子的筛选和分类. 高效地分离细胞和其他生物分子 ， 是生物学和医学研究中的關键技术. 例如 ， 特定细胞亚群的分离对癌症、自身免疫疾病和遗传疾病细胞疗法的进步至关重要， 而分离和研究肿瘤干细胞群是提高对疾病的认知和探求新的治疗方案的关键.

过去的研究中 ， 为了实现对粒子的不同操控 ， 人们研究了各种不同的光场对粒子的作用：具有环形结构特征的贝塞尔光束被用于探讨布朗动力学[16]中的有趣问题;较小的物体借助热运动才能迁移到光束中心 ， 而较大的物体由于能感受到光束的包络 ， 可以更快地迁移到光束中心 ， 这项研究成果可用于无源光学分选;涡旋光束[17]与粒子的相互作用则是研究光的自旋和轨道角动量， 以及自旋轨道相互作用的理想工具[18-23];倏逝波是近场 ， 只能穿透介质表面很小的距离 ， 从而可以利用倏逝波在样品内部不受干扰的情况下 ， 在其表面移动粒子 ， 由于近场光学不受衍射极限的限制 ， 相比传统光镊 ， 倏逝波可以在更小的空间范围内俘获粒子;光阱阵列可以用来研究存在相互作用的胶体粒子系统 ， 该系统常用于模拟热力学系统;光力作用下胶体粒子间的相互作用等效于热力学系统中的微观相互作用 ， 优点是现在有条件测量和控制该微观相互作用 ， 同时也能够模拟系综的宏观热力学行为[24-25].

研究以上提到的复杂光力场对粒子的操控 ， 都需要首先清楚粒子在其中的运动情况. 因此 ， 粒子运动轨迹的正确描绘是光学微操控研究中至关重要的问题. 例如， 了解光镊中被囚禁粒子的运动轨迹， 有利于理解原子、分子的冷却机制;而光分选中， 粒子轨迹直接决定了分选的结果.光场调制下胶体粒子系统的运动轨迹模拟对众多物理现象的理解非常有帮助 ， 如受到驱动的电荷密度波[26]、二维电子气中的电子能态[27]、原子在晶体表面的迁移[28]、化学反应动力学和 II 型超导体[29]中的通量流动等.

基于此， 本文选择微纳胶体粒子作为研究对象， 运用欧拉-理查森算法计算粒子在光力场中的运动轨迹 ， 并给出判定该算法所得结果收敛性的方法 ， 以保证计算结果是粒子通过光场时运动情况的稳定解.物体所受的合力恒定不变或变化规律可以用某一明确的解析表达式描述时 ， 无论粒子的运动轨迹是直线还是较复杂的曲线 ， 都可以通过严格的解析方法求解出来. 但对于力的解析表达式未知、力场不均匀等情形 ， 粒子运动情况就非常复杂了 ， 无法通过解析方法得到粒子的运动轨迹 ， 从而不得不发展相应的数值方法来求解.本文给出的方法不仅仅适用于光力场 ， 同时也适用于任意大小的粒子在任意力场中的运动情况.

1  理论方法

这里介绍一种任意尺寸的粒子在任意光力场中运动轨迹的计算方法 ， 并给出判断每一步长下收敛的条件.在任意入射光场照射下 ， 作用在粒子上的时间平均电磁力可以利用麦克斯韦应力张量的面积分表示 ， 公式为

其中， ⟨⟩是时间平均的麦克斯韦应力张量 ， 其表达式为

其中， ε0、µ0分别是背景介质的介电常数和磁导率; Re表示实部;上标*表示复共轭;电场和磁场都是总场 ， 即入射场的 i （ i）和散射场的 s （ s）的总和 ， 即

（3）

式（1）中， 计算电磁力的积分应在粒子的外表面 S 上进行， ！n 和dσ分别表示表面 S 上单位法向矢量和面积元. 如果背景介质是介电常数ε0和磁导率µ0的无吸收损耗的介质 ， 则积分可在包围粒子的任何表面上进行. 该方法可用于计算任意尺寸粒子在任意光场中受到的光力.

粒子受到的光力随空间位置变化而形成光力场. 光力场通常较复杂 ， 因而粒子在光力场中运动轨迹的严格计算是一个难题. 常用的计算粒子运动轨迹的方法是欧拉算法 ， 该方法假定速度和加速度在时间步长∆t 内不发生显著变化.因此， 为了得到需要的数值解， 时间步长∆t 必须足够小.但是如果∆t 太小 ， 会出现一些数值计算中的其他问题. 比如 ， 时间步长∆t 越小意味着计算中迭代的次数越多 ， 随着迭代次数的增加 ， 有限精度的浮点数的截断误差会逐渐累积 ， 最终数值结果会变得不准确;此外 ， 迭代次数越多 ， 计算机程序完成计算所需的时间就越长. 除了这些问题外 ， 欧拉算法对于许多系统都是不稳定的 ， 这意味着误差以指数级累积 ， 数值解会很快变得不准确.因此 ， 发展更精确和稳定的数值算法非常必要.本文将基于改进的欧拉算法—欧拉-理查森算法得到粒子运动轨迹.

简单的欧拉算法表述公式为

其中 ， vn 、sn 、an 分别为tn 时刻的速度、位置和加速度 ， vn+1、sn+1分别为tn+1即tn +∆t 时刻的速度和位置.公式（4）、公式（5）假定速度和加速度在时间步长∆t 内无显著变化 ， 即在计算位移时认为粒子在∆t 时间内做匀速直线运动 ， 而在计算速度时又认为粒子做匀加速直线运动 ， 且假定初始时刻（即 tn 时刻）的加速度和速度是∆t 时间内粒子的速度和加速度. 简单欧拉算法在振荡系统里误差非常明显 .后来人们又发展了欧拉-克罗默算法， 该算法可提高振荡系统里的计算精度.欧拉-克罗默算法和简单欧拉算法的区别是 ， 在∆t 时间内粒子的运动速度采用时间段结束时的速度vn+1 ， 而不是开始时的速度vn来计算新位置sn+1 ， 公式为

计算速度的方法与简单欧拉算法无异 ， 公式为

但对于力场随着速度变化的复杂问题 ， 以上两种方法都会带来明显的误差. 人们发现 ， 选择时间步长中间的速度而不是开始或结束时的速度 ， 计算精度能够得到明显提高.基于该思想发展的算法就是欧拉-理查森算法 ， 它对于解决力场依赖速度的问题非常有效;当然该算法也适用于不依赖于速度的力场问题 ， 且计算更加精确有效. 欧拉-理查森算法是利用欧拉算法先求出在中间时刻∆t/2 ， 即 tmid = tn +∆t/2时的位置 smid 和速度 vmid ;然后计算 t = tmid 时的力 F（smid， vmid， tmid）和加速度 amid ;最后将vmid 、amid 代入欧拉算法中求出 tn+1时刻 ， 即tn+1 = tn +∆t 时的位置 sn+1和速度vn+1. 欧拉-理查森算法的计算方法和步骤（参见文献[30]）数学上可以表示为

以及

尽管在得到每个时间步长的最后结果前 ， 需要两次使用欧拉算法 ， 但是因为该算法中的时间步长可以设置得比欧拉算法和欧拉-克罗默算法中的更长一点 ， 所以欧拉-理查森算法的计算速度仍然优于欧拉算法 ， 并且还具有比欧拉算法和欧拉-克罗默算法更高的精度.

以上的欧拉-理查森算法的数学表示针对的仅仅是一维情况 ， 对于二维或三维情况 ， 在每一维度上分别单独应用上述公式 ， 即可计算任何矢量力场中粒子的运动轨迹.但是选定步长之后 ， 该算法的计算结果是否收敛 ， 这点在数值计算过程中尤为重要.下面以二维情况（在平面直角坐标系中描述粒子位置）举例介绍判断收敛的方法. 判断的方法如图 1所示.

图1所示的流程图围绕步长∆t 的选择和确定展开. 利用欧拉-理查森算法进行计算时 ， 时间步长是非常关键的一个因素：步长太大时 ， 粒子运动的路径上某一点的力场若发生突变 ， 则计算结果会出现较大误差 ， 甚至不会收敛;而步长太小时 ， 计算步数太多 ， 计算量过大 ， 计算时间太长. 本文采取的办法是針对每一步长都进行判断 ， 只有满足设定的判断条件才认为该步长的选取是合适的 ， 或者说计算结果已经收敛 ， 程序才可以开始下一步长的计算. 具体判断收敛的过程如下.

第一步 ， 根据事先设定的位移条件 ， 初步确定一个时间步长. 本文是以入射光波长的1/64作为时间步长∆t 内粒子可移动的位移最大值. 给定一个初始时间步长 ， 由欧拉-理查森算法计算这个步长中间时刻即tmid = tn +∆t/2和tmid = tn+1时的速度、位置和加速度 ， 并计算出该步长内粒子运动的位移.此位移需小于预设的位移最大值;如果大于位移最大值 ， 表示选取的时间步长过大.此时 ， 将步长减半后重复上面的计算 ， 对新步长下算出的位移再进行判断 ， 重复该步骤直到满足位移条件为止 ， 从而初步得到一个时间步长.

第二步 ， 在第一步依据位移条件给出时间步长的基础上 ， 进一步判断该步长是否满足加速度条件 ， 由加速度条件最终确定合适的时间步长. 加速度条件需要区分不同的情形 ， 分别进行判断. 这里的加速度是指 tn 和 tn+1时刻处 （在判断中分别称之为初位置和末位置）的加速度. 首先 ， 如果 x 和 y 方向的初始和最终加速度均为0， 表示这个时间步长内粒子几乎不受力而做匀速直线运动 ， 算法给出的是严格的结果 ， 可以进入下一时间步长的计算. 此处要记录这一步长的末速度、末位置 ， 以备下一步计算使用. 其次 ， 若初位置的加速度均为0， 但末位置的加速度不为0， 则判断末位置相对初位置在 x 方向和 y 方向上的速度相对误差（即相对应方向的速度变化量除以初速度） ， 这两个方向上的相对误差均小于2‰， 则步长符合条件;若初位置的加速度在一个方向上不为0 而在另一个方向上为0， 则计算该方向上的加速度相对误差和另一方向上的速度相对误差 ， 这两个相对误差都小于2‰， 则步长符合条件;若初位置的加速度均不为0， 则计算 x 方向和 y 方向上的加速度相对误差 ， 这两个相对误差都小于1‰， 则步长符合条件. 上述几种情形下 ， 若不满足判断条件 ， 则将时间步长除以5 （可根据具体问题需要 ， 选择减小的倍数）作为新的时间步长 ， 重复上述计算和判断 ， 直到得到满足条件的时间步长. 最后 ， 将最终得到的时间步长代入欧拉-理查森算法 ， 计算出粒子运动的位移、速度和加速度 ， 再进入下一时间步长的计算.

2  周期性光力场中粒子运动轨迹的模拟

光学微操控中涉及粒子在光力场中运动的时候 ， 由于实际的光力场不一定有解析表达式 ， 粒子运动轨迹很难通过解析办法得到.此时 ， 本文前面介绍的数值方法就可以发挥独特的作用了.下面以二维周期性光力场中的粒子运动为例 ， 模拟粒子在该力场中的运动规律 ， 进而验证上述算法的有效性.

该二维光力场由4 束相同频率的线性极化平面波（每束平面波的强度为104 W/cm2）沿着两两互相垂直的方向入射形成 ， 如图2所示. 图2 中 ， 4束波的电场都在 x-y 面内振动. 为了保证这4束波全部相干 ， 要求它们在坐标原点的相位相同.在这样4束相干波产生的光场中 ， 无论多大的粒子受到的光力都只有梯度力而没有散射力[31] ， 即图2给出的是一个纯净的周期性梯度光力场. 虽然形式上这一光力场与传统光晶格没有差别 ， 但由于光晶格形成时对4束光的相位没有要求 ， 一般的光晶格不是纯净的梯度光力场 ， 仅仅是对瑞利粒子的一个近似的梯度光力场. 处在这样纯净的梯度光力场中的微纳粒子 ， 由于仅受到保守的梯度力 ， 将会被囚禁在某一个格子内. 图2给出了1.596µm ×1.064µm 的方形区域内光力场的空间分布图 ， 梯度力应是矢量.图2中显示的仅是梯度力的强度即jFg j ， 力的大小大约在皮牛（pN）量级 ， 图中的蓝色区域梯度力较弱 ， 红色区域较强.

本文选择密度为1.05 g/cm3的聚苯乙烯胶体球粒子作为研究对象 ， 且可根据需要 ， 通过改变它的半径大小来改变粒子尺寸.图2标示了半径为0.5µm 的粒子在该光力场中的初始位置（0， 0.45）， 即粒子处在图中正上方红色圆圈所包围的格子里.粒子被无初速释放后 ， 由于初始时刻粒子处在 x =0 这条线上 ， 粒子仅受到 y 方向的梯度力 ， 而不会受到 x 方向的作用力 ， 粒子将从初始位置 y =0.45开始沿 y 轴正方向运动 ， 到达约（0， 0.60）位置的时候速度降为0， 再开始反向运动 ， 到达初始位置后重复上面的振荡运动过程. 总地来看 ， 粒子将在这个格子里、在该区间内沿 y 轴做简谐振动.

接下来利用欧拉-理查森算法 ， 计算粒子释放后的运动轨迹 ， 看轨迹是否如分析的那样. 图3（a）–图3（d）模拟了粒子释放后不同时间段内的轨迹. 注意图3中各图横坐标 x 方向的数值仅有10–16 µm.为了清楚地显示轨迹 ， 图3中 x 方向的长度被放大显示了 ， x 方向上的数值非常非常小. 因此粒子在 x 方向上基本没有运动 ， 只沿着 y 方向做往复运动. 图3 各图也同时显示了该运动在 y 轴的最远点是（0， 0.60）.初始点和最远点都被标记在了图2中. 由图2 可以看到 ， 这两个点在该格子中基本对称 ， 符合粒子在对称格子里运动应该对称的特点.无论是往复运动还是运动的最远点 ， 都与本文在模拟前的物理分析吻合 ， 这意味着本文对运动轨迹的数值模拟方法是非常有效的.从图3（a）到图3（d）， 模拟的时间段越来越长 ， 但粒子做往复运动的规律并没有变化 ， 这说明本文的方法并没有时间积累误差 ， 完全可以用来模拟长时间内粒子的运动轨迹.

3  结论

本文在米氏散射理论的基础上 ， 利用欧拉-理查森算法发展了计算微纳粒子在光力场中运动轨迹的方法.该方法可用于计算原子、分子、胶体粒子、病毒细菌等生物大分子或生物细胞在光力场中的运动. 在文中 ， 我们以球形聚苯乙烯粒子在周期性保守光力场中的运动为例 ， 验证了该方法的有效性.虽然文中仅模拟了真空背景下粒子在光力场中的运动 ， 但是对于涉及阻尼和随机布朗运动 ， 甚至粒子同时受到力矩的作用等情况 ， 该方法仍然适用.在这些情形下 ， 计算粒子在各个方向所受的合力将变得复杂 ， 但一旦得到了所有合力 ， 每個步长下的位移就可以通过欧拉-理查森算法求出.本文介绍的欧拉-理查森算法适用于布朗运动、胶体科学和生物医学中的光学微操控等诸多研究领域.

[參考文献]

[1]LEBEDEW P. Untersuchungen über die druckkrfte des lichtes [J]. Annalen Der Physik， 1901， 311（11）：433-458.

[2]NICHOLS E F， HULL G F. A preliminary communication on the pressure of heat and light radiation [J]. Physical Review， 1901， 13（5）：307-320.

[3]ASHKIN A. Acceleration and trapping of particles by radiation pressure [J]. Physical Review Letters， 1970， 24（4）：156-159.

[4]ASHKIN A. Stability of optical levitation by radiation pressure [J]. Applied Physics Letters， 1974， 24（12）：586-588.

[5]GHISLAIN L P， SWITZ N A， WEBB W W. Measurement of small forces using an optical trap [J]. Review of Scientific Instruments， 1994， 65（9）：2762-2768.

[6]ROHRBACH A， STELZER E H K. Trapping forces， force constants， and potential depths for dielectric spheres in the presence of spherical aberrations [J]. Applied Optics， 2002， 41（13）：2494-2507.

[7]LITVINOV R I， SHUMAN H， BENNETT J S， et al. Binding strength and activation state of single fibrinogen-integrin pairs on living cells [J]. Proceedings of the National Academy of Sciences， 2002， 99（11）：7426-7431.

[8]GITTES F， SCHMIDT C F. Signals and noise in micromechanical measurements [J]. Methods in Cell Biology， 1998， 55（55）：129-156.

[9]GITTES F， SCHMIDT C F. Interference model for back-focal-plane displacement detection in optical tweezers [J]. Optics Letters， 1998， 23（1）：7-9.

[10]PRALLE A， PRUMMER M， FLORIN E L， et al. Three-dimensional high-resolution particle tracking for optical tweezers by forward scattered light [J]. Microscopy Research and Technique， 2015， 44（5）：378-386.

[11]SVOBODA K， BLOCK S M. Optical trapping of metallic Rayleigh particles [J]. Optics Letters， 1994， 19（13）：930-932.

[12]KE P C， GU M. Characterization of trapping force on metallic mie particles [J]. Applied Optics， 1999， 38（1）：160-167.

[13]ASHKIN A. History of optical trapping and manipulation of small-neutral particle， atoms， and molecules [J]. IEEE Journal of Selected Topics in Quantum Electronics， 2000， 6（6）：841-856.

[14]SVOBODA  K，  MITRA  P  P，  BLOCK  S  M. Fluctuation  analysis  of motor  protein  movement  and  single  enzyme  kinetics [J].Proceedings of the National Academy of Sciences of the United States of America， 1994， 91（25）：11782-11786.

[15] BUSTAMANTE C， SMITH S B， LIPHARDT J， et al. Single-molecule studies of DNA mechanics [J]. Current Opinion in StructuralBiology， 2000， 10（3）：279-285.

[16] TATARKOVA S A， SIBBETT W， DHOLAKIA K. Brownian particle in an optical potential of the washboard type [J]. PhysicalReview Letters， 2003， 91（3）：038101.

[17] GAHAGAN K T， SWARTZLANDER G A. Optical vortex trapping of particles [J]. Optics Letters， 1996， 21（11）：827-829.

[18] O’NEIL A T， PADGETT M J. Three-dimensional optical confinement of micron-sized metal particles and the decoupling of the spinand orbital angular momentum within an optical spanner [J]. Optics Communications， 2000， 185（1/2/3）：139-143.

[19] HE H， FRIESE M E J， HECKENBERG N R， et al. Direct observation of transfer of angular momentum to absorptive particles from alaser beam with a phase singularity [J]. Physical Review Letters， 1995， 75（5）：826-829.

[20] ALLEN L， PADGETT M J， BABIKER M. IV the orbital angular momentum of light [J]. Progress in Optics， 1999， 39：291-372.

[21] O’NEIL A T， MACVICAR I， ALLEN L， et al. Intrinsic and extrinsic nature of the orbital angular momentum of a light beam [J].Physical Review Letters， 2002， 88：053601.

[22] CURTIS J E， GRIER D G. Structure of optical vortices [J]. Physical Review Letters， 2003， 90（13）：133901.

[23] SIMPSON N B， DHOLAKIA K， ALLEN L， et al. Mechanical equivalence of spin and orbital angular momentum of light： An opticalspanner [J]. Optics Letters， 1997， 22（1）：52-54.

[24] KORDA P T， SPALDING G C， GRIER D G. Evolution of a colloidal critical state in an optical pinning potential landscape [J].Physical Review B， 2002， 66：024504.

[25] MANGOLD K， LEIDERER P， BECHINGER C. Phase transitions of colloidal monolayers in periodic pinning arrays [J]. PhysicalReview Letters， 2003， 90（15）：158302.

[26] BROWN S E， MOZURKEWICH G， GRÜNER G. Subharmonic shapiro steps and devil’s-staircase behavior in driven charge-density-wave systems [J]. Physical Review Letters， 1984， 52（25）：2277-2280.

[27] WIERSIG J， AHN K H. Devil’s staircase in magnetoresistance of a periodic array of scatterers [J]. Physical Review Letters， 2001，87（2）：026803.

[28] PIERRE-LOUIS O， HAFTEL M. Oscillatory driving of crystal surfaces： A route to controlled pattern formation [J]. Physical ReviewLetters， 2001， 87（4）：048701.

[29] REICHHARDT C， NORI F. Phase locking， devil’s staircases， farey trees， and arnold tongues in driven vortex lattices with periodicpinning [J]. Physical Review Letters， 1998， 82（2）：414-417.

[30] GOULD H， TOBOCHNIK J CHRISTIAN W.計算机模拟方法在物理学中的应用[M].影印版.3版.北京， 高等教育出版社， 2006.

[31] DU J J， YUEN C H， LI X， et al. Tailoring optical gradient force and optical scattering and absorption force [J]. Scientific Reports，2017， 7（1）：18042.

（责任编辑：李艺）