衛(wèi)壯志 薛文瑞 彭艷玲 程鑫 李昌勇

1)(山西大學物理電子工程學院,太原 030006)

2)(山西大學,量子光學與光量子器件國家重點實驗室,激光光譜研究所,太原 030006)

3)(山西大學,極端光學協(xié)同創(chuàng)新中心,太原 030006)

1 引 言

太赫茲(THz)波是頻率在0.1 THz到10 THz范圍內的電磁波,是當今世界的研究熱點之一[1?5].THz波在自由空間中的傳輸損耗較大,且遠距離傳輸時由于群速度色散所造成的時延差會導致信號失真[6].為了解決這些問題,用于THz波傳輸?shù)牟▽\而生[7?9].

自從2004年發(fā)現(xiàn)石墨烯以來,其獨特的性能引起了廣泛關注[10,11].由于在石墨烯的表面可以激發(fā)表面等離子激元(surface plasmon polaritons,SPPs)[12],其共振頻率落在THz波段至紅外波段,所以借助于石墨烯這種二維材料,可以構建THz波導,實現(xiàn)對THz波的傳輸[13?15].

基于石墨烯的納米帶波導是近年來用于傳輸THz波的波導.這種波導通常是將石墨烯納米帶鋪在硅基底上構成的[16].相較于獨立存在的單層石墨烯納米帶[17],最顯著的區(qū)別是等效折射率大幅度增加.較強的模式約束性通常會伴隨著較大的傳播損耗.為了降低損耗,人們在石墨烯和硅基底之間加入極薄的一層二氧化硅層作為緩沖層,既可以增加有效傳播距離,又可以提高波導的品質[18,19].在這種波導中,可以在石墨烯納米帶上激發(fā)邊緣模式和波導模式.研究表明,改變工作頻率或者納米帶的寬度,可以調節(jié)模式的數(shù)量、模式面積和傳播距離[20].

涂覆石墨烯的納米線波導是近年來發(fā)展起來的另外一種傳輸THz波的波導.這種波導通常是將單層的[21]或者多層的[22,23]石墨烯涂覆在單根[21]或者兩根圓柱形的[23,24]電介質納米線上構成的.由于沒有邊緣效應引起的損耗,這種波導受到了廣泛的關注.研究表明,涂覆單層石墨烯的單根納米線波導中石墨烯有助于提高模式的傳播長度[21].涂覆雙層石墨烯的單根納米線波導對模式具有較強的約束性,可以減小有效模式面積,增大傳播長度[22].在涂覆單層石墨烯的納米并行線之間的狹縫區(qū)域,場的增強效應非常明顯,可產生較高的梯度力[24].在涂覆雙層石墨烯的納米并行線波導的狹縫區(qū)域,可以得到更高的場增強[23].

迄今為止,尚未見解析分析過基于涂覆石墨烯的三根電介質納米線的THz波導.本文采用多極方法對這種THz波導的模式特性進行研究.通過改變工作頻率、中間納米線半徑、納米線之間的間距以及石墨烯的費米能,對模式的有效折射率的實部和傳播長度進行了詳細分析.這種波導在模分復用方面具有潛在的應用前景[25?27].

2 理論模型及計算方法

基于涂覆石墨烯的三根電介質納米線的THz波導的結構如圖1所示,該波導是由三根半徑分別為ρ0,ρ1和ρ2且軸心處于同一水平面上、涂覆了單層石墨烯的電介質納米線組成.假設結構左右對稱,且取典型值ρ0=ρ1=100 nm.相鄰兩根納米線的圓心間的距離均為a.納米線之間的距離為d.電介質納米線的相對介電常數(shù)為ε1,整個結構鑲嵌在相對介電常數(shù)為ε2的電介質中,并假設ε1=2,ε2=1.將石墨烯看成厚度為零的導體介質,其電導率σg=σintra+σinter,這里σintra和σinter分別是帶內和帶間電導率,它們可以由庫珀公式得到[28]:

其中e是電子的電荷量,kB是玻爾茲曼常量,T=300 K是環(huán)境溫度,~是約化的普朗克常量,ω=2πf是角頻率,Γ=2×1012rad/s是載流子散射率,EF是費米能級.由上式計算可知,石墨烯的電導率的大小主要取決于工作頻率f、費米能級EF以及溫度T.

圖1 基于涂覆石墨烯的三根電介質納米線的THz波導的橫截面示意圖,電介質納米線外側黑色的圓環(huán)為石墨烯Fig.1.Cross section of the THz waveguides based on three graphene-coated dielectric nanowires.The black rings on the outside of the dielectric nanowires are graphene.

本文采用多級方法[29?34]對圖1所示的結構支持的模式進行分析.假設波導結構中的模式在z方向上傳播,以三根電介質納米線的軸心為原點,分別建立極坐標系(r0,φ0),(r1,φ1)和(r2,φ2). 在各自的坐標系中,寫出z方向的電場和磁場分量如下:

納米線0的內部電磁場分量Ez01和Hz01的分布

納米線0的外部電磁場分量Ez02和Hz02的分布

納米線1的內部電磁場分量Ez11和Hz11的分布

納米線1的外部電磁場分量Ez12和Hz12的分布

納米線2的內部電磁場分量Ez21和Hz21的分布

納米線2的外部電磁場分量Ez22和Hz22的分布

其中,In和Kn是修正的貝塞爾函數(shù);

為待定系數(shù);,其中β是傳播常數(shù),μ0是真空的磁導率.利用加法定理[29],將納米線1與納米線2的外電磁場Ez12,Ez22,Hz12和Hz22分別變換到納米線0的坐標系中,則可以得到納米線0的外電磁場Ez|0外和Hz|0外.同理可以得到納米線1的外電磁場Ez|1外和Hz|1外,納米線2的外電磁場Ez|2外和Hz|2外.根據(jù)Maxwell方程組可由納米線外的z方向的電磁場計算出其他的場分量.把單層石墨烯看成厚度為零的導體邊界,利用切向邊界條件

其中i=0,1和2分別代表納米線0、納米線1和納米線2,就可以建立一個齊次線性代數(shù)方程組,

其中[X]為(3)—(8)式中的待定系數(shù)組成的一個列矢量.假設M為(3)—(8)式及其他場分量中n的上限值,也就是說,把無窮項求和截斷為M項求和,則[A]為一個方陣,且這個方陣由24×24個子方陣a(i,j)組成,每個子方陣a(i,j)中包含M×M個矩陣元am,n(i,j),其中i=0,1,···,23;j=0,···,23;m=0,1,···,M;n=0,1,···,M. 例如:

其中δmn為克羅尼克符號.

根據(jù)線性代數(shù)理論,齊次線性代數(shù)方程組(10)有解的充分必要條件是方陣[A]的行列式為0,即

通過求解方程(14),就可以得到模式的有效折射率的實部Re(neff)和有效折射率的虛部Im(neff)以及對應的模式的場分布.進一步可以得到傳播長度

其中c是光速,f是工作頻率.

3 結果與討論

首先,確定最低階模式.對于如圖1所示的結構,當頻率f=35 THz,半徑ρ0=ρ1=100 nm,ρ2=50 nm,間距d=25 nm,費米能EF=0.5 eV時,利用多級方法可以得到5個最低階模式,依次命名為模式1、模式2、模式3、模式4和模式5(Mode 1,Mode 2,Mode 3,Mode 4和Mode 5),如圖2所示.這5個最低階模式的形成可以歸因于三根納米線單獨存在時所支持的兩個最低階模式(見圖3所示的#0模式和#1模式)之間的5種組合,見圖2(a)—圖2(d).為了表述方便,采用符號“+”和“?”分別表示正的和負的表面電荷[24].可以看出,模式1是由#0,#0和#0組合而成;模式2是由#0,#1和#0組合而成;模式3是由#1,#0和#1組合而成;模式4和模式5均是由#1,#1和#1組合而成.

圖2 當工作頻率f為35 THz,納米線半徑ρ0=ρ1=100 nm,ρ2=50 nm,間距d=25 nm以及費米能EF=0.5 eV時,5種模式的場分布 (a)—(d)為5種組合方式;(e)—(i)為z方向的電場分布;(j)—(n)為電場強度分布Fig.2.The field distributions of the five modes when f=35 THz,ρ0= ρ1=100 nm,ρ2=50 nm,d=25 nm and EF=0.5 eV:(a)–(d)The five combinations;(e)–(i)the z-direction electric field distributions;(j)–(n)the electric field intensity distributions.

圖3 涂覆石墨烯的單根電介質納米線的THz波導所支持的兩個最低階模式的場分布 (a),(b)z方向電場分布;(c),(d)電場強度分布Fig.3.Field distribution of the two lowest order modes of the THz waveguides based on single graphenecoated dielectric nanowire:(a),(b)The z-direction electric field distributions;(c),(d)the electric field intensity distributions.

圖4 (a)有效折射率的實部和(b)傳播長度隨工作頻率f的變化關系,實線是采用FEM得到的數(shù)值解,數(shù)據(jù)點是采用多級方法得到的解析解Fig.4.The dependencies of(a)the real part of the effective refractive index and(b)the propagation length on the operating frequency f.The solid lines are the numerical solutions obtained by FEM,the data points are the analytical solutions obtained by using the multipole method.

其次,研究工作頻率f對波導的模式特性的影響. 圖4是當ρ0=ρ1=100 nm,ρ2=50 nm,d=25 nm,EF=0.5 eV時,各個模式的有效折射率的實部Re(neff)和傳播長度Lprop隨工作頻率f變化的關系.其中,實線是采用有限元方法(FEM)得到的數(shù)值解,數(shù)據(jù)點是采用多級方法得到的解析解,下文中均采用這種標注方法.由圖4可明顯看出,兩種方法所得到的結果十分符合.隨著頻率的增大,5種模式的有效折射率的實部都在增大,傳播長度都在減小.模式1和模式2的有效折射率的實部明顯大于其他模式的值.在頻率較低時,模式1與其他模式的傳播長度差值較大,隨著頻率增大,差值在逐漸減小.另外,模式3和模式4的有效折射率的實部和傳播長度基本是一樣的.在變化的過程中,模式1的有效折射率的實部和傳播長度始終是最大的,高階模式的傳播長度會出現(xiàn)交叉現(xiàn)象.

對于圖4(b)中高階模式的傳播長度出現(xiàn)的交叉現(xiàn)象,可以根據(jù)場分布得到解釋. 以Mode 2和Mode 5為例,圖5(a)—圖5(c)是當頻率分別為31,36.5和39 THz時Mode 2的電場分布圖,圖5(d)—圖5(f)是當頻率分別為31,36.5和39 THz時Mode 5的電場分布圖.36.5 THz是Mode 2和Mode 5達到相同傳播長度所對應的頻率.當頻率由31 THz變化到36.5 THz時,對比圖5(a)和圖5(b)可以看出,Mode 2對應的場在石墨烯表面的強度增強,場與石墨烯之間的相互作用增強,波導對模式的束縛性增強,導致傳播長度減小.但對于Mode 5,對比圖5(d)和圖5(e)可以看出,波導對模式的束縛性更強,使得傳播長度的減小量更大.當頻率由36.5 THz變化到39 THz時,對比圖5(b)和圖5(c)可以看出,Mode 2對應的場在石墨烯表面的強度增強,波導對模式的束縛性增強,導致傳播長度減小.但對于Mode 5,對比圖5(e)和圖5(f)可以看出,波導對模式的束縛性更強,使得傳播長度減小量更大.因此出現(xiàn)了圖4(b)中高階模式的傳播長度的交叉現(xiàn)象.

再次,研究納米線之間的間距d對波導的模式特性的影響. 圖6是當ρ0=ρ1=100 nm,ρ2=50 nm,f=35 THz,EF=0.5 eV時,各個模式的有效折射率的實部Re(neff)和傳播長度Lprop隨納米線之間的間距d變化的關系.可以看出,點與線十分符合.圖6顯示,隨著間距的增大,5種模式的有效折射率的實部都在減小,傳播長度都在增大并最終都趨于穩(wěn)定.隨著間距從10 nm增大到50 nm時,模式1和模式2的變化程度比其他模式要明顯,并且二者的有效折射率的實部明顯大于其他模式的值.相比較而言,模式5的變化比較平緩;而模式3和模式4的有效折射率的實部和傳播長度基本不受間距變化的影響,并且這兩個模式的值也是基本一樣;當d>16 nm時,模式1的傳播長度是所有模式中最大的,高階模式的傳播長度會出現(xiàn)交叉現(xiàn)象.

圖5 當ρ0= ρ1=100 nm,ρ2=50 nm,d=25 nm,EF=0.5 eV,(a)—(c)頻率分別為31,36.5和39 THz時Mode 2的電場分布,(d)—(f)頻率分別為31,36.5和39 THz時Mode 5的電場分布Fig.5.The distributions of the electric field of Mode 2 when the operating frequencies f are 31 THz(a),36.5 THz(b)and 39 THz(c),and the distributions of the electric field of Mode 5 when the operating frequencies f are 31 THz(d),36.5 THz(e)and 39 THz(f)with ρ0= ρ1=100 nm,ρ2=50 nm,d=25 nm and EF=0.5 eV.

圖6 (a)有效折射率的實部和(b)傳播長度隨間距d的變化關系,實線是采用FEM得到的數(shù)值解,數(shù)據(jù)點是采用多級方法得到的解析解Fig.6.The dependencies of(a)the real part of the effective refractive index and(b)the propagation length on the spacing d.The solid lines are numerical solutions obtained by FEM,the data points are analytical solutions obtained by using the multipole method.

對于圖6(b)中高階模式的傳播長度出現(xiàn)的交叉現(xiàn)象,可以根據(jù)場分布得到解釋.以Mode 2和Mode 5為例,圖7(a)—圖7(c)是當間距分別為15,28和45 nm時Mode 2的電場分布,圖7(d)—圖7(f)是當間距分別為15,28和45 nm時Mode 5的電場分布.28 nm是Mode 2和Mode 5達到相同傳播長度所對應的間距.當間距從15 nm變化到28 nm,再變化到45 nm時,對比圖7(a)、圖7(b)和圖7(c)可以看出,場在石墨烯表面的強度明顯減弱,場與石墨烯之間的的相互作用也明顯減弱,波導對模式的束縛性也明顯減弱,導致傳播長度明顯增大.對比圖7(d)、圖7(e)和圖7(f)可以看出,Mode 5對間距變化的依賴性小,傳播長度雖然也增大,但增量比較小.因此出現(xiàn)了圖6(b)中高階模式的傳播長度的交叉現(xiàn)象.

然后,研究中間納米線的半徑ρ2對波導的模式特性的影響.圖8是當ρ0=ρ1=100 nm,d=25 nm,f=35 THz,EF=0.5 eV時,各個模式的有效折射率的實部Re(neff)和傳播長度Lprop隨中間納米線的半徑ρ2變化的關系.可以看出,點與線十分符合.從圖8(a)可以看出,隨著半徑的增大,模式1,模式2和模式5的有效折射率的實部都在增大,模式2的變化最為明顯.從圖8(b)可以看出,當ρ2從25 nm 增大到75 nm時,模式1的傳播長度單調遞減,而模式2的傳播長度先減小后增大,在半徑約為40 nm處存在極小值現(xiàn)象.對于模式5,傳播長度隨著半徑的增大而緩慢增大.另外,模式3和模式4的有效折射率的實部和傳播長度基本不受半徑變化的影響,并且這兩個模式的值也基本相同.在變化的過程中,模式1的有效折射率的實部和傳播長度始終是最大的.

圖7 當ρ0=ρ1=100 nm,ρ2=50 nm,f=35 THz,EF=0.5 eV,(a)—(c)當間距分別為15,28和45 nm時Mode 2的電場分布,(d)—(f)是當間距分別為15,28和45 nm時Mode 5的電場分布Fig.7.The distributions of the electric field of Mode 2 when the spacing d are 15 nm(a),28 nm(b)and 45 nm(c),and the distributions of the electric field of Mode 5 when the spacing d are 15 nm(d),28 nm(e)and 45 nm(f)with ρ0= ρ1=100 nm,ρ2=50 nm,f=35 THz and EF=0.5 eV.

圖8 (a)有效折射率的實部和(b)傳播長度隨半徑ρ2的變化關系,實線是采用FEM得到的數(shù)值解,數(shù)據(jù)點是采用多級方法得到的解析解Fig.8.The dependencies of(a)the real part of the effective refractive index and(b)the propagation length on the radius ρ2.The solid lines are numerical solutions obtained by FEM,the data points are analytical solutions obtained by using the multipole method.

在圖8(b)中,當中間納米線半徑增大時,Mode 2的傳播長度存在極小值現(xiàn)象.這一現(xiàn)象可以根據(jù)場分布得到解釋.圖9給出了當中間納米線半徑分別為(a)25 nm,(b)40 nm和(c)70 nm時Mode 2的電場分布圖.可以看出,當中間納米線半徑由25 nm增加到40 nm時,場在中間納米線的表面上的強度不斷增強,場與石墨烯的相互作用增強,傳播長度減小;當中間納米線半徑由40 nm增加到70 nm時,場在中間納米線的表面上分布開始擴散,對場的束縛性減弱,因此傳播損耗也會降低,使得傳播長度變大.因此在40 nm處,出現(xiàn)了極小值現(xiàn)象.

最后,研究石墨烯的費米能EF對波導的模式特性的影響.基于石墨烯和貴重金屬的表面等離子波導之間的最大區(qū)別是:前者可以通過化學摻雜或施加在石墨烯表面的偏置電壓來改變石墨烯的費米能EF,并最終改變其表面電導率從而可以調控其支持的模式的傳輸特性.而后者不可以進行這些調節(jié).圖10是當ρ0=ρ1=100 nm,ρ2=50 nm,d=25 nm,f=35 THz時,各個模式的有效折射率的實部Re(neff)和傳播長度Lprop隨石墨烯的費米能EF的變化關系.可以看出,點與線十分符合.從圖10(a)可見,隨著費米能從0.4 eV增加到1.2 eV時,5個模式的有效折射率的實部都在減小,即模式束縛性減弱,但是模式1的值始終是最大的.從圖10(b)可見,除了模式4和模式5的傳播長度是先增大后減小,其他三個模式的傳播長度都隨著費米能的增大而增大.在費米能較低時,各個模式的傳播長度間的差異都較小,隨著費米能增大,這些差異在逐漸增大.

圖9 當ρ0=ρ1=100 nm,d=25 nm,f=35 THz,EF=0.5 eV,中間納米線半徑分別為25 nm(a),40 nm(b)和70 nm(c)時,模式2的電場分布Fig.9.The distributions of the electric field of Mode 2 when the radius ρ2are 25 nm(a),40 nm(b)and 70 nm(c)with ρ0= ρ1=100 nm,d=25 nm,f=35 THz and EF=0.5 eV.

圖10 (a)有效折射率的實部和(b)傳播長度隨費米能EF的變化關系,實線是采用FEM得到的數(shù)值解,數(shù)據(jù)點是采用多級方法得到的解析解Fig.10.The dependencies of(a)the real part of the effective refractive index and(b)the propagation length on the Fermi energy EF.The solid lines are numerical solutions obtained by FEM,the data points are analytical solutions obtained by using the multipole method.

圖11 當ρ0=ρ1=100 nm,d=25 nm,f=35 THz,ρ2=50 nm,費米能分別為0.4 eV(a),0.9 eV(b)和1.2 eV(c)時,模式5的電場分布Fig.11.The distributions of the electric field of Mode 5 when the Fermi energy EFare 0.4 eV(a),0.9 eV(b)and 1.2 eV(c)with ρ0= ρ1=100 nm,d=25 nm,f=35 THz and ρ2=50 nm.

當費米能從小到大變化時,Mode 5的傳播長度存在極大值現(xiàn)象.這一現(xiàn)象可以根據(jù)場分布得到解釋.圖11給出了當費米能分別為(a)0.4 eV,(b)0.9 eV和(c)1.2 eV時Mode 5的電場分布圖.對比圖11(a)和圖11(b)可以看出當費米能為0.9 eV時,場的分布范圍增大,場與石墨烯之間的相互作用減弱,結構對模式的束縛性減弱,導致傳播長度增大.對比圖11(b)和圖11(c)可以看出當費米能為1.2 eV時,三根納米線上的場的耦合增強,邊上的兩根納米線上的場相互吸引,導致中間納米線上的場收縮,場的強度增強,場與石墨烯的相互作用增強,損耗增大,傳播距離下降.因此在0.9 eV處出現(xiàn)了極大值現(xiàn)象.

4 結 論

采用多極方法對涂覆石墨烯的三根電介質納米線構成的THz波導的模式特性進行了解析分析.通過改變工作頻率、納米線之間的間距、中間納米線的半徑以及石墨烯的費米能,對5種低階模式的有效折射率的實部和傳播長度進行了研究.結果表明:隨著工作頻率的增大,5種模式的有效折射率的實部都在增大,傳播長度都在減小,并且變化的過程中會出現(xiàn)交叉現(xiàn)象.當納米線之間的間距不斷增大時,除了模式3和模式4基本不受影響以外,模式1、模式2和模式5的有效折射率的實部都減小,傳播長度都增大并最終趨于穩(wěn)定,并且變化的過程中會出現(xiàn)交叉現(xiàn)象.當中間納米線半徑不斷增大時,模式1和模式2的有效折射率的實部增大,模式1的傳播長度單調遞減,模式2的傳播長度先減小后增大,模式3和模式4的有效折射率的實部和傳播長度基本不受影響.對于模式5,有效折射率的實部和傳播長度在緩慢增大.當改變費米能時,每個模式的有效折射率的實部都隨著費米能的增大而減小.除了模式4和模式5的傳播長度是先增大后減小,其他三個模式的傳播長度都隨著費米能的增大而增大.計算表明,多極方法所得結果與有限元法所得結果完全一致.本文研究的基于涂覆石墨烯的三根電介質圓柱形的THz波導有望在模分復用方面得到應用.

[1]Siegel P H 2002IEEE Trans.Microw.Theory50 910

[2]Wang S H,Ferguson B,Zhang C L,Zhang X C 2003Acta Phys.Sin.52 120(in Chinese)[王少宏，B.Ferguson，張存林，張希成2003物理學報52 120]

[3]Chen Q,Tani M,Jiang Z P,Zhang X C 2001J.Opt.Soc.Am.B18 823

[4]Han H,Park H,Cho M,Kim J 2002Appl.Phys.Lett.80 2634

[5]Redo-Sanchez A,Zhang X C 2008IEEE J.Sel.Top.Quant.14 260

[6]Gallot G,Jamison S P,McGowan R W,Grischkowsky D 2000J.Opt.Soc.Am.B17 851

[7]Kawase K,Mizuno M,Sohma S,Takahashi T,Taniuchi T,Urata Y,Wada S,Tashiro H,Ito H 1999Opt.Lett.24 1065

[8]Quema A,Takahashi H,Sakai M,Goto M,Ono S,Sarukura N,Shioda R,Yamada N 2003Jpn.J.Appl.Phys.42 L932

[9]Chen L J,Chen H W,Kao T F,Lu J Y,Sun C K 2006Opt.Lett.31 308

[10]Novoselov K S,Geim A K,Morozov S V,Jiang D,Katsnelson M I,Grigorieva I V,Dubonos S V,Firsov A A 2005Nature438 197

[11]Ju L,Geng B S,Horng J,Girit C,Martin M,Hao Z,Bechtel H A,Liang X G,Zettl A,Shen Y R,Wang F 2011Nature Nanotechnol.6 630

[12]Wang J C,Song C,Hang J,Hu Z D,Zhang F 2017Opt.Express25 23880

[13]Jablan M,Buljan H,Solja?i? M 2009Phys.Rev.B80 245435

[14]He X Y,Kim S 2013J.Opt.Soc.Am.B30 2461

[15]Wang J C,Wang X S,Shao H Y,Hu Z D,Zheng G G,Zhang F 2017Nanoscale Res.Lett.12 9

[16]Donnelly C,Tan D T H 2014Opt.Express22 22820

[17]Christensen J,Manjavacas A,Thongrattanasiri S,Koppens F H L,Abajo F J G 2012ACS Nano6 431

[18]Hajati M,Hajati Y 2016Appl.Opt.55 1878

[19]Wang X S,Chen C,Pan L,Wang J C 2016Sci.Rep.UK6 32616

[20]He S L,Zhang X Z,He Y R 2013Opt.Express21 30664

[21]Gao Y X,Ren G B,Zhu B F,Wang J,Jian S S 2014Opt.Lett.39 5909

[22]Yang J F,Yang J J,Deng W,Mao F C,Huang M 2015Opt.Express23 32289

[23]Xing R,Jian S S 2016IEEE Photon.Tech.L.28 2779

[24]Zhu B F,Ren G B,Yang Y,Gao Y X,Wu B L,Lian Y D,Wang J,Jian S S 2015Plasmonics10 839

[25]Luo L W,Ophir N,Chen C P,Gabrielli L H,Poitras C B,Bergmen K,Lipson M 2014Nat.Commun.5 3069

[26]Yang H B,Qiu M,Li Q 2016Laser Photon.Rev.10 278

[27]Wu X R,Huang C R,Xu K,Shu C,Tsang H K 2017J.Lightwave Technol.35 3223

[28]Nikitin A Y,Guinea F,García-Vidal F J,Martín-Moreno L 2011Phys.Rev.B84 195446

[29]Wijngaard W 1973J.Opt.Soc.Am.63 944

[30]Wijngaard W 1974J.Opt.Soc.Am.64 1136

[31]Huang H S,Chang H C 1990J.Lightwave Technol.8 945

[32]Lo K M,McPhedran R C,Bassett I M,Milton G W 1994J.Lightwave Technol.12 396

[33]White T P,Kuhlmey B T,McPhedran R C,Maystre D,Renversez G,Sterke C M,Botten L C 2002J.Opt.Soc.Am.B19 2322

[34]Kuhlmey B T,White T P,Renversez G,Maystre D,Botten L C,Sterke C M,McPhedran R C 2002J.Opt.Soc.Am.B19 2331