Voltage Tunable Plasmon Propagation in Dual Gated Bilayer Graphene

TL;DR

This study theoretically explores how energy asymmetry in dual-gated bilayer graphene influences terahertz plasmon propagation, revealing tunable propagation length and confinement, with implications for reconfigurable plasmonic devices.

Contribution

It provides a theoretical framework for controlling plasmon propagation in bilayer graphene via dual gating and energy asymmetry, highlighting differences between AA and AB stacking.

Findings

Energy asymmetry increases TM plasmon propagation length.

AA stacking shows larger propagation length increase than AB.

Dual gating enables independent modulation of plasmon properties.

Abstract

In this paper, we theoretically investigate plasmon propagation characteristics in AB and AA stacked bilayer graphene (BLG) in the presence of energy asymmetry due to an electrostatic field oriented perpendicularly to the plane of the graphene sheet. We first derive the optical conductivity of BLG using the Kubo formalism incorporating energy asymmetry and finite electron scattering. All results are obtained for room temperature (300K) operation. By solving Maxwell's equations in a dual gate device setup, we obtain the wavevector of propagating plasmon modes in the transverse electric (TE) and transverse magnetic (TM) directions at terahertz frequencies. The plasmon wavevector allows us to compare the compression factor, propagation length, and the mode confinement of TE and TM plasmon modes in bilayer and monolayer graphene sheets and also study the impact of material parameters on plasmon characteristics. Our results show that the energy asymmetry can be harnessed to increase the propagation length of TM plasmons in BLG. AA stacked BLG shows a larger increase in propagation length than AB stacked BLG; conversely, it is very insensitive to the Fermi level variations. Additionally, the dual gate structure allows independent modulation of the energy asymmetry and the Fermi level in BLG, which is advantageous for reconfiguring plasmon characteristics post device fabrication.

Equations

$$\begin{array}{rl}& \{\begin{array}{l}\mathbf{E}=(E_{y_{+}}\hat{y})e^{-j(k_{z}z+k_{x}x)}\\ \mathbf{H}=(H_{x_{+}}\hat{x}+H_{z_{+}}\hat{z})e^{-j(k_{z}z+k_{x}x)}\end{array}x>0,\\ & \{\begin{array}{l}\mathbf{E}=(E_{y_{-}}\hat{y})e^{-j(k_{z}z-k_{x}x)}\\ \mathbf{H}=(H_{x_{-}}\hat{x}+H_{z_{-}}\hat{z})e^{-j(k_{z}z-k_{x}x)}\end{array}x<0,\end{array}$$

$$\begin{array}{rl}& \{\begin{array}{l}\mathbf{H}=(H_{y_{+}}\hat{y})e^{-j(k_{z}z+k_{x}x)},\\ \mathbf{E}=(E_{x_{+}}\hat{x}+E_{z_{+}}\hat{z})e^{-j(k_{z}z+k_{x}x)}\end{array}x>0,\\ & \{\begin{array}{l}\mathbf{H}=(H_{y_{-}}\hat{y})e^{-j(k_{z}z-k_{x}x)},\\ \mathbf{E}=(E_{x_{-}}\hat{x}+E_{z_{-}}\hat{z})e^{-j(k_{z}z-k_{x}x)}\end{array}x<0,\end{array}$$

$$1+\frac{j{\mu }_0\omega \sigma (\omega )}{2\sqrt{k_z^2-{\epsilon }_r{\omega }^2/c^2}}=0,\text{TE Mode},$$

$$1-\frac{j\sigma (\omega )\sqrt{k_z^2-{\epsilon }_r{\omega }^2/c^2}}{2{\epsilon }_r{\epsilon }_0\omega }=0,\text{TM Mode},$$

$$k_{z}=k_{\text{light}}\bigg{(}1+\frac{\eta^{2}}{4}(\sigma^{\prime\prime}(\omega)^{2}-\sigma^{\prime}(\omega)^{2})+j\frac{\eta^{2}}{2}\sigma^{\prime}(\omega)\sigma^{\prime\prime}(\omega)\bigg{)}^{1/2}$$

$$k_{z}=k_{\text{light}}\bigg{(}1+\frac{4}{\eta^{2}}\frac{\sigma^{\prime\prime}(\omega)^{2}-\sigma^{\prime}(\omega)^{2}}{|\sigma(\omega)|^{4}}-j\frac{8}{\eta^{2}}\frac{\sigma^{\prime}(\omega)\sigma^{\prime\prime}(\omega)}{|\sigma(\omega)|^{4}}\bigg{)}^{1/2}$$

$$\hat{H}=\left[\begin{array}{cccc}-\Delta /2& 0& 0& {\varphi }^{\ast }(\mathbf{k})\\ 0& \Delta /2& \varphi (\mathbf{k})& 0\\ 0& {\varphi }^{\ast }(\mathbf{k})& \Delta /2& {\gamma }_1\\ \varphi (\mathbf{k})& 0& {\gamma }_1& -\Delta /2\end{array}\right]\cdot$$

$$\begin{array}{rl}& {ϵ}_{\alpha }(\mathbf{k})=\\ & \pm \sqrt{\frac{{\gamma }_1^2}{2}+\frac{{\Delta }^2}{4}+|\varphi (\mathbf{k}){|}^2+(-1{)}^{\alpha }\sqrt{\frac{{\gamma }_1^4}{4}+|\varphi (\mathbf{k}){|}^2({\gamma }_1^2+{\Delta }^2)}}\cdot \end{array}$$

$$\hat{H}=\left[\begin{array}{cccc}-\Delta /2& 0& {\gamma }_1& \varphi (\mathbf{k})\\ 0& \Delta /2& {\varphi }^{\ast }(\mathbf{k})& {\gamma }_1\\ {\gamma }_1& \varphi (\mathbf{k})& \Delta /2& 0\\ {\varphi }^{\ast }(\mathbf{k})& {\gamma }_1& 0& -\Delta /2\end{array}\right]\cdot$$

$$\epsilon_{\alpha}({\bm{k}})=\pm\bigg{(}\sqrt{\gamma_{1}^{2}+\frac{\Delta^{2}}{4}}+(-1)^{\alpha}|\phi({\bm{k}})|\bigg{)}\cdot$$

$$\begin{split}\sigma^{\prime}_{\text{AB}}(\omega)=&\frac{2e^{2}}{\omega}\int_{-\infty}^{\infty}\frac{dE}{2\pi}[f(E)-f(E+\hbar\omega)]\\ &\int\frac{d^{2}{\bm{k}}}{(2\pi)^{2}}|v_{\bm{k}}|^{2}\bigg{[}A_{11}(E,\Delta)A_{44}(E+\hbar\omega,\Delta)\\ &+A_{44}(E,\Delta)A_{11}(E+\hbar\omega,\Delta)\\ &+A_{11}(E,-\Delta)A_{44}(E+\hbar\omega,-\Delta)\\ &+A_{44}(E,-\Delta)A_{11}(E+\hbar\omega,-\Delta)\\ &+2(A_{13}(E,\Delta)A_{13}^{*}(E+\hbar\omega,-\Delta)\\ &+A_{13}^{*}(E,-\Delta)A_{13}(E+\hbar\omega,\Delta))\bigg{]}\cdot\end{split}$$

$$G_{ij}(z)={\int }_{-\infty }^{\infty }\frac{dw}{2\pi }\frac{A_{ij}(w)}{z-w},$$

$$A_{ij}(w,\Delta )=\sum _{\alpha =1,2}[a_{ij}(\alpha ,\Delta )\delta (w-{ϵ}_{\alpha })+a_{ij}(\alpha ,-\Delta )\delta (w+{ϵ}_{\alpha })],$$

$$\begin{split}&a_{11}(\alpha,\Delta)=\frac{(-1)^{\alpha+1}\pi}{\epsilon_{2}^{2}-\epsilon_{1}^{2}}\bigg{(}[(\frac{\Delta}{2})^{2}+\gamma_{1}^{2}-\epsilon_{\alpha}^{2}](1-\frac{\Delta}{2\epsilon_{\alpha}})\\ &+|\phi({\bm{k}})|^{2}(1+\frac{\Delta}{2\epsilon_{\alpha}})\bigg{)},\\ &a_{44}(\alpha,\Delta)=\frac{(-1)^{\alpha+1}\pi}{\epsilon_{2}^{2}-\epsilon_{1}^{2}}\bigg{(}[(\frac{\Delta}{2})^{2}-\epsilon_{\alpha}^{2}](1-\frac{\Delta}{2\epsilon_{\alpha}})\\ &+|\phi({\bm{k}})|^{2}(1+\frac{\Delta}{2\epsilon_{\alpha}})\bigg{)},\\ &a_{13}(\alpha,\Delta)=\frac{(-1)^{\alpha}\pi}{\epsilon_{2}^{2}-\epsilon_{1}^{2}}(1-\frac{\Delta}{2\epsilon_{\alpha}})\gamma_{1}\phi({\bm{k}})\cdot\end{split}$$

$$\begin{split}\sigma^{\prime}_{\text{AA}}(\omega)=&\frac{2e^{2}}{\omega}\int_{-\infty}^{\infty}\frac{dE}{2\pi}[f(E)-f(E+\hbar\omega)]\\ &\int\frac{d^{2}{\bm{k}}}{(2\pi)^{2}}4|v_{\bm{k}}|^{2}\bigg{[}A_{11}(E+\hbar\omega)A_{11}(E)\\ &+A_{13}(E+\hbar\omega)A_{13}(E)\bigg{]},\end{split}$$

$$A_{11}(w,\Delta )=\sum _{\alpha =1,2}\frac{\pi }{2}[\delta (w-{ϵ}_{\alpha })+\delta (w+{ϵ}_{\alpha })],$$

$$A_{13}(w,\Delta )=\sum _{\alpha =1,2}\frac{\pi }{2}[\delta (w-{ϵ}_{\alpha })-\delta (w+{ϵ}_{\alpha })]\cdot$$

$$\delta (w)=\frac{1}{\pi }\frac{\Gamma }{w^2+{\Gamma }^2}\cdot$$

$${\sigma }^{\prime \prime }(\omega )=-\frac{2\omega }{\pi }\mathcal{P}{\int }_0^{\infty }\frac{{\sigma }_r({\omega }^{\prime })}{{\omega }^{\prime 2}-{\omega }^2}d{\omega }^{\prime },$$

$$n=n_t+n_b=-\frac{{\epsilon }_r{\epsilon }_0}{q}\delta E=\frac{{\epsilon }_r{\epsilon }_0}{q{t}_t}{V}_t+\frac{{\epsilon }_r{\epsilon }_0}{q{t}_b}{V}_b\cdot$$

$$\begin{array}{rl}n& =\sum _{\alpha }{n}_{\text{electrons}}-n_{\text{holes}}\\ & =\sum _{\alpha }{\int }_0^{\infty }{\rho }_{\alpha }(E)f(E)dE-{\int }_{-\infty }^0{\rho }_{\alpha }(E)[1-f(E)]dE,\end{array}$$

$${\rho }_{\alpha }(E)=4\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\delta (E-{ϵ}_{\alpha }(\mathbf{k}))\cdot$$

$$\begin{array}{r}\delta n=n_1-n_2=\frac{2}{\pi }\int \sum _{\alpha }(|{\psi }_{A1}^{\alpha }(k){|}^2+|{\psi }_{B1}^{\alpha }(k){|}^2\\ -|{\psi }_{A2}^{\alpha }(k){|}^2-|{\psi }_{B1}^{\alpha }(k){|}^2)kdk,\end{array}$$

$$\Delta =q(\overline{E}+E_{int})d=(V_t-V_b)qd/2t_t-\delta n(\Delta )\frac{q^{2}d}{2{\epsilon }_r{\epsilon }_0},$$

$$⟨\tau ⟩=\frac{{\sum }_{s,\alpha }\int {\rho }_{s,\alpha }(E){\tau }_{s,\alpha }(E)(-\frac{df(E)}{dE})dE}{{\sum }_{s,\alpha }\int {\rho }_{s,\alpha }(E)(-\frac{df(E)}{dE})dE}\cdot$$

$$\begin{array}{r}\frac{1}{{\tau }_{s,\alpha }(\mathbf{k})}=\frac{2\pi n_0}{\hslash }\int \frac{d^2{\mathbf{k}}^{\prime }}{(2\pi {)}^2}|⟨V_{s,\alpha }(\mathbf{k},{\mathbf{k}}^{\prime })⟩{|}^2\times \\ (1-\cos ({\theta }_{\mathbf{k},{\mathbf{k}}^{\prime }}))\delta (ϵ(\mathbf{k})-ϵ({\mathbf{k}}^{\prime })),\end{array}$$

$$⟨V_{s,\alpha }(\mathbf{k},{\mathbf{k}}^{\prime })⟩={\psi }_{s,\alpha }^{†}(\mathbf{k}){\psi }_{s,\alpha }({\mathbf{k}}^{\prime })V(\mathbf{q}),$$

Full text

Bilayer Graphene, Electrostatic gating, Energy asymmetry, Optical conductivity, Plasmon propagation, Terahertz plasmonics

I introduction

The two-dimensional (2D) material graphene offers strong light-matter interaction over a broad electromagnetic spectrum, ranging from several hundreds of GHz to optical frequencies Grigorenko et al. (2012). An interesting and useful attribute of graphene is its ability to support the propagation of surface plasmon polaritons (SPPs). SPPs (referred to as plasmons, hereafter) are collective oscillations of electrons, typically occurring at the interface between a metal and a dielectric, that can be coupled with electromagnetic waves and enable the manipulation of light below its diffraction limitJablan et al. (2013). Plasmon oscillations can occur in various metallic and semiconducting material systems including 2D electron gases Allen Jr et al. (1977); Liu et al. (2008); Theis (1980). Typically, plasmons have a smaller wavelength compared to the wavelength of light at the same frequency allowing significant mode compression.

Plasmons supported on a graphene sheet have a tighter confinement and a higher propagation length relative to those of surface plasmons in metals Koppens et al. (2011); Zhou et al. (2012). Moreover, graphene can also support the propagation of transverse electric (TE) modes Mikhailov and Ziegler (2007), which are absent in the case of a metal/dielectric system. The characteristics of graphene plasmons can also be tuned via electrostatic gating or chemical doping over a wide frequency range. In particular, the electrostatic tuning of plasmons enables reconfigurable plasmonic devices with graphene – a feature that is unique to graphene and not afforded by metal plasmonics.

Recently, plasmon oscillations in monolayer and multilayer graphene heterostructures were experimentally imaged Fei et al. (2015); Menabde et al. (2016); Chen et al. (2012); Woessner et al. (2015); Ju et al. (2011). These experiments pave the way to realize tunable plasmonic devices that are currently being researched and developed including waveguidesVakil and Engheta (2011), nanoantennasJornet and Akyildiz (2013), photodetectorsVicarelli et al. (2012), ocsillatorsRana (2008), and modulatorsLiu et al. (2011). From a theoretical standpoint, plasmons have been studied extensively in gated Rakheja and Sengupta (2016) and ungated Jablan et al. (2009); Hanson (2008) monolayer graphene (MLG) and also in ungated Jablan et al. (2011) bilayer graphene (BLG). However, not much attention has been given to plasmons in BLG structures in which the bandstructure and, therefore, optical properties can be readily tuned via electrostatic gating.

In this paper, we develop theoretical models to study the propagation of TE and TM (transverse magnetic) plasmons in gated BLG structures for various material parameters of graphene and the dielectrics surrounding the graphene sheet. We also highlight the distinguishing features of plasmon propagation in BLG and MLG. While both MLG and BLG exhibit tunable optical properties, the tunability is significantly enhanced in BLG. This enhancement is a result of the electric-field-induced energy asymmetry in BLG structures. Essentially, a vertical electric field in a metal-oxide-graphene geometry, creates an energy asymmetry between the two parallel layers of BLG Castro et al. (2008); McCann (2006); McCann et al. (2007). The energy asymmetry changes the bandstructure of BLG, which in turn alters its optical properties significantly Ohta et al. (2006); Zhang et al. (2009). In addition to the energy asymmetry, we also consider other phenomena, such as the electron scattering rate and the stacking forms, that affect the optical properties of BLG.

The remainder of this paper is organized as follows. In Section II, we derive the plasmon wavevector for both TE and TM polarized modes in a graphene sheet using the solutions of Maxwell’s equations. We use a device setup in which the graphene sheet is modeled as an impedance boundary condition with a finite dynamical surface conductivity. To obtain the dynamical conductivity, we first derive the bandstructure of BLG using tight binding model in Section III. The bandstructure is obtained in terms of the energy asymmetry between the layers and for two common stacking orders, that are AB and AA stacked BLG. The bandstructure is used to solve for the dynamical conductivity of BLG within the Kubo formalism. The effects of gate voltage and electrostatic screening on the carrier concentration and energy asymmetry are included using a self-consistent Hartree approximation. We present the results of our analysis in Section IV in which we quantify the impact of physical parameters of BLG on the plasmon characteristics such as compression factor, propagation length, and mode confinement. The key findings of the paper are summarized in Section V.

II Plasmon Waves in Graphene

To obtain plasmon wavevector in graphene heterostructures, we consider the device setup shown in Fig. 1 in which an infinite 2D sheet of graphene is immersed in a uniform, linear dielectric medium with a relative dielectric constant of $\varepsilon_{r}$. We note that graphene in this device setup could be either a monolayer or a bilayer sheet. The propagating electric and magnetic fields in the $z$-direction for TE-polarized plasmon mode are

[TABLE]

and for TM-polarized plasmon mode are

[TABLE]

where time dependence term $e^{j\omega t}$ is dropped and $k_{z}(k_{x})$ is the wavevector in $z$-($x$-) direction as illustrated in Fig. 1.

If we plug Eqs. (1) and (2) into Helmholtz equation ($\nabla^{2}{\bm{E}}+k^{2}_{\text{light}}{\bm{E}}=0$), we get $k_{x}^{2}+k_{z}^{2}=k^{2}_{\text{light}}$, where $k_{\text{light}}=\omega\sqrt{\varepsilon_{r}\varepsilon_{0}\mu_{0}}$ is the wavevector of light in the medium. Here, $\omega$ is the radial frequency, and $\mu_{0}$ and $\varepsilon_{0}$ are the permeability and permittivity, respectively, of free space.

By solving Maxwell’s equations in the device setup of Fig. 1 and applying appropriate boundary conditions (impedance boundary with a finite surface conductivity for graphene), we obtain the following dispersion relations Mikhailov and Ziegler (2007); Hanson (2008)

[TABLE]

where $\sigma(\omega)$ is a complex quantity that represents surface conductivity of the graphene sheet, and $c=1/\sqrt{\epsilon_{0}\mu_{0}}$ is the velocity of light in free space. By replacing $\sigma(\omega)=\sigma^{\prime}(\omega)-j\sigma^{\prime\prime}(\omega)$ and $\eta=\sqrt{\mu_{0}/\varepsilon_{r}\varepsilon_{0}}$ in Eq. (3) and rearranging the equations, we obtain the plasmon wavevectors

[TABLE]

for TE mode and

[TABLE]

for TM mode. From Eqs. (4) and (5), it is evident that the propagation of TE (TM) mode requires $\sigma^{\prime\prime}(\omega)$ to be negative (positive). This is due to the fact that for propagating plasmon modes, the imaginary part of $k_{z}$, which represents the propagation loss, must be positive.

While plasmons propagate in both z-direction (longitudinal) and x-direction (transverse), they are expected to decay more rapidly in the x-direction. The decay of plasmons can be represented using a length scale called the propagation length. The propagation length is the distance at which the intensity of the wave degrades to $1/e$ of its initial value, where $e=2.718$ is the base of natural logarithm. For given wavevectors $k_{z}=k^{\prime}_{z}-jk^{\prime\prime}_{z}$ and $k_{x}=k^{\prime}_{x}-jk^{\prime\prime}_{x}$, the propagation length is obtained as $1/2k^{\prime\prime}_{z}$ and $1/2k^{\prime\prime}_{x}$ for longitudinal and transverse directions, respectively. The measure of confinement of the plasmon wave to the graphene sheet can be represented by its propagation length in the transverse direction. A lower transverse propagation length means higher mode confinement on the graphene surface. Another important characteristic of a plasmon wave is its compression factor which determines how compressed the plasmon wavelength is relative to the wavelength of light calculated in the absence of the graphene sheet. The compression factor is mathematically given as the ratio of the longitudinal phase constant to the phase constant of light $k^{\prime}_{z}/k_{\text{light}}$.

An important ingredient in obtaining the plasmon dispersion is the dynamical conductivity, $\sigma(\omega)$, which depends on the bandstructure of graphene, electron scattering rate, Fermi level, and the operating frequency. In the next section, we present the model for $\sigma(\omega)$ for AB and AA stacked BLG. For brevity, we do not include the model for $\sigma(\omega)$ in MLG. Interested readers are referred to Refs. Falkovsky and Varlamov, 2007; Falkovsky and Pershoguba, 2007; Falkovsky, 2008; Hanson, 2008; Horng et al., 2011 for detailed discussion pertaining to MLG conductivity.

III Model of Bilayer Graphene

III.1 Bandstructure and Stacking forms

To study the optical properties of graphene, in particular its dynamical conductivity, we start with the calculation of bandstructure of BLG. Here, we investigate two of the most common stacking forms namely AB (Bernal stacked) and AA stacked BLG. The natural stacking form of BLG is AB, which is also found in graphite. Although AA stacked BLG is likely to be metastable Rozhkov et al. (2016), it is still interesting to study it from a theoretical point of view. Figure 2 illustrates the configuration of both the stacking forms. Each graphene layer consists of a triangular Bravais lattice with a basis of two atoms which we denote as ‘A’ and ‘B’. In an AB staked BLG, ‘A’ atoms of one layer align with ‘B’ atoms of the other layer and all other atoms lie at the center of the others layer’s hexagons. In an AA stacked BLG, all the atoms in one layer align exactly with other layer’s atoms correspondingly.

By utilizing tight binding model, we can write the Hamiltonian of an AB stacked BLG that is exposed to a vertical electrostatic field. The electric field introduces an energy asymmetry between the two layers, which is denoted as $\Delta$. It lowers the energy of the first layer by $\Delta/2$ and raises the energy of the second layer by the same amount. The Hamiltonian for an AB stacked BLG is written in the following form Nicol and Carbotte (2008)

[TABLE]

Here, $\gamma_{1}=0.4$ eV is the interlayer hopping energy between $B_{1}$ and $A_{2}$ sitesZhang et al. (2008), and $|\phi({\bm{k}})|=|\hbar v_{F}ke^{i\theta}|$ is the energy dispersion of MLG in the vicinity of the K point of the Brillouin Zone within a radius of $3$ eV. Also $v_{F}=\sqrt{3}\gamma_{0}b/2\hbar$ denotes the Fermi velocity, $\gamma_{0}=3$ eV is nearest neighbor hopping energy within each layer, and $b=0.246$ nm is graphene’s lattice constant. The eigenvalues of the Hamiltonian describe the bandstructure, which in the case of AB stacked BLG is given as

[TABLE]

Here, $\alpha=1,2$ is the subband index. We can see that for $\gamma_{1}=0$ and $\Delta=0$ this energy dispersion reduces to that of MLG. Figure 3a shows the bandstructure of AB stacked BLG for two different values of $\Delta$. Contrary to MLG that exhibits a gapless linear energy dispersion relationship, AB stacked BLG shows a nonlinear energy dispersion with a finite energy bandgap. We note that the energy bandgap ($E_{g}=\Delta\gamma_{1}/\sqrt{\Delta^{2}+\gamma_{1}^{2}}$) occurs at $k\not=0$ and is smaller than $\Delta$. Similarly, we can write the Hamiltonian of AA stacked BLG for a finite electrostatic bias as followsTabert and Nicol (2012)

[TABLE]

The eigenvalues of this matrix describe the bandstructure of AA stacked BLG, which is given as

[TABLE]

The above expression shows that the energy asymmetry does not open a bandgap in the case of AA stacked BLG, and $\Delta$ can be treated as a correcting factor to $\gamma_{1}\rightarrow\sqrt{\gamma_{1}^{2}+\Delta^{2}/4}$. Figure 3b shows the bandstructure of AA stacked BLG for two different values of $\Delta$.

III.2 Optical Conductivity

We can write the real part of frequency dependent conductivity $\sigma^{\prime}(\omega)$ using Kubo formalism. Following procedure described in Ref. Nicol and Carbotte, 2008, for an AB stacked BLG we have

[TABLE]

Here, $\mu$ is the Fermi level, $f(E)=1/(1+\exp((E-\mu)/k_{B}T))$ is the Fermi-Dirac distribution function with $k_{B}$ as the Boltzmann’s constant and $T$ as the temperature, and $|v_{\bm{k}}|$ is the velocity of electrons and is approximated by $v_{F}$. The inner integral is being evaluated over the K point in Brillouin zone. The $A_{ij}$ functions are the spectral representations of Green’s functions. They are related to the Green’s functions through following equation:

[TABLE]

and can be written in the form of

[TABLE]

in which $\alpha$ denotes the $\alpha^{\text{th}}$ subband, and $\epsilon_{\alpha}$ is the subband energy. The coefficients $a_{ij}(\alpha,\Delta)$ are derived from Eqs. (11), (12), and the Hamiltonian through $\hat{G}^{-1}=z\hat{I}-\hat{H}$.

[TABLE]

In a similar fashion, we can write the conductivity of AA stacked BLG which is described in Ref. Tabert and Nicol, 2012. The real part of conductivity is written as follows

[TABLE]

in which the spectral functions are

[TABLE]

To incorporate the effect of electron scatterings in graphene, the delta functions in Eqs. (12) and (15) are broadened using Lorentzians defined as

[TABLE]

The electron scattering rate is represented by $\Gamma$ which is generally a function of energy and the bandstructure. In the present work, we use a constant value that represents the energy-dependent scattering rate averaged over energy. The impact of bandstructure and energy asymmetry on the scattering rate is discussed in Appendix A. For all simulations in this work, $\Gamma$ is chosen as $20$ meV, unless otherwise stated, which is a typical experimental value of scattering rate in a graphene sheetMonteverde et al. (2010). The imaginary part of conductivity is derived from its real part using Kramers-Kronig relation

[TABLE]

where $\mathcal{P}$ denotes the principal value of the integral.

Figure 4 shows the optical conductivity of AB and AA stacked BLG for different values of $\mu$ and $\Delta$ using an energy-independent value of electron scattering rate $\Gamma=20$ meV at room temperature $T=300$ K. The plots are normalized to the universal conductivity of graphene $\sigma_{0}=q^{2}/4\hbar$. In Fig. 4 (a, c), the real parts of the optical conductivity ($\sigma^{\prime}(\omega)$) of AB and AA stacked BLG are illustrated for different values of the Fermi level, $\mu$ at $\Delta=0$ eV. We can see that for a fixed energy asymmetry, the optical conductivity can be tuned by the Fermi level. For instance in the case of AA stacked BLG, the interband transition is shifted in frequency as the Fermi level changes. In the case of AB stacked BLG, there is a strong absorption peak at frequencies around $\gamma_{1}$ due to interband transitions. Figures 4 (e, g) illustrate $\sigma^{\prime}(\omega)$ for different values of $\Delta$ at a fixed Fermi level $\mu=0.3$ eV. From these plots we can see that the impact of energy asymmetry on the optical conductivity of BLG is quite significant specially in the case of AB stacked BLG where the absorption peak shifts toward higher frequencies as $\Delta$ increases. This phenomenon can be seen as an additional degree of freedom for tuning the optical conductivity of BLG that is absent in the case of MLG. At frequencies close to the visible region, the real part of conductivity converges to $2\sigma_{0}$ as expectedNair et al. (2008). The imaginary parts of the optical conductivity $\sigma^{\prime\prime}(\omega)$ are plotted in Figs. 4 (b, d, f, h). These figures show the existence of different frequency bands where TE plasmons ($\sigma^{\prime\prime}(\omega)<0$) and TM plasmons ($\sigma^{\prime\prime}(\omega)>0$) can be supported and tuned in AA and AB stacked BLG structures.

III.3 Gate Voltage and Electrostatic Screening

Experimental observations Zhang et al. (2009); Oostinga et al. (2008); Ohta et al. (2006); Castro et al. (2007) as well as theoretical studies McCann (2006); McCann et al. (2007); Castro et al. (2008); Nilsson et al. (2007); Min et al. (2007); Avetisyan et al. (2009) have demonstrated that applying an electrostatic field perpendicular to the BLG plane forms an energy asymmetry between the two layers of graphene and consequently alters its bandstructure. In Sec. III.1 and III.2, we discussed that the presence of energy asymmetry creates a bandgap in the bandstructure of AB stacked BLG, while it modifies the interlayer hopping energy $\gamma_{1}$ of AA stacked BLG. These phenomena introduce an additional degree of freedom (along with the Fermi level) for tuning the optical conductivity of BLG. This eletctric-field-induced energy asymmetry can be realized by gating the graphene layers as shown in Fig. 5. We use the dual gate structure shown in this figure to obtain the energy asymmetry, $\Delta$, in terms of top and bottom gate voltages ($V_{t}$ and $V_{b}$ respectively).

In this metal-oxide-graphene geometry, the gate voltage applied at the metal electrode, modulates the energy asymmetry and the carrier density in BLG via the electrostatic field in the oxide. The physical thickness of the oxide is denoted as $t_{t}$ for the top oxide and $t_{b}$ for the bottom oxide. The parameter $d=0.35$ nm is the interlayer distance of BLG and is equal to the thickness of each layer. The use of dual gates allows us to independently control the Fermi level and energy asymmetry in the BLG. The fields of top and bottom gates are $E_{t}$ and $E_{b}$, respectively. It has been shown that the difference between these two fields $\delta E=E_{b}-E_{t}$ determines the Fermi level, and their average $\overline{E}=(E_{b}+E_{t})/2$ introduces the energy asymmetryZhang et al. (2009). For given gate voltages, $V_{t}>0$ and $V_{b}<0$, by utilizing Gauss’s law, we derive the total carrier density induced on the graphene layers as

[TABLE]

Here $n_{t}$ and $n_{b}$ are the carrier densities induced by the the top and bottom gates respectively and $q$ is the magnitude of the electron charge. We note that for the negative charged carriers, $n$ is positive. The Fermi level is related to the total carrier density by the following equation

[TABLE]

where $\rho_{\alpha}(E)$ is the density of states of the $\alpha^{\text{th}}$ subband and is defined as follows Ashcroft and Mermin (1976)

[TABLE]

We have included the valley and spin degeneracies using a factor of four in the above equation. Using Eqs. (18) and (19), we find the Fermi level as a function of the top and bottom gate voltages. Figure 6 shows how the Fermi level of BLG varies with respect to $V_{t}+V_{b}=t_{t}\delta E$ for a fixed $\Delta$. Here we have assumed that top and bottom oxides are made of SiO2 with permittivity $\varepsilon_{r}=3.9$ and that they have equal thicknesses $t_{t}=t_{b}=300$ nm. As we can see from the plots, the Fermi level of monolayer and AB stacked BLG show a higher sensitivity to the gate voltage than the Fermi level of AA stacked BLG. Also, at high gate voltages, different energy asymmetries lead to different curves. The inset illustrates the curves for gate voltages close to zero. The average electric field induced by the gate contacts ($\overline{E}=(V_{t}-V_{b})/2t_{t}$) polarizes the graphene layers, which leads to an internal electric field given as $E_{int}=-q\delta n/2\varepsilon_{r,g}\varepsilon_{0}$ with $\delta n$ as the charge imbalance between the two layers and $\varepsilon_{r,g}=1$ as the relative permittivity of BLG. The internal electric field acts to cancel out the external electric field. This effect, known as electrostatic screening, plays a major role in realistic multilayer graphene devices and is carefully modeled in our analysis. Due to the finite electrostatic screening in BLG, the electron density on the top layer $n_{1}=(n+\delta n)/2$ is not equal to the electron density on the bottom layer $n_{2}=(n-\delta n)/2$. We note that the charge conservation is satisfied, $n_{1}+n_{2}=n_{t}+n_{b}$. The charge imbalance $\delta n=n_{1}-n_{2}$ is given as the relative weights of electron wavefunctions as follows McCann et al. (2007)

[TABLE]

where \big{(}\psi^{\alpha}_{A1}(k),\psi^{\alpha}_{B2}(k),\psi^{\alpha}_{A2}(k),\psi^{\alpha}_{B1}(k)\big{)} represents the $\alpha^{\text{th}}$ eigenvector of the Hamiltonian in Eqs. (6) and (8), and the integral is evaluated over all occupied states. Therefore, we can write the energy asymmetry as

[TABLE]

which is the same as the self-consistent Hartree approximation in Ref. Castro et al., 2008. By setting $V_{t}+V_{b}$ to a constant value and applying a differential voltage $\delta V_{t,b}$ to the gate contacts ($V_{t}$ becomes $V_{t}+\delta V_{t,b}/2$ and $V_{b}$ becomes $V_{b}-\delta V_{t,b}/2$), we fix the Fermi level and can independently study the impact of gate voltage on the energy asymmetry. Figure 6 illustrates the energy asymmetry as a function of differential gate voltage, $\delta V_{t,b}$, for a variety of fixed Fermi levels. As we can see the energy asymmetry has a linear dependence on $\delta V_{t,b}$. Moreover, for the given gate voltages, AA stacked BLG shows a higher energy asymmetry than AB stacked BLG. We note that the curves are also valid for larger differential gate voltages as long as the oxide field is below the dielectric breakdown strength of SiO2, which is larger than $1$ V/nm.

IV Results

In this section, we discuss the features of propagating TM and TE plasmons, such as compression factor, mode confinement, and propagation length, in various graphene structures. From the perspective of graphene plasmonic devices (antennas and waveguides), it is desirable to simultaneously achieve confined (low transverse propagation length $1/2k^{\prime\prime}_{x}$) waves with a low loss (high longitudinal propagation length $1/2k^{\prime\prime}_{z}$) and a high compression factor $k^{\prime}_{z}/k_{\text{light}}$. We use the complex wavevector derived in Sec. II and combine it with the conductivity model from Sec. III.2 in the presence of gate effects (Sec. III.3) for analysis in this section.

In Fig. 7, we plot the following plasmon characteristics versus the operating frequency: (i) relative longitudinal propagation length $1/2\lambda k^{\prime\prime}_{z}$ normalized to the wavelength of light $\lambda=2\pi/k_{\text{light}}$, (ii) relative phase constant $k_{z}^{\prime}/k_{\text{light}}$ normalized to $k_{\text{light}}$, and (iii) relative transverse propagation constant $1/2\lambda k^{\prime\prime}_{x}$ normalized to $\lambda$. The results in this figure are obtained for AB and AA stacked BLG, and MLG in the terahertz and infrared regimes for $\mu=0.3$ eV and $\Delta=0$ eV. The conductivity of MLG is derived using Green’s functions approach described for BLG in Sec. III.2. It can also be derived using the procedure outlined in Ref. Falkovsky and Pershoguba, 2007. In the terahertz frequency band (left column of Fig. 7), all the graphene structures support TM mode propagation with similar propagation characteristics. As the frequency increases the plasmon waves become more compressed and more confined to the graphene sheet. While MLG has a relatively higher compression factor, its propagation length is smaller than that of BLG. In the infrared regime (right column of Fig. 7), there exists separate frequency bands that allow the propagation of both TM and TE modes in graphene structures. Generally, TE mode in both BLG and MLG has a phase constant ($k^{\prime}_{z}$) very close to the phase constant of light ($k_{light}$) and shows a very high propagation length in longitudinal direction but they are very loosely bounded to the graphene sheet. Conversely, TM mode mostly shows a high compression factor, but very low propagation length which is two orders of magnitude lower than the wavelength of light. However, they are more confined in the transverse direction than TE modes and so are better guided by the graphene structure. The trade-off between the highly confined TM mode and the highly propagating TE mode is a general characteristic of plasmons and has been reported previouslyHanson (2008); Jablan et al. (2009); Maier (2007). The reason is that high confinement to the graphene sheet implies that a large amount of plasmon wave energy is confined to the sheet which has a higher loss than the surrounding dielectric medium.

To study the impact of chemical potential, we plot the various plasmon characteristics in graphene structures as a function of the chemical potential $\mu$ for a constant frequency in the terahertz band and a constant frequency in the infrared band in Fig. 8.

The gate voltages corresponding to the Fermi levels are labeled on the top axis. We note that very high Fermi levels are not achievable in practice by mere gating since it would require electric fields higher than the break down field of SiO2. In our device setup, oxide voltages cannot exceed $300$ V, however to obtain high Fermi levels, one can utilize chemical doping which would basically act as a gate with a constant induced carrier density. For the $1$ THz frequency (left column of Fig. 8), it is evident that as the Fermi level increases, AB stacked BLG and MLG show an increase in both longitudinal and transverse propagation lengths and a decrease in compression factor. However, plasmons in AA stacked BLG show little sensitivity to changes in the Fermi level. Moreover, similar to the trade-off between TM and TE modes, we observe a trade-off between compression factor and propagation length of TM mode plasmons as we chcange the Fermi level. At infrared frequencies (right column of Fig. 8), by changing the Fermi level, we can switch the mode of propagation between the highly confined TM mode and the loosely bounded TE mode.

In Sec. III.2 we discussed that in addition to the Fermi level, the energy asymmetry in BLG can be utilized as an extra degree of freedom to tune its optical conductivity. To study the impact of energy asymmetry, we plot the propagation characteristics of plasmons in BLG versus frequency for two values of energy asymmetry ($\Delta=0$ and $0.8$ eV) and a fixed Fermi level in Fig. 9. The Fermi level needs to be sufficiently far from the charge neutrality point so that there are enough carriers to enable plasmon oscillations. Specifically, due to the induced bandgap in AB stacked BLG, there exists a plasmon cut-off region ($\mu<0.1$ eV) where the carrier concentration is not sufficient for plasmon oscillations.Fei et al. (2015); Chen et al. (2012); Woessner et al. (2015). Here, the Fermi level is fixed at $0.3$ eV. Each curve in Figs. 9a and 9b corresponds to a different value of scattering rate for both $\Delta=0$ eV (solid) and $\Delta=0.8$ eV (dashed). As we see from the figure, scattering rate has an enormous impact on plasmon waves decay. Hence, the ability to produce high quality samples with low impurity and defect concentrations becomes an essential challenge of building high performance graphene based plasmonic devices.

Interestingly, the propagation length of plasmons in AB and AA stacked BLG can be improved by increasing $\Delta$. Plasmons in AB stacked BLG can propagate 2 to 5 times longer, depending on the frequency band and the scattering rate. Similarly, plasmons in AA stacked BLG can propagate 1.5 to 100 times longer. The improvement in propagation length is mainly due to the reduction in absorption ($\sigma^{\prime}(\omega)$) caused by supression of interband transitions such as second valance subband to second conduction subband ($\epsilon_{-,2}\rightarrow\epsilon_{+,2}$) and first conduction subband to second conduction subband ($\epsilon_{+,1}\rightarrow\epsilon_{+,2}$). This a significant benefit of BLG as it will enable low energy and low leakage plasmonic antennas and waveguides for THz and low infrared bands communication. However, these improvements come at the expense of smaller compression factor (Fig. 9c, d) as mentioned earlier.

Finally we illustrate how the plasmon characteristics are affected by the medium dielectric. Figure 10 illustrates the plasmon characteristics as a function of the relative permittivity of dielectric medium for TM and TE modes of propagation. As we can see, while increasing dielectric constant, lowers propagation length of the TM mode, it increases propagation length of the TE mode. This effect can also be seen directly from Eqs. (4) and (5) where the dependence on dielectric constant appears in $\eta$. Moreover, higher compression factors can be achieved for TM mode by increasing dielectric constant.

V Conclusions

The focus of this paper is to quantify the impact of material parameters on the tunability of plasmon propagation characteristics in bilayer graphene (BLG) structures. Toward this end, we consider a dual gate metal-oxide-graphene geometry to independently tweak the Fermi level and energy asymmetry in BLG through the gate terminal voltages. We numerically obtain the relationship between energy asymmetry and the applied gate voltage using the self-consistent Hartree approximation. The optical conductivity of BLG is modeled as a function of Fermi level and energy asymmetry using the Kubo formalism. The propagation of surface plasmons is studied numerically to quantify their propagation length, mode confinement, and compression factor in various graphene structures including AB stacked and AA stacked BLG, and monolayer graphene (MLG). Important trade-offs between propagation length and mode confinement that must be considered for device design are discussed. Impact of scattering rate on the propagation length of plasmons is also quantified. Finally, we demonstrate that the propagation length of TM plasmons in BLG can be improved by increasing the energy asymmetry. Such an advantage is absent in MLG. Overall, the results presented in this work provide insight into the design and optimization of plasmonic structures utilizing graphene as their platform.

Appendix A Impact of energy asymmetry on scattering rate

In general, the scattering rate depends on the bandstructure and the Fermi velocity. A larger energy asymmetry would reduce the Fermi velocity of electrons. However, Fermi velocity is not the only factor that affects scattering rate and correspondingly the various figures of merit of plasmons. Factors including polarizability and screening, which depend on $\Delta$, are important to consider, particularly in the case of charged impurity scatterers. These factors could outweigh the reduction of Fermi velocity and even decrease the scattering rate. To study the nontrivial dependence of scattering rate on $\Delta$, we numerically calculate the average scattering rate, $\Gamma=1/2\braket{\tau}$, where $\braket{\tau}$ is the average carrier relaxation time and is given as Hwang and Sarma (2008a)

[TABLE]

Here, $\tau_{s,\alpha}(E)$ is the energy-dependent relaxation time corresponding to the band $s$ and subband $\alpha$. Density of states is denoted as $\rho_{s,\alpha}(E)$, and $f(E)$ is the Fermi-Dirac distribution function. The k-dependent relaxation time $\tau_{s,\alpha}({\bm{k}})=\tau_{s,\alpha}(E)|_{E=\epsilon_{s,\alpha}({\bm{k}})}$ is written using Boltzmann approximation as followsSarma et al. (2010)

[TABLE]

where $n_{0}$ is the density of scatterers, $\theta_{\bf k,k^{\prime}}$ is the angle between the incident wavevector ${\bf k}$ and the scattered wavevector ${\bf k^{\prime}}$, and $\braket{V_{s,\alpha}({\bf k,k^{\prime}})}$ is the matrix element of the disorder potential. $\braket{V_{s,\alpha}({\bf k,k^{\prime}})}$ depends on the scattering mechanism and the electron wavefunctions $\psi_{s,\alpha}({\bm{k}})$ and is written as

[TABLE]

Here, ${\bm{q}}={\bm{k}^{\prime}}-{\bm{k}}$ and $V({\bm{q}})=\int V({\bf r})e^{i{\bf q}\cdot{\bf r}}d{\bf r}$ is the Fourier transform of $V(\bm{r})$, which is the scattering potential.

We consider long-range Coulomb scattering and short-range defect scattering mechanisms. In the presence of screening, for charged impurities located at distance $d$ from the BLG sheet we have $V({\bm{q}})=v(q)\exp(-qd)/\varepsilon(q)$, where $v(q)=e^{2}/2\varepsilon_{r}\varepsilon_{0}q$ is the bare coulomb potential. The dielectric function is $\varepsilon(q)=1+v(q)\Pi(q)$, and $\Pi(q)$ represents the polarizability of BLG. The details of calculation of polarizability are provided in Refs. Hwang and Sarma, 2008b and Lv and Wan, 2010. For defect scattering $V({\bm{r}})$ is a delta function in space; therefore, its Fourier transform is constant and independent of the dielectric function and polarizability i.e. $V({\bm{q}})=V_{0}$. Here we assume that $V_{0}=10\,$eVnm2(Ref. Sarma et al., 2011). Figure A1 illustrates k-dependent scattering rates, evaluated numerically using Eq. 24, for both coulomb and defect scattering mechanisms, and for two values of energy asymmetry (i.e. $\Delta=0$ and $\Delta=0.8\,$eV). Material-specific parameters used to obtain the results reported in Fig. A1 are provided in the figure caption. The figure highlights the opposite behaviors of Coulomb and defect scattering in response to introducing a large $\Delta$. For Coulomb scattering, scattering rate decreases with an increase in momentumMonteverde et al. (2010) because of the increase in screening. Due to the absence of screening in the case of defect scattering, scattering rate increases with wavevector.

To study the impact of energy asymmetry on the average scattering rate, we use the results from Fig. A1 in Eq. 23. The average scattering rate for both Coulomb and defect scattering mechanisms versus $\Delta$ is plotted in Fig. A2. We see that as $\Delta$ increases from [math] to $0.8\,$eV, the average scattering rate increases by 70% for the case of defect scatterers, whereas it decreases by 30% for the case of Coulomb scatterers. The reason for the unexpected behavior of Coulomb scatterers is the nontrivial behavior of polarizability. Our calculations show that for $\mu=0.3\,$eV, introducing a large energy asymmetry will increase the polarizability of BLG and consequently its dielectric function. This is the effect that outweighs the reduction in Fermi velocity and leads to a smaller average scattering rate overall in a sample dominated by Coulomb scatterers.

Average scattering rate determines the intraband optical absorption and, therefore, the propagation length of plasmons in BLG. Depending on the dominant scattering mechanism, the propagation length could be improved or degraded. In some experimental works, defect scattering is believed to be dominantMonteverde et al. (2010); Sarma et al. (2010), in some other works, it is Coulomb scattering that limits carrier transportZhu et al. (2009); Xiao et al. (2010). For Coulomb-dominated samples, given the specific Fermi level of $0.3\,$eV, introducing a large $\Delta$ will reduce the average scattering rate and consequently increase the plasmon propagation length, while for defect-dominated samples the opposite effect will happen. To elucidate the situation, we plot the plasmon propagation length in Fig. A3 in the absence and presence of energy asymmetry for (i) constant scattering rate, (ii) Coulomb-dominated scattering, and (iii) defect-dominated scattering. For comparison, we choose disorder concentration in a way that for the gapless BLG, all the three cases have the same average scattering rate i.e. $\Gamma=20\,$meV (or equivalently $4.84\times 10^{12}\,$1/s). The disorder concentrations are $n_{i}=9.57\times 10^{16}\,$m*-2* for Coulomb-dominated scattering and $n_{d}=0.012\times 10^{16}\,$m*-2* for defect-dominated scattering. Introducing an energy asymmetry of $\Delta=0.8\,$eV will increase the average scattering rate of defect-dominated case to $\Gamma_{d}=34\,$meV which will degrade the propagation length relative to the constant scattering case. On the other hands, it will decrease the average scattering rate of Coulomb-dominated case to $\Gamma_{d}=14\,$meV which will improve the propagation length relative to the constant scattering case.

Acknowledgements.

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