光力場中粒子運動軌跡的計算

謝海磊 杜駿杰

摘要： 在采用米氏散射理論嚴格計算微納粒子所受光力的基礎(chǔ)上 ， 研究了基于歐拉-理查森算法計算粒子運動軌跡的問題. 相比歐拉算法和歐拉-克羅默算法 ， 歐拉-理查森算法精度更高且收斂速度更快 ， 是非常適合描繪粒子運動軌跡的方法.納米粒子在周期性保守光力場中的運動軌跡與物理分析完全吻合 ， 驗證了該方法的有效性和穩(wěn)定性.給出的計算方法 ， 可用于更高效地研究光學微操控中膠體粒子和生物大分子的囚禁、輸運、分類 ， 以及宏觀粒子的冷卻等.

關(guān)鍵詞：歐拉-理查森算法;? 光力場;? 粒子運動軌跡

中圖分類號： 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年前開普勒將彗星尾的形成解釋為光對物質(zhì)的機械效應(yīng). 隨著麥克斯韋電磁理論的發(fā)展 ， 20世紀初 Lebedew[1]、Nichols 和 Hull[2]建立了光壓的概念. 激光產(chǎn)生后 ， 光學微操控即借用光力俘獲和控制原子、分子、生物與膠體粒子等中性粒子成為可能 ， 其也成為光學研究中的重要方向. 最早的光學微操控實驗由 Ashkin 于1970年完成[3-4] ， 該項工作實現(xiàn)了利用單光束或雙光束牽引微粒子 ， 同時還得到了兩個反向傳播的光場俘獲粒子的結(jié)果 ， 這就是最早的用于俘獲粒子的光阱. 隨后 ， 1986年出現(xiàn)了單光束梯度阱 ， 它利用了強聚焦光束提供的梯度力[5-10] ， 捕獲和移動尺寸在幾十納米到幾十微米范圍內(nèi)的粒子[11-12] ， 該梯度阱后來被稱為“光鑷”. 光鑷目前已成為生物學[13-15]、物理、化學和軟凝聚態(tài)物理領(lǐng)域的重要研究工具 ， 它能夠在沒有物理接觸的情況下俘獲和移動微納粒子以及活體細胞.除此之外 ， 光力還可以實現(xiàn)對微粒子的篩選和分類. 高效地分離細胞和其他生物分子 ， 是生物學和醫(yī)學研究中的關(guān)鍵技術(shù). 例如 ， 特定細胞亞群的分離對癌癥、自身免疫疾病和遺傳疾病細胞療法的進步至關(guān)重要， 而分離和研究腫瘤干細胞群是提高對疾病的認知和探求新的治療方案的關(guān)鍵.

過去的研究中 ， 為了實現(xiàn)對粒子的不同操控 ， 人們研究了各種不同的光場對粒子的作用：具有環(huán)形結(jié)構(gòu)特征的貝塞爾光束被用于探討布朗動力學[16]中的有趣問題;較小的物體借助熱運動才能遷移到光束中心 ， 而較大的物體由于能感受到光束的包絡(luò) ， 可以更快地遷移到光束中心 ， 這項研究成果可用于無源光學分選;渦旋光束[17]與粒子的相互作用則是研究光的自旋和軌道角動量， 以及自旋軌道相互作用的理想工具[18-23];倏逝波是近場 ， 只能穿透介質(zhì)表面很小的距離 ， 從而可以利用倏逝波在樣品內(nèi)部不受干擾的情況下 ， 在其表面移動粒子 ， 由于近場光學不受衍射極限的限制 ， 相比傳統(tǒng)光鑷 ， 倏逝波可以在更小的空間范圍內(nèi)俘獲粒子;光阱陣列可以用來研究存在相互作用的膠體粒子系統(tǒng) ， 該系統(tǒng)常用于模擬熱力學系統(tǒng);光力作用下膠體粒子間的相互作用等效于熱力學系統(tǒng)中的微觀相互作用 ， 優(yōu)點是現(xiàn)在有條件測量和控制該微觀相互作用 ， 同時也能夠模擬系綜的宏觀熱力學行為[24-25].

研究以上提到的復雜光力場對粒子的操控 ， 都需要首先清楚粒子在其中的運動情況. 因此 ， 粒子運動軌跡的正確描繪是光學微操控研究中至關(guān)重要的問題. 例如， 了解光鑷中被囚禁粒子的運動軌跡， 有利于理解原子、分子的冷卻機制;而光分選中， 粒子軌跡直接決定了分選的結(jié)果.光場調(diào)制下膠體粒子系統(tǒng)的運動軌跡模擬對眾多物理現(xiàn)象的理解非常有幫助 ， 如受到驅(qū)動的電荷密度波[26]、二維電子氣中的電子能態(tài)[27]、原子在晶體表面的遷移[28]、化學反應(yīng)動力學和 II 型超導體[29]中的通量流動等.

基于此， 本文選擇微納膠體粒子作為研究對象， 運用歐拉-理查森算法計算粒子在光力場中的運動軌跡 ， 并給出判定該算法所得結(jié)果收斂性的方法 ， 以保證計算結(jié)果是粒子通過光場時運動情況的穩(wěn)定解.物體所受的合力恒定不變或變化規(guī)律可以用某一明確的解析表達式描述時 ， 無論粒子的運動軌跡是直線還是較復雜的曲線 ， 都可以通過嚴格的解析方法求解出來. 但對于力的解析表達式未知、力場不均勻等情形 ， 粒子運動情況就非常復雜了 ， 無法通過解析方法得到粒子的運動軌跡 ， 從而不得不發(fā)展相應(yīng)的數(shù)值方法來求解.本文給出的方法不僅僅適用于光力場 ， 同時也適用于任意大小的粒子在任意力場中的運動情況.

1? 理論方法

這里介紹一種任意尺寸的粒子在任意光力場中運動軌跡的計算方法 ， 并給出判斷每一步長下收斂的條件.在任意入射光場照射下 ， 作用在粒子上的時間平均電磁力可以利用麥克斯韋應(yīng)力張量的面積分表示 ， 公式為

其中， ??是時間平均的麥克斯韋應(yīng)力張量 ， 其表達式為

其中， ε0、μ0分別是背景介質(zhì)的介電常數(shù)和磁導率; Re表示實部;上標*表示復共軛;電場和磁場都是總場 ， 即入射場的 i （ i）和散射場的 s （ s）的總和 ， 即

（3）

式（1）中， 計算電磁力的積分應(yīng)在粒子的外表面 S 上進行， ！n 和dσ分別表示表面 S 上單位法向矢量和面積元. 如果背景介質(zhì)是介電常數(shù)ε0和磁導率μ0的無吸收損耗的介質(zhì) ， 則積分可在包圍粒子的任何表面上進行. 該方法可用于計算任意尺寸粒子在任意光場中受到的光力.

粒子受到的光力隨空間位置變化而形成光力場. 光力場通常較復雜 ， 因而粒子在光力場中運動軌跡的嚴格計算是一個難題. 常用的計算粒子運動軌跡的方法是歐拉算法 ， 該方法假定速度和加速度在時間步長?t 內(nèi)不發(fā)生顯著變化.因此， 為了得到需要的數(shù)值解， 時間步長?t 必須足夠小.但是如果?t 太小 ， 會出現(xiàn)一些數(shù)值計算中的其他問題. 比如 ， 時間步長?t 越小意味著計算中迭代的次數(shù)越多 ， 隨著迭代次數(shù)的增加 ， 有限精度的浮點數(shù)的截斷誤差會逐漸累積 ， 最終數(shù)值結(jié)果會變得不準確;此外 ， 迭代次數(shù)越多 ， 計算機程序完成計算所需的時間就越長. 除了這些問題外 ， 歐拉算法對于許多系統(tǒng)都是不穩(wěn)定的 ， 這意味著誤差以指數(shù)級累積 ， 數(shù)值解會很快變得不準確.因此 ， 發(fā)展更精確和穩(wěn)定的數(shù)值算法非常必要.本文將基于改進的歐拉算法—歐拉-理查森算法得到粒子運動軌跡.

簡單的歐拉算法表述公式為

其中 ， vn 、sn 、an 分別為tn 時刻的速度、位置和加速度 ， vn+1、sn+1分別為tn+1即tn +?t 時刻的速度和位置.公式（4）、公式（5）假定速度和加速度在時間步長?t 內(nèi)無顯著變化 ， 即在計算位移時認為粒子在?t 時間內(nèi)做勻速直線運動 ， 而在計算速度時又認為粒子做勻加速直線運動 ， 且假定初始時刻（即 tn 時刻）的加速度和速度是?t 時間內(nèi)粒子的速度和加速度. 簡單歐拉算法在振蕩系統(tǒng)里誤差非常明顯 .后來人們又發(fā)展了歐拉-克羅默算法， 該算法可提高振蕩系統(tǒng)里的計算精度.歐拉-克羅默算法和簡單歐拉算法的區(qū)別是 ， 在?t 時間內(nèi)粒子的運動速度采用時間段結(jié)束時的速度vn+1 ， 而不是開始時的速度vn來計算新位置sn+1 ， 公式為

計算速度的方法與簡單歐拉算法無異 ， 公式為

但對于力場隨著速度變化的復雜問題 ， 以上兩種方法都會帶來明顯的誤差. 人們發(fā)現(xiàn) ， 選擇時間步長中間的速度而不是開始或結(jié)束時的速度 ， 計算精度能夠得到明顯提高.基于該思想發(fā)展的算法就是歐拉-理查森算法 ， 它對于解決力場依賴速度的問題非常有效;當然該算法也適用于不依賴于速度的力場問題 ， 且計算更加精確有效. 歐拉-理查森算法是利用歐拉算法先求出在中間時刻?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]）數(shù)學上可以表示為

以及

盡管在得到每個時間步長的最后結(jié)果前 ， 需要兩次使用歐拉算法 ， 但是因為該算法中的時間步長可以設(shè)置得比歐拉算法和歐拉-克羅默算法中的更長一點 ， 所以歐拉-理查森算法的計算速度仍然優(yōu)于歐拉算法 ， 并且還具有比歐拉算法和歐拉-克羅默算法更高的精度.

以上的歐拉-理查森算法的數(shù)學表示針對的僅僅是一維情況 ， 對于二維或三維情況 ， 在每一維度上分別單獨應(yīng)用上述公式 ， 即可計算任何矢量力場中粒子的運動軌跡.但是選定步長之后 ， 該算法的計算結(jié)果是否收斂 ， 這點在數(shù)值計算過程中尤為重要.下面以二維情況（在平面直角坐標系中描述粒子位置）舉例介紹判斷收斂的方法. 判斷的方法如圖 1所示.

圖1所示的流程圖圍繞步長?t 的選擇和確定展開. 利用歐拉-理查森算法進行計算時 ， 時間步長是非常關(guān)鍵的一個因素：步長太大時 ， 粒子運動的路徑上某一點的力場若發(fā)生突變 ， 則計算結(jié)果會出現(xiàn)較大誤差 ， 甚至不會收斂;而步長太小時 ， 計算步數(shù)太多 ， 計算量過大 ， 計算時間太長. 本文采取的辦法是針對每一步長都進行判斷 ， 只有滿足設(shè)定的判斷條件才認為該步長的選取是合適的 ， 或者說計算結(jié)果已經(jīng)收斂 ， 程序才可以開始下一步長的計算. 具體判斷收斂的過程如下.

第一步 ， 根據(jù)事先設(shè)定的位移條件 ， 初步確定一個時間步長. 本文是以入射光波長的1/64作為時間步長?t 內(nèi)粒子可移動的位移最大值. 給定一個初始時間步長 ， 由歐拉-理查森算法計算這個步長中間時刻即tmid = tn +?t/2和tmid = tn+1時的速度、位置和加速度 ， 并計算出該步長內(nèi)粒子運動的位移.此位移需小于預設(shè)的位移最大值;如果大于位移最大值 ， 表示選取的時間步長過大.此時 ， 將步長減半后重復上面的計算 ， 對新步長下算出的位移再進行判斷 ， 重復該步驟直到滿足位移條件為止 ， 從而初步得到一個時間步長.

第二步 ， 在第一步依據(jù)位移條件給出時間步長的基礎(chǔ)上 ， 進一步判斷該步長是否滿足加速度條件 ， 由加速度條件最終確定合適的時間步長. 加速度條件需要區(qū)分不同的情形 ， 分別進行判斷. 這里的加速度是指 tn 和 tn+1時刻處 （在判斷中分別稱之為初位置和末位置）的加速度. 首先 ， 如果 x 和 y 方向的初始和最終加速度均為0， 表示這個時間步長內(nèi)粒子幾乎不受力而做勻速直線運動 ， 算法給出的是嚴格的結(jié)果 ， 可以進入下一時間步長的計算. 此處要記錄這一步長的末速度、末位置 ， 以備下一步計算使用. 其次 ， 若初位置的加速度均為0， 但末位置的加速度不為0， 則判斷末位置相對初位置在 x 方向和 y 方向上的速度相對誤差（即相對應(yīng)方向的速度變化量除以初速度） ， 這兩個方向上的相對誤差均小于2‰， 則步長符合條件;若初位置的加速度在一個方向上不為0 而在另一個方向上為0， 則計算該方向上的加速度相對誤差和另一方向上的速度相對誤差 ， 這兩個相對誤差都小于2‰， 則步長符合條件;若初位置的加速度均不為0， 則計算 x 方向和 y 方向上的加速度相對誤差 ， 這兩個相對誤差都小于1‰， 則步長符合條件. 上述幾種情形下 ， 若不滿足判斷條件 ， 則將時間步長除以5 （可根據(jù)具體問題需要 ， 選擇減小的倍數(shù)）作為新的時間步長 ， 重復上述計算和判斷 ， 直到得到滿足條件的時間步長. 最后 ， 將最終得到的時間步長代入歐拉-理查森算法 ， 計算出粒子運動的位移、速度和加速度 ， 再進入下一時間步長的計算.

2? 周期性光力場中粒子運動軌跡的模擬

光學微操控中涉及粒子在光力場中運動的時候 ， 由于實際的光力場不一定有解析表達式 ， 粒子運動軌跡很難通過解析辦法得到.此時 ， 本文前面介紹的數(shù)值方法就可以發(fā)揮獨特的作用了.下面以二維周期性光力場中的粒子運動為例 ， 模擬粒子在該力場中的運動規(guī)律 ， 進而驗證上述算法的有效性.

該二維光力場由4 束相同頻率的線性極化平面波（每束平面波的強度為104 W/cm2）沿著兩兩互相垂直的方向入射形成 ， 如圖2所示. 圖2 中 ， 4束波的電場都在 x-y 面內(nèi)振動. 為了保證這4束波全部相干 ， 要求它們在坐標原點的相位相同.在這樣4束相干波產(chǎn)生的光場中 ， 無論多大的粒子受到的光力都只有梯度力而沒有散射力[31] ， 即圖2給出的是一個純凈的周期性梯度光力場. 雖然形式上這一光力場與傳統(tǒng)光晶格沒有差別 ， 但由于光晶格形成時對4束光的相位沒有要求 ， 一般的光晶格不是純凈的梯度光力場 ， 僅僅是對瑞利粒子的一個近似的梯度光力場. 處在這樣純凈的梯度光力場中的微納粒子 ， 由于僅受到保守的梯度力 ， 將會被囚禁在某一個格子內(nèi). 圖2給出了1.596μm ×1.064μm 的方形區(qū)域內(nèi)光力場的空間分布圖 ， 梯度力應(yīng)是矢量.圖2中顯示的僅是梯度力的強度即jFg j ， 力的大小大約在皮牛（pN）量級 ， 圖中的藍色區(qū)域梯度力較弱 ， 紅色區(qū)域較強.

本文選擇密度為1.05 g/cm3的聚苯乙烯膠體球粒子作為研究對象 ， 且可根據(jù)需要 ， 通過改變它的半徑大小來改變粒子尺寸.圖2標示了半徑為0.5μm 的粒子在該光力場中的初始位置（0， 0.45）， 即粒子處在圖中正上方紅色圓圈所包圍的格子里.粒子被無初速釋放后 ， 由于初始時刻粒子處在 x =0 這條線上 ， 粒子僅受到 y 方向的梯度力 ， 而不會受到 x 方向的作用力 ， 粒子將從初始位置 y =0.45開始沿 y 軸正方向運動 ， 到達約（0， 0.60）位置的時候速度降為0， 再開始反向運動 ， 到達初始位置后重復上面的振蕩運動過程. 總地來看 ， 粒子將在這個格子里、在該區(qū)間內(nèi)沿 y 軸做簡諧振動.

接下來利用歐拉-理查森算法 ， 計算粒子釋放后的運動軌跡 ， 看軌跡是否如分析的那樣. 圖3（a）–圖3（d）模擬了粒子釋放后不同時間段內(nèi)的軌跡. 注意圖3中各圖橫坐標 x 方向的數(shù)值僅有10–16 μm.為了清楚地顯示軌跡 ， 圖3中 x 方向的長度被放大顯示了 ， x 方向上的數(shù)值非常非常小. 因此粒子在 x 方向上基本沒有運動 ， 只沿著 y 方向做往復運動. 圖3 各圖也同時顯示了該運動在 y 軸的最遠點是（0， 0.60）.初始點和最遠點都被標記在了圖2中. 由圖2 可以看到 ， 這兩個點在該格子中基本對稱 ， 符合粒子在對稱格子里運動應(yīng)該對稱的特點.無論是往復運動還是運動的最遠點 ， 都與本文在模擬前的物理分析吻合 ， 這意味著本文對運動軌跡的數(shù)值模擬方法是非常有效的.從圖3（a）到圖3（d）， 模擬的時間段越來越長 ， 但粒子做往復運動的規(guī)律并沒有變化 ， 這說明本文的方法并沒有時間積累誤差 ， 完全可以用來模擬長時間內(nèi)粒子的運動軌跡.

3? 結(jié)論

本文在米氏散射理論的基礎(chǔ)上 ， 利用歐拉-理查森算法發(fā)展了計算微納粒子在光力場中運動軌跡的方法.該方法可用于計算原子、分子、膠體粒子、病毒細菌等生物大分子或生物細胞在光力場中的運動. 在文中 ， 我們以球形聚苯乙烯粒子在周期性保守光力場中的運動為例 ， 驗證了該方法的有效性.雖然文中僅模擬了真空背景下粒子在光力場中的運動 ， 但是對于涉及阻尼和隨機布朗運動 ， 甚至粒子同時受到力矩的作用等情況 ， 該方法仍然適用.在這些情形下 ， 計算粒子在各個方向所受的合力將變得復雜 ， 但一旦得到了所有合力 ， 每個步長下的位移就可以通過歐拉-理查森算法求出.本文介紹的歐拉-理查森算法適用于布朗運動、膠體科學和生物醫(yī)學中的光學微操控等諸多研究領(lǐng)域.

[參考文獻]

[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.計算機模擬方法在物理學中的應(yīng)用[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.

（責任編輯：李藝）