【毕业学位论文】（Word原稿）基于平面波展开法计算光子晶体带结构-光信息科学与技术

毕 业 论 文 基于平面波展开法计算光子晶体 带结构 指导教师 刘景 锋 讲师 学院名称 理学院 专业名称 光信息科学与技术 论文提交日期 2011 年 4 月 论文答辩日期 年 月 答辩委员会主席 _ 评 阅 人 _ 摘 要 光子晶体可以有效的控制原子的自发辐射，从而为新型光量子器件的发展提供了先进的技术平台。它也是量子光学研究光与物质相互作用的研究热点。我们利用 平面波展开 法 计算光子晶体 能 带结构。 文中我们计算了 二维方格子 圆柱子 、二维三 角格子圆柱 子、含空气 空圆柱子的 三角 格 子、 三维 金刚石晶格和 木堆积晶格 等 结构 ，并探究 二维光子晶体具有比较大带隙时的形状特点。 光子晶体 就是 具有周期性分布的电介质 结构 。 本文 从电动力学中的麦克斯韦方程组 和电磁场的 波动方程 出发， 利用 固体物理 中的布洛 赫 定理，推导出关于 磁场强度 H 的本征方程 ，再 运用傅里叶变 换把 该本征中的物理 连续量进行离散化 ，把求解该本征方程 的特征解 转化为求解 矩阵方程的特征矩阵。 本文对光子晶体 带结构计算使用美国麻省理工学院（ 发的 以及周边的 计算程序包 ，并且对其 脚本化的使用方法进行深入学习探究 。 对 既定的光子晶体结构 经 算后，把数据综合整理，绘制成带结构图。 通过对平面波展开法进行详细的推导 以及 算程序的使用， 除了计算本文提及的光子晶体结构外，结合实验实践的可行性，可 从理论上寻找具有比较大范围带隙的光子晶体结构。这将 对光子晶体制备上具有很强大的指导意义。 关键词 光子晶体 能带 平面展开法 金刚石晶格 木堆积结构 目 录 1 前言 . 1 光子晶体提出 . 1 光子晶体简介 . 1 电子能带 . 1 光子带隙 . 2 第一块光子晶体 . 3 2 光子晶体计算 . 4 计算方法和依据 . 4 平面波展开法简介 . 4 时域有限差分法简介 . 5 平面波展开法推导 . 5 本征方程推导 . 5 厄密算符性质 . 8 布洛赫定理 . 9 矩阵表达形式 . 14 3 算程序包 . 14 4 结果及讨论 . 16 简单二维正方格子圆柱子 . 16 反二维正方格子圆柱子 . 17 简单二维三角格子圆柱子 . 19 圆柱子含空气孔的二维三角格子 . 20 二维脉络方格子 . 22 二维正方格子中圆柱子和脉络复合 . 23 金刚石球结构 . 25 木堆积结构 . 26 致谢 . 29 参考文献 . 30 . 31 毕 业论文成绩评定表 . 32 1 1 前言 光子晶体 提出 1987 年，美国物理评论快报 分别在三月和五月刊登了两篇 对光学和材料科学的发展具有重大 意义的论文。两位作者分别是普林斯顿大学（ 理系的 （目前在加拿大多伦多大学物理系）和美国贝尔通讯研究所的 （目前在加州大学洛杉矶分校电子电气工程系）。 当时， 授为了“减少激光器在原子自发辐射而浪费的能量”为研究出发点，重新研究新型激光器内部结构，意外地发现了某些晶体中在某特定频率电磁波范围内不能被“传输”，从而可以使原子不能以自发辐射的方式对外释放能量。而 子在含有缺陷的晶体中具有局限性”出发，类推出“光子在含有缺陷的晶体中受到局限”。 他借用了半导体晶体中 “ 电子带隙 ” 的表述，设想一种新型的光子微结构能与半导体晶体对电子 起的作用一样，对光子产生类似的约束。 光子晶体的这两点特性使得它可以有效地对光进行控制，因此光子晶体被认为是可以控制和操纵光的材料。 自然 写道：“只要能制作出使特定频率的电磁波不能通过的介电材料，各种各样几乎幻想的事都会称为可能。” 光子晶体简介 凡是具有周期性分布的电介质的结构都可以称为光子晶体。 以周期性 的空间拓扑结构方式 ，可以把光子晶体分为 三大类：一维光子晶体、二维光子晶体和三维光子晶体。一维光子晶体是介电常数不同的两种介质交互重叠而成，这种结构实际上已经被研究许多年 ，并且已经广泛应用在各种光学系统上，比如，布拉格 反 射镜，就是一种四分之一波长的一维光子晶体。 二维光子晶体是由介质柱周期排列而成，三维光子晶体由球形或其他晶胞按面心立方、体心立方，金刚石 结构 或其他方式排列而成。 若果考虑周边环境的影响或结构自身关于时间维度的变化特性，可以认为存在四维光子晶体 ，例如受到环境温度和声光效应影响的光子晶体 。 电子能带 而 在固体物理研究中发现，晶体中的周期性排列的原子所产生的周期性电势场对电子有一个特殊的约束作用。在这样的空间周期性电势场中的电子运动是由如下的薛定谔方程决定 的 ： 2 2 22 ( ( ) ) ( , ) 0m E V r r t ( ) ( )V r V r R (其中 () 它具有空间周期性 (求解方程 (以发现 ，电子的能量 E 只能取某些特殊值 ， 在某些能量区间内该方程无解 ， 也就是说电子的能量不可能落在 任意 的能量区间。 光子带隙 与普通晶体一样，光子晶体的周期排列具有能带结构，光子能带之间可能存在光子带隙（或光子禁带）。在这个频率范围里的电磁波不能在这个光子晶体传播。而频率位于能带里的电磁波则能在光子晶体里几乎无损地传播。带隙的宽度和位置与光子晶体的介电常数比值，周期排列的尺寸及排列规律都由关系。 光子 频率 带隙 ， 就是光子晶体本身结构具有阻挡某频率范围电磁波在内部传输 的特性，而存在的这个频率波段就称为“带隙” ，又称为“禁带” 。 光子晶体的带隙上分为两类，一类是完全带隙，一类是非完全带隙。完全带隙是指 ，在光子晶体内部各个方向都不能传播此频率范围的电磁波，而非完全带隙是指，在光子晶体内部的某些方向上不能传播此频率范围的电磁波而在某些方向上可以传播。 光子晶体最主要特征是具有光子禁带，落在光子禁带的频率的光在光子晶体中是禁止传播的。光子禁带的特征依赖于两种介质的介电常数之比，一般随着比值得增大而增大；同时还依赖于晶体的几何构性，如三维的金刚石结构的光子晶体中存在完全光子带隙。利用光子带隙这个特性可以实现对光子的极优良的滤波性能。这是由于光子晶体的滤波带宽可以做得比较大 ， 这种大范围的滤波作用利用传统的滤波器是 难以实现的。光子禁带也使得光子晶体可以制造出高效率低损耗反射镜等器件。 光子晶体中的光子态密度具有反常的分布特征，从而使得光子晶体展现出奇特的量子光学性质。如光子晶体可以抑制自发辐射。若光子晶体中的一个原子的自发辐射的光频率正好落在光子禁带中 , 由于自发辐射的几率与光子所在频率的态的密度成正比，而该频率光子的态的数目为零 , 因此自发辐射几率为零 , 自发辐射也就被抑制。 就像半导体中原子点阵可以控制电子传播一样，光子晶体中不同介电常数的排列可以控制一定频率的光的传播。如果介电常数的差异足够大，在电介质交界面 上发生布拉格散射，产生能量的禁带。在三维严格的光子晶体中，禁带内的电磁波无法向任何一个 3 方向传播。如果晶体中出现缺陷射出，在某个位置上电介质排列的周期性被打破了， 那么光波就可以从这个缺陷射出。如果这种缺陷是 晶体中 的线缺陷，那么光就沿着缺陷构成 的通道传播，从而实现对光的传输方向的控制，甚至可以让光转过很锐的拐 弯。 禁带 中心波长在微波波段或者在远红外波段的光子晶体相对而言 是比较容易实现，这就是为什么早期在实验中制造出来的光子晶体禁带波长多数落于这些波段。然而，要实现在可见光波段存在带隙的光子晶体，尤其是 在 三维光子 晶体，是比较难实现。 主要原因是这么小的晶格常数在工艺上难于实现 。 不过随着现代新的光刻工艺技术进步，越来越多高复杂度的光子晶体被制造出来。 第一块光子晶体 在 出光子带隙的不久之后，他就用一种十分巧妙的方法制作出了第一块具有全带隙 的光子晶体。该光子晶体是金刚石结构。他采用的制备方法简单描述如下：在一块介质上，放置一层具有三角形排列、带圆形洞的模板。在每个圆形洞处打孔，其打孔方向与模板平面的法方向呈 角，而三个空之间彼此的夹角为 120。这样就可以形成了三维金刚 石结构，证实了光子晶体禁带的存在。 这种光子晶体是在一块介质上打孔而得到，带隙的中心频率大约在 1大约是无线电通讯应用的频带范围 。而在此之前尽管 有不同方法的理论分析 预言，但是人们对光子晶体带隙的 是否 存在抱有怀疑的态度。 图 1 作光子晶体方法示意图 4 在此之后，人们在工艺上的努力招呼要是为了将带隙的中心频率提高，那就是向红外以及可见光范围推进，然而这遇到了极大的困难。可以这样说，关于光子晶体近几年来的理论研究以及实验探索非常活跃，可是研究员努力的结果或多或少有些不尽人意 。光子晶体的制作难易程度主要取决于光子晶体的维数和带隙所在的波段。随着光刻工艺技术的提高，对于 三维光子晶体 的实验探究有了进一步的成果 。 2 光子晶体计算 计算 方法 和依据 在实验上制备特定结构的光子材料在目前仍有很大的困难，而通常在理论上设计各种结构并对它进行分析研究则实际得多，并对实验有良好的指导作用。由于 麦克斯韦 方程组能在线性介质上得到精确求解，因此在理论上设计光子晶体结构并研究其性质，对光子晶体的实验制备和应用研究具有重要指导意义。在以往研究中，几种从经典电磁理论上发展起来的常用算法被运用于 各种不同的系统中。 理论研究主要是根据要计算的对象和目标来选择不同的方法。平面波展开法 （ 以下简称 是计算频带结构最常用的方法，应用布洛赫（ 原理，求出特定波矢对应的本征频率。由于它利用了晶体的周期性，从而适用于频带结构的计算。转移矩阵法通常用于计算有限大小光子晶体的透射系数和反射系数以及那些介电函数与频率相关或含虚部（有吸收）的结构。基于散射理论的一些方法则将光子晶体材料看成一组散射元，光子禁带被视为所有散射波干涉的结果。此外，时域有限 差分法 （ 能求出波动方程时域上的解，常用于腔品质因子的计算。 平面波展开法 简介 平面波展开法 是对于某一入射方向 k，以平面波的形式展开电磁场，将麦克斯韦方程组化成一个本征方程，求出 k 对应的一组频率本征值。 这种方法计算了第一个光子带隙结构，共用了几千个平面波，在理论上首次证实了完全带隙的存在。尽管实验中这个结构未能在光波段上建立，但仍在微波波段得到实现和测量。由于 麦克斯韦 方程组遵循 标度不变性，所以，若实验上能在 级上得到带隙，则随着结构周期常数的减小，必能在相应的波段得到更高中心频率的带隙。 这种方法的优点是 ， 没有引入假设条件，为频带结构的计算提供了一个稳定可靠的算法；编程简单，可以借助现有算法库中的傅立叶变换（ 矩阵对角 5 化等标准程序。缺点是 ， 当介电函数是频率的函数时，没有确定本征值方程，计算会复杂和困难得多；在实验上通常的测量量如折射率、透射率等，是与特定的频率相联系的，需要求出所有在该频率下的布洛赫波，于是，在计算中必须重新就某一频 率计算 k。 时域有限差分法 简介 而 时域有限差分法 是 将一个单位原胞划分成许多网格单元，直接用有限差分式代替麦克斯韦时域场旋度方程中的微分式，选取合适的边界条件，将麦克斯韦方程组化成矩阵形式的特征方程。这个矩阵是准对角化的，其中只有为数不多的非零矩阵元，明显地减少了计算量，节省了计算机内存。但是，有限差分法没有考虑晶格格点的形状，遇到具有特殊形状格点的光子晶体时，要求得精确解就比较困难。这种方法的计算精度较高，能得到时域解，通过傅立叶变换可以得到频域解，因此， 有一次时域计算代替 频域上逐点计算的潜力。但是 的计算需要大量的存储空间，计算时间长；同时，晶体含有金属时，收敛缓慢。 在电子能带结构的计算上，有两种经典的模型： “近自由电子模型 ”和 “紧束缚模型 ”。前者假设电子波函数近似于平面波，而后者则认为电子与原子核之间有强大的相互作用，以至于电子只能偶尔从一个原子附近跳到另一个原子附近。平面波展开法就类似于光子的近自由模型 ，而时域有限差分法是 类似紧束缚模型的计算方法 。 平面波展开法推导 本征方程推导 从 麦克斯韦 方程组出发，以平面波展开法推导出电磁波在光子晶 体中的本征方程 。在此方程的推到 中，我们假设介质是无源的、非磁性、 介观尺度的 和各向同性 ；研究对象是频率分布很小范围的电磁波，故可以忽略相对介电常数 的色散效应；介质是无吸收的，故是 ()r 是空间位移的实函数。则 程组可表示 微分形式 为： ( , ) = 0( , ) = 0( , ) = ( , )( , ) = ( , ) r B rH r D r(其中， D 为 电位移矢量， E 为电场强度， H 为磁场强度， B 为磁感应强度，上述各量均为时间 t 和空间位移 r 的函数。在上述假设情况下 ： 6 00( , ) = ( , ) , ( , ) = ( ) ( , )t t t t B r H r D r r E r(把 (式子分别代入 (第三 个 和第四个方 程 。其中，t是对时间变量 t 求偏导数，对常量和不含 t 变量的函数不其作用： 00( , ) = ( , )1 ( , ) = ( , )() E r H rH r E 对 (2方程分别两边求旋度，得到， 00( , ) = ( , )1 ( , ) = ( , )() E r H rH r E 旋度操作符 是对矢量作用，对时间 t 标量不作用， 00( , ) = ( , )1 ( , ) = ( , )() E r H rH r E 对 方程组 (方程等号右边的 ( , )t ( , )t (三第四式代换 ： 00( , ) = ( , )1 ( , ) = ( , )() E r D rH r B (用 (换 (的 B 和 D ， 整理得： 200 2200 21 ( , ) = ( , )()1 ( , ) = ( , )() E r E r H (其中，0是真空中的介电常数 ，0为真空中的磁导率 。由001c，即00 21c，代入 (2 7 22222211( , ) ( , )()11( , ) ( , )() E r E r H (假设方程组 (解 （ 即波动方程的本征解 ） 具有以下形式： ( , ) = ( , ) = (E r E r )H r H r )(将 (入 ( 22222211()11()i t i ti t i E r ) E r )rH r ) H r )r(根据 算符运算公式， F F F (有 以下化简 ， ( ( ( ( ( ( (i t i t i t i ti t i t i t i te e e ee e e e E r ) E r ) E r ) = E r )H r ) H r ) H r ) = H r )(把 (入 ( ， 22222211 ()11()i t i ti t i E r ) E r )rH r ) H r )r(进一步用 (简方程组 ( 得： 22221 ()1 () E r ) E r )rH r ) H r )r(磁场分量 H 满足： 21 ( ) ) ( )()( ) = 0c H r H (根据 ( 以求解 H ，进而通过 ( (得 E ： 8 0() ()i E r H ( r )r (这里 定义 21 ( ) ( ) ) ( )() c H r H r H ， 容易证明 为 线性算符 。 厄密 算符性质 关 于 厄密（ 符 。先 引入一个运算规则 ， 由两向量 () *3, dF G F r G r r (“ * ”是求共轭运算。定义 *,F G = G F ， 且 ,(是实数 、 非负，即使 F 是复数。如果 满足 ,( F G ) = ( F G )的算符 ，就称为 厄密算符 。 注意不是所有算符都是厄密算符，只有满足上述关系才是厄密算符。 对 进行 检验， *3* 3*31,()1()1() F G = F G F G F G r = F , 其中应用 F G F G F 可以知道， 是厄密算符 ， 其本征矢量正交完备。 另一方面， 电场分量 E 满足： 2201 ()( ) ( , ) = 0 E r ) E r ) r(根据 ( 以求解 E ，进而通过 ( (得 H ， 0() i H r E ( r ) (定义算符 21 ( ) ( ) ) ( )() c E r E r E ， 检测是否为 厄密算符 ， 9 *3*3*31,()1()1 () F G = F G F G F G r F , 算符 并不是厄密算符，求解 E 的问题并不能简单的归类为求解本征值问题，要复杂的多。所以，处理此问题我们采取的思路是： (1)找到 () (2)通过求解 程组来寻找 () 21 ( ) ) ( )()( ) = 0c H r H (0() ()i E r H ( r )r (布洛 赫 定理 (程组 的 第一个方程 为本征方程 ，若 ()r 具有周期性，即厄密算符 具有周期性，则根据 固体物理中 布洛赫 理论，本征矢量可以表示为： ( r ) = u r , u r + R u rk k k k (其 中 k 为波的传播方向。 将 (入方 程 (，得： 21 )()=0 k r k r u 利用 ( 性质，化简 ( 21()=0i i e e e e k r k r k rk r k u r u r u u r u rk k 进一步化简 ： 10 2()=0 u r u u 利用 (性质，化简 ( 21()=0ii kk k u r u u 其中波矢 k 是守恒的， 有限体积的原胞确定，故 是离散值(n k) ，本征方程为： 2()1( ) ( ) ( ) ( )() c kk k u r u rr k n k n(我们要得到频带关系也就是要知道 (n k)依赖于所有“ n ”和所有“ k ”的取值。要解决这个问题，我们需要考虑以下三部分： （ 1） k 的有限取值 ： 根据 布洛 赫 理论，本征方程的解是以 k 为周期的，这样就将整个 k 空间的问题 转化到第一布里渊区内，根据本征矢量的在 k 空间的对称性可将求解范围再简化到不可约布里渊区内 ，如图 2。 图 2 二维方格子第一布里渊区 （ 2） 有限自由度的限定 ： 在满足 ( ) ( ) = 0i k u ()元展开为1(Nk n m u r ) = b r )的形式，求解本征值问题 2 (k n k n u r ) = u r )转化有限 11 矩阵问 题： 2A h B h (其中 m l m bb，ml m 有限元法。其中平面波法即将 ()e r ) = h，并满足约束条件 ( ) 0 + k。这是依照均匀的“格子”来划分 ，满足是周期性的边条件。平面波法展开实际上是对周期场的傅立叶展开。 （ 3） 、本征值解法的效率 ： 在解决本征方程中，即计算矩阵 (，可以通过计算 A 和 B 来求得本征值，但是这需要大的内存和计算时间。另外一种方法是从某一个初始假设的本征矢量0迭代过程中进行优化。这种方法需要的内存小，时间快。迭代的方法有 许多方法，其中包括小 化方法。根据微扰理论，厄密矩阵的最小本征值0满足： 20 m h A h (可以通过计算最小值来求得本征值0，在计算最小值时，可以采用最速下降、共轭梯度、前承条件最速下降、前承条件共轭梯度算法，不同算法将直接影响计算精度。 二维周期性介质中平面波法求解光场本征方程 。 在周期性介质中，定义介质周期结构的原胞基矢 ()12a ,a，任意矢量都可以表示成： 1 1 2 2+ = r + R = r +(其中 R 是平移矢量，12,是整数。引入倒空间原胞基矢 ()12b ,b，满足 ， 2112123122 2a a kk a a， (类 似有平移倒格矢量： 12其中12,是整数。 12 H 分量的处理 。 根据 布洛赫 理论，光场的磁场分量具有式 (形式，由于 以傅立叶展开为： 1 , 2 e e (则磁场分量可以表示为： ()1 , 2 ( e e G + k H r )k(其中， 1,2)是波矢为（ G k）的平面波磁场分量的两个偏振方向的单位矢量，且12k + ,k + 构成右手坐标系；对应平面波磁场分量在 对本征方程 (2 行求解，根据 ( ) ( ) ( )f f f A A A， 有： ()1 , 2( ) ( )1 , 2()1 , 2()1 , 2( = = ( )j e eH e e e eH e e e G + k G + k r G + k G + k G + k H r )k + 周期性介质中， (r) 可以傅立叶展开为： ( G rGr ) G )(或者将 1() 1 () G )r (则有： 1 ( )1 , 21 ( ( ) ( )() e e G k G G H r ) G k + Gr k( 13 1 ( )1 , 21 ( )1 , 21 ( ) ( )1()( ) ( )= ( ) ( ) = ( ) ( ) ( ) ( )j e e i e e e e G k G G G k G G G k G r G k G r ) k + k + G k G k + G k + 1 , 2j G G (根据 F G G F G F F G F () k + G， ( ) ( ) ( ) ( ) ( ) 0j j j j je e e e e k + G k + G k + G k + G k + G (故： 1 ( )1 , 21 ( ( ) ( ) ( )() e e k + G + G G H r ) G k + G + G k + Gr k (用 ()G替式 (的 G ,上式等价于： 1 ( )1 , 21 ( ( ) ( ) ( )() e e k + G G H r ) G - G k + G k + Gr k (将 (、 (代入 (中，得到： 21 ( ) ( ) 21 , 2 1 , 2 ( ) ( ) ( ) j j e e H e k + G r G + k G ,G G G G - G k + G k + G (在 (两边分别乘以因子 ()G k 并在整个空间积分，可以得到 21 21 , 2 1 , 2 ( ) ( ) ( )j j j e H G , G , - G k + G k + G (在 (两边点乘据 ( ) ( )C A B B C A ，有 ， ( ) ( ) ( ) ( ) j j j je e e e k + G k + G k + G k + G (其中有差乘的性质 A B B A 。 14 则 (可以表示为： 2121 , 2 ( ) ( ) ( ) j j j e e G , G , - G k + G k + G (将 (写成算符表达式为： 2) ) ) ) 2j j j k( G , ( G ( G ( 其中 ： 1) ) ( ) ( ) ( ) j j j k( G , ( G G - G k + G k + G(矩阵表达形式 将 (转化为矩阵形式，首先，将 (中的算符中的部分表达式改写成： 3 3 ( ) ( ) j j j je e e e e e k + G k + G k + G k + G (由3 1 2 3 2 1 ,e e e e e e ，可得： 2 2 2 13 3 1 2 1 1 e e ee e e ee e e e (则 (中的算符可以表示为： 1) ) 2 2 2 111 2 1 1 ( ) ( ) ( ) () j j j e e ee e e e k( G , ( G G - G k + G k + - G k + G k + G (本征方程转化为： 2 2 2 2 1 1 112 1 2 1 1 2 () G G G G Ge e e e H He e e e H G G G G G 2G - G k + G k + G (当 G取1 2 3, , ,. . . . . . G G，算符展开为 2n 阶 复数厄密矩阵，其本征根为 2