WANG Shuo，HU Bin，LIU Juan

（School of Optics and Photonics, Beijing Institute of Technology, Beijing 100081, China）

* Corresponding author，E-mail: hubin@bit.edu.cn

Abstract: Compared with the commonly used simulation algorithms such as Finite Element Method (FEM)and Finite-Difference Time-Domain (FDTD) method, the Boundary Element Method (BEM) has the advantages of high accuracy, small memory consumption, and ability to deal with complex structures. In this paper,the basic principle of three-dimensional BEM is given, the corresponding program based on C++ language is written, and the Surface Plasmon Resonance (SPR) characteristics of a graphene nano-disk structure are studied. The Scattering Cross-Section (SCS) spectral lines of a graphene nano-disk under different chemical potentials, as well as the distributions of electromagnetic fields at the resonance wavelengths are calculated. The electromagnetic response of the graphene nano-disk in the infrared band is analyzed. In addition, considering the common corrugations of graphene materials caused by defects during processing, we study the influence of the geometric parameters of a convex structure in the center of the graphene nano-disk on the resonance intensity, wavelength and field distributions. A spring oscillator model of charge movement is used to explain the simulation results.

Key words: three-dimensional boundary element method; graphene; surface plasmon resonance; scattering characteristics

1 Introduction

When light is incident on the metal surface, the free charges in the metal will be driven by the incident light field to form collective oscillations with strong scattering, absorption or coupling. This characteristic is called Surface Plasmon Polaritons(SPPs) resonance. The SPPs have a strong local field enhancement effect. Because their transmission wavelength is smaller than that of incident light, they can be used in optical waveguide[1], lithography[2]and other fields. When the metal size is reduced to tens or hundreds of nanometers, the electromagnetic wave transmitted along the surface will be confined on the surface, generating Localized Surface Plasmonic Resonance (LSPR). Since the resonance wavelength, resonance intensity and scattering spectrum of LSPR can be modulated by the material, shape and other factors of nanostructures,LSPR has been widely studied and applied in the fields such as biosensing[3], energy[4]and information[5].

Graphene is a two-dimensional crystal made of a single layer of carbon atoms linked into honeycomb lattices[6], and is a good LSPR material with strong local field enhancement effect and low loss.In addition, the properties of graphene can be adjusted by the methods such as chemical doping and applied electric field, so graphene can be made into a new type of LSPR optical device with external modulation[6]. Graphene has a broad application prospect in optical waveguides[7-9], biosensors[10], metasurfaces[11], metamaterials[12,13], photoelectric devices[14]and other optical fields. Graphene nanodisk is a commonly used graphene structure. C.X.Cong et al.successfully prepared an ordered graphene nanodisk array by combining nanospheric lithography with reactive ion etching[15]. Sukosin Thongrattanasiri et al. investigated the application of graphene nanodisks in optical absorption and demonstrated a periodic graphene nanodisk array capable of achieving 100% optical absorption[16]. Zheyu Fang et al.proposed a tunable optical absorber based on graphene nanodisks and investigated the relationship between its absorption characteristics and the size, spacing and chemical potential of the nanodisks[17]. Hugen Yan et al. studied the coupling effect of graphene nanodisks and proved that different coupling modes would lead to different electromagnetic responses[18]. Jialong Peng et al. reported the study of a tunable terahertz half-wave plate based on coupled graphene nanodisks, using a reflective structure to achieve the function of half-wave plate[11]. Lauren Zundel et al. achieved infrared molecular sensing with subwavelength-level spatial resolution by utilizing the electrical tuning properties of the plasmon polaritons of graphene nanodisk arrays[19]. Vasilios D. Karanikolas et al. studied the role of graphene nanodisks between quantum emitters and demonstrated that the strong resonance of graphene nanodisks could increase the interaction distance between quantum emitters by an order of magnitude[20].

However, all of these studies focused on the ideal planar graphene. As a kind of flexible material,graphene can produce convex/concave structures,corrugations and other non-planar defects in the machining process.

Meanwhile, with the progress of modern micro-nano machining technique, the three-dimensional structures with surface corrugations can also be fabricated by artificial stretching, squeezing, mechanical vibration and other methods, and their dimensional parameters can be controlled[21-23]. These three-dimensional structures can also modulate the SPP resonance of graphene. Penghong Liu et al.studied the SPP resonance mode of wedge-shaped graphene waveguides and obtained the SPP propagation constants and local field distributions in different modes[8]. Slipchenko et al. studied the SPP reflection properties of corrugated graphene layers and quantitatively described the effect of corrugations on the reflection intensity[24]. Shengxuan Xia et al. demonstrated the tunable electromagnetically-induced transparency of sinusoidal graphene layers[25].Li Wang et al. investigated the LSPR properties of bent graphene nanoribbon arrays molded on soft substrates and discovered a new decoupling mechanism[26]. However, the above studies were all based on two-dimensional structures. There are few reports on the non-planar graphene with a three-dimensional defect structure.

In this paper, the three-dimensional boundary element method (3D-BEM or BEM) was used to study the SPR properties of flat and convex graphene nanodisks. Compared with other simulation algorithms in common use, the 3D-BEM can simplify the calculation of three-dimensional objects into the solution of the surface field of objects,that is, it can reduce the three-dimensional calculation to two-dimensional calculation. Therefore, this method can reduce the occupation of computer resources, and improve the calculation speed. We have written a 3D-BEM numerical calculation program based on C++ language. From the calculation results, it can be seen that the incident light in the mid-infrared band can stimulate the SPR phenomenon of graphene structure. By adjusting the electrochemical potential of graphene, we found that the resonance wavelength, the scattering intensity, and the FWHM (Full Width at Half Maximum) of the spectral lines all changed regularly. By adding a convex structure to the graphene structure, the scattering spectrum can also be shifted. We have made necessary explanation of the change law given in this paper by using the law of charge motion and the spring oscillator model. Our research will help understand the LSPR properties of graphene, and also expand the application of BEM algorithm.

2 3D-BEM theory

The schematic of scattering in three-dimensional space is shown in Fig. 1. In an isotropic free space, there is a three-dimensional scatterer, which is hit by infinite incident light. The scatterer surfaceSdivides the scatterer space into two parts, namely the outer regionV2and the inner regionV1. The outer normal vector of the scatterer is →n. The refractive index, relative dielectric constant and permeability of the regionV1aren1,ε1andμ1, respectively. The refractive index, relative dielectric constant and permeability of the regionV2aren2,ε2andμ2, respectively.

Fig. 1 Schematic of three dimensional scattering圖 1 三維空間散射示意圖

Fig. 2 Graphene model and surface-element division method (a) Graphene nanodisk model; (b) convex graphene nanodisk model; (c) mesh generation of graphene nanodisk model; (d) outer normal vectors of the meshes generated in graphene nanodisk(shown by blue arrows)圖 2 石墨烯模型以及面元?jiǎng)澐址椒?a)石墨烯圓盤模型；(b)凸起石墨烯圓盤模型；(c)石墨烯圓盤建模的網(wǎng)格劃分；(d)石墨烯圓盤所建網(wǎng)格的外法向矢量(藍(lán)色箭頭所示，見網(wǎng)絡(luò)彩圖)

Fig. 3 Flow chart of boundary element method圖 3 邊界元算法流程圖

By combining the classical Maxwell's equations with Green's function of vector potential, the following equations can be derived to show the distribution of electric and magnetic fields in the whole space[27]:

Next, we present the numerical method ofFirstly, the vector equations (3)and (4) are decomposed inx,yandzdirections to obtain 6 scalar equations[29], as shown in the appendix. Secondly, the surface of the scatterer is discretized and then divided intoMsurface elements.By using the values of electric and magnetic fields at the center of a surface element to approximate the electric and magnetic field distributions of the whole surface element, we can express the 6 scalar integral equations in the form of summation. Thus we can obtain the linear equation set of 6Munknowns[29], namely A6M×6MX6M=B6M. The 6Munknowns are the three components ofandin every surface elements, respectively, as shown in the appendix. Finally, by solving the linear equations, we can obtain the electric and magnetic field values of each surface element. It is worth noting that the calculated electric field value is only the tangential component of the actual electric field.

After the distributions of electric and magnetic fields on the surface of the scatterer are obtained,the magnetic field at any point in space can be obtained by further discretizing the equation (1) and equation (2). Then, by using Maxwell's equationthe electric field at this point can be calculated, and the distributions of electric and magnetic fields in the whole space can be obtained.After that, we can calculate other physical quantities, such as Scattering Cross Section (SCS), Absorption Cross Section (ACS), etc.

The SCS and ACS can be calculated by using the equations (5) and (6)[30]:

whereWsandWaare the scattering power and the absorption power. They can be calculated by using a closed surfaceГsurrounding the scatterer and the energy flux density (Poynting vector):

3 Modeling of the graphene nonadisk

We first consider the graphene nanodisk structure, with a diameter ofd=60 nm, located in thex-yplane, as shown in Figure 2(a). To study the graphene structure with a defect, we assume that there is a bulge at the center of the graphene nanodisk. This convex structure can be obtained by the mechanical vibration of a resonator[23]. The height of the bulge ish, as shown in Figure 2(b). The width of the bulge, marked as “w”, is defined as the length when the height is changed intoh/exp(1). Then the coordinates of the wrinkled nanodisk in thezdirection can be described as:

Next, we need to divide the surface elements of the graphene nonodisk structure through mesh generation. In the modeling process, the graphene thickness is set as 1 nm[32]. We consider the graphene nanodisk as a cylindrical model with very small height. Since the thickness of two-dimensional graphene material is very small relative to the size of the whole structure, the calculation instability can be caused easily. Therefore, an appropriate mesh generation method should be selected. The mesh generation here mainly follows the following two principles. Firstly, the aspect ratio of each mesh must not be too large. The maximum aspect ratio in this paper is 1.5:1. The structural change in the edge area is more complex than that in the central area, so the edge meshes should be denser. According to these two principles, meshes are generated in the disk, as shown in Figure 2(c). The upper and lower surfaces are divided into 11 rings. The width of each ring and the number of surface elements are shown in Table 1.

The flank is divided into two layers, each of which has 116 surface elements. Therefore, there are 1224 surface elements in total. The position of the center point of each face element and its outer normal vector are shown in Figure 2(d) (Color online). The mesh generation in a convex structure is realized by projecting the meshes of flat disk structure onto the convex structure along thezaxis.

Tab. 1 Parameters of mesh generation on the upper and lower surfaces of graphene nanodisk表 1 石墨烯納米圓盤上下表面網(wǎng)格劃分參數(shù)

The optical parameters of graphene are calculated using the method in the reference [32]. Its dielectric constant can be calculated by the following equation:

whereε0is the vacuum dielectric constant,ωis the angular frequency of the incident light,σis the conductivity of graphene, i is the imaginary unit, andtis the thickness of graphene. The conductivity of graphene in the far-infrared and terahertz bands can be expressed as[7,33]:

whereeis the charge per unit,EFis the chemical potential of graphene, andτis the carrier relaxation time.

Our numerical calculation program is written in C++ language. Under the sub-surface element condition mentioned above, the program occupies about 1.6 G of memory. It takes about 5 minutes for an ordinary personal computer (configuration: i7-8550U,4.0 GHz, 8G RAM) to calculate a wavelength point.The calculation process of the whole algorithm is shown in Figure 3.

4 Simulation results and analysis

4.1 Comparison of the results from BEM and FDTD

In order to verify the correctness of the threedimensional boundary element program written, we first calculated the SCS spectrum of graphene nanodisk in the wavelength range of 5~11 μm and the electric and magnetic field distributions at the resonant position, and then compared the results obtained by BEM with those obtained by the commercial software Lumerical FDTD Solutions based on Finite-Difference Time Domain (FDTD) method. The graphene disk under calculation has a diameter ofd=60 nm and a chemical potential ofEF=0.25 eV.The incident light perpendicular to the disk and incident along the ?zaxis is linearly polarized light,and the polarization direction of electric field is along theyaxis. The normalized SCS spectra of the two algorithms are shown in Figure 4(a) and Figure 4(d) respectively. It can be seen that, in the calculated band, only one scattering resonance peak is obtained in either of the two algorithms. The resonant peaks in BEM and FDTD are at 7.04 μm and 7.3 μm, respectively, with the relative error of 3.6%.

The electric and magnetic field distributions calculated by BEM and FDTD under the resonant wavelength are shown in Fig. 4(b)~4(c) and 4(e)~4(f), respectively. By comparison, it can be seen that the electric and magnetic field distributions obtained by the two methods are very similar.For the electric field, the corresponding light fields are mainly concentrated in the two areas along the polarization direction of the incident light (ydirection). The magnetic field distribution is perpendicular to the electric field. This is mainly because the electric field component of the incident light forces the free electrons inside the graphene to oscillate collectively in the direction of electric field to form the plasmon resonance. Therefore, the resonant peaks in Figures 4(a) and 4(d) are the enhanced scattering phenomena caused by the SPR of graphene.

Fig. 4 Comparison of the results from BEM and FDTD. (a) SCS spectrum obtained by BEM; (b) magnetic field distribution obtained by BEM under the resonance wavelength; (c) electric field distribution obtained by BEM under the resonance wavelength; (d) SCS spectrum obtained by FDTD; (e) magnetic field distribution obtained by FDTD under the resonance wavelength; (f) electric field distribution obtained by FDTD under the resonance wavelength圖 4 邊界元算法與有限時(shí)域差分算法結(jié)果對(duì)比圖。（a）邊界元獲得的SCS譜線；（b）共振波長(zhǎng)下邊界元獲得的磁場(chǎng)分布；（c）共振波長(zhǎng)下邊界元獲得的電場(chǎng)分布；（d）有限時(shí)域差分獲得的SCS譜線；（e）共振波長(zhǎng)下有限時(shí)域差分獲得的磁場(chǎng)分布；（f）共振波長(zhǎng)下有限時(shí)域差分獲得的電場(chǎng)分布

4.2 Dependence of SCS spectrum on the chemical potential of graphene

One advantage of graphene over other materials is that its chemical potential can be dynamically adjusted by an applied electric field and other methods. Therefore, we calculated the influence of the change of chemical potential on the scattering characteristics of a graphene disk, without changing other parameter settings as shown in Figure 4. The chemical potentialEFwas set as 0.1, 0.15, 0.2 and 0.25 eV respectively. The calculation results are shown in Figure 5(a). It can be seen from the calculation results that, with the increase of chemical potential, the resonance wavelength is blue-shifted, the scattering intensity gradually increases, and the resonance peak width gradually decreases.

According to the reference [17], the relationship between the plasmon resonance frequency and electrochemical potential of graphene is ωp～(EF/D)1/2. Therefore, as the chemical potential of graphene increases, the resonant angular frequency will increase and the resonant wavelength will be blue-shifted. The number of charge carriers and the density of free charges gathered at both ends of the disk will increase with the chemical potential of graphene, leading to stronger scattering field intensity. The spring oscillator model presented in[32, 34] can explain many SPR phenomena. The scattering spectrum and absorption spectrum in the models meet the following linear law:

whereris the parameter related to structure and material, and Γaand Γsare the absorption coefficient and scattering coefficient respectively.

Since the SCS spectrum obtained by calculation is similar to the ACS spectrum, we only use the SCS spectrum to study the properties of graphene structure. By fitting the equations (12) and (13) with the ACS spectrum and SCS spectrum, we can obtain the scattering coefficientГsand the absorption coefficientГa. It is worth noting that since only the equation (13) is used to fit the scattering spectrum to obtain three relevant parameters, we should pay attention to the way of fitting and combine equation(13) with equation (12) in order to obtain accurate results. Ifω=ω0, both the SCS and ACS will be peaks, and the results will be as follows:

From the above two equations,andcan be obtained.Cabs_maxandCsca_maxcan be obtained from the scattering and absorption spectra. Therefore, we can fit the one parameterГsonly through equation (13), and obtain accurate results.

In Fig. 5(b), we show the relationship between the scattering coefficientГsand the chemical potential. It can be seen thatГsincreases with the chemical potential, which further explains the change rule of the scattering intensity. Also as shown in [32], for the narrow-band model structure of graphene in our paper, the scattering coefficient is directly proportional to the square of the charge involved in resonance and inversely proportional to the square of the resonance wavelength. This can also explain the change rule of the scattering coefficient. According to Reference [32], the half peak width of the scattering spectrum is directly proportional to the extinction coefficientГand the square of the wavelengthλ, that is, Δλ～?！う?, where the extinction coefficientГis the sum of the scattering coefficientГsand the absorption coefficientГa. The relationship between the physical quantity ?！う?and the chemical potential is given in Figure 5(c). The physical quantity will decrease with the increase of chemical potential, thus explaining the change rule of half peak width.

Fig. 6 Influence of convex shape on SCS spectrum. (a) Relationship between SCS spectrum and convex height; (b) relationship between SCS spectrum and convex width圖 6 凸起的形狀對(duì)SCS譜線的影響。（a）SCS譜線與凸起高度的關(guān)系；（b）SCS譜線與凸起寬度的關(guān)系

4.3 Scattering spectrum changing with convex shape

In practical use, graphene may exhibit structural deformation due to the reasons mentioned earlier.Therefore, we also considered the influence of nonplanar convex structure on the LSPR of graphene.We calculated a disk-shaped planar graphene structure (as shown in Figure 4) with a bulge in the center and a chemical potential ofEF=0.2 eV. At first,the total width of the bulge was set as 20 nm, and its height was set as 10, 20, 30, and 40 nm respectively.The SCS spectrum results were obtained, as shown in Figure 6(a). We found that, with the increase of the bulge height, the resonance peak moved to the long-wave direction, resulting in a red shift. At the same time, the resonance peak intensity increased.The red shift of the resonance wavelength was mainly caused by two reasons. First, with the increase of the bulge height, the path of charge movement from one end of the graphene structure to the other end became longer[32]. Second, as the charge passed through the bulge, it could only be driven by the component of the electric field of the incident light in the direction of the path. Therefore, the increase of the bulge height led to the decrease of the charge-driving force, thereby slowing down the response of a charge to the incident light field[26]. That is to say, the bulge hindered the movement of the charge. Both reasons led to a longer cyclical movement period of the free charge, a lower resonance frequency and a red shift of the resonance wavelength. On the other hand, the increase of scattering intensity was possibly caused by the increase of the bulge height, which might lead to the increase of charge carriers inside the graphene and the increase of the charges gathering at both ends of the graphene structure to form a stronger scattering field. From this, we inferred that if the height of the bulge was fixed and its width was changed, a similar changing trend could also be generated to cause a resonance red shift and enhanced scattering. Then we did the calculation and obtained the results shown in Figure 6(b). We set the height ash=30 nm and the width changing from 10 nm to 20.5 nm. The calculation results were similar to what we expected, except that the SCS spectrum had no significant change until the bulge width was more than 18 nm.

5 Conclusion

In this paper, the SPR properties of the graphene nanodisk were studied by using the boundary element method. We established the numerical model of the geometrical structure of a graphene nanodisk, wrote a 3D-BEM algorithm based on C++language, and calculated the scattering spectrum of the nanodisk and its electric and magnetic field distributions at resonant wavelength.At first, we compared the results obtained by BEM with those obtained by the commercial software based on FDTD and verified the correctness of our program. Then we calculated the scattering spectra of the graphene disk under different chemical potentials. We found that with the increase of chemical potential, the resonance wavelength was blue-shifted, the scattering intensity increased, and the resonance peak width decreased gradually. Finally, we added a convex structure to the original graphene nanodisk model,and considered the influence of the convex structure on the scattering properties of graphene. Our research will contribute to the understanding of the physical law of graphene LSPR and the influence of surface corrugations on the LSPR phenomenon.

——中文對(duì)照版——

1 引 言

當(dāng)光照射到金屬表面，金屬內(nèi)的自由電荷在入射光場(chǎng)的驅(qū)動(dòng)下形成集體震蕩，表現(xiàn)出強(qiáng)烈的散射、吸收或者耦合，這個(gè)特性稱為表面等離子體激元共振（Surface Plasmon Polaritons, SPPs）。SPPs具有很強(qiáng)的局域場(chǎng)增強(qiáng)效應(yīng)，由于其傳導(dǎo)波長(zhǎng)比入射光小，因此可以應(yīng)用于光波導(dǎo)[1]、光刻[2]等領(lǐng)域。當(dāng)金屬尺寸縮小至幾十納米或者幾百納米時(shí)，沿表面?zhèn)鲗?dǎo)的電磁波將會(huì)被束縛在表面，形成局域表面等離子體共振（Localized Surface Plasmonic Resonance, LSPR）。由于LSPR的共振波長(zhǎng)、共振強(qiáng)度、散射譜線等能夠通過(guò)改變納米結(jié)構(gòu)的材料、形狀等因素調(diào)節(jié)，所以在生物傳感[3]、能源[4]、信息[5]等領(lǐng)域被廣泛研究和應(yīng)用。

石墨烯是由碳原子按蜂窩晶格鏈接而成的單原子層二維晶體[6]，是一種良好的LSPR材料，可以實(shí)現(xiàn)超強(qiáng)的局域場(chǎng)增強(qiáng)效應(yīng)，同時(shí)損耗較低，此外，石墨烯的性質(zhì)能夠通過(guò)化學(xué)摻雜或者外加電場(chǎng)等方法進(jìn)行調(diào)節(jié)，因而石墨烯可以用于實(shí)現(xiàn)外部調(diào)制的新型LSPR光學(xué)器件[6]，在光波導(dǎo)[7–9]、生物傳感[10]、超表面[11]、超材料[12,13]、光電器件[14]等光學(xué)領(lǐng)域具有廣闊的應(yīng)用前景。石墨烯納米圓盤是一種常用的石墨烯結(jié)構(gòu)，C.X.Cong等人通過(guò)納米球光刻技術(shù)與反應(yīng)離子刻蝕技術(shù)相結(jié)合，成功地制備了有序石墨烯納米盤陣列[15]。Sukosin Thongrattanasiri等人研究了石墨烯納米圓盤在光學(xué)吸收中的應(yīng)用，展示了一種能夠?qū)崿F(xiàn)100%光學(xué)吸收的周期性石墨烯納米圓盤陣列[16]。Zheyu Fang等人提出了基于石墨烯納米圓盤的可調(diào)節(jié)光學(xué)吸收器，研究了其吸收特性與納米盤尺寸、間距以及化學(xué)勢(shì)的關(guān)系[17]。Hugen Yan等人研究了石墨烯納米圓盤的耦合效應(yīng)，證明了不同的耦合方式會(huì)帶來(lái)不同的電磁響應(yīng)[18]。Jialong Peng等人報(bào)道了基于耦合石墨烯納米盤的可調(diào)太赫茲半波片，采用反射式結(jié)構(gòu)實(shí)現(xiàn)了半波片的功能[11]。Lauren Zundel等人利用石墨烯納米圓盤陣列的等離激元的電調(diào)諧特性，實(shí)現(xiàn)了亞波長(zhǎng)空間分辨率的紅外分子傳感[19]。Vasilios D. Karanikolas等人研究了石墨烯納米盤在量子發(fā)射極之間的作用，研究表明石墨烯納米盤的強(qiáng)烈共振能夠使量子發(fā)射器之間的相互作用距離提高一個(gè)數(shù)量級(jí)[20]。

然而，以上研究全部針對(duì)于理想的平面石墨烯。石墨烯作為一種柔性材料，在加工過(guò)程中會(huì)產(chǎn)生凸起、凹陷、褶皺等非平面的缺陷結(jié)構(gòu)，同時(shí)，隨著現(xiàn)代微納加工技術(shù)的進(jìn)步，也可以通過(guò)人為拉伸、擠壓或者機(jī)械振動(dòng)等方法，來(lái)制造具有表面褶皺的三維結(jié)構(gòu)，并可控制它們的尺寸參數(shù)[21-23]。這些三維結(jié)構(gòu)對(duì)石墨烯的表面等離激元共振也會(huì)產(chǎn)生調(diào)制作用。Penghong Liu等人研究了楔槽形狀的石墨烯波導(dǎo)的SPPs模式，獲得了不同模式下的表面等離子體激元傳播常數(shù)以及局域場(chǎng)分布[8]。Slipchenko等人研究了帶波紋的石墨烯層的SPPs反射特性，定量描述了褶皺對(duì)反射強(qiáng)度的影響[24]。Shengxuan Xia等人通過(guò)正弦型石墨烯層證明了可調(diào)節(jié)的電磁誘導(dǎo)透明特性[25]。Li Wang等人研究了在軟基底上塑造的彎折石墨烯納米帶陣列的LSPR特性，并發(fā)現(xiàn)了新的解耦合機(jī)制[26]。然而，以上工作全部是二維結(jié)構(gòu)，對(duì)于存在三維缺陷結(jié)構(gòu)的非平面石墨烯的研究還鮮有報(bào)道。

This sounded like a reasonable idea to all of us kids, so we kept on going with the stories. My mom knew the true story, though. Bobby s mom was a single parent, and she suspected that they just couldn t afford the Easter Bunny.

本論文采用三維邊界元算法（three dimensional Boundary Element Method，3D-BEM or BEM），研究了平坦及帶凸起石墨烯納米圓盤的表面等離子體共振特性。相對(duì)于其他常用仿真算法，3D-BEM能夠?qū)⑷S物體的計(jì)算簡(jiǎn)化為物體表面場(chǎng)的求解，即將三維計(jì)算降低到了二維，能夠減小計(jì)算機(jī)資源的占用，從而提高計(jì)算速度。編寫了基于C++語(yǔ)言的三維邊界元數(shù)值計(jì)算程序，從計(jì)算結(jié)果可以看出，中紅外波段的入射光能夠激發(fā)石墨烯結(jié)構(gòu)的表面等離子體共振現(xiàn)象。通過(guò)調(diào)節(jié)石墨烯的電化學(xué)勢(shì)發(fā)現(xiàn)，共振波長(zhǎng)，散射強(qiáng)度，譜線的半高寬都會(huì)發(fā)生規(guī)律性變化。通過(guò)對(duì)石墨烯結(jié)構(gòu)加入凸起結(jié)構(gòu)，也會(huì)使散射譜線發(fā)生移動(dòng)。本文從電荷運(yùn)動(dòng)規(guī)律，以及借助彈簧振子模型，對(duì)文中的變化規(guī)律進(jìn)行了必要的解釋。本文研究將有助于理解石墨烯的LSPR特性，同時(shí)也推廣了BEM算法的應(yīng)用范圍。

2 三維邊界元算法理論

本文考慮的三維空間散射問題如圖1所示。在各向同性的自由空間中，存在一個(gè)三維的散射物體。無(wú)窮遠(yuǎn)的入射光照射到散射體上。散射體表面S將空間分成兩部分，散射體外部區(qū)域V2和內(nèi)部區(qū)域V1。散射體的表面外法向量為區(qū)域V1和區(qū)域V2的折射率、相對(duì)介電常數(shù)、磁導(dǎo)率分別為n1、ε1、μ1，n2、ε2、μ2。

從經(jīng)典麥克斯韋方程組出發(fā)，結(jié)合矢量格林定理，可以推導(dǎo)出以下方程，用以表示整個(gè)空間的電場(chǎng)磁場(chǎng)分布情況[27]：

獲得了散射體表面上的電場(chǎng)磁場(chǎng)分布之后，再進(jìn)一步離散化方程(1)和(2)，便可以獲得空間中任意一點(diǎn)的磁場(chǎng)。然后利用麥克斯韋方程可以計(jì)算得到該點(diǎn)的電場(chǎng)，進(jìn)而獲得整個(gè)空間的電場(chǎng)磁場(chǎng)分布。之后，可以計(jì)算其他物理量，比如散射截面（Scattering Cross Section，SCS）、吸收截面（Absorption Cross Section，ACS）等。

散射截面和吸收截面可以由式（5）和式（6）獲得[30]：

其中Ws與Wa為散射功率和吸收功率，其可以選擇一個(gè)包圍散射體的閉合表面Г，借助能流密度（坡印廷矢量）來(lái)計(jì)算獲得：

3 石墨烯納米圓盤的建模

首先考慮石墨烯納米圓盤形的結(jié)構(gòu)，其直徑為d=60 nm，位于x-y平面內(nèi)，如圖2(a)所示。針對(duì)有缺陷的石墨烯結(jié)構(gòu)，考慮石墨烯納米圓盤中心處有一個(gè)凸起，這種凸起型結(jié)構(gòu)可以通過(guò)諧振器的機(jī)械振動(dòng)來(lái)實(shí)現(xiàn)[23]，凸起高度為h，如圖2(b)所示。定義凸起寬度為高度變?yōu)閔/exp(1)時(shí)的長(zhǎng)度，標(biāo)記為w，則該褶皺石墨烯納米圓盤在z方向上的坐標(biāo)可以描述為：

側(cè)面分為兩層，每一層分別為116個(gè)面元。因此所有網(wǎng)格面元合計(jì)共1224個(gè)。面元中心點(diǎn)位置及外法向向量如圖2(d)(彩圖見期刊電子版)所示。針對(duì)凸起結(jié)構(gòu)的網(wǎng)格劃分是在平整圓盤結(jié)構(gòu)的基礎(chǔ)上，沿z軸往凸起結(jié)構(gòu)做投影產(chǎn)生的。

石墨烯的光學(xué)參數(shù)計(jì)算采用文獻(xiàn)[32]中的方法。其介電常數(shù)可以由下式獲得：

其中ε0為真空介電常數(shù)，ω為入射光的角頻率，σ為石墨烯的電導(dǎo)率，i為虛數(shù)單位，t為石墨烯厚度。石墨烯電導(dǎo)率在遠(yuǎn)紅外及太赫茲波段可以表示為[7, 33]：其中e為單位電荷電量，EF為石墨烯化學(xué)勢(shì)，τ為載流子弛豫時(shí)間。數(shù)值計(jì)算程序采用C++語(yǔ)言編寫，在上文所述的分面元條件下，程序占用內(nèi)存大約為1.6 G，普通的個(gè)人電腦（配置：i7-8550U，4.0 GHz，8G RAM）計(jì)算一個(gè)波長(zhǎng)點(diǎn)大約需要5 min。整個(gè)算法的計(jì)算流程如圖3所示。

4 仿真結(jié)果與分析

4.1 邊界元算法與有限時(shí)域差分算法的結(jié)果對(duì)比

為了驗(yàn)證編寫的三維邊界元程序的正確性，首先計(jì)算了石墨烯納米圓盤在5~11 μm波長(zhǎng)范圍內(nèi)的SCS譜線和共振峰位置處的電磁場(chǎng)分布，并將邊界元算法獲得的結(jié)果與基于有限時(shí)域差分算法(Finite-Difference Time Domain，F(xiàn)DTD)的商業(yè)軟 件Lumerical FDTD Solutions的 結(jié) 果 做 了 對(duì)比。所計(jì)算的石墨烯圓盤直徑為d=60 nm，化學(xué)勢(shì)為EF=0.25 eV。入射光垂直于圓盤沿?z軸方向正入射，為線偏光，電場(chǎng)偏振方向沿y軸方向。兩種算法的歸一化SCS譜線分別如圖4(a)和4(d)所示?？梢钥闯觯谒?jì)算波段內(nèi)，兩種算法都只得到了一個(gè)散射共振峰。邊界元算法的共振峰位置在7.04 μm，而有限時(shí)域差分算法的共振峰位置在7.3 μm，相對(duì)誤差為3.6%。

圖4 (b)~4(c)、4(e)~4(f)分別給出了在共振峰波長(zhǎng)下，邊界元算法和有限時(shí)域差分算法計(jì)算出的電場(chǎng)和磁場(chǎng)分布。通過(guò)對(duì)比可以看到，兩種方法獲的電場(chǎng)和磁場(chǎng)分布十分相似。對(duì)電場(chǎng)而言，光場(chǎng)主要聚集在沿入射光偏振方向（y方向）的兩個(gè)區(qū)域。而磁場(chǎng)分布則與電場(chǎng)垂直。這主要是由于入射光的電場(chǎng)分量使得石墨烯內(nèi)部的自由電子在電場(chǎng)方向產(chǎn)生集體振蕩，形成等離子體共振。因此，圖4(a)、4(d)中的共振峰是由于石墨烯的表面等離子體共振造成的散射增強(qiáng)。

4.2 散射譜對(duì)石墨烯化學(xué)勢(shì)的依賴關(guān)系

石墨烯材料相對(duì)于其他材料的一個(gè)優(yōu)勢(shì)在于，其化學(xué)勢(shì)可以由外加電場(chǎng)等加以動(dòng)態(tài)調(diào)節(jié)。因此本文計(jì)算了化學(xué)勢(shì)的改變對(duì)于石墨烯圓盤散射特性的影響，其他參數(shù)設(shè)置與圖4保持一致?；瘜W(xué)勢(shì)EF分別選取0.1、0.15、0.2、0.25 eV，其結(jié)果如圖5(a)所示。由結(jié)果可以看出，隨著化學(xué)勢(shì)的增加，共振波長(zhǎng)發(fā)生藍(lán)移，散射強(qiáng)度逐漸增強(qiáng)，且共振峰的寬度逐漸減小。

通過(guò)文獻(xiàn)[17]可知，因?yàn)槭┑牡入x子體共振頻率和電化學(xué)勢(shì)滿足以下關(guān)系：ωp～(EF/D)1/2。所以，當(dāng)石墨烯的化學(xué)勢(shì)增加時(shí)，共振角頻率增加，共振波長(zhǎng)發(fā)生了藍(lán)移。當(dāng)石墨烯化學(xué)勢(shì)增加時(shí)，載流子的數(shù)量增多，聚集在圓片兩端的自由電荷密度增加，由此產(chǎn)生的散射場(chǎng)強(qiáng)度會(huì)更大，進(jìn)而增強(qiáng)了散射強(qiáng)度。文獻(xiàn)[32, 34]所展示的彈簧振子模型，可以解釋很多表面等離子體共振的現(xiàn)象。其中散射譜線與吸收譜線滿足以下線型規(guī)律：

其中r為與結(jié)構(gòu)、材料等有關(guān)的參量，Γa與Γs別為吸收系數(shù)與散射系數(shù)。

由于計(jì)算得到的SCS譜線與ACS譜線相似，故本文只采用SCS譜線來(lái)研究石墨烯結(jié)構(gòu)性質(zhì)，通過(guò)ACS譜線與SCS譜線來(lái)擬合式（12）和（13），可以得到散射系數(shù)Γs與 吸收系數(shù)Γa。值得注意的是，由于僅使用式（13）擬合散射譜線可獲得3個(gè)相關(guān)參數(shù)，需要注意擬合的方式并且結(jié)合式（12），才能獲得精確的結(jié)果。當(dāng)ω=ω0時(shí)，散射、吸收截面均為峰值，結(jié)果如下：

圖5 (b)展示了散射系數(shù)Гs與化學(xué)勢(shì)的關(guān)系，可以看出Гs隨著化學(xué)勢(shì)的提高而提高，這也進(jìn)一步解釋了散射強(qiáng)度的變化規(guī)律。同時(shí)文獻(xiàn)[32]也表明，對(duì)于石墨烯這種窄帶模型結(jié)構(gòu)，散射系數(shù)同時(shí)正比于參與共振的電荷量的平方而反比于共振波長(zhǎng)的平方，由此也可以解釋散射系數(shù)的變化規(guī)律。在文獻(xiàn)[32]中，散射譜線的半峰寬度Δλ～?！う?，半峰寬度正比于消光系數(shù)Г與波長(zhǎng)λ的平方，其中消光系數(shù)Г是散射系數(shù)Гs與吸收系數(shù)Гa之和，圖5(c)中給出了物理量 Γ ·λ2與化學(xué)勢(shì)的關(guān)系，化學(xué)勢(shì)的提高會(huì)引起此物理量的降低，這也解釋了半峰寬度的變化規(guī)律。

4.3 散射譜線隨凸起形狀的變化

在石墨烯材料的實(shí)際使用過(guò)程中，由于上文所述的原因，石墨烯會(huì)產(chǎn)生結(jié)構(gòu)上的形變，因此本文也考慮了非平面的凸起結(jié)構(gòu)對(duì)石墨烯局域等離子體共振的影響。在圖4圓盤型平面結(jié)構(gòu)的基礎(chǔ)上，計(jì)算了中心帶有凸起的石墨烯結(jié)構(gòu)，其化學(xué)勢(shì)為EF=0.2 eV。首先設(shè)置凸起的總寬度為20 nm，其高度分別為10、20、30、40 nm，其SCS譜線結(jié)果如圖6(a)所示。可以發(fā)現(xiàn)，隨著凸起高度的增加，共振峰的位置往長(zhǎng)波方向移動(dòng)，產(chǎn)生紅移。同時(shí)共振峰的強(qiáng)度隨之提高。共振波長(zhǎng)紅移主要是由兩個(gè)原因?qū)е碌模湟皇请S著凸起高度的增加，電荷從石墨烯結(jié)構(gòu)的一端運(yùn)動(dòng)到另一端的路徑變長(zhǎng)[32]；第二，當(dāng)電荷經(jīng)過(guò)凸起的時(shí)候，入射光的電場(chǎng)只有在路徑方向上的分力才會(huì)對(duì)電荷產(chǎn)生驅(qū)動(dòng)作用。因此當(dāng)凸起結(jié)構(gòu)的高度增加時(shí)，導(dǎo)致電荷驅(qū)動(dòng)力降低，使得電荷對(duì)入射光場(chǎng)的響應(yīng)減緩[26]。也就是說(shuō)，凸起對(duì)電荷運(yùn)動(dòng)起到了阻礙作用，這兩者同時(shí)導(dǎo)致自由電荷作周期性運(yùn)動(dòng)的周期變長(zhǎng)，共振頻率降低，共振波長(zhǎng)紅移。至于散射強(qiáng)度的增加，分析認(rèn)為是所設(shè)置的凸起變高之后，可能使石墨烯內(nèi)部的載流子數(shù)量增加，聚集在石墨烯結(jié)構(gòu)兩端的電荷會(huì)更多，散射場(chǎng)會(huì)更強(qiáng)，這就導(dǎo)致散射強(qiáng)度的增加。由此也可以預(yù)想，如果固定石墨烯凸起的高度不變，而改變石墨烯凸起的寬度，也會(huì)產(chǎn)生相似的變化規(guī)律，共振紅移，散射增強(qiáng)。相關(guān)計(jì)算結(jié)果如圖6(b)所示。固定高度h=30 nm，寬度從10 nm變化到20.5 nm，結(jié)果與預(yù)期相似，不同的是，寬度在18 nm以下時(shí)，SCS譜線變化緩慢，而超過(guò)18 nm之后，SCS譜線才有明顯變化。

5 結(jié) 論

本文采用邊界元算法研究了石墨烯納米圓盤的表面等離子體共振特性。針對(duì)石墨烯納米圓盤的幾何結(jié)構(gòu)進(jìn)行數(shù)值建模，編寫了基于C++語(yǔ)言的三維邊界元算法，并計(jì)算了納米圓盤的散射譜和共振波長(zhǎng)下的電磁場(chǎng)分布。首先將邊界元與基于有限時(shí)域差分算法的商業(yè)軟件做了結(jié)果對(duì)比，驗(yàn)證程序的正確性。然后計(jì)算了石墨烯圓盤在不同化學(xué)勢(shì)下的散射譜線。發(fā)現(xiàn)隨著化學(xué)勢(shì)的增加，共振波長(zhǎng)產(chǎn)生藍(lán)移，散射強(qiáng)度增強(qiáng)，并且共振峰的寬度在逐漸減小。最后，在原石墨烯納米圓盤模型的基礎(chǔ)上，加入了凸起結(jié)構(gòu)，研究了凸起對(duì)石墨烯散射性質(zhì)的影響。本文研究將有助于理解石墨烯表面等離子體局域現(xiàn)象的物理規(guī)律以及表面褶皺對(duì)該局域現(xiàn)象的影響。

附錄：

標(biāo)量化后的方程如(16)?(21)所示：

其中P1，P2，Q1，Q2分別為：

離散化后的方程分別如(22)?(27)所示：