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.

## I. INTRODUCTION

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 slab^{25} 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 $\chi eff(3)$ in graphene is typically in the range of 10^{–}^{17}–10^{–}^{19} m^{2}/V^{2} 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 dielectrics^{44} 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.

## II. DISPERSION RELATION FOR GUIDED PLASMONS

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=\xb1d/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 as^{37}

where $\u03f511$ 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\u2212$ 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 as^{46}

In the above set of equations, $E\tau $ is the tangential electric field in graphene, *β* is the plasmon propagation constant in the $z\u2212$ direction, $q22=(\beta 2\u2212\u03f52k02)$ is the attenuation factor in the linear dielectric medium, and $k0=\omega /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 as

^{48}

where *e* is the elementary charge, *μ* is the chemical potential, *k _{BT}* is the thermal energy, $\u210f$ 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 as

^{24}

where $\Omega =\u210f\omega /\mu ,\u2009VF\u2248106$ 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\u2212$ direction. For the harmonic time-dependence of the incident radiation ($e\u2212i\omega t$), we can write the electric field of the propagating modes in the *j ^{th}* medium (

*j*= 1, 2) according to $Ej=(Exjx\u0302+iEzjz\u0302)ei\beta z$. Here, $Ex(x)$ and $Ez(x)$ must satisfy Maxwell's equations

^{25}

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

where $(q1\u2212)2=\beta 2\u2212\u03f51\u2212k02,\u2009\u03f51\u2212=\u03f511+\alpha |E1(d/2)|2$, and $|E1(d2)|2=|Ex1(d2)|2+|Ez1(d2)|2$. The $z\u2212$ component of the electric field in the linear dielectric for $x\u2265d/2$ can be written as

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 $\u03f51(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

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

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

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, $\alpha T\u22640.1$.

^{24,49}Hence,

*α*is ignored in this paper, which is in accordance with prior approaches.

_{T}^{33,46}For negligible

*α*, the real part of $\sigma 3(\Omega )$ 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

_{T}For finding the relationship between $Ez1(Ex1)$ and $\u03f51(Ex1)$, we solve the conservation law stated as^{25,35,36}

where *C* is a function of the electric field intensity (*E*_{0}) 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 $\u03f51(Ex1)$ according to

where $X(Ex1)=(2\alpha Ex12)2+C(1+2\chi \alpha Ex12)$ and $\chi =(k0/\beta )2$. From Eqs. (5) and (14), we can derive the plasmon dispersion relationship as follows:^{25}

The unknowns in Eq. (15) include $Ex1(0)$, $Ex1(d/2)$, and $Ez1(Ex1)$. $Ex1(0)$ is the $x\u2212$ 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 *E*_{0} 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 $\u03f51\u2212$, which can be found from the conservation law,^{25} which is stated as

Equations (15) and (16) build a nonlinear system with two unknowns, frequency and $\u03f51\u2212$, 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

### A. Plasmon modes in absence of nonlinearity

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

where symbols have been defined previously. Each of the above dispersion relation leads to two plasmon branches. Ignoring the nonlinear optical response of graphene ($\sigma 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}

### B. Nonlinear optical response of multilayer graphene

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 *E* – *k* 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(\sigma g+\sigma NL|E\tau |2)E\tau $, where *σ _{g}* and

*σ*are defined for monolayer graphene and

^{NL}*N*is the number of layers in the graphene stack.

## III. RESULTS AND DISCUSSION

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 SiO_{2} with a dielectric permittivity of 3.9. The core material is a Kerr-type dielectric with $\u03f511=2.5$ and $\alpha =10\u221214$ m^{2}/V^{2}. 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, $\sigma NL|E\tau |2$ must be comparable to *σ _{g}* in Eq. (2a). This implies that the field strength must exceed $\u223c106$ 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

*α*($\u223c10\u221219$ m

^{2}/V

^{2}) as in Ref. 37 will require electric field strengths exceeding the breakdown threshold of graphene. In this work, we choose $\alpha =10\u221214$ m

^{2}/V

^{2}, 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 works^{24,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 $\u03f5ave=(\u03f5SiO2+\u03f51\u2212)/2$. Since $\u03f51\u2212$ 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 × 10^{16} m^{–2} and defect density of 10^{14} m^{–}^{2}, we obtained the scattering rate of 1.86 × 10^{13} 1/s for the symmetric mode at $\beta d=2$ for *d* = 30 nm and the scattering rate of 3.9 × 10^{12} 1/s for the anti-symmetric mode at $\beta d=7$ for *d* = 150 nm. Details are presented in the Appendix.

### A. Mode propagation in the nonlinear waveguide

Figure 2 shows the dispersion of the guided plasmonic modes for two different thicknesses *d* = 30 nm and *d* = 150 nm for $\alpha E02=0.7$ ($E0=8.36\xd7106$ 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 ($\u2248$38 THz).^{45}

As Fig. 2 shows, the waveguide supports four branches of guided plasmon modes, which we denote as *Sa*_{1}, *Sa*_{2}, *Sb*_{1}, and *Sb*_{2} for the symmetric case and *Aa*_{1}, *Aa*_{2}, *Ab*_{1}, and *Ab*_{2} for the anti-symmetric case. These modes correspond to the roots of Eq. (17). The *Sa*_{2} mode does not exist at $\beta d>$ 0.2 and $f>$ 32.2 THz for *d* = 150 nm, and $\beta d>$ 0.04 and $f>$ 32.2 THz for *d* = 30 nm; therefore, *Sa*_{2} is not shown in Fig. 2. Additionally, modes marked by $[\alpha =0,\sigma NL\u22600]$ 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 $[\alpha =0,\sigma 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 ($\alpha \u22600$, $\sigma 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 ($\alpha =6.4\xd710\u221212$ m^{2}/V^{2}) and using low electric field intensity ($E0\u22643.3\xd7105$ 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 *b*_{1} and *b*_{2} display similar dispersion features in low frequencies, these modes have a physically distinct origin. Mode *b*_{2} corresponds to a negative value of the parameter Δ in Eq. (17). Further, the $z\u2212$ component of the electric field at the location of graphene-Kerr interface ($x=d/2$) has opposite sign for *b*_{1} and *b*_{2} modes. In fact, mode *b*_{2} occurs as an extra branch of mode *b*_{1} in the presence of graphene optical nonlinearity. To verify our hypothesis, we consider the results in Ref. 25, in which mode *a*_{2} is referred to as an extra branch of mode *a*_{1} resulting from the optical nonlinearity of the core dielectric. The value of Δ for the *a*_{2} mode is negative of that for the *a*_{1} mode. Additionally, the amplitude of the electric field at $x=d/2$ is the same for both *a*_{1} and *a*_{2} modes, but the value of *E _{z}* for the extra branch

*a*

_{2}is the negative of that for the

*a*

_{1}mode. Therefore, based on the terminology adopted in Ref. 25, we conclude that the mode

*b*

_{2}indeed appears as an extra branch of mode

*b*

_{1}.

There are a few interesting differences between the various plasmon branches for both symmetric and anti-symmetric modes that we discuss below. First, modes *a*_{2}, *b*_{1}, and *b*_{2} have positive group velocity for all frequencies considered in this work. However, the group velocity of mode *a*_{1} 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

*b*

_{1}and

*b*

_{2}is two orders of magnitude greater than that for modes

*a*

_{1}and

*a*

_{2}(see Fig. 3). These results imply that modes

*b*

_{1}and

*b*

_{2}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

*Ab*

_{2}(

*Sb*

_{2}) propagates for $\beta 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.

### B. Effect of dielectric and graphene nonlinearity

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 *a*_{1} acquires a negative group velocity is higher than the corresponding value of *βd* for which mode *a*_{1} 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 $\alpha =0,\sigma NL\u22600$ for $\beta d>$ 6 and *d* = 30 nm.

For a perfectly linear waveguide $[\sigma g=\sigma NL=0]$, there exists only one plasmon mode that exhibits a positive group velocity for all frequencies considered. This indicates that modes *b*_{1} and *b*_{2} appear due to the nonlinear optical response of graphene with a Kerr-type core dielectric.

### C. Guided modes in multilayer graphene

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, $\alpha 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 *Sa*_{1} can exist is $\u2248$ 10.6 THz. However, with *N* = 5, we find that mode *Sa*_{1} can be supported on the waveguide at frequencies above 35 THz. Contrary to mode *a*_{1}, the frequency at which modes *b*_{1} and *b*_{2} exist is reduced in the case of multilayer graphene with *N* = 5. Note that *Sa*_{2} is not plotted in Fig. 4(a).

### D. Propagation velocity and localization length

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 ($\beta (\omega )$), phase velocity is calculated as $vp=\omega /\beta $, while group velocity is given as $vg=\u2202\omega /\u2202\beta $. In the case of dispersion-less propagation, *v _{p}* =

*v*. Therefore, the deviation in

_{g}*v*and

_{p}*v*characterizes the distortion in a Gaussian signal propagating through the waveguide.

_{g}^{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 *E*_{0} = 8.36 × 10^{6} V/m are shown. Two values of core optical nonlinearity are chosen: *α* = 0 and *α* = 10^{–14} m^{2}/V^{2} to highlight the key differences in propagation velocity due to dielectric nonlinearity. The results in this figure show that the phase velocity of *Sa*_{1} exceeds that of *Sb*_{1} and *Sb*_{2} implying that for fast communication, the *Sa*_{1} mode must be utilized in the waveguide. However, the *Sa*_{1} mode tends to have more dispersion compared to *Sb*_{1} and *Sb*_{2}. Therefore, a Gaussian signal excited in mode *Sa*_{1} 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 *Sb*_{1} and *Sb*_{2} modes of the waveguide with the Kerr type nonlinear dielectric core.

The localization length of various plasmon modes in the parallel-plate waveguide is shown in Fig. 6 for *d* = 30 nm, *N* = 1, *E*_{0} = 8.36 × 10^{6} V/m for two values of the optical nonlinearity of the core dielectric. Since the localization length of mode *b*_{2} is similar to that of mode *b*_{1}, we choose to focus on the results corresponding to mode *b*_{1}. 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 *Sa*_{1} 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.

A summary of key results obtained in this paper is presented in Table I for four different structures: (i) *α* = 0 and $\sigma NL=0$, (ii) $\alpha \u22600$ and $\sigma NL=0$, (iii) *α* = 0 and $\sigma NL\u22600$, and (iv) $\alpha \u22600$ and $\sigma NL\u22600$.

α | σ _{NL} | Results |

A: anti-symmetric S: symmetric | ||

b: modes due to graphene optical nonlinearity | ||

subscript “2”: extra modes corresponding to the main modes | ||

× | × | Aa_{1}, Sa_{1} |

$vg>0$ (Ref. 58) | ||

✓ | × | Aa_{1}, Aa_{2} |

Sa_{1}, Sa_{2} | ||

$vg(a1)<0$ for high light intensity (Ref. 25) | ||

× | ✓ | Aa_{1}, Ab_{1} |

Sa_{1}, Sb_{1} | ||

$vg(a1)<0$ for high light intensity | ||

✓ | ✓ | Aa_{1}, Aa_{2}, Ab_{1}, Ab_{2} |

Sa_{1}, Sa_{2}, Sb_{1}, Sb_{2} most localized than α = 0, $\sigma NL\u22600vg(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 | ||

× | × | Aa_{1}, Sa_{1} |

$vg>0$ (Ref. 58) | ||

✓ | × | Aa_{1}, Aa_{2} |

Sa_{1}, Sa_{2} | ||

$vg(a1)<0$ for high light intensity (Ref. 25) | ||

× | ✓ | Aa_{1}, Ab_{1} |

Sa_{1}, Sb_{1} | ||

$vg(a1)<0$ for high light intensity | ||

✓ | ✓ | Aa_{1}, Aa_{2}, Ab_{1}, Ab_{2} |

Sa_{1}, Sa_{2}, Sb_{1}, Sb_{2} most localized than α = 0, $\sigma NL\u22600vg(a1)<0$ for high light intensity b is less dispersive than a |

## IV. CONCLUSION

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.

## ACKNOWLEDGMENTS

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

### APPENDIX: OPTICAL CONDUCTIVITY OF GRAPHENE

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(\beta )=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 *a*_{1} and *a*_{2}, the imaginary part of graphene conductivity is significantly larger than its real part, justifying our assumption. On the other hand, for modes *b*_{1} and *b*_{2}, 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 *b*_{1} and *b*_{2} follows a reasonable trend, which indicates only a small error in evaluating the mode characteristics.

Figure 8 shows $\u03f51\u2212$ versus *βd* for both symmetric and anti-symmetric modes. As expected, $\u03f51\u2212$ changes nearly by two orders of magnitude as $\beta 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 $\Gamma =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 $\alpha 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 *a*_{1} (solid line).