Modified rigorous coupled-wave analysis for grating-based plasmonic structures with delta-thin conductive channel. Far- and near-field study
Yurii M. Lyaschuk, Serhii M. Kukhtaruk, Vytautas Janonis, Vadym V. Korotyeyev
MModified rigorous coupled-wave analysis for grating-based plasmonic structures withdelta-thin conductive channel. Far- and near-field study.
Yurii M. Lyaschuk, Serhii M. Kukhtaruk,
1, 2
Vytautas Janonis, and Vadym V. Korotyeyev ∗ Department of Theoretical Physics, Institute of Semiconductor Physics NAS of Ukraine, pr. Nauki 45, 03028 Kyiv, Ukraine Experimentelle Physik 2, Technische Universit¨at Dortmund, Otto-Hahn-Str. 4a, 44227 Dortmund, Germany Center for Physical Sciences and Technology, Saul˙etekio al. 3, LT-10257 Vilnius, Lithuania
The modified rigorous coupled-wave analysis technique is developed to describe the optical char-acteristics of the plasmonic structures with the grating-gated delta-thin conductive channel in thefar- and near-field zones of electromagnetic waves. The technique was applied for analysis of theresonant properties of AlGaN/GaN heterostructures combined with deeply-subwavelength metallicgrating which facilitates the excitation of the two-dimensional plasmons in the THz frequency range.The convergence of the calculations at the frequencies near the plasmon resonances is discussed. Theimpact of the grating’s parameters, including filling factor and thickness of the grating, on resonantabsorption of the structure was investigated in details. The spatial distributions of electromagneticfield in a near-field zone were used for the evaluation of total absorption of the plasmonic structuresseparating contributions of the grating-gated two-dimensional electron gas and the grating coupler.
I. INTRODUCTION
Nowadays, structures with spatially-periodical lateralstructurization/metasurfaces are in a focus of studies asthe key elements of the many opto- and optoelectron-ics devices with broad application areas, including spec-troscopy, imaging, holography, optical lithography, bio-chemical sensing and military application (for more in-formation, see recent review in Ref. [1]). For example,diffractive gratings are widely utilized in different spec-tral ranges from microwaves to deep ultraviolet as dis-persive optical component [2], antenna elements [3], po-larizers [4], etc.Recently, great attention has been paid to exploitationof the subwavelength metasurfaced structures in THzand Far-infrared spectral ranges. Particularly, the semi-conductor structures with surface-relief gratings [5, 6],quantum well (QW) heterostructures [7, 8] or graphene-based structures [9] incorporated with metallic gratingsare widely discussed as potential efficient emitters anddetectors of the THz radiation [10]. In such structures,the gratings play a role of the coupler, facilitating the res-onant interaction between incident electromagnetic ( em )waves and charge density waves such as surface plasmon-polaritons or 2D plasmons. Moreover, detailed investiga-tions of the resonant properties of grating-based plas-monic structures can serve as additional tool for ba-sic characterization of the 2D electron gas (2DEG) inQWs [11, 12], graphene [13–17] and other novel 2D lay-ered materials [18].Mentioned research is faced with the problem of therigorous electrodynamic simulations of optical character-istics of grating-based structures containing very thin oreven atomically thin conductive layers. In the simula-tions, these layers should be treated as delta-thin withtwo-dimensional parameters such as the sheet conduc- ∗ [email protected] tance of the electron channel. For such structures, thetechnique of integral equations (IE) is often applied forthe solutions of the Maxwell’s equations. This techniqueuses Green function formalism and is based on an reduc-tion of the Maxwell’s system of equations to the linearintegral equations. Latter can be solved, for example, us-ing Galerkin schemes with guaranteed convergence. TheIE technique has been exploited for the different problemsincluding investigations of 2D plasmon instabilities underthe grating [8, 19, 20], detection of THz radiation [21–23]and interaction of THz radiation with conductive-stripgratings [24, 25]. In spite of apparent advantages withrespect to fast convergence, this method is typically for-mulated for modeling structures with simple geometriesand grating is treated as delta-thin.Another powerful method for solution of grating-related electrodynamic problem is known as rigorouscoupled-wave analysis (RCWA) [26–28]. This methodbelongs to the matrix-type methods and operates withsystems of algebraic equations which are formulated forcoefficients of the Fourier expansion of the actual com-ponents of the em field. The RCWA can be applied forarbitrary complexity of the grating-based structures withany geometry of gratings. Typically, the realization of theRCWA relates to the structures where each layer is de-scribed by the bulk parameters. The RCWA method forstructures with delta-thin conductive channels requiresessential modifications which were recently discussed inRefs [29, 30] for plasmonic structures with graphene. Theaim of this paper is to present such modifications in moredetail, focusing on study of convergence of the proposedmethod on example of calculations of optical character-istics of metallic grating-based resonant plasmonic struc-ture with QW in THz frequency range. Effect of thegrating depth on the plasmon resonance, analysis of thenear-field pattern, and comparison to other methods arediscussed. A total absorption of the plasmonic structuresseparating contributions of the grating-gated 2DEG andmetallic grating couplers is also simulated using spatialdistributions of the em field in the near-field region. a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Mathematical formalism of modified RCWA methodis presented in Section II. RCWA will be formulated forplanar diffraction problem (plane of incidence is perpen-dicular to grating strips) including the cases of TM andTE polarizations. The investigations of the convergenceof proposed method and comparison with IE techniquewill be illustrated in Section III on example of the calcu-lations of transmission, reflection and absorption spectraof QW plasmonic structure with deeply subwavelengthmetallic grating. The effect of the finite thickness of themetallic grating will be studied in details. In Section IV,we will perform the analysis of the near-field patternsunder conditions of the plasmon resonances. The mainresults will be summarized in Section V.
II. MATHEMATICAL FORMALISM
Let us assume that multilayered structure with thegrating is illuminated by plane em wave of TM polar-ization with frequency, ω , and angle of incidence, θ (seeFig.1). The grating of period a g is formed by the in-finitely long in y − direction rectangular bars with width, w g and height, h g . The structure consists of N layerswith thicknesses, d j = z j − z j − , j = 1 ...N . Each j − layer,including the grating region, is described by the own di-electric permittivity, (cid:15) ω,j ( x ). The delta-thin conductivechannel is placed between j − j layers and describedby high-frequency conduction current, (cid:126)J D ( x, t ) δ ( z − z j ).The whole structure is set between two non-absorbinghalf-spaces with (cid:15) (at z <
0) and (cid:15) N +1 (at z > D ,where D is the total thickness of the structure).Assuming that all components of em field oscillate intime as exp( − iωt ), the Maxwell’s equations written foramplitudes take the form:rot (cid:126)H ω = − ik (cid:15) ω ( x, z ) (cid:126)E ω + 4 πc (cid:126)J Dω ( x ) δ ( z − z j ) , rot (cid:126)E ω = ik (cid:126)H ω . (1)where k = ω/c , (cid:15) ω ( x, z ) = (cid:15) ω,j ( x ) at z ∈ [ z j − , z j ] andeach non-zero components of the vectors (cid:126)E and (cid:126)H arethe functions of x and z coordinates. For simplicity, weconsider non-magnetic structure with magnetic permit-tivity is equal to 1. For the case of planar diffraction andTM polarization, the non-zero components are E x , E z and H y . Then, Eqs. (1) can be rewritten as: ∂H ω,y ∂z = ik (cid:15) ω ( x, z ) E ω,x − πc J Dω,x ( x ) δ ( z − z j ) ,∂H ω,y ∂x = − ik (cid:15) ω ( x, z ) E ω,z ,∂E ω,x ∂z − ∂E ω,z ∂x = ik H ω,y (2)According the Floquet theorem, we can search for the a g {E N+1 ,H N+1 } z N-1 z N z z .... z xz0 2DEGk E H
E(E x ,0,E z )H(0,H y ,0){E ,H }{E ,H }{E ,H }{E N ,H N } for TM: a g w g h g Figure 1. A schematic sketch of the geometry of multilayeredplasmonic structure with 2DEG. solutions in the form of Fourier expansion: H ω,y ( x, z ) E ω, { x,z } ( x, z ) J Dω,x ( x ) = M (cid:88) m = − M H ω,m,y ( z ) E ω,m, { x,z } ( z ) J Dω,m,x exp( iβ m x ) , (3)where β m = k √ (cid:15) sin θ + q m with q m = 2 πm/a g . Thetruncation rank M (actual number of Fourier harmon-ics) is selected in such way to provide the convergenceof the solution with a given accuracy. Using the expan-sion (3), system (2) written for each j − layer in Fourierrepresentation reads as: ∂ H y,j ∂z = ik ˆ E − ,j E x,j , ˆ β H y,j = − k ˆ E j E z,j ,∂ E x,j ∂z = ik H y,j + i ˆ β E z,j . (4)Here, we introduced the matrix notations: where Fourier-vectors H y,j , E { x,z } ,j contain 2 M + 1 correspondingFourier components, H ω,m,y , E ω,m, { x,z } ; ˆ β is the diagonalmatrix formed by the elements β m δ m,m (cid:48) (here δ denotesthe Kronecker delta symbol). Elements of the matrix ˆ E j and ˆ E inv ,j are expressed through the spatial profile of thedielectric permittivity (cid:15) ω,j ( x )[ ˆ E j ] m,m (cid:48) = (cid:90) (cid:15) ω,j (¯ x ) exp( − πi [ m − m (cid:48) ]¯ x ) d ¯ x, [ ˆ E inv ,j ] m,m (cid:48) = (cid:90) (cid:15) − ω,j (¯ x ) exp( − πi [ m − m (cid:48) ]¯ x ) d ¯ x, (5)respectively. Here ¯ x = x/a g . Note, that the emergenceof the ˆ E − ,j matrix in the first equation of system (4)relates to the Fourier factorization rule [31, 32] of a prod-uct of two discontinues functions (cid:15) ω ( x ) and E ω,x ( x ) withconcurrent jump in the region of the grating (product (cid:15) ω E ω,x is the electrical induction which is continuous inthe x − direction). In the spatially uniform layers, matri-ces ˆ E − ,j and ˆ E j are diagonal and identical. The applica-tion of this rule considerably improves the convergence ofthe results comparing with the old formulation of RCWAmethod[26, 28], where classical Laurent rule, representingthe product (cid:15) ω E ω,x as a convolution-type sum in conven-tional form ( ˆ E j was used instead of ˆ E − ,j ).System (4) can be rewritten in terms of the vector H y,j , ∂ H y,j ∂z = k ˆA j H y,j , where (6) ˆA j = ˆ E − ,j (cid:34) ˆ β ˆ E − j ˆ βk − ˆ I (cid:35) , and ˆ I is the identity matrix. Having H y,j , we can findFourier-vectors of the electric field components: E x,j = − ik ˆ E inv ,j ∂ H y,j ∂z , E z,j = − k ˆ E − j ˆ β H y,j . (7)The system (6) should be solved consequently for each j − layer with appropriate boundary conditions on the j − interfaces: E x,j ( z j ) = E x,j +1 ( z j ) , H y,j ( z j ) − H y,j +1 ( z j ) = 4 πc J Dx,j , (8)where J Dx,j is the Fourier-vector formed by J Dω,m,x com-ponents. The first equation in (8) expresses the conti-nuity of the tangential component of electric field andsecond equation describes the discontinuity of magneticfield component due to the presence of the conductive 2Dlayer [33, 34]. In the frames of the linear response theory, J Dω,m,x = σ Dω,m E ω,m,x , where high-frequency sheet con-ductivity, σ Dω,m , takes into account both frequency andspatial dispersion of the 2DEG. Our approach is valid forarbitrary form of the conductivity of 2DEG. Particularexamples of σ Dω,m for the 2DEG with parabolic spectrumcan be found in Refs. [8], [19]. For the electrons withDirac spectrum, as in the graphene, σ Dω,m can be ob-tained using Kubo formalism (see Refs. [35–37]). In thecase of doped graphene in steady-state applied electricfield, σ Dω,m can be found in Ref. [38].Thus, the second boundary condition in Eqs.(8) canbe rewritten in the form: H y,j ( z j ) − ˆΓ Dj E x,j ( z j ) = H y,j +1 ( z j ) , (9)where diagonal matrix ˆΓ Dj is formed by the elements4 π/c × σ Dω,m δ m,m (cid:48) .The Eqs. (6) compose the system of ordinary secondorder differential equations with constant coefficients.The solution of this system can be expressed in termsof the eigen values, λ m,j , and eigen vectors, (cid:126)w m,j of thematrix ˆA : H y,j ( z ) = M +1 (cid:88) ν =1 (cid:126)w ν,j (cid:2) C + ν,j exp( − k ¯ λ ν,j ( z − z j − ))+ C − ν,j exp( k ¯ λ ν,j ( z − z j )) ] , (10) where ¯ λ ν,j is the square root (with positive real part) ofthe eigenvalues λ ν,j . Two terms in the square bracketsdescribe two waves: transmitted(+) and reflected ( − )into j-layer. The expression for Fourier vector, E x,j ( z )can be obtained by means the first equation of (7) andwritten as follows E x,j ( z ) = M +1 (cid:88) ν =1 (cid:126)v ν,j (cid:2) C + ν,j exp( − k ¯ λ ν,j ( z − z j − )) − C − ν,j exp( k ¯ λ ν,j ( z − z j )) ] , (11)where vector (cid:126)v ν,j = i ˆ E inv ,j ¯ λ ν,j (cid:126)w ν,j .Matching magnetic and electric fields on j − interfaceaccording to the boundary conditions (8) and (9), wecome to the following recurrence relationship betweenconstants of integration, C ± ν,j and C ± ν,j +1 : M +1 (cid:88) ν =1 (cid:126)w − ν,j exp( − k ¯ λ ν,j d j ) C + ν,j + (cid:126)w + ν,j C − ν,j = M +1 (cid:88) ν =1 (cid:126)w ν,j +1 (cid:2) C + ν,j +1 + C − ν,j +1 exp( − k ¯ λ ν,j +1 d j +1 ) (cid:3) ; M +1 (cid:88) ν =1 (cid:126)v ν,j (cid:2) exp( − k ¯ λ ν,j d j ) C + ν,j − C − ν,j (cid:3) = M +1 (cid:88) ν =1 (cid:126)v ν,j +1 (cid:2) C + ν,j +1 − C − ν,j +1 exp( − k ¯ λ ν,j +1 d j +1 ) (cid:3) , (12)where (cid:126)w ± ν,j = (cid:126)w ν,j ± ˆΓ Dj (cid:126)v ν,j contains parameters of2DEG entered into matrix ˆΓ Dj . The Eqs.(12) can bewritten in the compact matrix form: (cid:18) ˆW − j ˆΛ j , ˆW + j ˆV j ˆΛ j , − ˆV j (cid:19)(cid:32) (cid:126)C + j (cid:126)C − j (cid:33) = (cid:18) ˆW j +1 , ˆW j +1 ˆΛ j +1 ˆV j +1 , − ˆV j +1 ˆΛ j +1 (cid:19)(cid:32) (cid:126)C + j +1 (cid:126)C − j +1 (cid:33) (13)where each vector, (cid:126)C ± j , contains 2 M + 1 integrationconstants, the matrixes ˆW j and ˆV j are formed by el-ements of the corresponding eigen vectors, (cid:126)w ν,j and (cid:126)v ν,j ,defined for all eigen values, λ ν,j . Matrixes ˆW ± j = ˆW j ± ˆΓ Dj ˆV j , ˆΛ j is the diagonal matrix with the ele-ments, exp( − k ¯ λ ν,j d j ) δ ν,ν (cid:48) .Relationship (13) allows us to couple amplitudes ofthe reflected wave (in region z < z ) and transmittedwave (in the region z > z N ), and consequently to calcu-late transmission and reflection coefficients for differentdiffraction orders. Indeed, components of the Fourier-vectors of y − magnetic and x − electric fields in the region z < j = 0) are: H ω,m,y = δ m, exp( − ¯ λ m, k z ) + r ω,m exp(¯ λ m, k z ) , (14) E ω,m,x = cos θ √ (cid:15) δ m, exp( − ¯ λ m, k z ) − i ¯ λ m, (cid:15) r ω,m exp(¯ λ m, k z ) , and in the region z > D ( j = N + 1) are: H ω,m,y = t ω,m exp( − ¯ λ m,N +1 k ( z − z N )) , (15) E ω,m,x = i ¯ λ m,N +1 (cid:15) N +1 t ω,m exp( − ¯ λ m,N +1 k ( z − z N )) , where ¯ λ m, { ,N +1 } = (cid:113) β m − (cid:15) { ,N +1 } k /k if β m > √ (cid:15) { ,N +1 } k and ¯ λ m, { ,N +1 } = − i (cid:113) (cid:15) { ,N +1 } k − β m /k if otherwise. The quantities r ω,m and t ω,m are the nor-malized magnetic-field amplitudes of the m-th backward-diffracted (reflected) and forward diffracted (transmit-ted) waves, respectively. These amplitudes form theFourier-vectors R and T .Using Eqs. (14), (15), boundary conditions (8) and(9), and relationship (13), we found that Fourier-vectors R and T are coupled through following matrix equations: (cid:32) ˆI − cos θ √ (cid:15) ˆΓ D θ √ (cid:15) ˆI (cid:33) (cid:126)δ + (cid:18) ˆI − ˆ Z I ˆΓ D ˆ Z I (cid:19) R = N (cid:89) j =1 (cid:18) ˆW j , ˆW j ˆΛ j ˆV j , − ˆV j ˆΛ j (cid:19)(cid:18) ˆW − j ˆΛ j , ˆW + j ˆV j ˆΛ j , − ˆV j (cid:19) − (cid:18) ˆI ˆ Z II (cid:19) T , (16)where matrices ˆ Z I and ˆ Z II are diagonal with elements − i ¯ λ m, /(cid:15) δ m,m (cid:48) and i ¯ λ m,N +1 /(cid:15) N +1 δ m,m (cid:48) , respectively, (cid:126)δ is the vector with elements δ m, . The (16) is the mastersystem of equations of the modified RCWA method pro-viding the way for calculation of transmission and reflec-tion coefficients for different diffraction orders. In con-trast to the previous formulation of the RCWA method[27, 28], account of 2DEG leads to the nontrivial modi-fications. Particularly, the master system (16) containsthe matrices ˆW ± j in the right-hand side and ˆΓ D in theleft-hand side. Latter term describes the possible exis-tence of the 2DEG on the top of the grating.It should be noted that usage of this system of equa-tions, written in the present form, applying to the situa-tion of deep surface grating and optically dense materialscan face with problem of the computational instability.This instability is associated with the procedure of thenumerical inversion of the second matrix in the r-h-s ofthe (16) when exponentially small terms, exp( − k ¯ λ ν,j d j )(standing in the ˆΛ j ) becomes smaller than machine preci-sion. To avoid this obstacle, authors in Ref.[28] proposedto use the following decomposition of the badly invertedmatrix: (cid:18) ˆW − j ˆΛ j , ˆW + j ˆV j ˆΛ j , − ˆV j (cid:19) − = (cid:18) ˆΛ j , ˆ0ˆ0 , ˆI (cid:19) − (cid:18) ˆW − j , ˆW + j ˆV j , − ˆV j (cid:19) − (17)Note, that second inversion matrix in the r-h-s of Eq. (17)is the regular and can be numerically inverted withoutany difficulties. Let’s the product of the last N − termsin Eq. (16) (cid:18) ˆW N , ˆW N ˆΛ N ˆV N , − ˆV N ˆΛ N (cid:19) (cid:18) ˆΛ N , ˆ0ˆ0 , ˆI (cid:19) − (cid:18) ˆX N ˆY N (cid:19) T ≡ (cid:18) ˆf N ˆ g N (cid:19) T (18) where we introduced the following designation (cid:18) ˆX N ˆY N (cid:19) ≡ (cid:18) ˆW − N , ˆW + N ˆV N , − ˆV N (cid:19) − (cid:18) ˆI ˆ Z II (cid:19) . Making substitution, T = ˆX − N ˆΛ N T N , term (cid:18) ˆΛ N , ˆ0ˆ0 , ˆI (cid:19) − (cid:18) ˆX N ˆY N (cid:19) T = (cid:18) ˆΛ N , ˆ0ˆ0 , ˆI (cid:19) − (cid:18) ˆΛ N ˆY N ˆX − N ˆΛ N (cid:19) × T N = (cid:18) ˆΛ N , ˆ0ˆ0 , ˆI (cid:19) − (cid:18) ˆΛ N , ˆ0ˆ0 , ˆI (cid:19)(cid:18) ˆIˆY N ˆX − N ˆΛ N (cid:19) T N = (cid:18) ˆIˆY N ˆX − N ˆΛ N (cid:19) T N , and we can obtain that (cid:18) ˆf N ˆ g N (cid:19) T = ˆW N (cid:104) ˆI + ˆΛ N ˆY N ˆX − N ˆΛ N (cid:105) ˆV N (cid:104) ˆI − ˆΛ N ˆY N ˆX − N ˆΛ N (cid:105) T N (19)Substituting relationship (19) into Eqs.(16) and sequen-tially performing above mentioned transformations foreach j -th term in the product we can rewrite master sys-tem of the equations (16) in the computationally stableform: (cid:32) ˆI − cos θ √ (cid:15) ˆΓ D θ √ (cid:15) ˆI (cid:33) (cid:126)δ + (cid:32) ˆI − cos θ √ (cid:15) ˆΓ D ˆ Z I (cid:33) R = (cid:18) ˆf ˆg (cid:19) T . (20)The matrices ˆf and ˆg can be found from the followingrecurrence relationship: (cid:18) ˆf j − ˆg j − (cid:19) = ˆW j − (cid:104) ˆI + ˆΛ j − ˆY j − ˆX − j − ˆΛ j − (cid:105) ˆV j − (cid:104) ˆI − ˆΛ j − ˆY j − ˆX − j − ˆΛ j − (cid:105) , (21)where (cid:18) ˆX j − ˆY j − (cid:19) = (cid:18) ˆW − j − , ˆW + j − ˆV j − , − ˆV j − (cid:19) − (cid:18) ˆf j ˆg j (cid:19) (22)and j is varied from N to 2. The vector T relates toFourier-vector of transmission coefficient, T , as follows T = (cid:89) j = N ˆX − j ˆΛ j T . (23)Components of the vectors T and R provide the trans-mission, T m , reflection, R m , coefficients for any m-thdiffraction order as well as total absorption, L . Particu-larly, T m = √ (cid:15) (cid:112) (cid:15) N +1 k − β m k (cid:15) N +1 cos θ | T [ m ] | ,R m = (cid:112) (cid:15) k − β m k √ (cid:15) cos θ | R [ m ] | , L = 1 − (cid:88) m T m + R m . (24)Here, the summation is taken over all numbers of visi-ble diffraction orders, i.e, over such values of m whichkeep the positive expressions under the square roots inthe nominators. In the case of subwavelength gratings,only zero diffraction order ( m = 0) occurs and we havethat T = | T [0] | , R = | R [0] | (if (cid:15) = (cid:15) N +1 ), and L ≡ L = 1 − T − R . The Eqs. (20-23) togetherwith Eqs. (24) finalize the computationally stable realiza-tion of the modified RCWA method for characterizationof the multi-layered plasmonic structures with grating-gated delta-thin conductive channels in the case of TM-polarization.This method can be extended for the case of TE po-larization of the incident radiation. Now, actual com-ponents of the em waves are: E y , H x and H z . Matrixequation (6) is formulated for the Fourier-vectors E y,j ,where ˆA j = ˆ β /k − ˆ E j . (25)The Fourier-vectors of the magnetic field components aregiven as follows: H x,j = ik ∂ E y,j ∂z , H z,j = ˆ βk E y,j and boundary conditions reads as E y,j ( z j ) = E y,j +1 ( z j ) , H x,j ( z j )+ˆΓ Dj E y,j ( z j ) = H x,j +1 ( z j ) . Finally, master system of equations (16) takes the form: (cid:18) ˆI ˆΓ D − √ (cid:15) cos θ ˆI (cid:19) (cid:126)δ + (cid:18) ˆI ˆΓ D + ˆ Z I (cid:19) R = N (cid:89) j =1 (cid:18) ˆW j , ˆW j ˆΛ j ˆV j , ˆV j ˆΛ j (cid:19)(cid:18) ˆW j ˆΛ j , ˆW j ˆV + j ˆΛ j , ˆV − j (cid:19) − (cid:18) ˆI ˆ Z II (cid:19) T , (26)where ˆW j is formed by eigen vectors, (cid:126)w ν,j obtainedfor the each λ ν,j eigen values of the matrix (25). Ma-trix ˆV j contains vectors, (cid:126)v ν,j = − i ¯ λ ν,j (cid:126)w ν,j and ˆV ± j = ˆV j ± ˆΓ Dj ˆW j . The matrices ˆ Z I and ˆ Z II are the diagonalwith elements i ¯ λ m, δ m,m (cid:48) and − i ¯ λ m,N +1 δ m,m (cid:48) , respec-tively. Applying procedure (see above) for stable compu-tation of system of equations (26), we can find vectors T and R , and calculate transmission, reflection coefficientsfor any m-th diffraction order as well as total absorptionfor the case of TE-polarized incident radiation: T m = (cid:112) (cid:15) N +1 k − β m k √ (cid:15) cos θ | T [ m ] | , R m = (cid:112) (cid:15) k − β m k √ (cid:15) cos θ | R [ m ] | ,L = 1 − (cid:88) m T m + R m . (27)The proposed modified RCWA method is appliedfor investigation of particular plasmonic structures withgrating-gated 2DEG channel. Such structures possess theresonant properties in the THz frequency range (wave-length of order of 100 µ m) for TM-polarized incident radiation due to excitations of plasmons in conductivechannel of 2DEG [7, 11, 19, 39]. Below, we will studyspectral characteristics of the AlGaN/GaN-based plas-monic structure with deeply subwavelength (micron pe-riod) metallic grating, calculating transmission ( T ), re-flection ( R ) and absorption ( L ) coefficients, their con-vergence vs number of the Fourier harmonics and thenear-field mapping. Also, we will pay attention to thedependence of the plasmon resonances vs geometry ofthe grating. III. FAR-FIELD CHARACTERISTICS ANDTHEIR CONVERGENCE
Here and below, we will study the case of TM-polarizedincident em wave with incidence angle, θ = 0. The struc-ture under test is formed by N = 3 media, embeddedinto the air, including the rectangular metallic gratingwith (cid:15) ω, ( x ) = (cid:15) ω,M Θ( w g − x ) + (cid:15) Θ( x − w g ) (whereΘ( x ) stands the Heaviside step function, x ∈ [0 , a g ], (cid:15) ω,M = 1 + 4 πiσ M /ω and σ M = 4 × s − that cor-responds to the gold), AlGaN barrier and GaN bufferlayers with constant dielectric permittivities (cid:15) = 9 . (cid:15) = 8 .
9, respectively. The 2D conductive channelis formed in the plane z = z . The matrix [Γ Dj ] m,m (cid:48) =4 π/c × σ Dω,m δ m,m (cid:48) δ j, , where for description of the high-frequency properties of 2DEG we used Drude-Lorentzmodel σ Dω,m = e n D τ D /m ∗ (1 − iωτ D ) with electroneffective mass, m ∗ = 0 . × m e , concentration of 2DEG, n D = 6 × cm − and effective scattering time, τ D = 0 . . . . em wave withthe plasmons in the channel of the 2DEG. At the res-onance, plasmon excitation of 2DEG with wavevectorsdetermined by the grating period can effectively absorbenergy of the incident em wave. Physics of 2D plasmonsand their resonant interaction with em radiation are well-described in Refs. [7, 19, 39–41].The considered plasmonic structure can provide con-siderable absorption of the THz radiation. For example,at 1.7 THz the absorption coefficient L ∼ h g =5.0 m L /2 , THz (e) /2 =1.7 THz T M M M M /2 =3.5 THz T M(d)
Figure 2. (color online): Spectra of the transmission (a),reflection (b) and absorption (c) coefficients for zero diffrac-tion order calculated at three values of the grating depth: h g ≡ d = 0 . , , µ m. Grating period, a g = 1 µ m andwidth of grating bars, w g = 0 . µ m. Thickness of AlGaNbarrier, d = 0 . µ m and GaN buffer layer, d = 1 µ m.Black dashed lines are the results of the IE method assum-ing delta-thin grating with 2D conductivity, σ Dg = 2 × cm/s. Dash-dotted lines are the results for the structure with-out 2DEG. All spectra are obtained at M = 100. Panels (d)and (e): Dependencies of the transmission coefficient, T vsnumber of Fourier harmonics, M , for two selected frequencies.Insets: relative errors, δ M , are plotted in logarithmic scale asfunction of M . ing with h g = 0 . µ m almost coincide. An increase ofthe grating depth only leads to the additional dispersionof transmission and reflection coefficients in the higherfrequency range. This dispersion is also observed in mod-eling structure without 2DEG (see dashed-dotted lines).Very weak dependence of the absorption of the plasmonicstructure vs grating depth (even for deep grating with h g = 5 µ m) indicates that subwavelength highly con-ductive grating plays the role of almost-lossless waveg-uide for incident TM-polarized wave of THz frequencies.For example, absorption of the deep grating in the struc-ture without 2DEG does not exceed 2% in the consideredspectral range. The calculations of the partial losses as-sociated with grating and 2DEG gas in the grating-2DEGplasmonic structure can be performed using the pattern h g =5.0 m R L /2 , THz (a)(b)(c)(e) T (d) /2 =1.38 THz0 20 40 60 80 100 120 1400.30.60.9 M M /2 =3.78 THz T M Figure 3. (color online): The same as in Fig.2 at w g = 0 . µ m. of the near-field and one will be done below in the Sec.IV.The convergence of the RCWA method vs number ofFourier harmonics, M, is illustrated in Fig.2(d) and (e)on example of T coefficient. The proposed method pro-vides fast convergence and results with reasonable accu-racy can be already obtained using M ≈ −
50 for allconsidered cases in Fig.2(d) and (e). In order to quantifythe convergence of the RCWA method, we introduce rel-ative error defined as follows δ M = | T ( M )0 − T (200)0 | /T (200)0 (see insets in Fig.2(d) and (e)). For example, for shallowgrating (black circles), δ = 0 .
2% for resonant frequency1.7 THz and δ = 0 .
01% for the frequency 3.5 THz. Con-vergence becomes worse for the deep gratings: at reso-nant frequency 1.7 THz δ = 0 .
5% (for h g = 1 µ m) and δ = 1 .
3% (for h g = 5 µ m); at non-resonant frequency3.5 THz δ = 0 .
1% (for h g = 1 µ m) and δ = 0 .
2% (for h g = 5 µ m). Thus, estimations show that convergenceof the RCWA method exhibits dependence on gratingdepth and frequency of the incident radiation. The casesof the deep gratings and resonant frequencies of the plas-mon excitation require account of the larger numbers ofFourier harmonics.In the case of the plasmonic structure with narrow-slitgrating, w g = 0 . µ m and a g = 1 µ m, our calculationspredict much more pronounced features in the opticalcharacteristics including intensity of the plasmon reso- Figure 4. Spatial distributions of amplitudes of H y (a), E x (b) and E z (c) components of em field in the units of amplitude ofincident wave, E ins , for the structure with a g = 1 µ m, w g = 0 . µ m, h g = 0 . µ m at frequency ω/ π = 1 . M = 450. nances vs grating depth (see Figs.3). Moreover, narrow-slit grating provides more efficient coupling between in-cident radiation and plasmon excitations that leads toan emergence of well-pronounced multiple plasmon res-onances which are red-shifted in comparison to the pre-vious case. The red-shift of the resonant frequency isthe result of a larger contribution of the gated region of2DEG where phase velocity of the plasmons is smallerthan in the ungated region of 2DEG [42].The first plasmon resonance occurs at frequency of 1 . ∼ L at resonant frequency of 3 .
78 THz is decreasedfrom 28% for shallow grating ( h g = 0 . µ m) to 15% forthe deepest grating ( h g = 5 µ m). Apparently, this effectrelates to an essential increase of the reflectivity of theplasmonic structures with thicker gratings as shown inFig.3(b). It means that for the deeper gratings, a smallerportion of the em energy is concentrated in 2DEG as itwill be further illustrated in Section IV.It should be noted that application of the RCWA meth-ods for accurate calculations of far-field spectral charac-teristics of plasmonic structure with narrower-slit gratingrequires larger number of Fourier harmonics (see Fig.3(d)and (e)). Now, the relative errors δ for resonant fre-quency 1 .
38 THz are equal to 0 . , .
2% and 2 .
5% for h g = 0 . , , µ m, respectively. The relative errors lessthan 1% for deep gratings is achievable at M >
60. Simi-larly to the previous case, the convergence of the RCWAmethod is improved at higher frequencies. So, at fre-quency of 3 .
78 THz accuracy of computation with rela-tive errors δ M <
1% is achieved at
M >
30. All spectrashown in Figs. 2 and 3 are obtained at M = 100. IV. NEAR-FIELD STUDY
Together with calculations of optical characteristics re-lating to the far-field, RCWA method allows us to studygeometry of the near-field. Especially, we will pay at-tention to the spatial distributions of the E x ( x, z ) and E z ( x, z )-components of the em fields. Absolute values ofthese components determine the local absorption of the em wave and can be used for extraction of partial lossesin metallic grating and 2DEG. Having RCWA data onFourier vectors T (taking from the solutions of mastersystem (16)) we can find vectors of the integration con-stants (cid:126)C ± j in the each j-layer of the structure: (cid:32) (cid:126)C + j (cid:126)C − j (cid:33) = (cid:18) ˆIˆY j ˆX − j ˆΛ j (cid:19) T j , (28)where T j can be found recurrently T j +1 = ˆX − j ˆΛ j T j .Substituting found constants into Eq. (10) (with known ˆW j and ˆΛ j matrices), we can calculate Fourier vec-tors H y,j ( z ) and reconstruct a spatial distribution of the H y -component in the each { x, z } point inside plasmonicstructures using Eq. (3). The Eqs. (14) and (15) are usedfor reconstruction of the near-field distribution of H y and E x components outside the plasmonic structures.The spatial distribution of | H y ( x, z ) | for particular caseof the plasmonic structure with shallow grating is shownin Fig.4(a). This component is tangential to the gratingsides and exhibits smooth behavior with partial pene-tration into the grating bar. Calculations give that skindepth, δ = c/ √ πσ M ω , of the gold at frequency of 1 . µ m which is comparable with heightof the grating bar. The cold zone of H y component oc-cupies middle region of the bar near the bottom face. Inthe plane of 2DEG, | H y ( x, z ) | is a discontinuous quantityaccording boundary conditions (8).The electric components E x (panel(b)) and E z (panel(c)) show more interesting behavior with highlynon-uniform distributions. The E z ( x, z ) can be directlyobtained from second relationship in Eqs. (7) (using al-ready found Fourier vectors, H y,j ( z )) and Eq. (3). For Figure 5. The same as in Fig.4 for h g = 1 µ m. the correct reconstruction of the E x component, we fol-low the method discussed in Refs.[43, 44]. In the re-gion of the grating, z ∈ [0 , h g ], x-component of theelectric field is the normal to the grating bar’s sidesand one has a discontinuity. It is more effective to re-construct a continuous quantity, the component of thedisplacement field, D x ( x, z ), which can be easily calcu-lated from the derivative of the H y component with re-spect to z-coordinate (see first equation in (2)). Then E x ( x, y ) = D x ( x, z ) /(cid:15) ( x, z ), where dielectric permittiv-ity (cid:15) ( x, z ) is the known discontinuous function. Suchmethod allows us partially avoid an emergence of theunphysical spurious oscillations, known as Gibb’s phe-nomenon. Nevertheless, reconstruction of the near-fieldpatterns requires account of the much more Fourier har-monics than for calculations of the far-field characteris-tics. This circumstance was discussed in Ref. [44].As seen, both E x and E z components demonstrate thefield concentration effect. The energy of the em field ismainly concentrated near the ridges (hot zone I) of themetallic bars and in the region between grating bars and2DEG (hot zone II). In the hot zones (I) and (II) bothelectric components are essentially enhanced. In the hotzone (II), E z component predominantly dominates. Thespecific formation of the cold zone for E z component atthe vertical axis x = w g / x = ( a g + w g ) / E x component mainlypenetrates to the grating’s bars from the upper and backfaces and E z from the side faces as a tangential ones forcorresponding faces.The similar geometry of the near-fields is realized forthe case of the deep grating (see Figs. 5). The map-pings of the E x and E z components show that incidentwave passes through the subwavelength grating in theform of TEM mode, i.e in the grating slit, the wave havepredominantly polarization along x-direction with almostconstant amplitude.Additionally, we used COMSOL Multiphysics (cid:114) [45] to validate independently the obtained results by finite ele-ment method. The Wave Optics module [45] is used tosolve Maxwell equations for the system, which is shown inFig. 1. The 2DEG was introduced as the surface currentdensity at the interface of AlGaN and GaN. The uniformquadratic mesh with 5 nm size is used to resolve near-field components (see Fig. 6). The COMSOL’s results ofthe electric components distributions for the case of theshallow grating with h g = 0 . µ m are shown in Fig. 6.Excellent agreement between the modified RCWA andfinite element methods is demostrated in both far- andnear field studies. Figure 6. The same as in Fig.4 obtained by usage of COMSOLMultiphysics.
The spatial distributions of the | E x | ( x, z ) and | E z | ( x, z ) can be used for calculation of the partial lossesin the grating, L gr and 2DEG, L DEG : L gr = 4 πσ M c (cid:82) w g (cid:82) h g dxdz | E x | + | E z | (cid:82) a g dx | E ins | (29) h g =5 m (b) | E x ( x , z ) | (a) gated region0.00 0.25 0.50 0.75 1.000123 /2 =3.5 THz x/a g /2 =1.7 THz Figure 7. (color online): Distribution of | E x ( x, z ) | on onespatial period of the plasmonic structure with w g = 0 . µ mand a g = 1 µ m at two resonant frequencies. and L DEG = 4 π Re[ σ Dω ] c (cid:82) a g dx | E x ( x, z ) | (cid:82) a g dx | E ins | . (30)where E ins is the amplitude of the incident wave.For the case in Fig.4, we obtained that L gr = 0 . . E x -componentand 0 .
11% contribution of E z -component. Losses in2DEG are considerably larger, L DEG = 37 . L = L DEG + L gr =38 . L from far-field characteris-tics gives the almost same value 38 . L gr = 0 . . .
49% contributions for E x and E z components,respectively) with almost same value of absorption in2DEG L DEG = 37 . L = 38 . .
58% obtainedfrom far-field characteristics.Calculations of the partial losses at the frequency ofthe 1-st order plasmon resonance indicate that incident em wave is mainly absorbed by 2DEG and this absorp-tion weakly depends on thickness of grating bars. Thisfact is illustrated by the spatial distribution of the ampli-tude of the x − component of the electric field, | E x ( x, z ) | ,in the plane of 2DEG, calculated at the frequencies of the1-st (Fig. 7(a)) and 2-nd (Fig. 7(b)) plasmon resonancesat three values of the grating depth. As seen, all three(grey, red, green) curves for lower frequency ( ω/ π = 1 . em energyis absorbed in the gated region i.e. under metallic strip.For the higher frequency ( ω/ π = 3 . | E x ( x, z ) | quantity (Fig.7(b)) ac-quires more complicated form with several spatial oscil-lations in the gated region. At this, the effect of thegrating thickness becomes visible, i.e, the deep gratingsstarts to screen the interaction of em wave with 2DEG. /2 =3.78 THz/2 =1.38 THz0.00 0.25 0.50 0.75 1.0002468 (b) | E x ( x , z ) | x/a g (a) Figure 8. (color online): The same as in Fig.7 for w g = 0 . µ m. The emergence of the several spatial oscillations in distri-bution of | E x ( x, z ) | leads to suppression of absorptivityof the plasmonic structures at higher order plasmon reso-nances. Also, according Drude model, at higher frequen-cies response of electron gas on em wave becomes weakerwhich leads to a decrease of the prefactor standing inEq.30. This prefactor, 4 π Re[ σ Dω ] /c , for two consideredfrequencies ω/ π = 1 . ω/ π = 3 . .
049 and 0 . h g = 0 . , , µ m, L DEG = 4 . , . . L = 4 . , . , . em radi-ation. The distributions of | E x ( x, z ) | calculated for thestructure with w g = 0 . µ m at two resonant frequencies ω/ π = 1 .
38 THz (1-st order plasmon resonance) and ω/ π = 3 .
78 THz (3-rd order plasmon resonance) areshown in Figs. 8. As seen, the geometry of the distri-butions obtained for the frequency of 1-st order plasmonresonance (Fig.8(a)) is similar to the previous case de-picted in Fig.7(a). However, the wider gated region of2DEG integrally provides larger contribution to the ab-sorption of em wave by 2DEG. The corresponding val-ues of L DEG are following: 47 .
5% (for h g = 0 . µ m),46 .
6% (for h g = 1 µ m) and 41 .
2% (for h g = 5 µ m). Atthis, L = 48 . , . , . h g =5 µ m essentially suppresses theplasmon absorption of em wave. The corresponding val-ues of L DEG are following: 27 .
2% (for h g = 0 . µ m),24 .
1% (for h g = 1 µ m) and 11 .
8% (for h g = 5 µ m). Atthis, L = 28 . , . , . V. SUMMARY
We have developed computationally stable RCWAmethod for solution of Maxwell’s equations in the caseof the multi-layered plasmonic structures with delta-thingrating-gated conductive channel. The method was for-mulated for planar diffraction problem for both TM andTE polarization of incident wave. Method was imple-mented for investigation of far- and near-field charac-teristics of the particular plasmonic structures based onAlGaN/GaN heterostructure with deeply subwavelengthmetallic grating coupler.The calculations of the far-field characteristics includ-ing transmission, reflection and absorption coefficientsfor zero diffraction order were performed in THz fre-quency range where considered structure has multipleresonances related to the excitations of 2D plasmons inconductive channel of AlGaN/GaN heterostructure. Thedependence of these characteristics vs grating parametersand their convergence vs number of the Fourier harmon-ics were analyzed.We found that spectra of transmission and reflectioncoefficients in the lower frequency range, 0 .. h g /a g = 0 . µ m/1 µ m) and deep ( h g /a g = 1,5 µ m/1 µ m) gold grating are almost identical and co-incide with the results of IE method where grating istreated as delta-thin. In higher frequency range, 3 .. h g becomes more pronounced fornarrower-slit grating with w g /a g = 0 . µ m/1 µ m thanfor wide-slit grating with w g /a g = 0 . µ m/1 µ m. Weshowed that convergence of the calculations depends ongeometrical parameters of the grating and frequencies.Better convergence is achieved for shallow and wide-slitgrating with relative errors of ∼ . − .
2% (in depen-dence on frequency) with 30 Fourier harmonics. For deepand narrow-slit grating, relative errors of ∼ . ..
1% isachieved at M ∼ em wave’s components inthe near field-zone reveals the main physical peculiaritiesof the interaction of the plasmonic structure with inci-dent radiation. It was shown that subwavelength metal-lic grating plays a role of perfect waveguide for incidentwave, concentrator of the em energy and polarizationrotator. In the region of the grating slit em wave haspredominantly lateral polarization with amplitude closeto amplitude of incident wave. The hot zone is formed inregion between grating bars and 2DEG where em wavehas predominantly vertical polarization with amplitudesthat can in 100 times exceed the amplitude of incidentwave.The pattern of the near-field also was used for the cal- culations of the partial losses related to the grating and2DEG. It was shown that at the frequencies of the plas-mon resonances the structure can efficiently absorb THzradiation with absorption coefficient values in the orderof 20 −
50 % (in dependence of the grating filling fac-tor and order of plasmon resonance). We found thatcontribution of 2DEG to the total losses is dominant atlow-frequency plasmon resonances with weak dependenceon the grating depth. At high-frequency plasmon reso-nances, the effect of the grating depth becomes essential.The deep gratings can effectively screen interaction of the em waves with plasmon oscillation in 2DEG that leadsto an decrease of the total absorption of THz radiation.Also, it should be noted that the proposed modifiedRCWA has several advantages over conventional vol-umetric RCWA. First one, our realization of RCWAmethod allows us to avoid additional numerical manipu-lation with matrices that can reduce the computationaltime. For considered structure, we have 25% in termof computation time savings in comparison with conven-tional RCWA at the volumetric treatment of conductivelayer. This value can be increased in simulation of struc-tures with the stack of 2D conductive layers. Second one,we operate with one parameter, two-dimensional concen-tration, n D , instead of two independent parameters ofbulk concentration, n D , and thickness of the layer, d . Itcan be convenient for metrology of the structures at theprocessing of the experimental data.We suggest that proposed modified RCWA algorithmcan be effectively used for the modeling of the opticalcharacteristics of various kinds of plasmonic structureswith 2D conductive channels, including quantum wells-or graphene-based structures and results of the paperprovide deeper insight on physics of the interaction ofTHz radiation with grating-gated plasmonic structures. VI. FUNDING
The work was supported by the Research Coun-cil of Lithuania (Lietuvos mokslo taryba) under the”KOTERA-PLAZA” Project (Grant No. DOTSUT-247)funded by the European Regional Development Fundaccording to the supported activity ”Research ProjectsImplemented by World-class Researcher Groups” un-der Measure No. 01.2.2-LMT-K-718-0047. SMK wassupported by the Bundesministerium f¨ur Bildung undForschung through the project VIP+ ”Nanomagnetron”.
VII. ACKNOWLEDGMENTS
Authors thanks to Prof. V. A. Kochelap (ISP NASU,Ukraine) and Dr. I. Kaˇsalynas (FTMC, Lithuania) forfruitful discussions of the various aspects of this work.
Disclosures.
The authors declare no conflicts of inter-est.1 [1] E. Popov ed. “Gratings: Theory and Numeric Applica-tions”, First Edition, Presses universitaires de Provence(PUP), (2012).[2] W. Neumann, “Fundamentals of Dispersive Optical Spec-troscopy Systems”, SPIE Press, Bellingham, Washing-ton, USA (2014).[3] H. F. Hammad, Y. M. M. Antar, A. P. Freundorfer andM. Sayer, “A new dielectric grating antenna at millimeterwave frequency,” IEEE Transactions on Antennas andPropagation, , 36-44, (2004).[4] S. Shena, Y. Yuana, Z. Ruana, and H. Tan, “Optimizingthe design of an embedded grating polarizer for infraredpolarization light field imaging,” Results in Physics, ,21-31 (2019).[5] G. A. Melentev, V. A. Shalygin, L. E. Vorobjev, V. Yu.Panevin, D. A. Firsov, L. Riuttanen, S. Suihkonen, V. V.Korotyeyev, Yu. M. Lyaschuk, V. A. Kochelap, and V.N. Poroshin, “Interaction of surface plasmon polaritonsin heavily doped GaN microstructures with terahertz ra-diation,” J. Appl. Phys. , 093104 (2016).[6] V. Janonis, S. Tumenas, P. Prystawko, J. Kacperski,and I. Kaˇsalynas, “Investigation of n-type gallium nitridegrating for applications in coherent thermal sources,”Appl. Phys. Lett. , 112103 (2020).[7] V. V. Popov, D. V. Fateev, O. V. Polischuk, and M. S.Shur, “Enhanced electromagnetic coupling between ter-ahertz radiation and plasmons in a grating-gate transis-tor structure on membrane substrate,” Opt. Express ,16771 (2010).[8] V. V. Korotyeyev, V. A. Kochelap, S. Danylyuk, and L.Varani, “Spatial dispersion of the high-frequency conduc-tivity of two-dimensional electron gas subjected to a highelectric field: Collisionless case,” Appl. Phys. Lett. ,041102 (2018).[9] V. Ryzhii, T. Otsuji, and M. Shur, “Graphene basedplasma-wave devices for terahertz applications”, Appl.Phys. Lett. , 140501 (2020).[10] T. Otsuji and M. Shur, “Terahertz Plasmonics: Good Re-sults and Great Expectations,” IEEE Microw. Mag. ,43-50 (2014).[11] D. Pashnev, T. Kaplas, V. Korotyeyev, V. Janonis, A.Urbanowicz, J. Jorudas and I. Kaˇsalynas, “Terahertztime-domain spectroscopy of two-dimensional plasmonsin AlGaN/GaN heterostructures,” Appl. Phys. Lett. ,051105 (2020).[12] D. Pashnev, V. Korotyeyev, V. Janonis, J. Jorudas, T.Kaplas, A. Urbanowicz, and I. Kaˇsalynas, “Experimen-tal evidence of temperature dependent effective mass inAlGaN/GaN heterostructures observed via THz spec-troscopy of 2D plasmons” Appl. Phys. Lett. , 162101(2020).[13] Bo Yan, Jingyue Fang, Shiqiao Qin, Yongtao Liu,Yingqiu Zhou, Renbing Li and Xue-Ao Zhang “Exper-imental study of plasmon in a grating coupled graphenedevice with a resonant cavity” Appl. Phys. Lett. ,191905 (2015)[14] M. M. Jadidi, A. B. Sushkov, R. L. Myers-Ward, A. K.Boyd, K. M. Daniels, D.K. Gaskill, M. S. Fuhrer, H. Den-nis Drew and T. E. Murphy “Tunable Terahertz HybridMetal - Graphene Plasmons” Nano Lett.
5, 7099 -7104(2015). [15] Bo Zhao and Zhuomin M. Zhang “Strong Plasmonic Cou-pling between Graphene Ribbon Array and Metal Grat-ings” ACS Photonics , 1611 -1618 (2015).[16] Hua Lu, Jianlin Zhao and Min Gu “Nanowires-assistedexcitation and propagation of mid-infrared surface plas-mon polaritons in graphene” J. Appl. Phys. , 163106(2016).[17] S.M. Kukhtaruk, V.V. Korotyeyev, V.A. Kochelap andL. Varani, “Interaction of THz Radiation with PlasmonicGrating Structures Based on Graphene,” Proceedings ofInternational Conference on Mathematical Methods inElectromagnetic Theory, Lviv, pp. 196-199 (2016).[18] T. Low, A. Chaves, J. D. Caldwell, A. Kumar, N.X. Fang, Ph. Avouris, T. F. Heinz, F. Guinea, L.Martin-Moreno, F. Koppens, “Polaritons in layered two-dimensional materials,” Nature Materials , 182-194(2017).[19] S.A. Mikhailov, “Plasma instability and amplification ofelectromagnetic waves in low-dimensional electron sys-tems,” Phys. Rev. B , 1517 (1998).[20] V. V. Korotyeyev and V. A. Kochelap, “Plasma waveoscillations in a nonequilibrium two-dimensional electrongas: Electric field induced plasmon instability in the tera-hertz frequency range,” Phys. Rev. B , 235420 (2020).[21] D. V. Fateev, V. V. Popov, and M. S. Shur, “Plasmonspectra transformation in grating-gate transistor struc-ture with spatially modulated two-dimensional electronchannel”, Semiconductors , 1455 (2010) [Fiz. Tekh.Poluprovodn. (St. Petersburg) , 1455-1462 (2010)].[22] V. V. Popov, D. V. Fateev, T. Otsuji, Y. M. Meziani,D. Coquillat, and W. Knap, “Plasmonic terahertz detec-tion by a double-grating-gate field-effect transistor struc-ture with an asymmetric unit cell,” Appl. Phys. Lett. ,243504 (2011).[23] Y. M. Lyaschuk and V. V. Korotyeyev, “Theory of de-tection of terahertz radiation in hybrid plasmonic struc-tures with drifting electron gas,” Ukr. J. Phys. (10),889 (2017).[24] Olga V. Shapoval, Juan Sebastian Gomez-Diaz, JulienPerruisseau-Carrier, Juan R. Mosig and Alexander I.Nosich, IEEE Transactions on Terahertz Science andTechnology “Integral Equation Analysis of Plane WaveScattering by Coplanar Graphene - Strip Gratings in theTHz Range” (5), 666-674 (2013).[25] O. V. Shapovala and A. I. Nosich “Finite gratings ofmany thin silver nano strips: Optical resonances and roleof periodicity,” AIP Advances , 042120 (2013).[26] M. G. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc.Am (7), 811-818 (1981).[27] M. G. Moharam, E. B. Grann, D. A. Pommet and T.K. Gaylord, “Formulation for stable and efficient imple-mentation of the rigorous coupled-wave analysis of binarygratings,” J. Opt. Soc. Am. A (5) 1068-1076 (1995).[28] M. G. Moharam, Drew A. Pommet, Eric B. Grann andT. K. Gaylord, “Stable implementation of the rigor-ous coupled-wave analysis for surface-relief gratings: en-hanced transmittance matrix approach,” J. Opt. Soc.Am. A (5), 1077-1086 (1995).[29] S. Inampudi and H. Mosallaei “Tunable wideband-directive thermal emission from SiC surface using bun- dled graphene sheets”, Phys. Rev B , 125407 (2017).[30] S. Inampudi, V. Toutam and S. Tadigadapa “Robust visi-bility of graphene monolayer on patterned plasmonic sub-strates”, Nanotechnology , 015202 (2019).[31] L. Li, “Use of Fourier series in the analysis of discontinu-ous periodic structures,” J. Opt. Soc. Am. A (9), 1870-1876 (1996).[32] E. Popov, M. Nevi`ere, and N. Bonod, “Factorization ofproducts of discontinuous functions applied to Fourier-Bessel basis,” J. Opt. Soc. Am. A (1), 46-52 (2004).[33] E. Bleszynski, M. Bleszynski and T. Jaroszewicz“Surface-Integral Equations for Electromagnetic Scatter-ing from Impenetrable and Penetrable Sheets,” IEEE An-tennas and Propagation Magazine (6), 14-25 (1993).[34] T. L. Zinenko and A. I. Nosich “Plane Wave Scatteringand Absorption by Resistive-Strip and Dielectric-StripPeriodic Gratings,” IEEE Transactions on Antennas andPropagation (10), 1498-1505 (1998).[35] V. P. Gusynin, S. G. Sharapov, and J. P. Carbotte, “Sumrules for the optical and Hall conductivity in graphene,”Phys. Rev. B
5, 165407 (2007).[36] L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multi-layer,” Phys. Rev B , 153410 (2007).[37] M. V. Balaban, O. V. Shapoval, and A. I. Nosich “THzwave scattering by a graphene strip and a disk in thefree space: integral equation analysis and surface plasmonresonances,” J. Opt.
5, 114007 (2013). [38] S.M. Kukhtaruk, V.A. Kochelap, V.N. Sokolov, K.W.Kim, “Spatially dispersive dynamical response of hot car-riers in doped graphene,” Physica E: Low-dimensionalSystems and Nanostructures, , 26-37, (2016).[39] Y. M. Lyaschuk and V. V. Korotyeyev, “Interaction of aTerahertz electromagnetic wave with the plasmonic sys-tem ”grating-2D-gas”. Analysis of features of the nearfield,” Ukr. J. Phys. (5), 495-504 (2014).[40] A. V. Chaplik, Surf. Sci. Rep. “Absorption and emissionof electromagnetic waves by two-dimensional plasmons” , 289 (1985).[41] M. Dyakonov and M. Shur, “Shallow water analogy fora ballistic fild effct transistor: New mechanism of plasmawave generation by dc current” Phys. Rev. Lett., ,2465-2467 (1993).[42] Parametrical studies of the plasmon resonance at dif-ferent configurations of grating coupler is reported inRef.[39] and recent experimental paper Ref.[11].[43] K.-H. Brenner, “Aspects for calculating local absorptionwith the rigorous coupled-wave method,” Optics Express (10), 10369 (2010).[44] M. Weismann, Dominic F G Gallagher, and Nicolae CPanoiu, “Accurate near-field calculation in the rigor-ous coupled-wave analysis method,” Journal of Optics, (12), 125612 (2015).[45] COMSOL Multiphysics ®®