3.1.

Metagate-Induced Plasmonic Crystal

The dispersion relation $\mathbf{q}(\omega )$ of the GPPs propagating in homogeneous graphene with a uniform chemical potential $E_F$ surrounded by an embedding medium with a dielectric constant $\epsilon$ is given by

Equation (1) implies that a periodic modulation in $E_F$ generates a periodic texture in the effective refractive index for the propagating GPPs on graphene, thereby creating a photonic crystal. A wide range of periodic modulation of $E_F$ can be achieved using a metagate;³¹ see Fig. 1(a) for the metagate geometry employed in the rest of this work. The metagate is a thin slab of metal (perfect electric conductor or PEC) perforated by a periodic array of unit cells containing a circular hole of radius $r_1$ at its center and two identical holes of radius $r_2$ at the periphery. When graphene is capacitively doped by this structured gate, the area above the empty holes attracts a smaller carrier density than the area above the unperforated regions. A periodic landscape of $E_F$ produced by metagate-induced doping is shown in Fig. 1(b) (a specific case of $r_1=r_2$ is depicted). With an eye on future experimental realizations, our electrostatic simulation assumes that graphene is hBN-encapsulated and accounts for graphene’s quantum capacitance self-consistently under electrochemical equilibrium in graphene.^31,46

For the $r_1=r_2$ case shown in Fig. 1(b), the metagate is patterned as a honeycomb lattice with a smaller unit cell consisting of a single hole, as shown by the dashed line in Fig. 1(b). The photonic band structure (PBS) of the terahertz (THz)-frequency GPPs is shown in Fig. 1(d); it is calculated along the irreducible boundary of the Brillouin zone (BZ) defined from the smaller unit cell [dashed line in Fig. 1(c)]. As expected, the honeycomb photonic lattice produces a twofold degenerate conical dispersion at the two equivalent ($\mathbf{k}$ and ${\mathbf{k}}^{\prime }$) valleys. The presence of a metallic gate in close proximity of the graphene layer is captured by the linear dispersion relation of the so-called acoustic GPPs⁹ at small momenta $|\mathbf{q}|\ll a^{-1}$ (i.e., near the $\gamma$-point of the BZ). Also, it is previously reported that the screening through this periodically structured gate, along with $E_F$ landscape, contributes to stronger modulation of PBS.³¹ We also note the subwavelength nature of the metagate (and, correspondingly, of the imprinted photonic crystal): the unit cell periodicity $a=600\mathrm{nm}$ is much smaller than the free-space wavelength of the THz radiation. For example, at the $\mathbf{k}({\mathbf{k}}^{\prime })$ point, we obtain an ultrahigh confinement ratio ${\lambda }_{\text{vacuum}}/{\lambda }_{\mathrm{GPP}}\sim 140$.

The folded band structure in Fig. 1(e) is identical to the one in Fig. 1(d), except that it is drawn along the smaller BZ defined by the original unit cell for generic ($r_1\ne r_2$) cases. The double degeneracy of the Dirac cones at $\mathbf{\Gamma }$ originates from the trivial band folding from the “true” (minimal) first BZ to an artificially shrunken BZ [i.e., from the dashed hexagon to the solid hexagon in Fig. 1(c)]. The four-fold degenerate modes at $\mathbf{\Gamma }$ shown here are the fundamental building blocks^26,27 for the topological crystalline phases discussed below.

3.2.

Kekulé-Textured Honeycomb Lattice

To highlight the underlying physics, we provide a well-known tight-binding (TB) model for a Kekulé-textured graphene-like system.²⁶ Figure 2(a) illustrates the setting of the TB model and its mapping onto the metagate-tuned graphene structure shown in Fig. 1. The hexagonal unit cell of our metagate can be divided into six equivalent segments [see the shaded triangle in Fig. 2(a)] supported by ${\mathcal{C}}^6$ point group symmetry. If $r_1>r_2$ ($r_1<r_2$), the distance between the center-of-mass positions of the two adjacent metagate segments is shorter (greater) when they are positioned across the unit cell boundary versus within the same unit cell. Qualitatively, if we consider each segment to be an effective TB site, then the effective hopping strength within a unit cell $t_{\text{in}}$ is weaker (stronger) than that across adjacent unit cells $t_{\text{out}}$ in the case of $r_1>r_2$ ($r_1<r_2$).

Fig. 2

(a) Schematic of the TB model on a Kekulé lattice. Shaded triangular segment of a metagate maps onto a lattice site of the TB model. (b) PBS calculated for $t_{\text{in}}=1.1t/t_{\text{out}}=0.9t$ (left) and $t_{\text{in}}=0.9t/t_{\text{out}}=1.1t$ (right). Color: ratio between the projection amplitudes onto the dipole-like orbitals $|p_{\mathrm{1,2}}⟩$ and onto the quadrupole-like orbitals $|d_{\mathrm{1,2}}⟩$. (c) PBS calculated for a structure interfacing the two configurations from (b) along the zigzag-oriented domain wall. Bandgap-crossing purple lines: dispersion curves for the edge states propagating along the domain wall.

Any eigenstates on six hexagonal vertices can be described by linear combinations of six orbital states: a monopole-like ($|s⟩=[1;1;1;1;1;1]$), two dipole-like ($|p_1⟩=[1;1;0;-1;-1;0]$, $|p_2⟩=[1;-1;-2;-1;1;2]$), two quadrupole-like ($|d_1⟩=[1;-1;0;1;-1;0]$, $|d_2⟩=[1;1;-2;1;1;-2]$), and an octopole-like ($|f⟩=[1;-1;1;-1;1;-1]$).²⁶ The monopole and octopole modes are at high energies, playing negligible roles for low-energy physics (high/low compared to $t_{\text{in}}$ and $t_{\text{out}}$). When two hopping strengths are the same, the dipole and quadrupole modes are degenerate at zero Bloch momentum just like in Fig. 1(e). This is because, if $t_{\mathrm{in}}=t_{\mathrm{out}}$, $|p_{\mathrm{1,2}}⟩$ and $|d_{\mathrm{1,2}}⟩$ become equivalent with respect to the redefinition of the unit hexagon (around the center hole or around a periphery hole). This fourfold degeneracy is lifted when $t_{\text{in}}\ne t_{\text{out}}$, as shown in Fig. 2(b). Specifically, the dipole (quadrupole) modes can be shown to reside below (above) the bandgap for $t_{\text{out}}<t_{\text{in}}$. Band inversion occurs between these four orbitals around the $\mathbf{\Gamma }$-point for the $t_{\text{out}}>t_{\text{in}}$ case, thus creating an insulating phase that is topologically distinct from the “topologically trivial” insulating phase realized in the case of $t_{\text{out}}<t_{\text{in}}$.

In electronic systems, such topologically nontrivial insulators are sometimes referred to as quantum orbital-Hall insulators²⁶ because their underlying cause—band inversion of orbital-locked states—is reminiscent of that in quantum spin-Hall insulators. Indeed, the low-energy effective four-band Hamiltonian around the $\mathbf{\Gamma }$ point is characterized by a ${\mathbb{Z}}_2$ topological invariant^16,26,27 ($Z_2=0$ for $t_{\text{out}}<t_{\text{in}}$ and $Z_2=1$ for $t_{\text{out}}>t_{\text{in}}$) as in quantum spin-Hall insulators (see the Supplementary Material). Thus following the bulk-edge correspondence,^11–13 linear edge states should arise at the interface between the topological and trivial domains. Indeed, we observe two expected domain wall states in the bandgap, as shown in Fig. 2(c). The forward (backward) propagating mode is linked to the orbital states such that the relative phase of the lattice sites within a unit cell is increasing in the counter-clockwise (clockwise) direction.

3.3.

Topological Crystalline Phase in Graphene Plasmons

Next, we confirm that our metagate-tuned GPPs indeed exhibit the two topologically distinct phases introduced above. This is done by comparing the band structures and examining the plasmonic mode profiles above and below the bandgap. Figure 3(a) provides the GPP PBS calculated for two different cases: $r_1=90\mathrm{nm}$ and $r_2=150\mathrm{nm}$ (topologically trivial case), and the topologically nontrivial case of $r_1=160\mathrm{nm}$ and $r_2=110\mathrm{nm}$. Here $r_1$ and $r_2$ are chosen such that any adjacent holes are clearly separated by at least a 40-nm gap. Plasmonic bandgaps of similar sizes are opening for both configurations, and the overall shape of the bands agrees with the band structures of the TB model shown in Fig. 2(b). The dynamic charge density in graphene associated with the propagating GPPs is plotted in Fig. 3(b) for a single unit cell.

Fig. 3

(a) Plasmonic band structures calculated for $r_1=90\mathrm{nm}/r_2=150\mathrm{nm}$ (left: topologically trivial domain) and $r_1=160\mathrm{nm}/r_2=110\mathrm{nm}$ (right: topologically nontrivial domain). Line color: dipole-to-quadrupole ratios (see text). (b) Color: AC surface charge density in graphene associated with the GPPs at the $\mathbf{\Gamma }$-point (below the bandgaps). Parameters: same as in Fig. 1.

Next, we color-coded [Fig. 3(a)] each eigenstate according to the ratio between the dipole and quadrupole strengths. By comparing the left and right plots in Fig. 3(a), we note that the band inversion predicted by the TB model is indeed observed for the GPPs propagating through the modulated two-dimensional Fermi energy landscape $E_F(x,y)$ in graphene. As depicted in Fig. 3(b), the plasmonic modes near the $\mathbf{\Gamma }$ point display dipole (quadrupole) charge distributions below the bandgap when $r_1$ is smaller (greater) than $r_2$. Although the visualized dipole-quadrupole ratio does not show a perfect agreement between the plasmonic band structure and the TB model, this is due to the fact that the actual plasmonic modes are not strongly bounded as in the idealized TB model. To account for weak localization more precisely, the effective TB description would need to include the second (or higher) nearest neighbor couplings. These minor discrepancies, however, can be neglected since the simplest TB model can capture all the topological features of our system.

3.4.

Topologically Protected Plasmonic Edge States

To observe the topologically protected edge states propagating along the domain wall between the trivial and nontrivial topological domain, it is necessary to carry out additional engineering of the interface between the two domains. This is needed because the microscopic field profiles of the dipole and quadrupole modes are not identical inside the two domains owing to the change of the unit cell structure. Therefore, when the two distinct metagate domains are interfaced, these abrupt changes across the domain wall can introduce a small bandgap between two edge states.^26,27 To remove this gapping of the edge states, we engineer the domain wall by setting the radius $r_0$ of the peripheral metagate holes (i.e., the holes at the domain wall) to be the average of the two $r_2$ values from each domain: $\overline{r}=(r_{2e}+r_{2s})/2$. Here the subscripts $s$ and $e$ stand for the shrunken (topologically trivial) and expanded (topologically nontrivial) configuration, respectively. This additional interface engineering mitigates the local ${\mathcal{C}}^6$ symmetry breaking at the interface.

Next, we calculate the PBS of the GPPs propagating along the zigzag-shaped domain wall oriented along the $x$ direction. The computational domain is comprised of 18 unit cells on each side of the domain wall. The results are shown in Fig. 4(a), where the black thin lines correspond to the projected bulk modes (i.e., having the same value of $k_x$), whereas the red (blue) solid lines correspond to the forward (backward) propagating chiral edge states. The edge states have been found to be exponentially localized at the domain wall, as shown in Figs. 4(c) and 4(d), while the bulk modes are delocalized (not shown). The effectiveness of the interface engineering described above is apparent from the fact that the gap between the forward- and backward-propagating edge modes is extremely small, as can be seen in Fig. 4(a). The existence of an extremely small but finite gapping of the edge modes due to symmetry breaking at the domain wall has been noted earlier.²⁷

Fig. 4

(a) Projection of the PBS of the GPPs supported by a graphene nanoribbon gated by a metagate comprised of the two domains shown in Fig. 3. Each domain contains 18 unit cells, zigzag-type domain wall runs parallel to the $x$-direction, the radius of peripheral holes at the domain wall: $\overline{r}=130\mathrm{nm}$. Black thin lines: $x$-projected bulk bands; red (blue) lines and circles: forward (backward)-propagating chiral edge states. (b) Metagate-induced Fermi energy profile on graphene near the domain boundary (white dotted line). (c) AC plasma charge density profiles for the forward- and backward-propagating edge states. 2-D color code: the magnitude and phase of the charge density. The gray arrows: guiding to the eye for the directions of the phase increase. (d) Reflection-free propagation of topological edge modes along the domain wall with sharp corners. The forward (backward) mode on the left (right) is excited at $f=4.05\mathrm{THz}$ with a surface current density (green arrows) ${\mathbf{J}}_{+}=\widehat{\mathbf{x}}+i\widehat{\mathbf{y}}$ (${\mathbf{J}}_{-}=\widehat{\mathbf{x}}-i\widehat{\mathbf{y}}$). Color plots are absolute squares of the out-of-plane electric field. Parameters: same as in Figs. 1 and 3.

As predicted by the TB model, the forward-propagating mode (red circle) has the AC surface charge density distribution with the phase increasing around the center holes in the counterclockwise direction, see Fig. 4(c). Due to the mirror symmetry in the nanoribbon, the backward-propagating mode (blue circle) is a mirror image of the forward-propagating mode. Therefore, these two modes form a chiral pair. We next investigated the topological robustness of these chiral edge states by modeling their propagation along a complex pathway with several types of sharp corners along the domain wall. For this simulation, each mode was selectively excited by circularly polarized external surface currents [green lines in Fig. 4(d): ${\mathbf{J}}_{\pm }=\widehat{\mathbf{x}}\pm i\widehat{\mathbf{y}}$]. The modeling results shown in Fig. 4(d) clearly indicate that the edge states are indeed immune to backscattering at all—60 deg, 90 deg, or 120 deg— lattice-preserving sharp turns.

3.5.

Midgap Topological Corner States

In addition to the ${\mathbb{Z}}_2$ invariant introduced above for characterizing the quantum orbital-Hall effect, the ${\mathcal{C}}^6$ crystalline insulator carries another binary invariant associated with topological corner charges:³⁵ $q_c\equiv 0$ (mod 1) for $t_{\text{out}}<t_{\text{in}}$ and $1/2$ for $t_{\text{out}}>t_{\text{in}}$ (see the Supplementary Material). The nontrivial half corner charge, given that the chiral symmetry is respected both inside the domain and at its termination, gives rise to a corner-localized mode in the middle of the bandgap.³⁴ This midgap corner state is a key-signature of second-order topological effects, where a band topology of a $d$-dimensional crystal is manifested in the $(d-2)$-dimensional observable. For 2-D crystalline PTIs, such observables are corner states that can be viewed as topologically robust and highly localized photonic cavities. Below we demonstrate that such nanocavities for GPPs can be developed using the metagate-graphene platform that enables us to embed a topologically nontrivial PTI inside a topologically trivial one.

One challenge in designing corner states in 2-D plasmonic systems is that, unlike the case of a waveguide-based platform^34,43 that is naturally described within the framework of a TB model, our plasmonic crystal does not possess strongly localized/bounded lattice sites. Therefore, any abrupt termination of a domain causes large shifts in effective on-site potentials, i.e., breaks the local chiral symmetry around a corner of the structure. To be specific, when graphene is set to be undoped or removed outside the domain boundary of our plasmonic crystal, the resulting structure is mapped to a TB model with severe redshifts in on-site potentials of boundary lattice sites (see the Supplementary Material). To avoid such unwanted perturbation along the boundary, we choose to imbed the topological domain within a trivial ${\mathcal{C}}^6$ crystalline insulator. In other words, while we are utilizing the same configuration as the one used for demonstrating bandgap-spanning chiral edge states in Fig. 4, we are now pursuing a different objective: to produce a sufficient bandgap between the two edge states in order to spectrally isolate the midgap state inside it. Therefore, we again used interface engineering to maximize the gapping of the edge states. As shown in Fig. 5(a), the size of the edge bandgap can be tuned from almost negligible [as was done in Figs. 4(a) and 4(b), where peripheral hole radii at the domain wall were chosen as $\overline{r}=(r_{2s}+r_{2e})/2$] to large. The latter objective is accomplished by staggering the radii of the peripheral holes between $r_{2s}$ and $r_{2e}$, as shown on the right side of Fig. 5(a). Such interface engineering amounts to edge roughening and enables us to increase the edge bandgap in such a way that the corner states appear near the middle of the bandgap.

Fig. 5

Designer interfaces between topologically distinct PTIs enable nanoscale waveguides and cavities. (a) Left: a uniform domain wall configuration used in Fig. 4 for supporting bandgap-spanning chiral edge modes. Right: staggered domain wall for supporting cavity-like corner modes. Small icons above the arrow denote the 1-D band structures for each domain wall type. (b) Left: gapping of the edge states supported by the staggered domain wall ($r_{1e}=165\mathrm{nm}$, $r_{2e}=120\mathrm{nm}$, $r_{1s}=100\mathrm{nm}$, and $r_{2s}=155\mathrm{nm}$) enables corner states. Right: DOS for the topological (expanded) domain ($10\times 10$ unit cells) imbedded into a trivial (shrunken) domain. Inset: the field profile of a midgap corner state. (c) A TB model qualitatively mimics the imbedded structure with a staggered domain wall. Small circles: sublattice sites, color-coded links: hopping amplitudes. (d) Same as (b), but for the TB model shown in (c).

Figure 5(b) shows the emergence of the expected topological corner state. An eigensimulation for a structure with 120-deg-angled corners reveals an eigenstate per corner that is well inside the edge bandgap and localized around the corner. Also, the mode profile of this corner state [Fig. 5(b) inset] features dominant amplitudes at two effective lattice sites that are the second corner-most, which exactly matches the reported result in an earlier work.³⁴ The density of states (DOS) was defined as $\mathrm{Im}[\sum _n\frac{1}{{\omega }_n-\omega }]$, where ${\omega }_n$ is the complex eigenfrequency (energy).

The proposed approach of embedding an SOTI into a topologically trivial structure and appropriately engineering the interface between them can be applied to any subwavelength platform. For example, we show in Figs. 5(c) and 5(d) that a generic TB model predicts the same results. Figure 5(c) describes the way we incorporate the edge roughening effect into the model. The relative strengths of hopping amplitudes are determined by considering the widths of metagate junctions, at which two adjacent effective lattice sites [Fig. 2(a)] meet together. Comparing to the setting drawn on the right side of Fig. 5(a), the hopping strengths are assigned to be proportionate to the junction widths. Then the following results in Fig. 5(d) confirm that the midgap corner state is guaranteed in this scheme. The reason for this is that our plasmonic crystal is not strictly chiral symmetric, as it is observed in Fig. 3(a) that the lower two bands are not exact mirror images of the upper two bands with respect to the midgap frequency. The exact emulation of the chiral symmetric Hamiltonian is not possible due to the weak localization of effective lattice sites and the presence of additional bands at higher frequencies.

To ensure that the DOS peak corresponding to a corner state is clearly standing out from the DOS due to the edge states, intrinsic graphene losses must be taken into account. At cryogenic temperatures ($T=20$ to 100 K), the intrinsic plasmonic loss in graphene is primarily due to the pseudomagnetic field induced by electron–phonon scattering:¹⁰ $\gamma /2\pi =[1.5\times (E_F\text{in}\mathrm{eV})\times (T\text{in}\mathrm{K})]\mathrm{GHz}$ (see the Supplementary Material). Thus, for the DOS calculation in Fig. 5(b) we used $\gamma /2\pi =20\mathrm{GHz}$ from $E_F=0.15\mathrm{eV}$ and $T=90\mathrm{K}$. Then the peak from the corner state contribution is clearly separated from other peaks, indicating that the observation of the corner state would be experimentally feasible at cryogenic settings.

3.6.

Nonlocal Optical Response of Dirac Electrons

So far, we have used the semiclassical Drude model [Eq. (2)] to describe graphene optical conductivity. Several experimental works^3,4,8,9 have, however, revealed that the quantum nonlocal response of Dirac electrons needs to be considered for the precise prediction of graphene plasmon dispersion. The nonlocal optical conductivity in homogeneously doped graphene is described accurately enough at the random phase approximation (RPA) level^2–4,31

Equation (2) is valid only when $\frac{v_F|\mathbf{q}|}{\omega }<1$ and $\hslash \omega \ll |E_F|$. Note that $f(z)$ is an increasing function in $z$ bounded below by 1. This momentum-dependent rise of conductivity stems from velocity matching between the Dirac electrons, traveling at $v_F$, and plasmons, traveling at $v_p=\omega /|\mathbf{q}|$; as $z={(v_F/v_p)}^2\to 1$. The electrons interact longer with each plasma wavefront and, therefore, exhibit an enhanced response.⁴

For nonlocal RPA conductivity under a spatially varying Fermi level, as in our metagate-tuned graphene, we adopt a semiphenomenological extension from the homogeneous case^31,47

Figure 6 shows how much the topological bandgap in our metagate-tuned plasmonic crystal gets shifted by the quantum nonlocal effect. We consider two cases: (a) $V_0=1\mathrm{V}$, as assumed so far, and (b) $V_0=0.05\mathrm{V}$, to demonstrate a case where nonlocal responses become significant. With $V_0=1\to 0.05\mathrm{V}$, GPP frequencies are roughly halved; $\omega \propto V_0^{1/4}$ [Eq. (1)]. As stated earlier, Eqs. (2) and (4) are valid only when $\rho \equiv \frac{v_F|\mathbf{q}|}{\omega }<1$. Here the dimensionless quantity $\rho$ is a measure of the importance of the nonlocal quantum effects.³¹ Bloch eigenstates in a periodic system contain many wavevector components, but here we take only the primary wavevector for simplicity in estimating $\rho$.¹⁸ The primary wavevector for dipole and quadrupole states around the plasmonic bandgap shown in Fig. 3(a) is $\mathrm{\Gamma }$ with the first band folding; that is, its magnitude is $|\mathbf{q}|=\frac{4\pi }{\sqrt{3}a}$. For case (a), $\omega /2\pi$ is around 4 THz, and thereby we get $\rho \approx 0.5$. As a result, the resonance shift shown in Fig. 6(a) is insignificant compared to the bandgap size, and the nonlocal effect converges rapidly at the first two orders. For case (b), we get $\rho \approx 1.0$ from $\omega /2\pi \approx 2\mathrm{THz}$. With $\rho$ approaching near unity, the nonlocal correction already exceeds the bandgap size at the first order and keeps increasing as the higher-order terms are added [see Fig. 6(b)]. Thus Eq. (4) is applicable only under high gate voltages and fails to be an effective description as $\rho \to 1$.

Fig. 6

Shifts of plasmonic resonance frequencies due to nonlocal responses; eigenfrequencies are evaluated at $\mathrm{\Gamma }$ above and below topological bandgap, using ${\sigma }_{\text{pheno}}^{[n]}$ in Eq. (4) for $n=0$ (local Drude limit) and $n=\mathrm{1,2}$ (nonlocal corrections up to the first or the second order). (a) With $V_0=1\mathrm{V}$, the bandgap is located around 4 THz and its size is about 0.3 THz. (b) With $V_0=0.05\mathrm{V}$, the bandgap is located around 2 THz and its size is about 0.15 THz.