Coupled TM surface plasmon features of graphene-metal layered structure in the sub-THz frequency range

TM surface plasmon (SP) characteristics of a four-layer structure, consisting of air as the superstrate, a monolayer of graphene, a dielectric buffer layer and metal as the substrate are analyzed at sub-THz frequencies. TM SPs in such case are represented by metal-like and graphene-like branches. For small frequencies the metal-like plasmon splits up into two branches depending on the graphene electron concentration; one of the branches goes into cutoff at the point where the branch features Brewster-type characteristics. Graphene-like plasmon modes are converted into short-range modes for small buffer thicknesses. Brewster-type SP modes can be effectively modulated in the vicinity of their cutoff thicknesses by means of the graphene electron concentration.


Introduction
Electronic, optical and mechanical properties of graphene were extensively studied after the appearance of the first publication devoted to experimental investigations of graphene [1]; results of these studies were widely published [2][3][4]. Free carriers in graphene can be created either by doping or by field effects. A graphene electron gas of high mobility forms an ideal two dimensional sheet [1]. The optical surface conductivity of a graphene monolayer was studied in a wide range of frequencies, from the terahertz to the visible range [2,[5][6][7][8]. Electromagnetic wave (EMW) propagation in structures containing graphene layers was investigated in a number of publications [9][10][11][12] and various perspective applications for graphene in photonics and optoelectronics have been discussed [13]. Special attention was paid to plasmons supported by graphene layers and which are tunable by a gate voltage. Perspective applications of graphene plasmons in nanophotonics due to their high confinement around the graphene layer have been discussed in many papers [14][15][16][17][18][19][20][21][22][23].
In spite of the fact that graphene plasmons were investigated in a number of structures containing a few graphene monolayers or graphene bilayers, the case of graphene plasmons coupled with other types of structure modes (surface metal-plasmon or waveguide modes) has not been thoroughly examined. The coupling of graphene plasmons with surface metal plasmons in a structure containing a graphene layer and metal substrate separated by an air gap was studied [24]: the solution of the coupled plasmon mode shows a linear dispersion behavior in a specific parameter range.
In this article we analyze the propagation characteristics of TM coupled modes in a structure containing a metal substrate, dielectric buffer layer, a monolayer of graphene, and air as the superstrate.
The article is divided into two parts. In the first part, the dispersion equation for TM eigenmodes and some of its approximate solutions are given on the basis of Maxwell's equations solution. In the second part, dispersion characteristics of coupled graphene-metal plasmons are analyzed on the basis of numerical solutions of the dispersion equation for sub-THz frequencies. The results are summarized in the conclusions.

Dispersion equation for the coupled graphenemetal surface plasmons (MSPs) in four-layered structure
We consider the propagation of TM-polarized EMWs along the Z direction in the structure shown in figure 1; magnetic field of TM-polarized EMW is directed along X axis and electric field vector of TM-polarized EMW is parallel to (Y, Z) plane, the magnetic and electric field component respectively can be written as: H = (H x , 0, 0), and E = (0, E y , E z ). The structure of figure 1 consists of a metal substrate with relative dielectric permittivity ε , 3 a buffer layer of thickness d 2 and relative dielectric permittivity ε , 2 a graphene monolayer created on the buffer layer surface y = 0 (Y-axis is orthogonal to the multi-layer structure) and air as the lossless superstrate medium behind the graphene layer with relative dielectric permittivity ε = 1 1 .
The standard electrodynamic boundary conditions represented by the continuity of the electric field tangential component at the interfaces y = 0, y = −d 2 and by the continuity of the magnetic field tangential component at the interface y = −d 2 are used to calculate the EMW propagation in the structure of figure 1. The boundary condition for the tangential component of the magnetic field at y = 0 is modified by the presence of the monolayer of graphene, which features a surface electron concentration n s and a surface conductivity σ. In the SI units, the boundary condition can be written down as follows: where H x, 1 and H x,2 are the tangential components of the magnetic field in the superstrate and the buffer layer, respectively; = E y ( 0) z is the amplitude of the tangential component of the electric field at y = 0.
Straightforward calculations for the TM eigenmodes lead to the spatial distribution of the magnetic field amplitude given by equations (2a) and (2c) and to the dispersion equation (3): y > 0: where q is the propagation constant in z direction, k 0 is the electromagnetic wavenumber in vacuum, ω is the EMW angular frequency, ε μ = c 1/ 0 0 where c is the EMW speed in vacuum, ε 0 and μ 0 are permittivity and permeability of vacuum.
Throughout the paper we ignore the spatial dispersion of graphene conductivity, suggesting that ω ≪ qv ,  [25], however temporal dispersion of graphene conductivity is taken into account therefore σ ω σ ω = q ( , ) ( ). We use equation (4a) for graphene's conductivity σ ω ( ) consisting of intraband σ ω ( ) intra and interband σ ω ( ) inter contributions [2]:  Dependence of Fermi level energy E F (from now on this energy is denoted as Fermi level) versus graphene surface electron concentration n s and temperature t is defined by the following integral equation: Interband transitions become more important for higher optical frequencies and lower Fermi levels. At room temperature, and low terahertz frequencies f < 30 THz, intraband transitions give the largest contribution to the optical conductivity of the graphene monolayer for E F > 0.05 eV (EMW frequency is denoted as f, in contrast to electron distribution function which is denoted as f(ε)).
Equation (3) was derived in [26] although the coupling of MSPs and surface carrier plasmons created by a two dimensional surface charge en s embedded in the interface y = 0 had not been studied.
For a very wide buffer layer thickness when → ∞ d , 2 equation (3) reduces to the dispersion equations of both decoupled TM MSPs (equation (6)) and TM graphene plasmons (equation (7), see also [7,12,26]): Rewriting equation (6) leads to the well-known dispersion relation for MSPs [27] given by equation (8): Equation (7) is described by the 4th power equation for the q 2 value and its general solutions are too cumbersome to write them down. However, the solution is considerably simplified for the case ε ≫ q k 0 1,2 when the dispersion relation for graphene plasmons can be written down as follows [7]: Taking into account only the intraband contribution to graphene's optical conductivity given by equation (4b) for τ → ∞, equation (3) is reduced to equation (10) for is the metal electron plasma frequency): ⎡ ⎣ ⎤ ⎦ is the graphene plasmon frequency defined by equation (9). Equation (10) is coincident with equation (22) derived in [24]; the equation describes frequencies of coupled metal-graphene plasmons for ε ε = = 1 1 2 when the plasmon absorption is totally neglected. However, this equation is valid only for large propagation constant values when the conditions ε ≫ q k j 0 are fulfilled and these inequalities are not fulfilled for long-range modes of low terahertz frequencies propagating in the structure of figure 1 with material parameters under study (refer to the graphs plotted in figure 2). In this case it is necessary to solve equation (3) to find out the dispersion of the coupled modes. The terms long/ short-range modes are used to define the amount of mode damping: for long-range or short-range modes the propagation length is respectively larger or smaller than the mode's wavelength.

Results of numerical calculations and discussion
Numerical calculations were performed for the structure of figure 1 with the following material parameters for the different layers: the relative dielectric permittivity of air ε = 1, 1 for the buffer layer we consider SiO 2 with ε = + i 4 0.04 , 2 valid for EMW frequency ⩽ f 3 THz) [28], and finally for the metal we used Au. Frequency dependent dielectric relative permittivity ε 3 of Au is described by Drude model with plasma and damping frequencies taken from [29]. The optical conductivity of graphene depends on the electron relaxation time τ, whereby the τ value is evaluated from the DC mobility μ according to the relation τ μ π = ℏ n ev .  values of the real and imaginary parts of the propagation constant with small (short range) or large (long range) propagation lengths. For further classification of the SP mode branches we first assume an infinite thick buffer, such that one can define purely decoupled but Fermi level dependent graphene surface plasmons (GSP) and MSP solutions. The MSP solution exist over the full considered Fermi level range (0-0.5 eV), however the GSP do not.
The dependences of the real parts of the propagation constant of the decoupled MSP and GSP modes on the plasmon frequencies for E F = 0.5 eV are depicted in figure 2 versus.
The pure MSP results from equation (8), which even approximates to ε = q k 0 2 for low THz frequencies. The pure GSP exists in the three-layer system air/graphene/buffer: equation (7) yields the exact solution and equation (9) for ε ≫ q k 0 1,2 yields an approximate solution. The approximate solution of equation (9) leads to considerable errors over the whole considered frequency band. It is also seen that the difference of decoupled MSP and GSP propagation constant vanishes for low THz and GHz frequencies. Hence we may expect stronger coupling effects in this frequency range. All subsequent numerical calculations are mainly performed in the strong coupling regime for f = 0.3 THz.
The dispersion dependence of decoupled GSP on graphene's electron Fermi energy level is shown in figure 3 for the EMW frequency f = 0.3 THz. From figure 3 one can observe that the decoupled GSP is represented by a longrange mode only for ⩾ E 0.32 eV, F when π < Im q Re q ( ) ( ) 2 and the GSP propagates at least one wavelength.
It is obvious that for a thinner buffer layer decoupled plasmon dispersion dependences will modify. In the extreme case of a very thin buffer layer when → d 0, 2 the only longrange solution of equation (3) exists and this is represented by the MSP for the two-layer superstrate/metal structure. This solution is very close to that given by equation (8)  although slightly modified due to the presence of the monolayer graphene. The long-range SP tends to become the surface metal plasmon of the three-layer structure superstrate/buffer/metal in the case of small electron concentrations in graphene. Increasing the graphene electron concentration leads to a considerable modification of the coupled SP spectrum.
We classify coupled SPs as metal-like surface plasmons (MLSP) and graphene-like surface (GLSP) plasmons depending on their behavior for → ∞ d 2 . At this juncture the dependences of coupled plasmon modes on the buffer thickness and the Fermi level E F in the graphene will be elaborated.
Dependences of MLSP propagation constant values versus buffer thickness for different Fermi levels are shown in For larger Fermi levels, the MLSP spectrum is considerably modified: as illustrated in figure 5.  The MLSP splits up into two mode branches. One of them exists within the whole range of buffer thicknesses and this solution is called the continuous MLSP branch, the other MLSP solution exists only for buffer thicknesses smaller than the cutoff thickness d 2c and this is called the cutoff MLSP mode. For → ∞ d , 2 the propagation constant solution converges to q , 3,2 and it converts from long-range solution to a short-range solution for buffer layer thinner than d . 1 also the imaginary part of the propagation constant vanishes at the cutoff thickness d 2c . Furthermore, the cutoff MLSP branch converts from the bound type mode for d 2 < d 2c to the growing type mode for d 2 > d 2c and this branch is represented by an incident plane wave without any reflection for d 2 = d 2c . These kinds of modes are called Brewster-type modes [30]. The whole incident wave is actually absorbed inside the structure [31,32]. The propagation constant values real parts q c at the cutoff conditions for TM modes can be found from the solution of equation (11): where = ± ± … N 0, 1, 2, . Buffer layer thicknesses at the cutoff conditions are found using the relation: The splitting of the MLSP is due to its interaction with an additional cutoff mode solution of dispersion equation (3), which behaves as a short-range mode for small buffer thickness d 2 . As illustration, the real and imaginary parts of the propagation constant values of the following modes have been plotted in figure 6 versus buffer thickness for f = 0.3 THz: for E F = 0.178 eV one  The continuous MLSP branch for E F = 0.179 eV turns out to be close to the additional cutoff mode at E F = 0.178 eV for buffer thickness d 2 < 126.5 μm and to the MLSP at E F = 0.178 eV for buffer thickness d 2 > 126.5 μm. On the other hand the cutoff MLSP branch for E F = 0.179 eV turns out to be close to the additional cutoff mode at E F = 0.178 eV for buffer thickness d 2 > 126.5 μm and to the MLSP at E F = 0.178 eV for buffer thickness d 2 < 126.5 μm.
It is seen that there is no crossing points for the MLSP branches dependences when E F = 0.179 eV unlike for MLSP mode and additional cutoff mode dependences at E F = 0.178 eV when these crossing points exist both for real and imaginary parts of the propagation constant. This is clearly seen in figure 6(a) and (b) graphs. Dependences of the GLSP propagation constant values on buffer thickness are shown in figure 7 for different Fermi levels E F . One can notice that for a thinner buffer layer, these modes become more damped and finally they are converted into short-range ones for small buffer thickness due to the stronger mode localization around the graphene layer. This effect is more pronounced for smaller Fermi level energies. GLSP for E F = 0.2 eV and E F = 0. 25  It is interesting to mention that GLSP dispersion dependences for E F = 0.35 eV and E F = 0.5 eV are markedly different (the same is valid for continuous MLSP branches shown in figure 5). This is explained by mode modification spectrum analogously to that resulting in the MLSP splitting. The point is that dispersion dependences for GLSP mode and continuous MLSP branch become very close for Fermi level about E F = 0.48 eV and buffer layer thickness about 750 μm (see figure 8). This leads to a modified mode spectrum. As a result, the dependences for the real parts of the mode propagation constant values are repulsed in this region of buffer thickness for E F = 0.5 eV and GLSP mode (continuous MLSP branch) propagation constant dependence for E F = 0.5 eV turns out to be close to continuous MLSP branch (GLSP mode) propagation constant dependence for E F = 0.48 eV in the region of buffer thickness μ < d 700 m 2 . If we replace GLSP dispersion dependence for E F = 0.5 eV in figure 7 by continuous MLSP branch dispersion dependences for E F = 0.5 eV taken from figure 5 we see that continuous MLSP branch dispersion dependences fits much better, the set of graphs plotted in figure 7 (the same is valid for the analogous replacement of graphs in figure 5). The modes are defined as graphene (metal)-like ones according to their behavior for very large buffer layer when d 2 tends to infinity, therefore GLSP branches can be replaced by the MLSP ones due to the mode modification. At higher frequencies, one can show that the difference between propagation constant values of decoupled metal and graphene plasmons (see figure 2) becomes larger as the graphene surface conductivity σ ω ( ) decreases due to the decreased  Both real and imaginary parts of the MLSP propagation constant can be effectively modulated by changing the Fermi level in graphene; modulation of the imaginary part of the propagation constant, in particular, can lead to a large increase of the MLSP damping and to the SP propagation blockage.
These observations are in agreement with what was found in figures 4 and 5 where the MLSP propagation constant was considerably modulated by changing the graphene electron concentration. This effect is especially strong near the regions when SP splits up into two branches.

Conclusions
The dispersion dependences of TM graphene-metal coupled SPs were analyzed versus the buffer layer thickness and Fermi level in graphene. Calculations were mainly performed for sub-THz frequencies.  It is shown that TM SPs are represented by metal-like and graphene-like branches depending on their behavior for very thick buffer layers. For frequencies smaller than 0.75 THz (for the concrete structure under study), the MLSP split up into two branches depending on the graphene electron concentration: one of the branches (continuous MLSP branch) exists in the whole range of the buffer thickness and being a short-range mode for small thicknesses, another one (cutoff MLSP branch) undergoes cutoff and exists only within a limited range of buffer thicknesses smaller than the cutoff thickness. It is found that cutoff MLSP branches are represented by Brewster-type modes at the cutoff thicknesses and cutoff MLSP branches are converted into unphysical growing modes for larger buffer thicknesses. The graphene-like plasmons are converted into short-range modes for small buffer thicknesses due to strong mode confinement around graphene.
It is demonstrated that propagation constant values of SP mode near cutoff can be effectively modulated by changing the graphene electron concentration in the vicinity of the modes cutoff thicknesses.