石墨烯表面等离激元的粒子模拟仿真方法与流程

本发明涉及计算机模拟仿真技术领域。

背景技术：

石墨烯(graphene)由于具有特殊的属性和巨大的潜在应用，已经成为现代科学技术中最具吸引力的研究领域。由于石墨烯中的等离子振荡频率介于1-50thz之间，属于太赫兹频段和红外线频段范围。因此，在光子与光电子应用研发领域，如超快激光、太阳能电磁、光学调节器、光电探测器、光发射装置，基于石墨烯的等离子的诸多功能扮演着重要的角色。理论和实验都以证明，当介电常数的实部为负时，石墨烯能够激励起表面等离激元。

石墨烯的表面等离激元的基本特性和激励方法类似于贵金属薄膜。例如，石墨烯表面等离子体也能够被入射的平面波、强非线性场和电子束激励。此外，由于石墨烯的电导率和介电常数能够通过调节加载电压和化学掺杂等方式而改变，因而石墨烯表面等离激元具有许多优于金属薄膜表面等离激元的特性。

基于对石墨烯表面等离激元研究的需求，开展石墨烯表面等离激元的粒子模拟仿真，尤其是通过仿真得到外部激励因素以及石墨烯本身的等离子体参数与激励起的表面等离激元的频率、场强及衰减速度之间的关系，将为采用石墨烯表面等离激元辐射机制，制备太赫兹源或红外光源的研究提供参考依据。

对金属薄膜中的表面等离激元进行粒子模拟时，通常采用drude模型对其电导率进行等效，这需考虑金属薄膜的厚度(通常为几十纳米)，并在厚度方向划分不小于10个网格，才能保证模拟的精度。

因石墨烯的厚度远小于金属薄膜，仅为零点几纳米。若在粒子模拟时采用drude模型对其中的等离子体振荡行为进行等效，在厚度方向划分10个或10个以上网格，其网格数将是金属薄膜模拟的数百倍。计算机将无法提供足够的内存对pic计算中网格上的电磁场等参数进行存储。

另外，由于粒子模拟中网格大小与时间步长呈正相关，且石墨烯中表面等离激元的周期(约为飞秒量级)远大于金属薄膜中表面等离激元的周期(约为皮秒量级)，所以其总时间步数将比金属薄膜大3-4个数量级。

粒子模拟的计算量＝(网格数+宏粒子数)х时间步数。因此，若采用对drude模型对石墨烯中的电导率进行等效，以此来实现表面等离激元的粒子模拟。在宏粒子数相同的情况下，其计算量比金属薄膜表面等离激元的粒子模拟大5-6个数量级。

总而言之，采用drude模型实现石墨烯表面等离激元的粒子模拟存在以下几个方面的不足：

(1)划分的网格数量巨大，计算机内存无法承受。

(2)模拟所需总时间步数极多，计算量极大，模拟耗时极长。

技术实现要素：

本发明所要解决的技术问题是，针对现有技术中因石墨烯的厚度过小，表面等离体激元周期相较时间步长过大等特性而造成的石墨烯表面等离激元粒子模拟仿真困难的问题，提供一种石墨烯表面等离激元的粒子模拟仿真方法，可以采用更大的网格和更大的时间步长，在减小内存消耗和计算量的同时，保证模拟的准确性。

本发明解决所述技术问题采用的技术方案是，石墨烯表面等离激元的粒子模拟仿真方法，包括下述步骤：

1)沿石墨烯厚度方向划分网格；

2)通过下式获得网格上的面电流密度在时间和空间上的迭代迭代关系，以及与石墨烯所在网格层电场的迭代关系：

其中，δt为模拟的时间步长，δl为垂直于石墨烯方向的模拟网格大小，e为电子的电量，kb为玻尔兹曼常量，t为环境温度，为普朗克常数，μc为石墨烯的化学势，τ为弛豫时间，et为平行于石墨烯方向的电场，(jf)t为石墨烯所在网格面上的面电流密度；

3)根据空间网格上第n步的电场、第n+1/2步的面电流密度以及第n+1/2步的磁场，由下式迭代求解得到第n+1步各个方向的电场强度ex、ey和ez：

上式中，ex,ey,ez分别代表x,y,z方向的电场，，i,j,k分别代表x,y,z方向的坐标索引值，ε为介质的介电常数，δt为模拟的时间步长，δx,δy,δz分别代表x,y,z方向的网格大小；

4)根据空间网格上第n+1/2步的磁场和第n+1步的电场，由下式计算第n+3/2步的磁场：

第n+3/2步的各个方向的磁场强度hx、hy和hz：

其中，ex,ey,ez分别代表x,y,z方向的电场，hx,hy,hz分别代表x,y,z方向的磁场，各电磁场分量的上标n+1/2和n+1和n+3/2表示场分量对应的时间步，i,j,k分别代表x,y,z方向的坐标索引值，μ0为真空的磁导率，δt为模拟的时间步长，δx,δy,δz分别代表x,y,z方向的网格大小。

本发明对二维石墨烯电导率模型提出的电导率公式进行差分，并将其与时域有限差分法相结合，这样在三维数值建模中可以采用一个二维面模拟石墨烯的电导率特性，从而在粒子模拟中可以采用更大的网格和更大的时间步长，在减小内存消耗和计算量的同时，保证模拟的准确性。采用本发明技术，模拟的结果能帮助人们更全面的了解石墨烯表面等离激元的物理特性，对采用石墨烯制备新型太赫兹辐射源及红外线辐射源的研究具有指导性作用。

附图说明

图1为本发明实施方式的流程图。

图2为石墨烯层两边的边界条件示意图。

图3为模拟模型示意图。该模型中，模型附着于一个相对介电常数为2.1的介质基底上，从阴极发射一个速度为0.04倍光速的电子团，并沿平行于石墨烯的方向运动。模型中，取石墨烯的化学势μc＝0.15ev，弛豫时间τ＝0.5ps，温度t＝300k。模型四周采用完全匹配层吸收边界作为模拟空间的边界条件。

具体实施方式

本发明的石墨烯表面等离激元的粒子模拟仿真方法包括：

步骤一、根据线电流密度的概念，将二维石墨烯模型提出的电导率公式与石墨烯所在网格的电场相乘，获得该网格上面电流密度的求解公式；(参见下文式7)

步骤二、采用时域有限差分方法，对面电流密度的求解公式进行差分，得到其在时间和空间上的迭代公式，及其与石墨烯所在网格层电场的迭代关系；(式10)

步骤三、设置模拟初始条件，采用时域有限差分法进行计算。在时域有限差分法中将电场e和磁场h在时间上相差半个步长交替计算，通常将e放在整数时间步，将h放在半整数时间步，另外电流密度jf也是放在半整数步进行计算。所以根据空间网格上第n步的电场、第n+1/2步的面电流密度以及第n+1/2步的磁场结合麦克斯韦磁场旋度差分方程，可求得第n+1步的电场；(式11-13)

步骤四、将石墨烯所在网格面上的第n+1/2步的面电流密度及步骤三求得的第n+1步的电场，代入步骤二所得的差分公式，求得石墨烯所在网格面上第n+3/2步的面电流密度。同时，根据空间网格上第n+1/2步的磁场和第n+1步的电场，求得第n+3/2步的磁场；(式10、14-16)

步骤五、按顺序循环执行步骤三和步骤四，即可实现整个模拟区间电磁场在时间上的迭代求解，最终实现对石墨烯表面等离激元的粒子模拟仿真。

本发明所述的第n+1/2步是指在第n步以后的半步，这是由于在时域有限差分法中，根据差分公式，电场和磁场需要交替求解，所以通常将电场放置于整数时间步求解，磁场和电流密度放置于半整数步求解，这样就能保证每个时间步内电场和磁场均求解了一次。例如，在第n个时间步内，电场在第n步实现求解更新，而磁场则要滞后半个时间步，即在n+1/2步完成求解更新。同理，在第n+1个时间步内，电场在第n+1步实现求解更新，磁场滞后半个时间步，即在n+3/2步完成求解更新。

本发明提供了一种石墨烯表面等离激元的粒子模拟仿真方法。由于石墨烯表面等离激元的产生，本质上是石墨烯中的等离子体受电磁场的影响产生的集体震荡现象。因此，首先对不同石墨烯的电导率模型进行分析，确定了采用二维石墨烯模型中的电导率公式，即将石墨烯等效成一个二维的平面后，其电导率如式(1)所示：

其中t是温度，kb是玻尔兹曼常量，τ是弛豫时间，μc是化学势。通常在数值计算中取τ＝0.5ps，t＝300k，根据掺杂浓度的不同化学势的取值范围为μc＝(0.08-0.15)ev。

如图2所示，石墨烯层介于相对介电常数分别为ε1和ε2的两层介质之间，那么石墨烯可以看作能激发表面等离激元的导电薄层。则石墨烯层两边的边界条件为：

e1t＝e2t(2)

d1n-d2n＝ρs(4)

b1n＝b2n(5)

其与石墨烯的相互作用主要是在石墨烯中激发起平行于石墨烯的电流，而石墨烯中的线电流密度可以表示为：

js＝σet(6)

由于上面所求得的是石墨烯中的线电流密度js，而在时域有限差分法求解电场的差分方程中用到的是面电流密度jf。因此，需要将js转换为垂直于石墨烯面的jf。根据它们之间的关系，jf可写为：

jf＝jsδz＝σetδl(7)

其中δl为垂直于石墨烯面的网格大小。

由(1)式可推得：

进而可得：

将上式差分可得

上式即为求解石墨烯所在网格面电流密度的差分迭代公式。其中，δt为模拟的时间步长，δl为垂直于石墨烯方向的模拟网格大小，e为电子的电量，kb为玻尔兹曼常量，t为环境温度，为普朗克常数，μc为石墨烯的化学势，τ为弛豫时间，et为平行于石墨烯方向的电场分量，(jf)t为石墨烯所在网格面上的面电流密度分量，各电场分量以及面电流密度分量的上标n+1/2和n+1和n+3/2表示场分量对应的时间步。

模拟时，若要求解第n+1步的电场强度、第n+3/2步的磁场强度以及第n+3/2步的面电流密度。可先采用以下时域有限差分迭代公式，求解第n+1步各个方向的电场强度ex、ey和ez：

上式中，ex,ey,ez分别代表x,y,z方向的电场，hx,hy,hz分别代表x,y,z方向的磁场，各电磁场分量以及面电流密度分量的上标n、n+1/2和n+1表示场分量对应的时间步，i,j,k分别代表x,y,z方向的坐标索引值，ε为介质的介电常数，δt为模拟的时间步长，δx,δy,δz分别代表x,y,z方向的网格大小。

而后，采用(10)式求解石墨烯所在平面网格的面电流密度(jf)t。

最后，采用以下时域有限差分迭代公式，求解第n+3/2步的各个方向的磁场强度hx、hy和hz：

上式中，ex,ey,ez分别代表x,y,z方向的电场，hx,hy,hz分别代表x,y,z方向的磁场，各电磁场分量的上标n+1/2和n+1和n+3/2表示场分量对应的时间步，i,j,k分别代表x,y,z方向的坐标索引值，μ0为真空的磁导率，δt为模拟的时间步长，δx,δy,δz分别代表x,y,z方向的网格大小。

这样就完成了一整个时间步的电磁场求解。不断对以上差分公式依次循环迭代，即可模拟石墨烯表面等离激元从产生到衰落的全过程。

自此，就实现了石墨烯表面等离激元的粒子模拟。

下面举一个实例对本发明进行说明。本实例的模拟将激励起频率为太赫兹波段的表面等离激元。其具体步骤包括：

1、确定pic模型的基本参数。

该模型采用一个平行掠过石墨烯表面的电子团进行激励，电子团的速度为0.04倍光速。石墨烯附着于一个相对介电常数为2.1的基底上，取石墨烯的，化学势μc＝0.15ev，弛豫时间τ＝0.5ps，温度t＝300k。

2、根据色散曲线公式理论计算得出，该电子团将在石墨烯中激励起频率为8.44thz的表面等离激元。

3、根据步骤1的有关参数，构建如图3所示的pic模型，启动模拟仿真。

4、查看模拟结果，分析模拟结果，获得该电子团在石墨烯中激励起的表面等离激元的频率，电场强度及衰落速度等参数。其中，模拟得到的频率为8.44thz与色散曲线公式计算结果一致，印证了模拟方法的正确性。

自此，就完成了一次石墨烯表面等离激元粒子模拟仿真。