Long-lived domain wall plasmons in gapped bilayer graphene

TL;DR

This paper predicts the existence of long-lived, tunable plasmons along topological domain walls in gapped bilayer graphene, offering a new platform for low-dissipation plasmonic applications.

Contribution

It introduces the concept of domain wall plasmons in bilayer graphene, highlighting their exceptional lifetimes and unique properties distinct from bulk plasmons.

Findings

Domain wall plasmons can propagate at zero bulk density.

These plasmons have lifetimes up to a few picoseconds at room temperature.

They exhibit tunable frequencies and sub-wavelength confinement.

Abstract

Topological domain walls in dual-gated gapped bilayer graphene host edge states that are gate- tunable and valley polarized. Here we predict that plasmonic collective modes can propagate along these topological domain walls even at zero bulk density, and possess a markedly different character from that of bulk plasmons. Strikingly, domain wall plasmons are extremely long-lived, with plasmon lifetimes that can be orders of magnitude larger than the transport scattering time in the bulk. While most pronounced at low temperatures, long domain wall plasmon lifetimes persist even at room temperature with values up to a few picoseconds. Domain wall plasmons possess a rich phenomenology including a wide range of frequencies (up to the mid-infrared), tunable sub-wavelength electro-magnetic confinement lengths, as well as a valley polarization for forward/backward propagating modes. Its unusual features render them a new tool for realizing low-dissipation plasmonics that transcend the restrictions of the bulk.

Figures

a in the main text). At high temperature, scattering from the bulk dominates the $\tau_{p}$ because $\gamma_{\rm tr}$ is about two orders of magnitude larger than $\gamma_{v}$. In Fig. 2c of main text, we displayed that an increase of $\kappa$ prolongs the DWP lifetime when $\omega$ is fixed. The dependence of lifetime on $\kappa$ can be discerned from Eq. ([S-18](#Sx5.E18)) by approximating $|q|\approx\omega/v_{0}\sqrt{1+\eta}$. This approximation was obtained from the wavevector of $T=0$ dispersion for DWP in Eq. ([8](#S0.E8)) of the main text (we have assumed that $\omega>\omega_{b}$). Plug

Equations

$${\rho }^{\nu }(\mathbf{r},t)$$

$${\mathbf{j}}^{\nu }(\mathbf{r},t)$$

$${\partial }_t{\rho }_b^{\nu }(\mathbf{r},t)+\nabla \cdot {\mathbf{j}}_b^{\nu }(\mathbf{r},t)=0,{\mathbf{j}}_b^{\nu }(\mathbf{r},t)={\mathbf{\sigma }}^{\nu }[-\nabla \varphi (\mathbf{r},t)]$$

$${\partial }_t{\rho }_e^{\nu }-\nu v_0{\partial }_y{\rho }_e^{\nu }+\mathcal{G}{\mathbf{j}}_b^{\nu }\cdot \hat{\mathbf{x}}=-{\gamma }_v(\delta {\rho }_e^{\nu }-\delta {\rho }_e^{-\nu }),$$

$$\varphi (\mathbf{r},t)=\int d{\mathbf{r}}^{\prime }U(\mathbf{r},{\mathbf{r}}^{\prime })\delta \rho ({\mathbf{r}}^{\prime },t),U(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{\kappa |\mathbf{r}-{\mathbf{r}}^{\prime }|},$$

$$\partial_{x}{\widetilde{\phi}}_{q}\big{|}_{0^{+}}-\partial_{x}{\widetilde{\phi}}_{q}\big{|}_{0^{-}}=\left(\partial_{x}U_{q}\big{|}_{0^{+}}-\partial_{x}U_{q}\big{|}_{0^{-}}\right)\delta{\widetilde{\rho}}_{q,e}.$$

$$\delta {\rho }_e=-\frac{({\mathcal{M}}^{K^{\prime }}+{\gamma }_v)\mathcal{G}{\mathbf{j}}_b^K\cdot \hat{\mathbf{x}}+({\mathcal{M}}^K+{\gamma }_v)\mathcal{G}{\mathbf{j}}_b^{K^{\prime }}\cdot \hat{\mathbf{x}}}{{\mathcal{M}}^K{\mathcal{M}}^{K^{\prime }}-{\gamma }_v^2},$$

$$({\partial }_x^2-2q^2){\tilde{\varphi }}_q(x)=\frac{-4\pi }{\kappa }|q|\delta {\tilde{\rho }}_q(x).$$

$$\omega =v_0|q|\sqrt{1+\eta },\eta =4\sqrt{2}\pi {\sigma }_H/v_0\kappa .$$

$${\tilde{\varphi }}_q(x)={\varphi }_0{e}^{-k_0|x|},k_0=\sqrt{2}|q|{\left(\frac{{\tilde{\omega }}^2-{\omega }_b^2+i{\gamma }_{\mathrm{tr}}\tilde{\omega }}{{\tilde{\omega }}^2-2{\omega }_b^2+i{\gamma }_{\mathrm{tr}}\tilde{\omega }}\right)}^{1/2},$$

$$k_0\left[{\tilde{\omega }}^2+2i\tilde{\omega }{\gamma }_v-(v_{0}q{)}^2+ϵ|q|(i\tilde{\omega }-2{\gamma }_v){\sigma }_{xx}\right]=ϵ{\sigma }_H{v}_0|q{|}^3,$$

$$\frac{1}{{\tau }_p}=\frac{1}{{\tau }_v}+\frac{1}{{\tau }_{\mathrm{tr}}}\left(\frac{2{\omega }_b^2}{{\omega }^2}\right).$$

$$\delta {\rho }_e^{\nu }(y,t)=\frac{1}{\omega +\nu v_{0}q}\left[\frac{2k_{0}D(\theta )}{\omega }-\nu 2{\sigma }_{H}q\right]{\varphi }_0{e}^{i(qy-\omega t)}.$$

$$-i\omega \delta {\tilde{\rho }}_{q,b}(x)+{\sigma }_{xx}(-{\partial }_x^2+q^2){\tilde{\varphi }}_q(x)=0,$$

$$({\partial }_x^2-2q^2){\tilde{\varphi }}_q(x)=\frac{-4\pi }{\kappa }|q|\delta {\tilde{\rho }}_{q,b}(x).$$

$$(k_0^2-2q^2){\tilde{\varphi }}_q(x)=\frac{4\pi }{i\omega \kappa }|q|{\sigma }_{xx}(k_0^2-q^2){\tilde{\varphi }}_q(x).$$

$$\left[k_0^2\left(1-\frac{2{\omega }_b^2}{i\omega {\gamma }_{\mathrm{tr}}+{\omega }^2}\right)-2q^2\left(1-\frac{{\omega }_b^2}{i\omega {\gamma }_{\mathrm{tr}}+{\omega }^2}\right)\right]{\tilde{\varphi }}_q=0.$$

$$k_0=\sqrt{2}|q|{\left(\frac{{\omega }^2-{\omega }_b^2+i{\gamma }_{\mathrm{tr}}\omega }{{\omega }^2-2{\omega }_b^2+i{\gamma }_{\mathrm{tr}}\omega }\right)}^{1/2},$$

$${\partial }_t{\rho }_e^K-v_0{\partial }_y{\rho }_e^K+\mathcal{G}{\mathbf{j}}_b^K\cdot \hat{\mathbf{x}}=$$

$${\partial }_t{\rho }_e^{K^{\prime }}+v_0{\partial }_y{\rho }_e^{K^{\prime }}+\mathcal{G}{\mathbf{j}}_b^{K^{\prime }}\cdot \hat{\mathbf{x}}=$$

$$\left(\begin{array}[]{c}\delta\rho_{e}^{K}\\ \delta\rho_{e}^{K^{\prime}}\end{array}\right)=\frac{-1}{\mathcal{M}^{K}\mathcal{M}^{K^{\prime}}-\gamma_{v}^{2}}\left(\begin{array}[]{cc}\mathcal{M}^{K^{\prime}}&\gamma_{v}\\ \gamma_{v}&\mathcal{M}^{K}\end{array}\right)\left(\begin{array}[]{c}{\mathcal{G}}{\bf j}_{b}^{K}\cdot\hat{{\bf x}}\\ {\mathcal{G}}{\bf j}_{b}^{K^{\prime}}\cdot\hat{{\bf x}}\end{array}\right),$$

$$\delta {\rho }_e=-\frac{({\mathcal{M}}^{K^{\prime }}+{\gamma }_v)\mathcal{G}{\mathbf{j}}_b^K\cdot \hat{\mathbf{x}}+({\mathcal{M}}^K+{\gamma }_v)\mathcal{G}{\mathbf{j}}_b^{K^{\prime }}\cdot \hat{\mathbf{x}}}{{\mathcal{M}}^K{\mathcal{M}}^{K^{\prime }}-{\gamma }_v^2},$$

$$\mathcal{G}{\mathbf{j}}_b^K\cdot \hat{\mathbf{x}}$$

$$\mathcal{G}{\mathbf{j}}_b^{K^{\prime }}\cdot \hat{\mathbf{x}}$$

$$\delta{\widetilde{\rho}}_{q,e}=\frac{\sigma_{xx}\left(\partial_{x}\phi^{>}\big{|}_{0^{+}}-\partial_{x}\phi^{<}\big{|}_{0^{-}}\right)2(i\omega-2\gamma_{v})+4\sigma_{H}v_{0}q^{2}\phi_{0}}{\omega^{2}+2i\omega\gamma_{v}-(v_{0}q)^{2}},$$

$$D=N{e}^2\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{v}^2(k)\left(-\frac{\partial f(ϵ)}{\partial ϵ}\right),$$

$$D=\frac{2N{e}^2}{\pi {\hslash }^2}{\int }_{\Delta }^{\infty }d{ϵ}_{+}\frac{{ϵ}_{+}^2-{\Delta }^2}{{ϵ}_{+}}\left(-\frac{\partial f({ϵ}_{+})}{\partial {ϵ}_{+}}\right),$$

$$D(\theta )=\frac{2N{e}^2\Delta }{\pi {\hslash }^2}\mathcal{F}(\theta ),\mathcal{F}(\theta )={\int }_1^{\infty }d\xi \left(1+\frac{1}{{\xi }^2}\right)\frac{1}{1+e^{\xi /\theta }},$$

$$k_0\left[{\tilde{\omega }}^2+2i\tilde{\omega }{\gamma }_v-(v_{0}q{)}^2+ϵ|q|(i\tilde{\omega }-2{\gamma }_v){\sigma }_{xx}\right]=ϵ{\sigma }_H{v}_0|q{|}^3,$$

$$f(\tilde{\omega })=\frac{ϵ}{\sqrt{2}}{\sigma }_H{v}_0{q}^2+(v_{0}q{)}^2,$$

$$f(\tilde{\omega })={\tilde{\omega }}^2+2i\tilde{\omega }{\gamma }_v+ϵ|q|(i\tilde{\omega }-2{\gamma }_v)\frac{D}{{\gamma }_{\mathrm{tr}}-i\tilde{\omega }}.$$

$$\frac{-2\omega }{{\tau }_p}+2{\gamma }_v\omega +ϵD|q|\omega \frac{({\gamma }_{\mathrm{tr}}-{\gamma }_v)}{({\gamma }_{\mathrm{tr}}-{\tau }_p^{-1}{)}^2+{\omega }^2}=0.$$

Full text

Edge states are a hallmark of the peculiar twisting of crystal wavefunctions in topological materials kane05 ; bernevig06 ; konig07 , and host a fermiology that departs from that of its parent bulk murakami06 ; bernevig06a ; hughes08 . Domain wall edge states (DWS) in gapped bilayer graphene are a particularly interesting example. Arising when the sign of the local gap in gapped bilayer graphene flips in real space martin08 ; jeiljung11 ; fanzhang13 ; eunah13 ; fengwang14 ; helin16 ; junzhu16 , DWS manifest in a number of different settings, e.g., at stacking faults (AB- and BA-) fanzhang13 ; eunah13 ; fengwang14 ; helin16 or in a split dual-gate geometry wherein perpendicular applied electric field in adjacent regions have opposite signs martin08 ; fanzhang13 ; junzhu16 . Domain walls host counter-propagating one dimensional (1D) edge states (DWS) living in separate $K$ and $K^{\prime}$ valleys martin08 ; fanzhang13 ; eunah13 , with valley filtered currents that are robust to disorder jeiljung11 . In contrast to helical edge states in intrinsic topological insulators kane05 ; bernevig06 ; konig07 , DWS in gapped bilayer graphene enjoy large and tunable bulk gaps up to $200\,{\rm meV}$ castro07 allowing their unusual behavior to manifest even at room temperature.

Here we show that the collective motion of carriers in DWS manifest unusual plasmon modes — domain wall edge plasmons (DWPs) — whose characteristics are distinct from conventional bulk plasmons (Fig. 1). Arising from collective charge density oscillations of carriers in the domain wall edge states (Fig. 1a), DWPs can exist even at zero bulk charge density (no doping) with a tunable frequency from the terahertz up to the mid-infrared ($\sim 200$ meV) (Fig. 1b,c) and disperse linearly in contrast to that expected from conventional 2D bulk plasmons.

Importantly, DWPs are long-lived and possess an insensitivity to bulk long-range disorder. While conventional plasmon lifetimes are limited by bulk transport scattering soljacic09 ; sorger12 ; principi13b , DWPs at low temperature transcend the restrictions of bulk transport scattering exhibiting DWP lifetimes orders of magnitude larger than the bulk transport scattering time (Fig. 2). As we argue below, these long lifetimes persist to high temperatures and can reach values of a few picoseconds at room temperature (for corresponding bulk transport scattering times of $\sim 0.5\,{\rm ps}$).

The topological edge states that host DWP are intimately locked to the difference of valley Chern number on either side of the domain wall martin08 ; fanzhang13 ; eunah13 ; DWPs possess valley polarization with backward/forward modes predominantly propagating in K/K’ valleys (Fig. 3). As we explain below, in addition to currents in the domain walls, DWP propagation also induces bulk undergap valley current flow, which renormalize the frequency of collective oscillations in the domain wall states. Control of the latter (e.g, via screening from a dielectric background) grants an unconventional knob to tune a myriad of DWP characteristics that range from its velocity and confinement, to the degree of DWP valley polarization.

We expect DWPs to manifest in experimentally available gapped bilayer graphene systems fanzhang13 ; eunah13 ; martin08 ; jeiljung11 ; fengwang14 ; helin16 ; junzhu16 such as along AB/BA stacking faults in globally gapped bilayer graphene, as well as electrostatically defined domain walls in split-dual-gate geometries. Indeed, both these methods have been recently employed to study topological domain walls experimentally fengwang14 ; helin16 ; junzhu16 . DWPs also feature subwavelength confinement of light, and can be probed by a variety of techniques that include gratings, and scanning near-field optical microscopy fei12 ; koppens12 .

Domain wall states and collective dynamics — We begin by considering domain walls in gapped bilayer graphene. These domain walls can be created in a number of ways, for e.g., (i) defined electrostatically where split-dual gates in bilayer graphene are biased to yield adjacent regions with layer potential of opposite signs martin08 ; fanzhang13 ; junzhu16 , and (ii) at AB-BA stacking faults where the bilayer graphene is globally gapped fanzhang13 ; eunah13 ; fengwang14 ; helin16 .

We account for both these types of domain walls phenomenologically by describing gapped bilayer graphene with a spatially varying band gap: ${\widetilde{\Delta}}(x)=\pm 2\Delta$ on either side of $x=0$. Reversing its sign at $x=0$, the domain walls at the zero node of ${\widetilde{\Delta}}(x)$ host DWS (Fig. 1a). We note, parenthetically, that the qualitative form of DWS is insensitive to the specific ${\widetilde{\Delta}}(x)$ profile used since DWS, arising from band inversion, is locked to its zero nodes. For electric field defined domain walls, ${\widetilde{\Delta}}(x)$ directly correlates with the layer potential difference. For domain walls at stacking faults, however, the physical band gap (layer potential difference) does not flip in real space. Instead, the chirality (in each valley) in AB and BA stacking regions are opposite, leading to opposite signs of valley specific Berry curvature and Chern number martin08 ; eunah13 ; fanzhang13 . We absorb this (chirality) sign into an effective ${\widetilde{\Delta}}(x)$. DWS are valley-helical states located inside the bulk band gap with backward (forward) moving DWS locked to the valley index $K(K^{\prime})$ (Fig. 1a) martin08 ; eunah13 ; fanzhang13 . For each valley, there are four edge states (DWS) with the same helicity propagating along $\hat{{\bf y}}$ stemming from layer and spin degrees of freedom martin08 ; fanzhang13 ; eunah13 ; junzhu16 .

In order to describe the dynamics of carriers in DWS, it is useful to separate out the density into bulk, $\rho_{b}({\bf r})$ (in $x<0$ and $x>0$ regions) as well as edge state density $\rho_{e}({\bf r})$ (situated at $x=0$) via

[TABLE]

where $\Theta(x)$ is the Heaviside function, and $\nu=\pm 1$ denote $K(K^{\prime})$ valleys. We note that the edge current ${\bf j}_{e}^{\nu}$ in each of the valleys arises from the chirality of the edge states: ${\bf j}_{e}^{\nu}({\bf r})=-\nu v_{0}\rho_{e}^{\nu}({\bf r})\hat{{\bf y}}$ where the edge states in valley $K$ and $K^{\prime}$ possess effective chiral velocity $-v_{0}\hat{{\bf y}}$ and $v_{0}\hat{{\bf y}}$, respectively supplement .

Bulk charge density evolves dynamically as

[TABLE]

where $-\nabla\phi({\bf r},t)$ is the electric field, and $\boldsymbol{\sigma}^{\nu}$ is the bulk conductivity tensor. $\boldsymbol{\sigma}^{\nu}$ contains both diagonal, $\sigma_{xx}$, as well as off-diagonal components, $\sigma_{xy}^{\nu}$. In gapped bilayer graphene, the latter arises from valley Hall currents xiao07 ; tarucha15 ; yuanbozhang15 and as we will see below, plays an integral role in DWP dynamics. Valley dependent Hall motion is characterized by the sign of the gap as well as the valley index; here we model $\sigma_{xy}^{\nu}(x)=\nu\ {\rm sign(x)}\sigma_{H},\ \sigma_{H}=4e^{2}/h$ where the factor $4$ corresponds with the number of DWS in each valley 111Here we have noted that there are two branches of DWS per valley/spin fanzhang13 ..

Similarly, the dynamics of the edge charge density can be discerned by applying the continuity relation to Eq. (1) and matching $\delta$-functions. We obtain

[TABLE]

where {\mathcal{G}}{\bf j}_{b}^{\nu}={\bf j}_{b,>}^{\nu}\big{|}_{0^{+}}-{\bf j}_{b,<}^{\nu}\big{|}_{0^{-}} and we have used $\partial_{x}\Theta(\pm x)=\pm\delta(x)$. While the second term describes dynamics arising from edge current flow within the DWS, the third term arises from bulk currents impinging into the DWS. The latter contribution include both valley Hall $\sigma_{H}$ as well as longitudinal $\sigma_{xx}$ currents. Valley relaxation is accounted for via a phenomenological inter-valley scattering rate $\gamma_{v}$.

Collective modes of the domain wall states emerge as self-sustained density oscillations of Eq. (1-3), and electric potential obeying

[TABLE]

in the non-retarded limit. Here $U({\bf r},{\bf r}^{\prime})$ is the Coulomb interaction, and $\delta\rho({\bf r},t)=\rho({\bf r},t)-\rho^{(0)}$ where $\rho^{(0)}$ is the equilibrium charge density. Since the system is translationally invariant along the edge ($y$ direction), DWPs propagate as waves of form $\phi({\bf r},t)={\widetilde{\phi}}_{q}(x,z)e^{i(qy-\omega t)}$ and $\delta\rho({\bf r},t)=\delta{\widetilde{\rho}}_{q}(x)\delta(z)e^{i(qy-\omega t)}$. Hereafter, we concentrate on the fields $\phi,\delta\rho$ at $z=0$.

In what follows, we will describe collective modes along the domain wall compactly in terms of $\phi$, by eliminating $\delta\rho$ from the dynamical equations. To do so, we first note that charge density localized on the domain wall, $\delta\rho_{q,e}$, produces a jump in the electric field as

[TABLE]

where $U_{q}(x)=\int dk\ e^{ikx}(q^{2}+k^{2})^{-1/2}/\kappa$ is the effective one dimensional (1D) Coulomb kernel. In obtaining Eq. (5) we have taken the derivative of Eq. (4), using the plane-wave forms of $\delta\rho,\phi$ and Eq. (1) above. Importantly, $\delta{\rho}_{e}=\delta\rho_{e}^{K}+\delta\rho_{e}^{K^{\prime}}$ in Eq. (5) can be directly related to the electric potential by inverting Eq. (3):

[TABLE]

where $\mathcal{M}^{\nu}=\partial_{t}+\gamma_{v}-\nu v_{0}\partial_{y}$ is an operator that acts on $\phi({\bf r},t)$ and $\nu=\pm 1$ for $K$ ($K^{\prime}$) valley; note that ${\bf j}_{b}^{K,K^{\prime}}$ depends on $\phi$ directly through Eq. (2).

In addition to continuity of ${\widetilde{\phi}}_{q}(x)$ and jump in electric field discussed above, electric potential of the plasmon, $\phi({\bf r},t)$, also satisfies Eq. (4); this yields $\phi({\bf r},t)$ as a solution to a non-local integro-differential problem. Instead, here we adopt a simplified Coulomb kernel ${\widetilde{U}}_{q}(x)=~{}\frac{1}{\kappa}\int dk\ 2qe^{ikx}/(2q^{2}+k^{2})$ 222$U_{q}(x)$ yields highly non-local integro-differential equation Eq. (4). We note there are other methods to (numerically) analyze integro-differential problem, for example, by the Wiener-Hopf method volkov88 ; katsnelson16 or multipole expansion xia94 ; kinaret11 . To illustrate the essential features of DWPs, here we adopt a simplified Coulomb kernel ${\widetilde{U}}_{q}(x)$ whose Fourier transform matches that of $U_{q}(x)$ up to leading order in $k/q$. which captures the essential long wavelength features of $U_{q}(x)$ fetter85 ; zabolotnykh16 . Using simplified ${\widetilde{U}}_{q}(x)$, we find ${\widetilde{\phi}}_{q}(x)$ follows

[TABLE]

Since Eq. (7) is local, $\phi_{q}(x)$ profile can be obtained in a straight-forward fashion as described below.

We first discuss the dispersive features of DWPs, focussing on the case $\gamma_{v}=0$ and Fermi energy inside the gap and $T=0$ so that no bulk carriers are excited; see below for a detailed discussion of the role of $\gamma_{v}$ and $\sigma_{xx}$. This yields $\sigma_{xx}=0$, $\delta\rho_{b}=0$ in the bulk, and a solution of Eq. (7) as ${\widetilde{\phi}}_{q}(x)=\phi_{0}e^{-\sqrt{2}|qx|}$. Plugging this ${\widetilde{\phi}}_{q}(x)$ profile into Eq. (5) and (6), we obtain the DWP dispersion (Fig. 1c):

[TABLE]

The first term inside the square root comes from the velocity of edge state carriers, whereas $\eta$ captures collective bulk valley Hall motion that moves along the DWP. We note that for $\hbar\omega\geq 2\Delta$, DWP enters the single particle continuum (SPC) (shaded region Fig. 1c) where particle-hole excitations damp the plasmon and destroy its coherence. When $q/q_{0}<1$, the SPC boundary, delineated by $\hbar\omega=2\Delta$, is nearly constant.

Strikingly, bulk valley Hall currents renormalize the collective mode velocity of DWS in Eq. (8). Estimating $v_{0}$ from Ref. martin08 , we obtain $v_{0}=4v_{F}(\sqrt{2}\Delta/t_{1})^{1/2}/3$ at zero Fermi energy, where $v_{F}=10^{6}\ {\rm m/s}$ is the monolayer Fermi velocity, $t_{1}=0.3\ {\rm eV}$ is the interlayer hopping parameter. Choosing $\Delta=0.1$ eV, $\sigma_{H}=1.5\times 10^{6}\ {\rm m/s}$ and $\kappa=1$, we estimate that the valley Hall contribution can be 27 times larger than the single particle edge state contribution (see orange bar vs purple bar in Fig. 1c). As a result, DWP group velocity can be five times larger than $v_{0}$.

Low plasmon velocities yield tight confinement of light when the plasmon is hybridized to form plasmon-polaritons. Indeed, taking $\Delta=0.1\ \rm eV$ we find a plasmon confinement of about 60 times smaller than free-space wavelength. For example, for $\hbar\omega=0.1\ \rm eV$ [below the single particle continuum (SPC) where $\hbar\omega\geq 2\Delta$, shaded area, right panel of Fig. 1c], this gives a confinement length as small as $200\ \rm nm$ (c.f. free-space wavelength for the same frequency of $12\ \mu\rm m$). Importantly, since $\eta$ depends strongly on background $\kappa$, screening can dramatically reduce $\eta$ and DWP velocity, further enhancing the confinement of DWP (dashed line in Fig. 1c; here $\kappa=20$). For very large $\kappa$ and small $\Delta=10\ \rm meV$, velocities dramatically slow down, giving a confinement that can be squeezed up to 3 orders of magnitude shorter than the free-space wavelength.

**Domain wall plasmon lifetime – ** The dynamics of (thermally activated) bulk charge as well as inter-valley scattering can contribute to the decay and damping of DWP. Employing Eq. (2) we find bulk charge dynamics: $-i{\widetilde{\omega}}\ \delta{\widetilde{\rho}}_{q,b}+\sigma_{xx}\left(-\partial_{x}^{2}+q^{2}\right){\widetilde{\phi}}_{q}=0$, where we model the bulk conductivity via a Drude model: $\sigma_{xx}=D(\theta)/(\gamma_{\rm tr}-i{\widetilde{\omega}})$, $\gamma_{\rm tr}=1/\tau_{\rm tr}$ is the carrier scattering rate and $D(\theta)$ is the Drude weight supplement that depends on $\theta=k_{B}T/\Delta$. Here we have used complex ${\widetilde{\omega}}$ to capture both plasmon oscillations (${\rm Re}({\widetilde{\omega}})$ denotes the plasmon frequency) as well as decay dynamics (${\rm Im}({\widetilde{\omega}})$ denotes its inverse lifetime). Using these, ${\widetilde{\phi}}_{q}(x)$ in the bulk takes the form supplement

[TABLE]

where $\omega_{b}^{2}=2\pi D(\theta)|q|/\kappa$ is the bulk plasmon frequency. We will first treat the case $|{\widetilde{\omega}}|^{2}\gg\omega_{b}^{2}$.

Substituting $\delta{\widetilde{\rho}}_{q,e}$ from Eq. (6) into Eq. (5) we obtain a complex DWP ${\widetilde{\omega}}$ obeying:

[TABLE]

where $\epsilon=8\pi/\kappa$. The plasmon dispersion and lifetime can be discerned from Eq. (10) by writing ${\widetilde{\omega}}(q)=\omega(q)-i/\tau_{p}(q)$, where $\tau_{p}(q)$ is the DWP lifetime. Solving Eq. (10) numerically, we plot the plasmon frequency ($\omega$) and lifetime ($\tau_{p}$) respectively in Figs. 2a and b; in these, we have used parameters $\Delta=0.1\ \rm eV$, and $\kappa=1$ as well as (disorder-limited) transport scattering time $\tau_{\rm tr}=0.5\ \rm ps$ which corresponds to a relatively high mobility $50,000\ {\rm cm^{2}/Vs}$ that can be realized in hBN-encapsulated bilayer graphene samples yuanbozhang15 ; junzhu16 . Here we use a temperature independent $\tau_{\rm tr}$. We have also used the intervalley scattering lifetime $\tau_{v}=1/\gamma_{v}=100\ \rm ps$ as estimated Ref. jeiljung11 .

In Fig. 2a, we show DWP frequency $\omega$ at low temperature (black solid line) and room temperature (red solid line) in comparison with their bulk plasmon frequencies $\omega_{b}$ (dashed lines). We note that DWP dispersion remains largely linear and exhibits little difference between room temperature ($300\,{\rm K}$) vs low temperature ($50\,{\rm K}$) due to the slow increase of the Drude weight with temperature supplement . The bulk plasmon frequency $\omega_{b}$ is negligible at low temperature. However, $\omega_{b}$ becomes comparable to $\omega$ at small $q$ and large temperatures. When $\omega=\alpha\omega_{b}$, kinematics allow DWPs to rapidly decay into bulk plasmons, when $\alpha$ is of order unity. While a detailed analysis of DWP to bulk plasmon emission is beyond the scope of this work, we delineate this regime in Fig. 2b with regions $\omega\lesssim\alpha\omega_{b}$ shown in white. As an illustration, we set $\alpha=2.5$. See supplement for a detailed comparison of $\omega$ to $\omega_{b}$.

Importantly, as shown in Fig. 2(b,c), DWPs [obtained numerically from Eq. (10)] can exhibit very long lifetimes $\sim 1.5\ \rm ps$ even at room temperature exceeding reported plasmon lifetime ($\sim 0.5\ \rm ps$) in hBN-encapsulated graphene woessner15 . Strikingly, $\tau_{p}$ exceeds the bulk transport scattering time of $\tau_{\rm tr}=0.5\ \rm ps$ (dotted black line, Fig. 2c), and clearly demonstrates how DWP $\tau_{p}$ can transcend the conventional limit set by bulk transport scattering woessner15 ; principi13b . We note that expected phonon-limited mobility at room temperature for bilayer graphene can reach values of $200,000\ \rm cm^{2}/Vs$ morozov08 ; castro10 ; mobilities of $125,000\ \rm cm^{2}/Vs$ wees11 at room temperature have been reported in hBN-encapsulated graphene. With those values of mobility, DWP lifetime may reach $\sim 6\ \rm ps$ at room temperature.

Enhanced lifetimes arises due to a suppression of bulk carrier density that provides a pathway for DWPs to decay. To illustrate this, we estimate DWP lifetime $\tau_{p}$ from Eq. (10) by taking the limit $\omega_{b}\ll\omega$ and $\tau_{v,p,{\rm tr}}^{-1}\ll\omega$. In this limit, $\tau_{p}$ takes on the simple form supplement :

[TABLE]

In obtaining Eq. (11) we have additionally assumed $\tau_{v}^{-1}\ll\tau_{\rm tr}^{-1}$. In Fig. 2c, we show $\tau_{p}$ as a function of temperature for a fixed $\hbar\omega=0.1\ \rm eV$ obtained from Eq. (11) (solid lines) and compare with numerical results (dashed lines) showing excellent agreement. Crucially, Eq. (11) shows explicitly how low bulk carrier density (encoded in the bulk plasmon frequency, $\omega_{b}$) quenches the role of bulk transport scattering in DWP lifetime. Indeed, neglecting $\tau_{v}^{-1}$ and for $\omega>\omega_{b}$, $\tau_{p}$ is enhanced by a factor of $\sim\omega^{2}/2\omega_{b}^{2}$ over $\tau_{\rm tr}$; in this regime, $\tau_{p}$ scales (approximately) linearly with $\tau_{\rm tr}$. Interestingly, the dependence of $\omega_{b}$ and $\omega$ on $\kappa$ in Eq. (11) indicate that DWP lifetime at room temperature can be further boosted by screening as shown in Fig. 2c [$\kappa=1$ (red) and $\kappa=20$ (blue)].

$\tau_{p}$ exhibits a distinct temperature dependence (see Eq. (11), Fig. 2c). At high temperature, since bulk Drude weight is thermally activated, $\tau_{p}$ similarly displays an exponential temperature dependence (Fig. 2c) sharply increasing as temperature drops. However, at low temperature $\omega_{b}$ vanishes (black dashed lines in Fig. 2a, see also Eq. (11)). As a result, intervalley scattering dominates DWP lifetime cutting the exponential rise of DWP lifetime $\tau_{p}\to\tau_{v}$ (see Fig. 2c).

Due to the valley-helical nature of DWS, $\tau_{v}$ can in-principle be very large. Indeed, Ref. jeiljung11 reported that 1D channel is insensitive to backscattering and long range disorder giving a mean free path of about $100\ \mu\rm m$ corresponding to $\tau_{v}$ as high as $100\ \rm ps$. We note that recent transport experiments along both electric field and stacking fault domain walls report shorter $\tau_{v}$ of about a few hundred $\ \rm fs$ fengwang14 ; junzhu16 . Shorter $\tau_{v}$ in electric field domain walls may arise from short-ranged disorder that can scatter between valleys such as grain boundaries fengwang14 , as well as wide electrostatic profile used to create ${\widetilde{\Delta}}(x)$ in electric field defined domain walls. In the latter, the electrostatic profile is characterized by a finite effective width $L_{0}$ of the domain wall junzhu16 . Although the electronic structure of DWS are relatively independent of $L_{0}$ – domain walls arise whenever the ${\widetilde{\Delta}}(x)$ flips sign – broad $L_{0}$ allow additional non-chiral (non-topological) states that can mediate scattering between DWS in separate valleys and consequently reduce $\tau_{v}$ junzhu16 . We note that the typical width used by Ref. junzhu16 was about $L_{0}\approx 100\ \rm nm$. With reduced $L_{0}$ and smooth potential profile, $\tau_{v}$ may reach long ballistic timescales characteristic of topological edge states.

Valley polarization – Single-particle carrier transport within DWS are completely filtered by valley index fanzhang13 : at $K(K^{\prime})$ valley, carriers in DWS propagate in the $-y$ ($+y$) direction. In contrast, the collective modes of DWP experience a mixture of both valley contributions since Coulomb interactions are long-ranged and do not discriminate between valleys. In order to quantify how much each valley contributes to the collective motion of DWP, we analyze Eq. (3). For simplicity, we specialize to the limit $\gamma_{\rm tr}=0,\gamma_{v}=0$. Using $\phi({\bf r},t)$ profiles from Eqs. (7) and (9), we obtain an oscillating charge density in each valley as

[TABLE]

exhibiting finite amplitudes of charge density for both valleys in both directions.

The distinction between valley contributions for opposing directions are particularly clear for very small $T$, where $D(\theta)=0$. In this limit, we find that valley polarization ${\widetilde{\rho}}_{e}^{\nu}(y,t)=\delta\rho_{e}^{\nu}/\delta\rho_{e}$ are $\pi$ out-of-phase with each other (i.e. for every $(q,y,t)$, ${\widetilde{\rho}}_{e}^{\nu}$ have opposite signs) and have different amplitudes (Fig. 3a,b); departing from perfect valley polarization regime (dashed lines of Fig. 3a). Non-zero amplitude ${\widetilde{\rho}}_{e}^{\nu}$ in both valleys for DWPs (and partial valley polarization) is a direct consequence of collective motion of bulk valley Hall currents. Indeed, DWP frequency renormalization in Eq. (8) originates from mixing of the two valleys. This contrasts with the $\eta=0$ case in Eq. (8) where DWPs traveling along $q<0$ ($q>0$) are fully $K$ ($K^{\prime}$) valley polarized. In Fig 3a, we have set $\kappa=20$ and obtain about $80\%$ vs $20\%$ mixture of $|{\widetilde{\rho}}_{e}^{\nu}|$. We note that at smaller $\kappa$, ${\widetilde{\rho}}_{e}^{\nu}$ will deviate even further from perfect filtering (dashed lines of Fig. 3a) due to a stronger Coulomb potential.

DWPs are long-lived and possess decay times that surpass conventional plasmon decay restrictions wherein plasmon lifetime is limited by the bulk’s transport scattering time. This property is unusual and stems from the distinct origin of DWPs: collective oscillations of carriers in the edge states. Indeed, the edge states enable large quality factors for DWP oscillations which can range from about $10^{2}$ (at room temperature) up to $10^{4}$ (at low temperatures) supplement . Surprisingly, DWPs’ long lifetime and high quality manifest without sacrificing subwavelength electro-magnetic confinement. A tantalizing prospect for utilizing DWPs are gate-defining topological domain walls in gapped bilayer graphene for DWP plasmonic waveguides. Together with high quality and long-lived DWPs, gate-defined domain walls provide a means for patterning low-dissipation (and valley polarized) plasmonic circuits ebbesen08 ; sorger12 .

Acknowledgements — We are grateful for useful conversations with Frank Koppens, Niels Hesp, and Mike Schecter. This work was supported by the Singapore National Research Foundation (NRF) under NRF fellowship award NRF-NRFF2016-05.

I Supplementary Information for “Long-lived domain wall plasmons in gapped bilayer graphene”

Sign of edge current

Edge current ${\bf j}_{e}^{\nu}({\bf r},t)$ in Eq. (1) of the main text arises from single particle motion in the domain wall edge states (DWS). We determine the direction ${\bf j}_{e}^{\nu}({\bf r},t)$ through bulk-edge correspondence: valley-helical edge ${\bf j}_{e}^{\nu}({\bf r},t)$ propagates in the same direction as the bulk undergap valley Hall current ${\bf j}_{b}^{\nu}({\bf r},t)$ close to the edge.

For clarity, we focus on electric field domain walls where ${\widetilde{\Delta}}(x)$ reflect the layer potential difference. We first note that near an edge, a confining potential $\mathcal{V}(x)$ creates electric field $-\partial_{x}\mathcal{V}\hat{{\bf x}}$ and the Hall current ${\bf j}_{b}^{\nu}({\bf r},t)=-\partial_{x}\mathcal{V}\hat{{\bf x}}\times\sigma_{xy}^{\nu}\hat{{\bf z}}$. Similarly, close to a domain wall where ${\widetilde{\Delta}}(x)$ flips sign, a layer dependent potential $\mathcal{V}_{\rm b,t}(x)$ exhibits a profile near the domain wall edge (see solid lines of Fig. 4), acting on carriers to produce a bulk $j_{b}^{\nu}({\bf r},t)$ in the valence band; for $k_{B}T\ll\Delta$ valence band carriers dominate the anomalous Hall current in each of the valleys. Note that since ${\widetilde{\Delta}}(x)$ flips sign, valence band carriers reside in different layers on either side of $x=0$ martin08 ; fanzhang13 and experience different $\mathcal{V}_{\rm b,t}$ potential profiles. In $x<0$ ($x>0$) region, $-\partial_{x}{\mathcal{V}}_{\rm b}\hat{{\bf x}}$ ($-\partial_{x}{\mathcal{V}}_{\rm t}\hat{{\bf x}}$) is pointing along $-x$ ($+x$) as shown in red (blue) arrow. Noting $\sigma_{xy}^{\nu}(x)=\nu\ {\rm sign(x)}\sigma_{H}$, ${\bf j}_{b}^{\nu}$ in both regions are directed in the $-\nu\hat{{\bf y}}$ direction where $\nu=\pm 1$ for $K$ or $K^{\prime}$, respectively. Matching the directions of ${\bf j}_{b}^{\nu}$ and ${\bf j}_{e}^{\nu}$, we write ${\bf j}_{e}^{\nu}=-\nu v_{0}\rho_{e}^{\nu}\hat{{\bf y}}$. We expect a similar reasoning also applies for domain walls at stacking faults where the opposite chirality in AB and BA stacking regions flips the sign of effective ${\widetilde{\Delta}}(x)$.

Inverse lateral DWP length, $k_{0}$

The self-induced potential around $x=0$ due to DWPs takes the form ${\widetilde{\phi}}_{q}(x)=\phi_{0}e^{-k_{0}|x|}$. To obtain the plasmon inverse lateral $k_{0}$ length in Eq. (9) of the main text we analyze the dynamics of bulk charge density in Eq. (2) of the main text:

[TABLE]

together with the simplified Coulomb kernel in Eq. (7) of the main text for $x\neq 0$:

[TABLE]

Substituting $\delta{\widetilde{\rho}}_{q,b}$ in Eq. (S-2) with Eq. (S-1) and replacing $\partial_{x}\to-k_{0}$, we obtain:

[TABLE]

Recalling $\sigma_{xx}=D(\theta)/(\gamma_{\rm tr}-i\omega)$ and bulk plasmon frequency $\omega_{b}^{2}=2\pi D(\theta)|q|/\kappa$, we can rewrite Eq. (S-3) in terms of $\omega_{b}$:

[TABLE]

As a result, non-trivial solutions to Eq. (S-4) occur when

[TABLE]

where we have taken only the positive root to ensure that the potential profile in Eq. (9) stays finite for all $x$. Note that in the limit of zero bulk density, $\sigma_{xx}\to 0$, and Eq. (S-5) reduces to $k_{0}=\sqrt{2}|q|$.

Dynamics of edge charge density

Dynamics of charge carriers within the domain wall states along the topological domain walls, together with a self-induced electric potential, yield DWP collective modes. We note that the dynamics of domain wall charge density arise from a number of contributions that include the chiral flow of valley charges in each of the edge states, bulk undergap valley currents that impinge into the edge states, as well as intervalley scattering (which we include via a phenomenological intervalley scattering rate $\gamma_{v}=1/\tau_{v}$). These processes are captured by continuity equation of edge charge density as shown in Eq. (2) of the main text, reproduced here for convenience:

[TABLE]

where {\mathcal{G}}{\bf j}_{b}^{K}={\bf j}_{b,>}^{K}\big{|}_{0^{+}}-{\bf j}_{b,<}^{K}\big{|}_{0^{-}}. To obtain the charge dynamics in terms of the current dynamics, we invert Eq. (S-6), to find

[TABLE]

where the operators $\mathcal{M}^{K}=\partial_{t}+\gamma_{v}-v_{0}\partial_{y}$ and $\mathcal{M}^{K^{\prime}}=\partial_{t}+\gamma_{v}+v_{0}\partial_{y}$. Summing both contributions, the total edge charge density $\rho_{e}=\rho_{e}^{K}+\rho_{e}^{K^{\prime}}$ is

[TABLE]

as shown in Eq. (6) of the main text.

To analyze DWP, we will describe its motion compactly in terms of electric potential, $\phi$, by eliminating $\delta\rho$ from the dynamical equations. To do so, we first note that the current flow in each of the valleys is directly related to the electric potential via Ohm’s law [see Eq. (2) of the main text]. Writing this out explicitly gives

[TABLE]

where \phi_{0}=\phi\big{|}_{x=0} [see Eq. (9) of main text] and we have noted that the opposite signs of $\sigma_{xy}^{\nu}$ on either side of the domain wall add when $\mathcal{G}$ acts on ${\bf j}_{b}^{\nu}$. Substituting the plane wave form $\delta\rho_{e}({\bf y},t)=\delta{\widetilde{\rho}}_{q,e}e^{i(qy-\omega t)}$ into Eq. (S-8) and using Eq. (S-9), produces a direct relation between $\delta{\widetilde{\rho}}_{q,e}$ and $\phi$:

[TABLE]

Finally, we recall that $\phi({\bf r},t)$ satisfy boundary conditions at $x=0$: $\phi({\bf r},t)$ is continuous at $x=0$, $\partial_{x}\phi$ may exhibit a jump as in Eq. (5) of the main text. Applying the form of $\phi_{q}(x)$ in Eq. (9) of the main text to Eq. (S-10) and the boundary conditions above, we obtain the plasmon dispersion (for complex ${\widetilde{\omega}}$) shown in Eq. (10) of the main text. Note that for $\sigma_{xx},\gamma_{v}\to 0$, Eq. (10) of the main text reduces to Eq. (8) in the main text as expected.

Drude weight for gapped bilayer graphene

The Drude weight for gapped bilayer graphene can be obtained semiclassically via

[TABLE]

where $N=4$ accounts for spin and valley degeneracy, $e$ is electron’s charge, $v=\partial\epsilon({\bf k})/\partial(\hbar{\bf k})$ is electron’s group velocity, and $f(\epsilon)$ is the Fermi-Dirac distribution. We adopt a simple two-band model of gapped bilayer graphene $\epsilon_{\pm}(k)=\pm\Delta[1+(k/q_{0})^{4}]^{1/2}$ fanzhang13 ; martin08 , where $\Delta$ is half gap, $q_{0}=\sqrt{\Delta t_{1}}/\hbar v_{F}$, $t_{1}$ is the interlayer hopping and $v_{F}$ is the Fermi velocity of monolayer graphene. Using the form of $\epsilon_{\pm}(k)$ above, we change integration variables in Eq. (S-11) from $k$ to $\epsilon_{+}$ yielding

[TABLE]

where the factor of 2 accounts for equal contributions of electrons in the conduction band and holes in the valence band. Integrating by parts, recalling $f(\infty)=0$, and making the integrand dimensionless, $\xi=\epsilon_{+}/\Delta$, we obtain

[TABLE]

where $\theta=k_{B}T/\Delta$, $k_{B}$ is the Boltzmann constant and we have written $f(\epsilon_{+})$ explicitly. For all plots in the main text and the supplement, Eq. (S-13) was integrated numerically.

DWP lifetime and $\omega/\omega_{b}$ ratio

We can estimate DWP lifetime $\tau_{p}$ from Eq. (10) in the main text, reproduced here for convenience:

[TABLE]

where $\epsilon=8\pi/\kappa$, $k_{0}=\sqrt{2}|q|[({\widetilde{\omega}}^{2}-\omega_{b}+i\gamma_{\rm tr}{\widetilde{\omega}})/({\widetilde{\omega}}^{2}-2\omega_{b}+i\gamma_{\rm tr}{\widetilde{\omega}})]^{1/2}$, $\omega_{b}^{2}=2\pi D|q|/\kappa$ and $\sigma_{xx}=D/(\gamma_{\rm tr}-i{\widetilde{\omega}})$. Specializing to the case $\omega_{b}\ll\omega$, yields $k_{0}=\sqrt{2}|q|$. Rearranging Eq. (S-14), we isolate terms containing ${\widetilde{\omega}}$ to the left hand side (LHS) so that Eq. (S-14) reads as

[TABLE]

where

[TABLE]

It is useful to note that Eq. (S-16) contains ${\widetilde{\omega}}$ and is a function of complex values, while the RHS is purely real. As a result, ${\rm Im}\,f({\widetilde{\omega}})=0$. Writing ${\widetilde{\omega}}=\omega-i/\tau_{p}$, this condition can be expressed as

[TABLE]

Taking the limits, $\gamma_{\rm tr},\tau_{p}^{-1}\ll\omega$ and $\gamma_{\rm tr}\gg\gamma_{v}$, produces a simple relation for the plasmon lifetime as shown in Eq. (11) of the main text, reproduced here for convenience

[TABLE]

where we have used relations $\gamma_{\rm tr}=\tau_{\rm tr}^{-1}$ and $\gamma_{v}=\tau_{v}^{-1}$.

As shown in Eq. (S-18), the small ratio between bulk plasmon frequency $\omega_{b}$ and DWP frequency $\omega$ suppresses the role of $\tau_{\rm tr}$ on DWP lifetime. In Fig. 5a, we compare DWP frequency to the frequency of bulk plasmons $\omega/\omega_{b}$ as a function of $q$ and $T$ using Eq. (10) of the main text. Strikingly, at low temperature, DWP frequency is about six orders of magnitude higher than $\omega_{b}$ since there are very few carriers in the bulk, resulting a very soft bulk plasmon mode. Since $\omega_{b}\ll\omega$ at low temperatures, decay dynamics from bulk transport scattering is completely quenched. Instead, damping from intervalley scattering dominates, and $\tau_{p}\to\tau_{v}$ [see Eq. (S-18) and Fig. 2c in the main text].

As temperature increases, the Drude weight in Eq. (S-13) also increases making the bulk plasmon stiffer. As a result, the difference between $\omega$ and $\omega_{b}$ shrinks and the ratio of $\omega/\omega_{b}$ drops. Nevertheless, even at room temperature, DWP can still exceed the bulk plasmon frequency by several times (see red solid and dashed lines of Fig. 2a in the main text). At high temperature, scattering from the bulk dominates the $\tau_{p}$ because $\gamma_{\rm tr}$ is about two orders of magnitude larger than $\gamma_{v}$. In Fig. 2c of main text, we displayed that an increase of $\kappa$ prolongs the DWP lifetime when $\omega$ is fixed. The dependence of lifetime on $\kappa$ can be discerned from Eq. (S-18) by approximating $|q|\approx\omega/v_{0}\sqrt{1+\eta}$. This approximation was obtained from the wavevector of $T=0$ dispersion for DWP in Eq. (8) of the main text (we have assumed that $\omega>\omega_{b}$). Plugging this estimate into $\omega_{b}$ of Eq. (S-18) yields a scaling

[TABLE]

where we have neglected effect of intervalley scattering at room temperature.

Recalling $\eta=4\sqrt{2}\pi\sigma_{H}/v_{0}\kappa$ and assuming $\eta\gg 1$, we obtain $\tau_{p}\propto\sqrt{\kappa}$. Indeed as shown in Fig. 2c in the main text, $\tau_{p}$ increases by about $\sqrt{20}$ from $\tau_{p}=1.5\ \rm ps$ at $\kappa=1$ to become $\tau_{p}=6.5\ \rm ps$ at $\kappa=20$ at $300$ K.

Quality factor

The quality factor or inverse loss function, $Q={\rm Re}\,{\widetilde{\omega}}/{\rm Im}\,{\widetilde{\omega}}=\omega\tau_{p}$, is a dimensionless quantity describing the number of plasmon oscillations performed before decay. We plot $Q$ for DWPs in Fig. 5b using a numerical solution of Eq. (10) in the main text showing large Q factors. At room temperature, $Q$ for DWPs can be several hundreds and increases exponentially to about $10^{4}$ at low temperatures (Fig. 5b). This size of $Q$ is large and particularly arresting when compared to conventional bulk plasmon $Q$ factors in graphene that have been experimentally observed ($\sim 20$) woessner15 and theoretically predicted ($\sim 100$) soljacic09 .

For ordinary bulk plasmon, dielectrics tend to reduce the $Q$ factor as they introduce additional scattering pathways for plasmon damping. However, $Q$ of DWP is surprisingly enhanced by dielectric background at room temperature and a fixed $\omega$ as $Q\propto\sqrt{\kappa}$ owing to Eq. (S-19). The enhanced $\kappa$ reduces the bulk contribution to DWP which consequently enhances $Q$ factor.