We theoretically model the dispersion characteristics of surface plasmons in a graphene-based parallel-plate waveguide geometry using nonlinear Kerr-type core (inter-plate) dielectric. The optical nonlinearity of graphene in the terahertz band under high light intensity is specifically included in the analysis. By solving Maxwell's equations and applying appropriate boundary conditions, we show that the waveguide supports four guided plasmon modes, each of which can be categorized as either symmetric or anti-symmetric based on the electric field distribution in the structure. Of the four guided modes, two modes are similar in characteristics to the modes obtained in the structure with linear graphene coating, while the two new modes have distinct characteristics as a result of the nonlinearity of graphene. We note that the group velocity of one of the plasmon modes acquires a negative value under high light intensity. Additionally, the optical nonlinearity of the core dielectric leads to a significant enhancement in the localization length of various plasmon modes. The description of the intra-band optical conductivity of graphene incorporates effects of carrier scatterings due to charged impurities, resonant scatterers, and acoustic phonons at 300 K. The proposed structure offers flexibility to tune the waveguide characteristics and the mode index by changing light intensity and electrochemical potential in graphene for reconfigurable plasmonic devices.

The two-dimensional (2D) material graphene has emerged as a tunable, dynamic platform to implement ultra-fast electronic and optical devices.1–6 Unlike the vast majority of 2D materials, graphene exhibits a strong light-matter interaction over a broad frequency spectrum ranging from the gigahertz to the ultraviolet.7–9 Additionally, the optoelectronic properties of graphene can be configured post-fabrication via electrostatic gating, which has been exploited in recent works to develop concepts for tunable optical devices.10,11 Graphene, along with other 2D materials, offers the flexibility to create new and cost-effective van der Waal heterostructures, thereby eliminating fabrication constraints resulting from lattice matching in three-dimensional materials.1,12 Various optoelectronic applications of 2D heterostructures, such as optical modulators, photodetectors, solar cells, and transparent electrodes, are currently being investigated.13–15 Graphene has a crucial role in advancing these applications due to its unique linear energy-dispersion relationship and superior carrier transport properties.16,17

In the terahertz (THz) band, graphene nanostructures have been shown to support highly localized surface plasmons that are polarized in the transverse magnetic (TM) direction.18,19 Plasmons are collective charge oscillations that exist at the metal-dielectric interface when excited by incident electromagnetic wave at the appropriate wavelength.20 Plasmons in 2D electron systems can be accessed via different methods including optical measurements, electron energy-loss spectroscopy,21 inelastic light scattering, angle-resolved photo-emission spectroscopy (ARPES),22 and scanning tunnelling spectroscopy.16,23 A graphene sheet essentially behaves as a thin metal sheet with a positive imaginary part of optical conductivity in the THz regime, which enables the generation of TM plasmons.24,25 The two main advantages of TM plasmon modes in graphene are the considerable enhancement of field localization and low propagation loss.24 Due to field localization, the wavelength of plasmons in graphene is (40–100)× smaller than the free-space wavelength in the THz regime.26,27 Therefore, graphene plasmons confine light below its diffraction limit, which can facilitate the heterogeneous integration of electronic and optical devices.28 Graphene can also support TE plasmons at high frequencies when the imaginary component of the optical conductivity of graphene becomes negative.29,30 Unlike TM plasmons, TE plasmons are weakly bound to the graphene interface, but they have remarkably low propagation loss.24 Surface plasmons are also observed in graphene-based double-layer or parallel-plate waveguide geometries embedded in dielectric media.25,31 In this case, the plasmons split into symmetric (even) and anti-symmetric (odd) modes.25 The odd plasmon modes are quasi-transverse electromagnetic in nature with the electric field localized strongly between the parallel plates. The odd plasmon modes propagate in the waveguide nearly dispersion-less and, therefore, support extremely high bandwidths in interconnect applications.32 

Despite the large body of work in the area of graphene plasmonics, very few research groups have studied the impact of optical nonlinearity of waveguide materials on plasmon propagation.33 References 34, 35, and 36 investigate plasmon characteristics in a metallic layer surrounded by Kerr-type nonlinear dielectric media. More recently, surface plasmons in the graphene sheet on a Kerr-type substrate,37 graphene parallel-plate waveguide bounded by Kerr media,19 and graphene-coated Kerr slab25 have been studied. While prior works consider the optical nonlinearity of the dielectric material surrounding graphene, the effect of optical nonlinearity of graphene on plasmon modes has been largely ignored. That is, prior works assume that there exists a linear relationship between the induced surface current in graphene and the amplitude of incident electromagnetic radiation. This assumption is strictly true under low light intensities and could be inapplicable in the case of graphene, which is known to produce a strong third-order nonlinear optical response.38 In recent works, the optical nonlinearity of graphene has been extensively investigated.29,39–43 Extracted effective bulk susceptibility χeff(3) in graphene is typically in the range of 1017–1019 m2/V2 depending on the measurement technique, the frequency of incident radiation, and the fabrication method.42 The nonlinearity of graphene is multiple orders of magnitude larger than that of bulk dielectrics44 and is also comparable to the nonlinearity of other strongly nonlinear materials including carbon nanotubes.24,40 In light of recent results, we recognize that ignoring the nonlinear contribution of graphene is unjustified because of the strong localization of the electromagnetic field for various plasmon modes.

In this work, we address the limitations of prior works by modeling plasmon propagation in a parallel-plate waveguide comprising a nonlinear graphene sheet within an inhomogeneous dielectric media (linear substrate and Kerr-type core dielectric). The device setup used in this work is illustrated in Fig. 1. We use the method described in Refs. 24, 25, and 36 to find the plasmon dispersion relationship in the device setup and extend the approach to incorporate the nonlinear optical response of graphene. Moreover, we include the effect of carrier scatterings due to the substrate and cladding to describe the optical conductivity of graphene.45 We elucidate the impact of material and geometrical parameters on plasmon characteristics using the models developed in this work. We show that by considering the nonlinear optical response of graphene, the parallel-plate waveguide supports up to four plasmon branches corresponding to the symmetric and anti-symmetric field profiles. This behavior is different from that reported in Ref. 25 in which only two plasmon modes are found. It is shown that there are two specific modes for which the third order susceptibility of graphene is significantly higher, which allows low dispersion in the mode propagation. The advantage of the proposed waveguide is that the plasmon characteristics, such as the localization length, propagation velocity, and mode index, can be configured easily by changing the dimensions, field strength, material type, and the number of layers in the multilayer graphene stack.

FIG. 1.

The structure under study consists of graphene-based parallel-plate waveguide with the Kerr-type inter-plate dielectric medium. The waveguide is embedded in a linear dielectric environment. The structure is symmetric about x = 0 plane. The plasmons propagate in the z direction.

FIG. 1.

The structure under study consists of graphene-based parallel-plate waveguide with the Kerr-type inter-plate dielectric medium. The waveguide is embedded in a linear dielectric environment. The structure is symmetric about x = 0 plane. The plasmons propagate in the z direction.

Close modal

In the device setup as shown in Fig. 1, we consider a parallel-plate waveguide geometry with nonlinear Kerr-type dielectric media between the two plates separated by a thickness of d. The parallel plates at the boundaries x=±d/2 are implemented using graphene. The parallel-plate geometry is enclosed in a linear dielectric environment. The linear dielectrics are considered semi-infinite on both sides. The relative permittivity of the Kerr-type core dielectric under optical nonlinearity is given as37 

ϵ1=ϵ11+α|E|2,
(1)

where ϵ11 and ϵ1 are the linear and the total dielectric permittivities, respectively, α is the Kerr coefficient of the dielectric, and E is the electric field in the dielectric. This structure supports TM plasmons propagating along the z direction. To include the optical nonlinearity of graphene in our analysis, we consider that the surface conducitivity of graphene is nonlinear. Therefore, the magnitude of the induced surface current in graphene is expressed as46 

J=σEτ=(σg+σNL|Eτ|2)Eτ,
(2a)
σNL=σ34(3+β2q22).
(2b)

In the above set of equations, Eτ is the tangential electric field in graphene, β is the plasmon propagation constant in the z direction, q22=(β2ϵ2k02) is the attenuation factor in the linear dielectric medium, and k0=ω/c is the free-space wave-vector. The linear part of graphene conductivity, σg, is obtained using the Drude model that is valid for intra-band scattering mechanisms.47 It is mathematically expressed as48 

σg=i2e2π2kBTω+iτ1ln[2cosh(|μ|2kBT)],
(3)

where e is the elementary charge, μ is the chemical potential, kBT is the thermal energy, is the reduced Planck's constant, and τ is the carrier relaxation time in graphene. In Eq. (2b), σ3 is the third order nonlinear conductivity of graphene, which is given as24 

σ3(Ω)=i332e2π(eVF)22μ4Ω3(1+iαT),
(4)

where Ω=ω/μ,VF106 m/s is the Fermi velocity of carriers in graphene and αT is the two-photon absorption in graphene.24 

We assume that the structure is illuminated by a TM-polarized wave for which the magnetic field is oriented in the y direction. For the harmonic time-dependence of the incident radiation (eiωt), we can write the electric field of the propagating modes in the jth medium (j = 1, 2) according to Ej=(Exjx̂+iEzjẑ)eiβz. Here, Ex(x) and Ez(x) must satisfy Maxwell's equations25 

(5)
dEzjdx=qj2βExj,
(5a)
dExjdx+Exjϵjdϵjdx=βEzj.
(5b)

Due to the symmetry of the parallel-plate structure, we focus on one interface at x=d/2 for purposes of theoretical calculations. From Eq. (5), we have

dEz1dx|x=d2=(q1)2βEx1(d2),
(6)

where (q1)2=β2ϵ1k02,ϵ1=ϵ11+α|E1(d/2)|2, and |E1(d2)|2=|Ex1(d2)|2+|Ez1(d2)|2. The z component of the electric field in the linear dielectric for xd/2 can be written as

Ez2(x)=Ez2(d2)eq2(xd2).
(7)

Note that in this analysis, the electromagnetic field distribution in the nonlinear medium is unknown. However, to obtain the plasmon dispersion relationship in the structure, we must know the variation of ϵ1 with Ex1 throughout the Kerr medium. Hence, we use two key concepts for finding the dispersion relationship: (i) application of appropriate boundary conditions at the interface x=d/2 and (ii) solving the conservation law to derive the closed-form expression of ϵ1(Ex1).25,35,36 In the following paragraphs, we discuss both concepts briefly.

By applying the standard boundary condition at x=d/2 for TM polarization, we obtain

(8)
Ez1(d2)=Ez2(d2),
(8a)
(ϵ2q22dEz2dxϵ1(q1)2dEz1dx)|x=d2=iωϵ0[σg+σNL|Ez1(d2)|2]Ez1(d2).
(8b)

Equations (6)–(8) lead to the following relationship between the x and z components of the electric field at x=d/2 in the nonlinear dielectric:

Ez1(d2)=ϵ1ΔβEx1(d2),
(9)
Δ=ϵ2q2iωϵ0[σg+σNL|Ez1(d2)|2].
(10)

Equation (10) holds for the general case of complex dynamical conductivity of graphene. Under the condition that Ohmic losses in graphene are negligible, we can assume that the conductivity of graphene is purely imaginary; therefore, Δ can be simplified to

Δ=ϵ2q2+1ωϵ0[σg+σNL|Ez1(d2)|2].
(11)

In the above derivation, we also assume that αT [defined in Eq. (4)] is equal to zero. Recent experiments show that the value of the two-photon absorption in graphene, αT0.1.24,49 Hence, αT is ignored in this paper, which is in accordance with prior approaches.33,46 For negligible αT, the real part of σ3(Ω) as defined in Eq. (4) is much smaller than its imaginary part. The validity of a purely imaginary conductivity of graphene is examined in the  Appendix. By invoking the relationship E12(d2)=Ex12(d2)+Ez12(d2) and using Eqs. (1) and (9), we can derive an expression for Ex1(d2), which is given as

αEx12(d2)=ϵ1ϵ111+(ϵ1)2/(βΔ)2.
(12)

For finding the relationship between Ez1(Ex1) and ϵ1(Ex1), we solve the conservation law stated as25,35,36

ϵ1[ϵ12(2χϵ1)αEx12]=C,
(13)

where C is a function of the electric field intensity (E0) at the center and can be calculated using Eqs. (5) and (6) in Ref. 36 for symmetric and anti-symmetric modes, respectively. Equation (13) can be solved to find the root ϵ1(Ex1) according to

ϵ1(Ex1)=2αEx12+X(Ex1)1+2χαEx12,
(14)

where X(Ex1)=(2αEx12)2+C(1+2χαEx12) and χ=(k0/β)2. From Eqs. (5) and (14), we can derive the plasmon dispersion relationship as follows:25 

βd2=Ex1(0)Ex1(d/2)ϵ1(Ex1)+2αEx12Ez1(Ex1)X(Ex1)dEx1.
(15)

The unknowns in Eq. (15) include Ex1(0), Ex1(d/2), and Ez1(Ex1). Ex1(0) is the x component of the electric field at the center of the structure and is strictly zero for the anti-symmetric mode, and it is equal to E0 for the symmetric mode, and Ez1(Ex1) can be found from Eq. (4) of Ref. 36. For calculating Ex1(d/2) from Eq. (12), we need to calculate ϵ1, which can be found from the conservation law,25 which is stated as

ϵ1Cϵ1=2(2χϵ1)(ϵ1ϵ11)1+(ϵ1)2/(βΔ)2.
(16)

Equations (15) and (16) build a nonlinear system with two unknowns, frequency and ϵ1, which can be solved by using any standard numerical technique. In this work, we use the particle swarm optimization method (PSO) owing to its simple structure, ability to handle vectorized inputs, and high performance.50 

The main difference between our work and Ref. 25 is that we consider the nonlinear optical response of graphene, which manifests itself in the term Δ defined in Eq. (17). As a result of the nonlinear optical response of graphene, Eq. (10) becomes a cubic equation in terms of Δ with three roots, each of which represents a unique plasmon mode. Mathematically, Δ is re-written as

Δ=ϵ2q2iωϵ0[σg+σNL×ϵ1ϵ11α|ϵ1|2|ϵ1|2+|Δβ|2].
(17)

To obtain the plasmon dispersion relationship in the case of a nonlinear graphene parallel-plate waveguide with a linear core medium (α0), we obtain symmetric and anti-symmetric modes as discussed in Ref. 36. For the symmetric mode,

(18)
Fs=σg+σNLE02q12β2sinh2(q1d2),
(18a)
coth(q1d2)=q1ϵ1ϵ2q2(1+iq2Fsωϵ0ϵ2).
(18b)
For the anti-symmetric mode,

Fas=σg+σNLE02cosh2(q1d2),
(19a)
tanh(q1d2)=q1ϵ1ϵ2q2(1+iq2Fasωϵ0ϵ2),
(19b)

where symbols have been defined previously. Each of the above dispersion relation leads to two plasmon branches. Ignoring the nonlinear optical response of graphene (σNL=0), Eqs. (18b) and (19b) reduce to the well-known dispersion relationship of TM surface plasmons of the typical symmetric graphene parallel-plate waveguide in the absence of any optical nonlinearities.51 

The optical conductivity of graphene scales with the number of layers in the graphene multilayer stack, which is considered in this work. We assume that the Ek relationship of multilayer graphene is linear. This is true in the case of multilayer graphene structures fabricated using chemical vapor deposition or epitaxial growth on silicon carbide. With such fabrication techniques, the layers in the graphene stack are decoupled and rotated relative to each other destroying the Bernal stacking order that is present in exfoliated graphene stacks.52,53 In the absence of electrostatic screening between the layers, the linear part of the optical conductivity of multilayer graphene becomes proportional to the number of layers.45 As discussed in Refs. 40, 44, 54, and 55, the optical nonlinearity of graphene also increases proportionally to the number of layers. To incorporate the effect of number of layers in our analysis, we modify Eq. (2a) as J=N(σg+σNL|Eτ|2)Eτ, where σg and σNL are defined for monolayer graphene and N is the number of layers in the graphene stack.

In this section, we study the impact of optical nonlinearity of the parallel-plate waveguide on the plasmon characteristics using the analytical model derived in Sec. II. For all analyses, we consider that the dielectric material surrounding the parallel-plate waveguide is SiO2 with a dielectric permittivity of 3.9. The core material is a Kerr-type dielectric with ϵ11=2.5 and α=1014 m2/V2. We take a lower value of α compared to Ref. 25 to have a stronger electric field and considerable nonlinear response of graphene at the THz frequency range. Note that to have a significant impact on plasmon mode characteristics due to graphene optical nonlinearity in THz regime, σNL|Eτ|2 must be comparable to σg in Eq. (2a). This implies that the field strength must exceed 106 V/m. By considering Kerr-type core dielectric with a considerably large value of α as in Ref. 25, it is found that nonlinear modes are excited only at sub-THz frequencies. Since this work focuses on frequencies in excess of 1 THz, we choose a Kerr-type dielectric core with a lower optical nonlinearity. On the other hand, a very low value of α (1019 m2/V2) as in Ref. 37 will require electric field strengths exceeding the breakdown threshold of graphene. In this work, we choose α=1014 m2/V2, which allows us to probe graphene nonlinearity at the desired frequency range and illustrate some critical mode characteristics in the presence of graphene nonlinearity that would otherwise not be possible.

A significant difference in this work compared to prior works24,25,46 is that we evaluate the carrier scattering rate as a function of the chemical potential, temperature, and the substrate and cladding materials of the waveguide. We model the optical conductivity of graphene using Eq. (3) in which the carrier scattering time due to charged impurities, resonant scatterers, and acoustic phonons is incorporated using the physics-based models developed in Ref. 45. We assume that the chemical potential in graphene is set to 0.3 eV, and the temperature is 300 K. For calculating carrier scatterings due to charged impurities, we model the average relative permittivity of the medium surrounding graphene according to ϵave=(ϵSiO2+ϵ1)/2. Since ϵ1 changes with the frequency of operation and the mode index, the approach adopted by prior works in assuming a constant carrier scattering rate is unjustified. We correct for this limitation by invoking appropriate models of carrier scattering rate in this work. For example, by using the charged impurities concentration of 2.2 × 1016 m–2 and defect density of 1014 m2, we obtained the scattering rate of 1.86 × 1013 1/s for the symmetric mode at βd=2 for d = 30 nm and the scattering rate of 3.9 × 1012 1/s for the anti-symmetric mode at βd=7 for d = 150 nm. Details are presented in the  Appendix.

Figure 2 shows the dispersion of the guided plasmonic modes for two different thicknesses d = 30 nm and d = 150 nm for αE02=0.7 (E0=8.36×106 V/m) for symmetric and anti-symmetric modes in the parallel-plate waveguide geometry. To validate the use of Drude optical conductivity model, we restrict our analysis to frequencies below the excitation frequency of optical phonons in graphene (38 THz).45 

FIG. 2.

Dispersion characteristics of symmetric (S*) and anti-symmetric (A*) plasmon modes for the structure under study: (a) symmetric mode with d = 30 nm, (b) anti-symmetric mode with d = 30 nm, (c) symmetric mode with d = 150 nm, and (d) anti-symmetric mode with d = 150 nm. d denotes the core dielectric thickness. Since Sa2 propagates at an extremely low value of βd, it is not plotted in the figure. α refers to the Kerr coefficient of the dielectric, while σNL refers to the optical nonlinearity of graphene. Results corresponding to α0,σNL=0 have the same behaviour as those reported in Ref. 25.

FIG. 2.

Dispersion characteristics of symmetric (S*) and anti-symmetric (A*) plasmon modes for the structure under study: (a) symmetric mode with d = 30 nm, (b) anti-symmetric mode with d = 30 nm, (c) symmetric mode with d = 150 nm, and (d) anti-symmetric mode with d = 150 nm. d denotes the core dielectric thickness. Since Sa2 propagates at an extremely low value of βd, it is not plotted in the figure. α refers to the Kerr coefficient of the dielectric, while σNL refers to the optical nonlinearity of graphene. Results corresponding to α0,σNL=0 have the same behaviour as those reported in Ref. 25.

Close modal

As Fig. 2 shows, the waveguide supports four branches of guided plasmon modes, which we denote as Sa1, Sa2, Sb1, and Sb2 for the symmetric case and Aa1, Aa2, Ab1, and Ab2 for the anti-symmetric case. These modes correspond to the roots of Eq. (17). The Sa2 mode does not exist at βd> 0.2 and f> 32.2 THz for d = 150 nm, and βd> 0.04 and f> 32.2 THz for d = 30 nm; therefore, Sa2 is not shown in Fig. 2. Additionally, modes marked by [α=0,σNL0] demonstrate the dispersion curves for a nonlinear graphene parallel-plate waveguide with a linear core medium obtained using Eqs. (18b) and (19b). Modes labeled by [α=0,σNL=0] represent the dispersion result of a typical graphene parallel-plate waveguide without considering any optical nonlinearity. To compare the results with prior work, we also evaluate plasmon modes by ignoring the nonlinear optical response of graphene (α0, σNL=0 in Fig. 2). The results in this case have the same behaviour as those reported in Ref. 25. In Ref. 25, surface plasmons in graphene-coated Kerr slab have been studied. However, by considering the large Kerr coefficient for the core dielectric (α=6.4×1012 m2/V2) and using low electric field intensity (E03.3×105 V/m), the effect of optical nonlinearity of graphene on plasmon modes is neglected. Additionally, the carrier scattering rate in graphene is considered as a constant number in Ref. 25. Hence, the effects of chemical potential, temperature, and waveguide materials on plasmons propagation are not studied.

While the nonlinear modes marked by b1 and b2 display similar dispersion features in low frequencies, these modes have a physically distinct origin. Mode b2 corresponds to a negative value of the parameter Δ in Eq. (17). Further, the z component of the electric field at the location of graphene-Kerr interface (x=d/2) has opposite sign for b1 and b2 modes. In fact, mode b2 occurs as an extra branch of mode b1 in the presence of graphene optical nonlinearity. To verify our hypothesis, we consider the results in Ref. 25, in which mode a2 is referred to as an extra branch of mode a1 resulting from the optical nonlinearity of the core dielectric. The value of Δ for the a2 mode is negative of that for the a1 mode. Additionally, the amplitude of the electric field at x=d/2 is the same for both a1 and a2 modes, but the value of Ez for the extra branch a2 is the negative of that for the a1 mode. Therefore, based on the terminology adopted in Ref. 25, we conclude that the mode b2 indeed appears as an extra branch of mode b1.

There are a few interesting differences between the various plasmon branches for both symmetric and anti-symmetric modes that we discuss below. First, modes a2, b1, and b2 have positive group velocity for all frequencies considered in this work. However, the group velocity of mode a1 changes the sign at higher values of βd. The negative group velocity of this mode is associated with a much stronger electric field and plasmon confinement at the graphene sheet. Further, simulations show that σNL for modes b1 and b2 is two orders of magnitude greater than that for modes a1 and a2 (see Fig. 3). These results imply that modes b1 and b2 are influenced more by optical nonlinearity than the other modes. We also note that there exists an upper cut-off value of βd for the propagation of various guided modes. In the case of d = 150 nm, only mode Ab2 (Sb2) propagates for βd> 4.5. Similar cut-off behavior is observed for a lower value of d = 30 nm. This is due to the fact that the parameter Δ only has one relevant root beyond a certain value of βd for both symmetric and anti-symmetric modes.

FIG. 3.

Nonlinear conductivity of graphene (σNL) for d = 30 nm and αE02=0.7 for Sa1 and Sb1 modes.

FIG. 3.

Nonlinear conductivity of graphene (σNL) for d = 30 nm and αE02=0.7 for Sa1 and Sb1 modes.

Close modal

To study the effect of the core dielectric nonlinearity on the plasmon dispersion characteristics, we examine the behavior of various modes for a nonlinear graphene parallel-plate geometry with linear core dielectric. For this structure, each of the symmetric and anti-symmetric plasmon modes exhibit two distinct branches. One of the two plasmon branches can propagate with either positive or negative group velocity depending on the mode index. This is similar to the mode propagation observed in the case with the nonlinear Kerr dielectric core. However, the value of βd for which the plasmon mode a1 acquires a negative group velocity is higher than the corresponding value of βd for which mode a1 in the nonlinear structure acquires a negative group velocity. This implies that negative group velocity originates due to the optical nonlinearity of the core dielectric, but it is efficiently modulated by the nonlinearity of the graphene sheet. Since the structure with a linear core displays only two main modes without the existence of extra branches, we confirm that extra mode branches result from the dielectric optical nonlinearity rather than from the graphene optical nonlinearity.

Also, there exists an upper cut-off value of βd for the propagation of guided modes in the nonlinear graphene parallel-plate geometry with linear core dielectric. Our simulation results reported in Figs. 2(a) and 2(b) demonstrate that there are no guided plasmon modes for α=0,σNL0 for βd> 6 and d = 30 nm.

For a perfectly linear waveguide [σg=σNL=0], there exists only one plasmon mode that exhibits a positive group velocity for all frequencies considered. This indicates that modes b1 and b2 appear due to the nonlinear optical response of graphene with a Kerr-type core dielectric.

To study the effect of number of layers in multilayer graphene on mode propagation, we obtain the dispersion relationship of the parallel-plate waveguide for d = 150 nm, αE02 = 0.7, and number of layers N = 5. Results are shown in Fig. 4. A comparison of the results plotted in Fig. 4 with the results in Figs. 2(c) and 2(d) (N = 1) shows that increasing the number of layers in the graphene stack expands the range of frequencies at which the waveguide can support guided plasmon modes. With N = 1, the maximum frequency at which mode Sa1 can exist is 10.6 THz. However, with N = 5, we find that mode Sa1 can be supported on the waveguide at frequencies above 35 THz. Contrary to mode a1, the frequency at which modes b1 and b2 exist is reduced in the case of multilayer graphene with N = 5. Note that Sa2 is not plotted in Fig. 4(a).

FIG. 4.

Dispersion of the symmetric and anti-symmetric modes in multilayer graphene with five layers and d = 150 nm: (a) symmetric mode and (b) antisymmetric mode. For each plot, solid line, dashed-dotted line, dashed line, and dotted line represent a1, a2, b1, and b2, respectively.

FIG. 4.

Dispersion of the symmetric and anti-symmetric modes in multilayer graphene with five layers and d = 150 nm: (a) symmetric mode and (b) antisymmetric mode. For each plot, solid line, dashed-dotted line, dashed line, and dotted line represent a1, a2, b1, and b2, respectively.

Close modal

It is important to characterize the properties of plasmon modes in terms of phase velocity, group velocity, and the localization length. For waveguiding applications using plasmon-based transmission, it is important to have a lower localization length for better confinement, while maximizing the speed at which information is communicated. From the dispersion relationship (β(ω)), phase velocity is calculated as vp=ω/β, while group velocity is given as vg=ω/β. In the case of dispersion-less propagation, vp = vg. Therefore, the deviation in vp and vg characterizes the distortion in a Gaussian signal propagating through the waveguide.56,57 The localization length (LL) of plasmon modes is mathematically given as LL=1/Re(q), where q is the plasmon wave-vector in the surrounding dielectric media.

In Fig. 5, group and phase velocities for symmetric modes with d = 30 nm, N = 1, and E0 = 8.36 × 106 V/m are shown. Two values of core optical nonlinearity are chosen: α  =  0 and α  =  10–14 m2/V2 to highlight the key differences in propagation velocity due to dielectric nonlinearity. The results in this figure show that the phase velocity of Sa1 exceeds that of Sb1 and Sb2 implying that for fast communication, the Sa1 mode must be utilized in the waveguide. However, the Sa1 mode tends to have more dispersion compared to Sb1 and Sb2. Therefore, a Gaussian signal excited in mode Sa1 will undergo more distortion on the waveguide. Hence, one can tradeoff distortion and speed for optimal signaling using the proposed nonlinear waveguide. Moreover, using linear dielectric as a core leads to higher phase velocity and fast communication, while its lower branch in Fig. 2 has significantly more dispersion compare to Sb1 and Sb2 modes of the waveguide with the Kerr type nonlinear dielectric core.

FIG. 5.

Phase and group velocities of the symmetric modes for one layer of graphene for d = 30 nm and E0 = 8.36 × 106 V/m: (a) Sa1 and upper branch of the structure with α  =  0, σNL0 in Fig. 2, (b) Sb1, Sb2, and lower branch of the structure with α  =  0, σNL0 in Fig. 2.

FIG. 5.

Phase and group velocities of the symmetric modes for one layer of graphene for d = 30 nm and E0 = 8.36 × 106 V/m: (a) Sa1 and upper branch of the structure with α  =  0, σNL0 in Fig. 2, (b) Sb1, Sb2, and lower branch of the structure with α  =  0, σNL0 in Fig. 2.

Close modal

The localization length of various plasmon modes in the parallel-plate waveguide is shown in Fig. 6 for d = 30 nm, N = 1, E0 = 8.36 × 106 V/m for two values of the optical nonlinearity of the core dielectric. Since the localization length of mode b2 is similar to that of mode b1, we choose to focus on the results corresponding to mode b1. We clearly see from the figure that the nonlinearity of the core dielectric enhances localization at all values of the mode index for both symmetric and anti-symmetric modes. As much as 11% of reduction in localization length of Sa1 mode is observed at 10 THz for finite α. Therefore, for high-frequency signaling purposes, it is favorable to use the Kerr-type core material to enhance field confinement between the waveguide plates and reduce the undesirable effects of cross-talk.

FIG. 6.

Localization length of (a) a1 and (b) b1 plasmon modes. Here, d = 30 nm, E0 = 8.36 × 106 V/m, and only one layer of grapnene is considered. Solid line and dotted line correspond to symmetric and anti-symmetric modes with α=1014 m2/V2, respectively. Dashed-dotted and dashed line correspond to symmetric and anti-symmetric modes with α  =  0 m2/V2.

FIG. 6.

Localization length of (a) a1 and (b) b1 plasmon modes. Here, d = 30 nm, E0 = 8.36 × 106 V/m, and only one layer of grapnene is considered. Solid line and dotted line correspond to symmetric and anti-symmetric modes with α=1014 m2/V2, respectively. Dashed-dotted and dashed line correspond to symmetric and anti-symmetric modes with α  =  0 m2/V2.

Close modal

A summary of key results obtained in this paper is presented in Table I for four different structures: (i) α  =  0 and σNL=0, (ii) α0 and σNL=0, (iii) α  =  0 and σNL0, and (iv) α0 and σNL0.

TABLE I.

The summary of characteristics of surface plasmons modes for different structures. vg denotes the group velocity of the plasmon mode.

α σNL Results 
A: anti-symmetric S: symmetric 
b: modes due to graphene optical nonlinearity 
subscript “2”: extra modes corresponding to the main modes 
× × Aa1, Sa1 
  vg>0 (Ref. 58
✓ × Aa1, Aa2 
Sa1, Sa2 
vg(a1)<0 for high light intensity (Ref. 25
× ✓ Aa1, Ab1 
Sa1, Sb1 
vg(a1)<0 for high light intensity 
✓ ✓ Aa1, Aa2, Ab1, Ab2 
Sa1, Sa2, Sb1, Sb2 most localized than α = 0, σNL0vg(a1)<0 for high light intensity b is less dispersive than a 
α σNL Results 
A: anti-symmetric S: symmetric 
b: modes due to graphene optical nonlinearity 
subscript “2”: extra modes corresponding to the main modes 
× × Aa1, Sa1 
  vg>0 (Ref. 58
✓ × Aa1, Aa2 
Sa1, Sa2 
vg(a1)<0 for high light intensity (Ref. 25
× ✓ Aa1, Ab1 
Sa1, Sb1 
vg(a1)<0 for high light intensity 
✓ ✓ Aa1, Aa2, Ab1, Ab2 
Sa1, Sa2, Sb1, Sb2 most localized than α = 0, σNL0vg(a1)<0 for high light intensity b is less dispersive than a 

In this paper, we perform theoretical calculations to study the characteristics of guided plasmon modes in a graphene-based parallel-plate waveguide geometry in the presence of optical nonlinearity. Our calculations show that due to the optical nonlinearity of the structure, four plasmon modes are created, each of which has either symmetric or anti-symmetric electric field distribution. Two of the four plasmon modes are interpreted as main modes, while the remainder are the extra branches of the main modes. We identified that the extra plasmon branches are a direct consequence of the optical nonlinearity of the Kerr-type core dielectric. Additionally, the impact of optical nonlinearity of graphene is strongly dependent on the mode index. As such, two plasmon modes in the structure that exist at lower THz frequencies display a higher third-order susceptibility and are also less dispersive with lower phase velocity. The group velocity of one of the less nonlinear plasmon modes acquires a negative value for high incident light intensity. The properties of various plasmon modes are also examined as functions of the geometrical and material parameters of the waveguide and the effect of various optical nonlinearities in the structure are delineated. We show that using nonlinear Kerr-type dielectric as a core leads to lower localization length and more confinement, which is a useful feature for waveguiding applications. In this work, the nonlinear optical conductivity of graphene is modeled by incorporating carrier scatterings due to charged impurities, resonant scatterers, and acoustic phonons. Therefore, the scattering rate and the optical conductivity are dependent on the mode index and frequency of operation. While we have assumed that the conductivity of graphene is purely imaginary, which validates the use of conservation laws, we are currently investigating numerical techniques to handle a complex-valued graphene conductivity.

The authors acknowledge the funding support of National Science Foundation through the grant no. CCF-1565656.

We discuss the validity of the assumption that the real part of optical conductivity of graphene is significantly lower than its imaginary part. This allows us to neglect absorption losses in graphene, implying zero energy dissipation in the system (Im(β)=0). Optical conductivity of graphene for d =30 nm for all symmetric and anti-symmetric modes is plotted in Fig. 7. The figure shows that for modes a1 and a2, the imaginary part of graphene conductivity is significantly larger than its real part, justifying our assumption. On the other hand, for modes b1 and b2, the assumption of a purely imaginary conductivity of graphene fails at lower values of βd. Yet, for small βd values, the dispersion plot for modes b1 and b2 follows a reasonable trend, which indicates only a small error in evaluating the mode characteristics.

FIG. 7.

Conductivity of graphene versus βd for d =30 nm: (a) symmetric mode Sa1, (b) symmetric mode Sb1 (solid line) and Sb2 (dashed line), (c) anti-symmetric mode Aa1 (Solid line) and Aa2 (dashed line), (d) anti-symmetric mode Sb1 (solid line), and Sb2 (dashed line).

FIG. 7.

Conductivity of graphene versus βd for d =30 nm: (a) symmetric mode Sa1, (b) symmetric mode Sb1 (solid line) and Sb2 (dashed line), (c) anti-symmetric mode Aa1 (Solid line) and Aa2 (dashed line), (d) anti-symmetric mode Sb1 (solid line), and Sb2 (dashed line).

Close modal

Figure 8 shows ϵ1 versus βd for both symmetric and anti-symmetric modes. As expected, ϵ1 changes nearly by two orders of magnitude as βd increases. Therefore, it would be unjustified to assume a frequency-independent carrier scattering rate while evaluating the dynamical conductivity of graphene. Based on the trend in Fig. 8, it is to be expected that the carrier scattering rate will strongly change with mode index of the plasmon. For calculating the linear part of the conductivity of graphene, method described in Ref. 45 must be used to properly account for the impact of Fermi level and temperature on various scattering channels in graphene. Considering a constant scattering rate Γ=1013 1/s in graphene, we obtain the dispersion relationship as shown in Fig. 9. Results are plotted for the symmetric mode in monolayer graphene with a chemical potential of 0.3 eV, d = 30 nm, and αE02=0.7. Curves indicated with * correspond to the results for constant Γ value. As expected, dispersion characteristics for constant Γ case are significantly different than the actual dispersion relationship for mode a1 (solid line).

FIG. 8.

ϵ1 versus βd for both symmetric and anti-symmetric modes for d = 30 and αE02=0.7. Solid line, dotted line, and dashed line correspond to a1, b1, and b2 nonlinear modes, respectively.

FIG. 8.

ϵ1 versus βd for both symmetric and anti-symmetric modes for d = 30 and αE02=0.7. Solid line, dotted line, and dashed line correspond to a1, b1, and b2 nonlinear modes, respectively.

Close modal
FIG. 9.

Dispersion of the symmetric mode for one layer of graphene with d = 30 nm, αE02=0.7, and μ = 0.3 eV. Solid line, dashed line, and dotted line represent a1, b1, and b2 nonlinear modes, respectively. Scattering rate for * curves is equal to 1013 1/s, while curves without * indicate actual dispersion relationship.

FIG. 9.

Dispersion of the symmetric mode for one layer of graphene with d = 30 nm, αE02=0.7, and μ = 0.3 eV. Solid line, dashed line, and dotted line represent a1, b1, and b2 nonlinear modes, respectively. Scattering rate for * curves is equal to 1013 1/s, while curves without * indicate actual dispersion relationship.

Close modal
1.
F.
Xia
,
H.
Wang
,
D.
Xiao
,
M.
Dubey
, and
A.
Ramasubramaniam
,
Nat. Photonics
8
,
899
(
2014
).
2.
S.
Rakheja
and
P.
Sengupta
,
J. Phys. D: Appl. Phys.
49
,
115106
(
2016
).
3.
S.
Rakheja
and
P.
Sengupta
,
IEEE Trans. Nanotechnol.
15
,
113
(
2016
).
4.
5.
F.
Bonaccorso
,
Z.
Sun
,
T.
Hasan
, and
A.
Ferrari
,
Nat. Photonics
4
,
611
(
2010
).
6.
T.
Mueller
,
F.
Xia
, and
P.
Avouris
,
Nat. Photonics
4
,
297
(
2010
).
7.
F.
Xia
,
H.
Yan
, and
P.
Avouris
,
Proc. IEEE
101
,
1717
(
2013
).
8.
H.
Choi
,
F.
Borondics
,
D. A.
Siegel
,
S. Y.
Zhou
,
M. C.
Martin
,
A.
Lanzara
, and
R. A.
Kaindl
,
Appl. Phys. Lett.
94
,
172102
(
2009
).
9.
Q.
Bao
and
K. P.
Loh
,
ACS Nano
6
,
3677
(
2012
).
10.
H.
Yang
,
J.
Heo
,
S.
Park
,
H. J.
Song
,
D. H.
Seo
,
K.-E.
Byun
,
P.
Kim
,
I.
Yoo
,
H.-J.
Chung
, and
K.
Kim
,
Science
336
,
1140
(
2012
).
11.
L.
Ju
,
B.
Geng
,
J.
Horng
,
C.
Girit
,
M.
Martin
,
Z.
Hao
,
H. A.
Bechtel
,
X.
Liang
,
A.
Zettl
,
Y. R.
Shen
 et al,
Nat. Nanotechnol.
6
,
630
(
2011
).
12.
A. K.
Geim
and
I. V.
Grigorieva
,
Nature
499
,
419
(
2013
).
13.
L.
Britnell
,
R.
Ribeiro
,
A.
Eckmann
,
R.
Jalil
,
B.
Belle
,
A.
Mishchenko
,
Y.-J.
Kim
,
R.
Gorbachev
,
T.
Georgiou
,
S.
Morozov
 et al,
Science
340
,
1311
(
2013
).
14.
W. J.
Yu
,
Y.
Liu
,
H.
Zhou
,
A.
Yin
,
Z.
Li
,
Y.
Huang
, and
X.
Duan
,
Nat. Nanotechnol.
8
,
952
(
2013
).
15.
K.
Roy
,
M.
Padmanabhan
,
S.
Goswami
,
T. P.
Sai
,
G.
Ramalingam
,
S.
Raghavan
, and
A.
Ghosh
,
Nat. Nanotechnol.
8
,
826
(
2013
).
16.
A.
Grigorenko
,
M.
Polini
, and
K.
Novoselov
,
Nat. Photonics
6
,
749
(
2012
).
17.
A. C.
Neto
,
F.
Guinea
,
N. M.
Peres
,
K. S.
Novoselov
, and
A. K.
Geim
,
Rev. Mod. Phys.
81
,
109
(
2009
).
18.
F. H.
Koppens
,
D. E.
Chang
, and
F. J.
García de Abajo
,
Nano Lett.
11
,
3370
(
2011
).
19.
H.
Hajian
,
A.
Soltani-Vala
,
M.
Kalafi
, and
P.
Leung
,
J. Appl. Phys.
115
,
083104
(
2014
).
20.
S. A.
Maier
,
Plasmonics: Fundamentals and Applications
(
Springer Science & Business Media
,
2007
).
21.
T.
Eberlein
,
U.
Bangert
,
R.
Nair
,
R.
Jones
,
M.
Gass
,
A.
Bleloch
,
K.
Novoselov
,
A.
Geim
, and
P.
Briddon
,
Phys. Rev. B
77
,
233406
(
2008
).
22.
A.
Bostwick
,
T.
Ohta
,
T.
Seyller
,
K.
Horn
, and
E.
Rotenberg
,
Nat. Physics
3
,
36
(
2007
).
23.
V. W.
Brar
,
S.
Wickenburg
,
M.
Panlasigui
,
C.-H.
Park
,
T. O.
Wehling
,
Y.
Zhang
,
R.
Decker
,
Ç.
Girit
,
A. V.
Balatsky
,
S. G.
Louie
 et al,
Phys. Rev. Lett.
104
,
036805
(
2010
).
24.
A.
Gorbach
,
Phys. Rev. A
87
,
013830
(
2013
).
25.
H.
Hajian
,
I. D.
Rukhlenko
,
P.
Leung
,
H.
Caglayan
, and
E.
Ozbay
,
Plasmonics
11
,
735
(
2016
).
26.
N. K.
Emani
,
D.
Wang
,
T.-F.
Chung
,
L. J.
Prokopeva
,
A. V.
Kildishev
,
V. M.
Shalaev
,
Y. P.
Chen
, and
A.
Boltasseva
,
Laser Photonics Rev.
9
,
650
(
2015
).
27.
V. W.
Brar
,
M. S.
Jang
,
M.
Sherrott
,
J. J.
Lopez
, and
H. A.
Atwater
,
Nano Lett.
13
,
2541
(
2013
).
28.
H.
Buljan
,
M.
Jablan
, and
M.
Soljačić
,
Nat. Photonics
7
,
346
(
2013
).
29.
S.
Mikhailov
and
K.
Ziegler
,
Phys. Rev. Lett.
99
,
016803
(
2007
).
30.
G. W.
Hanson
,
J. Appl. Phys.
103
,
064302
(
2008
).
31.
H.
Hajian
,
A.
Soltani-Vala
, and
M.
Kalafi
,
J. Appl. Phys.
114
,
033102
(
2013
).
32.
X.
Gao
,
L.
Zhou
, and
T. J.
Cui
,
Sci. Rep.
5
,
9250
(
2015
).
33.
X.
Jiang
,
J.
Bao
,
B.
Zhang
, and
X.
Sun
,
Nanoscale Res. Lett.
12
,
395
(
2017
).
34.
J.-H.
Huang
,
R.
Chang
,
P.-T.
Leung
, and
D. P.
Tsai
,
Opt. Commun.
282
,
1412
(
2009
).
35.
I. D.
Rukhlenko
,
A.
Pannipitiya
,
M.
Premaratne
, and
G. P.
Agrawal
,
Phys. Rev. B
84
,
113409
(
2011
).
36.
I. D.
Rukhlenko
,
A.
Pannipitiya
, and
M.
Premaratne
,
Opt. Lett.
36
,
3374
(
2011
).
37.
L.
Wang
,
W.
Cai
,
X.
Zhang
, and
J.
Xu
,
Opt. Lett.
37
,
2730
(
2012
).
38.
S.
Mikhailov
,
Europhys. Lett.
79
,
27002
(
2007
).
39.
S.
Mikhailov
and
K.
Ziegler
,
J. Phys.: Condens. Matter
20
,
384204
(
2008
).
40.
E.
Hendry
,
P. J.
Hale
,
J.
Moger
,
A.
Savchenko
, and
S.
Mikhailov
,
Phys. Rev. Lett.
105
,
097401
(
2010
).
41.
X.
Yao
and
A.
Belyanin
,
Phys. Rev. Lett.
108
,
255503
(
2012
).
42.
J.
Cheng
,
N.
Vermeulen
, and
J.
Sipe
,
New J. Phys.
16
,
053014
(
2014
).
43.
J. L.
Cheng
,
N.
Vermeulen
, and
J.
Sipe
,
Phys. Rev. B
91
,
235320
(
2015
).
44.
G.
Demetriou
,
H. T.
Bookey
,
F.
Biancalana
,
E.
Abraham
,
Y.
Wang
,
W.
Ji
, and
A. K.
Kar
,
Opt. Express
24
,
13033
(
2016
).
45.
S.
Rakheja
,
V.
Kumar
, and
A.
Naeemi
,
Proc. IEEE
101
,
1740
(
2013
).
46.
D. A.
Smirnova
,
A. V.
Gorbach
,
I. V.
Iorsh
,
I. V.
Shadrivov
, and
Y. S.
Kivshar
,
Phys. Rev. B
88
,
045443
(
2013
).
47.
L.
Falkovsky
and
A.
Varlamov
,
Eur. Phys. J. B
56
,
281
(
2007
).
48.
L.
Falkovsky
,
J. Phys.: Conf. Ser.
129
,
012004
(
2008
).
49.
T.
Gu
,
N.
Petrone
,
J. F.
McMillan
,
A.
van der Zande
,
M.
Yu
,
G.-Q.
Lo
,
D.-L.
Kwong
,
J.
Hone
, and
C. W.
Wong
,
Nat. Photonics
6
,
554
(
2012
).
50.
S.
Ebbesen
,
P.
Kiwitz
, and
L.
Guzzella
, in
American Control Conference (ACC), 2012
(IEEE,
2012
), pp.
1519
1524
.
51.
A.
Auditore
,
C.
De Angelis
,
A.
Locatelli
,
S.
Boscolo
,
M.
Midrio
,
M.
Romagnoli
,
A.-D.
Capobianco
, and
G.
Nalesso
,
Opt. Lett.
38
,
631
(
2013
).
52.
F.
Varchon
,
P.
Mallet
,
L.
Magaud
, and
J.-Y.
Veuillen
,
Phys. Rev. B
77
,
165415
(
2008
).
53.
A.
Reina
,
X.
Jia
,
J.
Ho
,
D.
Nezich
,
H.
Son
,
V.
Bulovic
,
M. S.
Dresselhaus
, and
J.
Kong
,
Nano Lett.
9
,
30
(
2009
).
54.
H.
Zhang
,
S.
Virally
,
Q.
Bao
,
L. K.
Ping
,
S.
Massar
,
N.
Godbout
, and
P.
Kockaert
,
Opt. Lett.
37
,
1856
(
2012
).
55.
J.
Khurgin
,
Appl. Phys. Lett.
104
,
161116
(
2014
).
56.
S.
Rakheja
,
IEEE Trans. Nanotechnol.
15
,
936
(
2016
).
57.
S.
Rakheja
, in
2017 18th International Symposium on Quality Electronic Design (ISQED)
(IEEE,
2017
), pp.
45
51
.
58.
Y.
Kurokawa
and
H. T.
Miyazaki
,
Phys. Rev. B
75
,
035411
(
2007
).