北京大学：《微纳光学》课程教学课件（超材料光学讲稿）第三章 介观光学的理论方法

第三章介观光学的理论方法由于介观体系的细微结构可以和波长相比拟，基至比波长小得多，在处理此类问题时，边界尤为复杂，可解析求解的两类特例已不能满足许多实际问题的需要。所以，目前发展出多种数值方法和理论方法，如：纯数值计算：有限时域差分有限元方法FEMFDTD电偶极近似的算法：离散偶极近似DDA耦合偶极子近似CDA周期性结构中常用算法：多重散射MMS和平面波展开PWE处理不规则纳米结构：格林函数方法云GFM转移矩阵方法TMM古英2012年4月ygu@pku.edu.cn

1 第三章 介观光学的理论方法 由于介观体系的细微结构可以和波长相比拟，甚至比波 长小得多，在处理此类问题时，边界尤为复杂，可解析 求解的两类特例已不能满足许多实际问题的需要。 所以，目前发展出多种数值方法和理论方法，如： 纯数值计算：有限时域差分 FDTD 有限元方法 FEM 电偶极近似的算法：离散偶极近似 DDA 耦合偶极子近似 CDA 周期性结构中常用算法：多重散射MMS和平面波展开PWE 处理不规则纳米结构：格林函数方法 GFM， 转移矩阵方法TMM 古英 ygu@pku.edu.cn 2012年4月

有限时域差分（FDTD）FiniteDifference Time Domain在电磁领域用的最广泛的纯数值算法-popular能算任意电磁波段，任意结构，任意介电常数-powerful研究物理问题的手段：实验、解析和数值计算FDTD是纯数值计算法，如格林函数和DDA都是半解析半数值的算法，所以FDTD的适用性更强。原理上说FDTD须用更多的计算资源和计算时间。自前快速发展的计算机技术和计算方法可解决这一问题。FDTD的计算结果需要清楚的头脑和物理思想去分析和提炼。22

2 2 在电磁领域用的最广泛的纯数值算法 －popular 能算任意电磁波段，任意结构，任意介电常数 －powerful 有限时域差分（FDTD） Finite Difference Time Domain 研究物理问题的手段：实验、解析和数值计算 FDTD是纯数值计算法，如格林函数和DDA都是半解析 半数值的算法，所以FDTD的适用性更强。原理上说， FDTD须用更多的计算资源和计算时间。目前快速发展 的计算机技术和计算方法可解决这一问题。 FDTD的计算结果需要清楚的头脑和物理思想去分析和 提炼

FDTD的实质就是数值解满足边界条件的maxwell方程。(1)从maxwell 方程出发：√×E=-aB/at(2)VxH=aD/at+JV.D=pV.B=0andD=&EB=uH加上初始和边条件的考虑后，只有（1），（2)起作用E和H可写成以下形式：V× Etotal = -uOHtoial / Ot -α* HiotalV× H iotal = OEtotal / Ot +oEtotal3Finite differencetime domain method for elctromagneticsKarlS.KunzandRaymondJ.Luebbers,1993CRCpress

3 FDTD 的实质就是数值解满足边界条件的maxwell 方程。 Finite difference time domain method for elctromagnetics, Karl S. Kunz and Raymond J. Luebbers,1993 CRC press. / / 0 and E Bt H DtJ D B D E B H ρ ε µ ∇ × = −∂ ∂ ∇× =∂ ∂ + ∇⋅ = ∇⋅ = = = 从maxwell 方程出发： 加上初始和边条件的考虑后，只有（1），(2) 起作用 （1） （2） E 和 H 可写成以下形式： / * / total total total total total total E HtH H E tE µ σ ε σ ∇× = − ∂ ∂ − ∇× = ∂ ∂ +

Etotal =Eincident +EscatteredHincident+HscatteredH total=将E和H代入后得：V×Einc =-μ.aHinc / tV×Hinc=C.oEinc/ataH scat9*9*HincH scatatuu和另一组方程(μ-lo) aHinc1scalatuuaEscataaEincEscatat8(-) Einc1at884Finite difference time domain method for elctromagneticsKarlS.KunzandRaymondJ.Luebbers,1993CRCpress

4 Finite difference time domain method for elctromagnetics, Karl S. Kunz and Raymond J. Luebbers,1993 CRC press. 将 E 和 H 代入后得： total incident scattered total incident scattered EE E HH H ≡ + ≡ + 0 0 / / inc inc inc inc E Ht H Et µ ε ∇× = − ∂ ∂ ∇× = ∂ ∂ 和另一组方程 0 0 * * ( ) 1 ( ) ( ) 1 ( ) scat scat inc inc scat scat scat inc inc scat H H H t H E t E E E t E H t σ σ µ µ µ µ µ µ σ σ ε ε ε ε ε ε ∂ =− − ∂ − ∂ − − ∇× ∂ ∂ =− − ∂ − ∂ − − ∇× ∂

上页得到的方程就是数值求解的基础。利用以下有限差分aff(x,t,)-f(x,t)f(x,t,)-f(x,t)limat△tA1-0△taff(x,t)-f(x,t)f(x,t)- f(xi,t)limaxAr-→0AxAx以及Escat和Hscat间的关系aH scatQEscataH scat2EscatxEscatatOzayOEscataH scataEscat1VHscatHscatatazayμo5Finite differencetimedomain method forelctromagneticsKarlS.KunzandRaymondJ.Luebbers,1993CRCpress

5 Finite difference time domain method for elctromagnetics, Karl S. Kunz and Raymond J. Luebbers,1993 CRC press. 上页得到的方程就是数值求解的基础。 21 21 0 21 21 0 (, ) (, ) (, ) (, ) lim ( ,) ( ,) ( ,) ( ,) lim t x f f xt f xt f xt f xt tt t f fxt fxt fxt fxt xx x ∆ → ∆ → ∂ − − = ≈ ∂∆ ∆ ∂ − − = ≈ ∂∆ ∆ 以及Escat和Hscat间的关系 利用以下有限差分 0 1 , , 1 , , scat scat scat x z y scat scat y y scat scat scat x z y scat scat y y E H H E E t yz H E E H H t yz ε µ ∂   ∂ ∂ = −   ∂ ∂∂   ∂   ∂ ∂ =− −   ∂ ∂∂  

可以得到scat2Escat,n-1Escat,n2AHAA△tAyAz6011scatn+scat.n22AEscat,nAEscat,HH11X△tAxAzuo△x << 2下面就可以编程计算了，但是要求：(△x)At <</3c可以看到：整个过程几乎没有任何公式推导，将麻烦交给算法和计算机。原则上，可以解决任何maxwell框架内的问题。自前已有几款商用及自由软件可选择：XFDTD,Solution,Comsol(FEM),EastFDTD自由软件：Meep6Finitedifferencetime domain methodforelctromagneticsKarlS.KunzandRaymondJ.Luebbers,1993CRCpress

6 Finite difference time domain method for elctromagnetics, Karl S. Kunz and Raymond J. Luebbers,1993 CRC press. 可以得到 1 1 , , , ,1 2 2 0 1 1 , , 2 2 , , 0 1 1 scat n scat n scat n scat n x x z y scat n scat n scat n scat n y y z x E E H H t yz H H E E t xz ε µ − − − + −   −   ∆ ∆ = − ∆ ∆∆   −   ∆ ∆ =   − ∆ ∆∆   下面就可以编程计算了，但是要求： ( ) 3 x x t c ∆ << λ ∆ ∆ << 可以看到：整个过程几乎没有任何公式推导，将麻烦交给算法 和计算机。原则上，可以解决任何maxwell框架内的问题。 目前已有几款商用及自由软件可选择: XFDTD, Solution, Comsol(FEM), EastFDTD 自由软件：Meep

CDA (Coupled dipoleapproximation)偶极耦合近似概述：将组成单元考虑成偶极子，组成单元之间的耦合也只考虑偶极作用。适用范围：可以看成偶极子的纳米小颗粒或由其组成的阵列。颗粒的尺度要比波长小得多，以使小球的多级共振不会造成大的影响。颗粒之间距离也比较大（一般大于半个波长即可）。组成单元多是球形或椭球形，也可以是其它形状。IneidentwaevectorElechiefieldpolariced

7 CDA (Coupled dipole approximation) 适用范围：可以看成偶极子的纳米小颗粒或由其组成的阵 列。颗粒的尺度要比波长小得多，以使小球的多级共振不 会造成大的影响。颗粒之间距离也比较大（一般大于半个 波长即可）。组成单元多是球形或椭球形，也可以是其它 形状。 概述：将组成单元考虑成偶极子，组成单元之间的耦合也 只考虑偶极作用。 偶极耦合近似

考虑一散射体系，包含多个组成单元，其位置和偶极极化率由r和α表示，则极化强度可以表示为P, = α,Eloc,i注意这里的Eloci是包括入射光场和其他小球散射场的影响。所以，对于一给定的波长。NEloc,i = Einc,i + Edipole,i = Eo exp(ikr)- Z A,Pjj=1,ji根据电动力学的知识，两个偶极子的耦合可以写成r,2P, -3r,(rP,)A,P, = k2eikg ×(y ×P,))ikrj(te353rij8Theoretical studies ofplasmon resonances in one-dimensional nanoparticle chains:narrowlineshapeswithtunablewidths,ShengliZouetal,Nanotechnology(2oo6)

8 根据电动力学的知识，两个偶极子的耦合可以写成 考虑一散射体系，包含多个组成单元，其位置和偶极极化率 由 和 表示，则极化强度可以表示为 注意这里的 是包括入射光场和其他小球散射场的影响。 所以，对于一给定的波长。 Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths, Shengli Zou et al, Nanotechnology (2006) ri αi P E i i loc i =α ,  Eloc i, , , ,0 1, exp( ) N loc i inc i dipole i i ij j j ji E E E E ik r A P = ≠ =+ = −   ∑ 2 2 , 3 5 ( ) 3( ) (1 ) ij ij ikr ij ij j ikr i j j ij ij j ij j ij ij ij r r P r P r rP A P k e e ikr r r × × −  = + −

如果散射体系有N个组成颗粒，那么就有N个类似上面的线性方程。我们通过解一个自洽的N阶线性方程组A'P=E就可以得到每个颗粒的极化强度P，接着就可以计算散射体系的消光，吸收截面等物理量N4元kZ Im(E)LPext2Eoj=1V4元k2/mP,(aj) p,-2kpabs2Eo9Theoreticalstudies ofplasmonresonancesin one-dimensional nanoparticle chains:narrowlineshapeswithtunablewidths,ShengliZouetal,Nanotechnology(2oo6)

9 就可以得到每个颗粒的极化强度 ，接着就可以计算散射 体系的消光，吸收截面等物理量 如果散射体系有N个组成颗粒，那么就有N个类似上面的线性 方程。我们通过解一个自洽的N阶线性方程组 AP E ′ =   Pi * 2 , 1 0 4 Im( ) N ext inc j j j k C π = = ∑ E P E  ( ) * 2 1* 3 2 1 0 4 2 Im 3 N abs j j j j j k C k π α− =     =   −       ∑ PP P E  Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths, Shengli Zou et al, Nanotechnology (2006)

利用CDA方法，并不要求组成颗粒一定是球形只要考虑偶极和偶极耦合反应体系的主要物理性质时，都可以使用。最简单的是小球，其偶极极化率可以由Mie理论中的展开系数直接求出来：3aadipole2h3而对于其他形状的偶极极化率，则需要建立近似模型甚至数值方法求出了10Theoreticalstudies ofplasmonresonancesin one-dimensional nanoparticlechains:narrowlineshapeswithtunablewidths,ShengliZouetal,Nanotechnology(2oo6)

10 而对于其他形状的偶极极化率，则需要建立近似模型甚至数 值方法求出了 利用CDA方法，并不要求组成颗粒一定是球形，只要考虑偶极 和偶极耦合反应体系的主要物理性质时，都可以使用。 最简单的是小球，其偶极极化率可以由Mie理论中的展开系数 直接求出来： 1 3 3 2 dipole a k α = Theoretical studies of plasmon resonances in one-dimensional nanoparticle chains: narrow lineshapes with tunable widths, Shengli Zou et al, Nanotechnology (2006)