李桂萍，王茂琰，李海龍，張小川，董宇亮，徐 軍

?

FDTD法模擬手征介質(zhì)柱的波傳播和洛倫茲力密度

李桂萍1，王茂琰1，李海龍1，張小川1，董宇亮2，徐 軍1

(1. 電子科技大學(xué)物理電子學(xué)院 成都 610054；2. 四川中測微格科技有限公司 成都 610021)

基于輔助差分方程時(shí)域有限差分法，模擬了色散手征介質(zhì)柱的電磁場和洛倫茲力密度分布。從本構(gòu)關(guān)系出發(fā)，給出了手征介質(zhì)中頻域電、磁極化強(qiáng)度與感應(yīng)電、磁極化強(qiáng)度和耦合電、磁極化強(qiáng)度之間的關(guān)系；并給出了波方程和電場的迭代公式。推導(dǎo)了手征介質(zhì)中含束縛電荷、電流和束縛磁荷、磁流密度的時(shí)均洛倫茲力密度表達(dá)式。與相關(guān)文獻(xiàn)結(jié)果進(jìn)行了對比，驗(yàn)證了輔助差分方程時(shí)域有限差分法和洛倫茲力密度方法的正確性。仿真了增益手征介質(zhì)柱的場和光力分布情況，討論了電磁流和電磁荷對洛倫茲力密度的貢獻(xiàn)，為手征介質(zhì)在光鑷和手征參數(shù)測量等工程應(yīng)用提供了理論指導(dǎo)。

手征介質(zhì); 色散; 時(shí)域有限差分方法; 力

自從文獻(xiàn)[1]發(fā)現(xiàn)單一波束產(chǎn)生的梯度力能夠吸引介質(zhì)微粒以來，激光對顆粒的捕獲效應(yīng)，即“光鑷”技術(shù)，由于其在物理、機(jī)械、化學(xué)和生命等學(xué)科的潛在應(yīng)用，引發(fā)了對光學(xué)操縱等的研究熱潮。為了得到負(fù)的光學(xué)力，文獻(xiàn)[2-3]通過設(shè)計(jì)特殊的貝塞爾光束或者通過如非磁各向異性介質(zhì)、梯度折射率介質(zhì)、增益介質(zhì)和手征介質(zhì)等特殊材料來實(shí)現(xiàn)。文獻(xiàn)[4]給出了排列成螺旋形狀的25個(gè)金屬球所受光力的解析解，文獻(xiàn)[5]討論了光場施加在介電系數(shù)虛部為負(fù)數(shù)時(shí)增益材料負(fù)的輻射光壓。而對于激光與增益手征介質(zhì)相互作用的復(fù)雜機(jī)理尚需要進(jìn)一步的研究。

增益材料包括固體、液體、氣體和半導(dǎo)體等。有些增益材料是非手征的，有些增益材料具有手征性，如手征異向介質(zhì)[5]、細(xì)胞中的綠色熒光蛋白[6]、碳納米管[7]等。目前手征介質(zhì)光學(xué)力的研究大部分采用Mie理論[8]或通過設(shè)計(jì)特殊的結(jié)構(gòu)光束[9]。Mie理論盡管具有準(zhǔn)確度高和速度快等優(yōu)點(diǎn)，但是不適合于求解非球形粒子，且Mie理論是基于Bohren的波分解技術(shù)，無法實(shí)時(shí)計(jì)算手征介質(zhì)的磁電耦合效應(yīng)。

與解析解相比[8]，時(shí)域有限差分方法(finite- difference time-domain, FDTD)作為一種計(jì)算簡單、表達(dá)直觀的數(shù)值方法[10-19]，具有廣泛的適用性，能夠模擬任意形狀手征介質(zhì)隨時(shí)間和空間變化的電磁場分布情況。與其他數(shù)值方法相比，F(xiàn)DTD方法可以模擬天然有機(jī)分子、人工玫瑰花型和等效色散的手征介質(zhì)[11, 13]。除了基于與Mie理論類似的波場分解技術(shù)的BI-FDTD方法[10]，色散的FDTD方法還可以直接處理本構(gòu)關(guān)系為磁電耦合的手征介質(zhì)[12]。麥克斯韋應(yīng)力張量(Maxwell’s stress tensor)和洛倫茲力(Lorentz force)是常用的兩種光學(xué)力計(jì)算方法，均基于微粒的電磁場分布。麥克斯韋應(yīng)力張量的方法，是在包含所計(jì)算結(jié)構(gòu)的任意閉合曲面S上對麥克斯應(yīng)力張量做面積分，能夠獲取寬頻段的光學(xué)力分布情況；而洛倫茲力方法可以基于所計(jì)算結(jié)構(gòu)的時(shí)諧電磁場分布，提取結(jié)構(gòu)中任意位置的光學(xué)力分布情況。文獻(xiàn)[3]基于FDTD方法，計(jì)算了兩個(gè)手征介質(zhì)板之間的作用力，無須考慮電荷和磁荷密度對手征介質(zhì)輻射光壓的影響；此外，與一維情況相比，二維FDTD方法能夠模擬電磁波斜入射時(shí)，復(fù)雜形狀手征介質(zhì)的電磁場和洛倫茲力分布情況。

本文基于FDTD方法模擬了電磁波在二維增益手征介質(zhì)柱的波傳播和洛倫茲力密度的分布情況。首先給出了基于輔助差分方程(auxiliary differential equation, ADE) FDTD方法中手征介質(zhì)的電極化和磁極化強(qiáng)度，推導(dǎo)了計(jì)算手征介質(zhì)的波方程和洛倫茲力密度。模擬了二維普通介質(zhì)板的場和洛倫茲力密度的分布，驗(yàn)證了本文方法和程序的正確性。最后分析了手征介質(zhì)柱的同極化和交叉極化力密度分布情況，討論其潛在工程應(yīng)用。

1 理論

1.1 本構(gòu)關(guān)系

頻域各向同性手征介質(zhì)中的磁電耦合本構(gòu)關(guān)系可表示為[18]：

式中，()、()、()和()分別是與頻率相關(guān)的介電系數(shù)、磁導(dǎo)系數(shù)、手征參數(shù)和非互易參數(shù)。本文僅討論純手征介質(zhì)的情況，即()=0。絕大部分天然和人工手征介質(zhì)的宏觀等效介質(zhì)參數(shù)[16-18]由其材料屬性、物理幾何結(jié)構(gòu)和電磁波入射角度等決定。洛倫茲模型一般用來表征手征介質(zhì)的介電系數(shù)和磁導(dǎo)系數(shù)，而Condon模型用來表示手征參數(shù)，即：

式中，ε、μ、∞和∞分別表示頻率為零和無窮時(shí)的相對介電系數(shù)和磁導(dǎo)系數(shù)；ω、ω和ω表示諧振角頻率；ξ、ξ和ξ表示阻尼系數(shù)；τ表示表征手征介質(zhì)旋光幅度的特征時(shí)間常數(shù)。人工和生物手征介質(zhì)的手征參數(shù)受其幾何尺寸，如螺旋半徑、寬度、厚度、輪廓長度和俯仰角等的影響。

手征介質(zhì)中感應(yīng)電流密度和磁流密度，耦合電流密度和磁流密度為：

如果考慮手征介質(zhì)的增益損耗特性，無源手征介質(zhì)的介電系數(shù)、磁導(dǎo)系數(shù)和手征參數(shù)滿足如下條件[20]：

如果介質(zhì)參數(shù)不滿足式(12)中的任意一條，手征介質(zhì)即變成有源增益材料。

利用式(6)～式(9)的時(shí)域表達(dá)式及麥克斯韋方程，可得出手征介質(zhì)中傳播模的方程為[3]：

二維情況下，各物理量與無關(guān)。采用ADE- FDTD法可推導(dǎo)出無源色散手征介質(zhì)中的迭代公式。鑒于篇幅原因，本文只給出TM極化波的E，J和K的迭表示達(dá)式為：

式中，有：

1.2 洛倫茲力密度

將式(23)代入麥克斯韋方程組中的磁感應(yīng)強(qiáng)度散度和電場旋度方程，有：

式(25)可重寫為：

根據(jù)束縛磁流密度定義并將式(26)帶入有：

式中，束縛磁流密度由磁極化強(qiáng)度和耦合電極化強(qiáng)度決定。

束縛磁荷密度可定義為：

采用相似的推導(dǎo)過程，可以得到手征介質(zhì)中束縛電荷密度e_bound和束縛電流密度e_bound的表達(dá)式為：

電磁場、電荷和電流密度都是時(shí)間和空間坐標(biāo)的函數(shù)。因此，施加在手征介質(zhì)的光力可通過計(jì)算洛倫茲力密度的時(shí)均值而求得。

通過對含點(diǎn)頻時(shí)諧電磁場、束縛電磁流和電磁荷的式(30)在一個(gè)周期內(nèi)進(jìn)行積分和求平均，可以得到手征介質(zhì)的時(shí)均洛倫茲力密度。

因?yàn)榇怪庇诮橘|(zhì)分界面的磁感應(yīng)強(qiáng)度⊥必須連續(xù)，即磁場⊥是不連續(xù)的，所以束縛磁荷僅存在于兩種相鄰介質(zhì)之間的表面上；束縛電荷同樣適用于類似的邊界條件。

電磁場施加在手征介質(zhì)上與時(shí)間和空間相關(guān)的洛倫茲力密度(,)[21]可表示為：

將式(31)進(jìn)一步用Yee元胞進(jìn)行數(shù)值離散，結(jié)合ADE-FDTD方法計(jì)算得到的電磁場，可求出單位周期內(nèi)的時(shí)均洛倫茲力密度。

2 FDTD模擬結(jié)果

2.1 驗(yàn)證算例

時(shí)間=1 100ΔTE極化波H的場分布圖如圖1所示，同時(shí)疊加了電場強(qiáng)度矢量(E,E)。激勵源為連續(xù)波高斯波束()=0.5sin(0)，波束的焦平面在圖1中總場左邊界處，隨空間變化為exp(-2/12)，其中1=0.9 μm，工作波長0=0.65 μm，介質(zhì)板折射率及相對介電常數(shù)ε分別為2和4。波的介質(zhì)波長λ=0.50，算例中距離和時(shí)間離散間隔分別取為=2λ/65和Δ=/2。

圖1 介質(zhì)板電磁場Hz (Ey, Ex)分布的FDTD數(shù)值結(jié)果

圖2 FDTD模擬的洛倫茲力密度Fx (Fy, Fx)分布圖

時(shí)間從=1 165Δ到=1 295ΔTE極化波入射時(shí)，介質(zhì)板時(shí)均力密度F分量的分布圖如圖2所示，同時(shí)疊加了矢量場(F,F)。高斯波束施加在圖2中半個(gè)介質(zhì)板上的單位面積沿軸積分的時(shí)均洛倫茲力密度為1.620 6 pN/m2，與文獻(xiàn)極化波結(jié)果1.657 7 pN/m2相比，相對誤差為2.2%，且圖2中數(shù)值仿真結(jié)果與文獻(xiàn)[18]中的極化結(jié)果幾乎完全一致，證明了本文ADE-FDTD和洛倫茲力密度方法和程序的正確性。

2.2 手征介質(zhì)柱洛倫茲力密度分布

圖3給出了TM極化波垂直入射，半徑為50 nm的手征介質(zhì)柱的電磁場以及電磁荷和電磁流所貢獻(xiàn)的洛倫茲力密度分布圖。入射波工作頻率選為0為461.54 THz。式(3)～式(5)中手征介質(zhì)板的介質(zhì)參數(shù)分別為ε=1.10，∞=0，μ=1.10，∞=0，ω=ω=π× 500 THz，ω=π×200 THz，ξ=1 500ω，ξ=1 500ω，ξ=3，τ=2.4×10-16。手征介質(zhì)板與之相應(yīng)的相對介電系數(shù)ε、磁導(dǎo)系數(shù)μ和手征參數(shù)κ在頻率0=461.54 THz時(shí)分別為1–j1.8×10-5、1–j1.8×10-5和-0.01-j0.016。計(jì)算頻率對應(yīng)的自由空間波長為0=650 nm。由于所選介質(zhì)板介電系數(shù)和磁導(dǎo)率與自由空間非常接近，因此ADE-FDTD方法的Yee元胞空間尺寸選為=0/130=5 nm，時(shí)間步長Δ=/2。

同極化TM和交叉極化TE波表面磁荷和電荷產(chǎn)生的力密度分別為(–2.75×10-5，2.1×10-8)PN/和(9.81×10-6，6.51×10-6) PN/m。平面波施加在手征介質(zhì)柱上總的洛倫茲力密度為(–2.25×10-5，–3.29× 10-6) PN/m?？梢钥闯銎矫嫒肷洳▽κ终鹘橘|(zhì)柱產(chǎn)生捕獲力，分析其根本原因是圖3模擬的手征介質(zhì)為增益材料，滿足Im2()>[Im()Im()/00]，該捕獲力主要由交叉極化波的電磁流密度貢獻(xiàn)。

a. 同極化TM波的|E|分布

b. 交叉極化TE波的|H|分布

c. TM波電磁流貢獻(xiàn)的F(F F)分布

d. TE波電磁流貢獻(xiàn)的F(F F)分布

e. TM波電磁流貢獻(xiàn)的F(F F)分布

f. TE波電磁流貢獻(xiàn)的F(F F)分布

g. TM波磁荷貢獻(xiàn)的F(F F) 分布

h. TE波電荷貢獻(xiàn)的F(F F)分布

i. TM波磁荷貢獻(xiàn)的F(F F)分布

j. TE波電荷貢獻(xiàn)的F(F F)分布

圖3 FDTD仿真手征介質(zhì)柱電磁場和洛倫茲力密度分布圖

3 結(jié)束語

本文仿真了平面波照射下，色散手征介質(zhì)柱的同極化和交叉極化電磁場和時(shí)均洛倫茲力密度分布。推導(dǎo)了模擬手征介質(zhì)的ADE-FDTD方法和洛倫茲力密度公式。與文獻(xiàn)中介質(zhì)板的電磁場和洛倫茲力密度計(jì)算結(jié)果對比，驗(yàn)證了本文算法和程序的正確性。數(shù)值模擬結(jié)果表明，平面波入射時(shí)，耦合的交叉極化波在手征介質(zhì)柱中產(chǎn)生了一個(gè)會聚點(diǎn)，從而捕獲該色散增益的手征介質(zhì)柱，該捕獲力主要?dú)w因于交叉極化波電磁流的貢獻(xiàn)。本文的研究工作可為光鑷在生物大分子中的應(yīng)用和手征參數(shù)的測量提供借鑒。寬頻段手征介質(zhì)的電磁波傳播、散射以及光學(xué)力分布情況是下一步研究工作的重點(diǎn)。

[1] ASHKIN A,DZIEDZIC J M, BJORKHOLM J E, et al. Observation of a single-beam gradient force optical trap for dielectric particles[J]. Optical Letter, 1986, 11(5): 288-290.

[2] GAO Dong-liang, NOVITSKY A, ZHANG Tian-hang, et al. Unveiling the correlation between non-diffracting tractor beam and its singularity in Poynting vector[J]. Laser& Photonics Reviews, 2015, 9(1): 75-82.

[3] WANG Mao-yan, LI Hai-long, GAO Dong-liang, et al. Radiation pressure of active dispersive chiral slabs[J]. Optics Express, 2015, 23(13): 16546-6553.

[4] DING Kun, NG J, ZHOU Lei, et al. Realization of optical pulling forces using chirality[J]. Physical Review A, 2014, 89(6): 063825.

[5] MIZRAHI A, FAINMAN Y.Negative radiation pressure on gain medium structures[J]. Optics Letters, 2010, 35 (20): 3405-3407.

[6] WANG Bing-nan, ZHOU Jiang-feng, KOSCHNY T, et al. Chiral metamaterials: Simulations and experiments[J]. Journal of Optics A: Pure and Applied Optics, 2009, 11(11): 114003-1-114003-10.

[7] GATHER M C, YUN S H. Single-cell biological lasers[J]. Nature Photonics, 2011, 5(7): 406-410.

[8] LIU Hua-ping, NISHIDE D, TANAKA T, et al. Large-scale single-chirality separation of single-wall carbon nanotubes by simple gel chromatography[J]. Nature Communications, 2011, 2(5): 309-1-309-8.

[9] SHANG Qing-chao, WU Zhen-sen, QU Tan, et al. Analysis of the radiation force and torque exerted on a chiral sphere by a Gaussian beam[J]. Optics Express, 2013, 21(7): 8677-8688.

[10] TKACHENKO G, BRASSELET E. Helicity-dependent three-dimensional optical trapping of chiral microparticles [J]. Natural Communication, 2014, 5(7): 4491-1-4491-8.

[11] AKYURTLU A, WERNER D H. BI-FDTD: a novel finite-difference time-domain formulation for modeling wave propagation in bi-isotropic media[J]. IEEE Trans Antennas and Propagation, 2004, 52(2): 416-425.

[12] ZHENG Kui-song, MU Zong-min, LUO Huan, et al. Electromagnetic properties from moving dielectric in high speed with Lorentz-FDTD[J]. IEEE Antennas and Wireless Propagations Letters, 2016, 15(3): 934-937.

[13] PEREDA J A, GRANDE A, GONZáLEZ O V et al. FDTD modeling of chiral media by using the mobius transformation technique[J]. IEEE Antennas and Wireless Propagation Letters, 2006, 5(3): 327-330.

[14] WANG Mao-yan, LI Hai-long, DONG Yu-liang, et al. Propagation matrix method study on THz waves propagation in a dusty plasma sheath[J]. IEEE Trans Antennas Propagations and Propagation, 2016, 64(1): 286-290.

[15] 劉廣東. Padé近似下模擬一般色散媒質(zhì)的FDTD改進(jìn)方案[J]. 電子科技大學(xué)學(xué)報(bào), 2015, 44(6): 845-850.

LIU Guang-dong. Simulation of general dispersive media in Padé approximation of the improved FDTD scheme[J]. Journal of University of Electronic Science and Technology of China, 2015, 44(6): 845-850.

[16] 賴生建, 王秉中, 黃廷祝. 共享內(nèi)存系統(tǒng)中高效并行FDTD計(jì)算方案[J]. 電子科技大學(xué)學(xué)報(bào), 2010, 39(5): 680-683.

LAI Sheng-jian, WANG Bing-zhong, HUANG Ting-zhu. Efficient parallel FDTD computing scheme in space shared memory system[J]. Journal of University of Electronic Science and Technology of China, 2010, 39(5): 680-683.

[17] 劉瑜, 梁正, 楊梓強(qiáng). 混合并行技術(shù)在FDTD計(jì)算中的應(yīng)用研究[J]. 電子科技大學(xué)學(xué)報(bào), 2009, 38(2): 222-226.

LIU Yu, LIANG Zheng, YANG Zi-qiang. Application of hybrid parallel technology in FDTD computing[J]. Journal of University of Electronic Science and Technology of China, 2009, 38(2): 222-226.

[18] ZAKHARIAN A R, MANSURIPUR M, MOLONEY J V. Radiation pressure and the distribution of electromagnetic force in dielectric media[J]. Optics Express, 2005, 13(7): 2321-2336.

[19] DEMIR V, ELSHERBENI A Z, ARVAS E. FDTD formulation for dispersive chiral media using the Z transform method[J]. IEEE Trans Antennas and Propagation, 2005, 53(10): 3374-3384.

[20] 葛德彪, 閆玉波. 電磁波時(shí)域有限差分方法[M]. 第二版.西安: 西安電子科技大學(xué)出版社, 2005.

GE De-biao, YAN Yu-bo. Finite-difference time-domain method for electromagnetic waves[M]. 2nd Edition. Xi’an: Xidian University Press, 2005.

[21] LINDELL V, SIHVOLA A H, TRETYAKOV S A, et al. Electromagnetic waves in chiral and bi-isotropic media[M]. Boston, USA: Artech House, 1994.

[22] MANSURIPUR M, ZAKHARIAN A R. Maxwell’s macroscopic equations, the energy-momentum postulates, and the Lorentz law of force[J]. Physical Review E, 2009, 79(2): 026608-1-026608-11

編 輯 黃 莘

Wave Propagation and the Lorentz Force Density of a Chiral Column Based on the FDTD Method

LI Gui-ping1, WANG Mao-yan1, LI Hai-long1, ZHANG Xiao-chuan1, DONG Yu-liang2, and XU Jun1

(1. School of Physical Electronics, University of Electronic Science and Technology of China Chengdu 610054;2. Sichuan Zhongce Microgrid Technology Co. Ltd. Chengdu 610021)

Based on the auxiliary differential equation (ADE)finite-difference time-domain (FDTD) method, distributions of electromagnetic fields and Lorentz force densities in a dispersive chiral column are simulated. Firstly, relationships between electromagnetic polarization densities and induced electromagnetic polarization densities, as well as coupled electromagnetic polarization densities of chiral media, are presented based on the constitutive relations. Wave equations and recurrence formula of electric are given. Secondly, the Lorentz force density in chiral media containing bound electric charge and electric current densities, as well as bound magnetic charge and magnetic current densities, is derived. Then, we verify the correctness of the ADE-FDTD method and the Lorentz force density method by comparing with literature’s results. Finally, distributions of fields and optical forces for an active chiral cylinder are simulated. The contribution of electromagnetic current and electromagnetic charge densities to the Lorentz force density is discussed. The work in this paper provides some theoretical guidance for chiral media’s potential engineering applications in optical tweezers and measurement of chiral parameter.

chiral media; dispersion; finite-difference time-domain method; force

O441

A

10.3969/j.issn.1001-0548.2017.06.010

2016-09-09；

2017-01-17

國家自然科學(xué)基金(41304119, 41104097)；中央高?；A(chǔ)研究基金(ZYGX2015J041, ZYGX2015J039)

李桂萍(1974-)，女，博士生，主要從事微波電路與系統(tǒng)、雙各向異性介質(zhì)和計(jì)算電磁學(xué)等方面的研究.