Nanoscale plasmonic excitations in graphene¹ have proven their versatile potential for nanophotonic device applications^2,3 and fundamental studies of strong light–matter interaction^4,5 owing to their field-effect tunability and ultrasubwavelength confinement.^6,7 In high-quality graphene encapsulated between hexagonal boron nitride (hBN) layers, these plasmon polaritons can propagate over long distances with minimal attenuation.^8–10 In the past two years, researchers have sought new ways of combining the established concepts of topological photonics^11–17 with this attractive nanophotonic platform in order to further enhance its robustness. The earliest proposals for topological graphene plasmonics were based on nanopatterned graphene with broken time-reversal symmetry (TRS) via external DC magnetic fields.^18,19 However, DC magnets are not compatible with overall device miniaturization and even present a hindrance to the plasmon detection tools such as scanning near-field microscopes.^1,3–10 Moreover, graphene patterning, even if done over hBN-encapsulated structure for maximal protection,²⁰ would considerably degrade the carrier mobility and introduce severe damping of the graphene plasmon polaritons (GPPs). Therefore, it is imperative to find an approach that preserves both the graphene quality and the TRS.

Although TRS-preserving photonic topological insulator (PTI) phases resembling spin-Hall electronic insulators have been realized by emulating spin-1/2 Kramers pairs using right/left circulating²¹ or transverse electric/magnetic²² modes, this approach is not applicable here because GPPs are essentially scalar-like excitations due to their quasielectrostatic nature. The remaining two options for TRS-invariant topological phases are then valley-Hall insulators^23–25 and topological crystalline insulators.^26–28 Earlier proposals^29,30 for realizing these topological phases with GPPs suggested specific spatial landscapes of graphene optical conductivity, and our recent work³¹ suggested a practical path to realizing a GPP-based valley-Hall insulator using metagate-tuned graphene. A metagate, a designer nanopatterned planar electrode, imprints its structural symmetry onto the static carrier density landscape in graphene through the field-effect gating. Modulation of the carrier density translates into modulation of the optical conductivity, which creates a photonic crystal structure for GPPs. A similar method, which imparts a pattern to the dielectric spacer between graphene and backgate, has been experimentally shown to create a photonic crystal structure for GPPs while maintaining their long propagation distance.³² These approaches are not deleterious to the graphene quality because all patterning is done on the dielectric (or metallic) substrate while leaving the graphene itself intact.

In this paper, we present a metagate-graphene geometry that emulates topological crystalline insulators with GPPs. The employed photonic structure corresponds to a honeycomb lattice with a Kekulé distortion,³³ which is known to produce two topologically distinct phases depending on the distortion directions.^26,27 Orbital-locked chiral states arise at the boundary between two distinct domains, and these confined edge modes are shown to propagate freely around sharp corners with no reflection. Remarkably, a slightly modified structure can be used to emulate topological midgap corner states^34,35 of second-order topological insulators (SOTIs), d-dimensional crystals supporting topological states of a dimension less than d−1.^36–43 We show that our metagate geometry successfully realizes these phenomena by providing the first-principle calculations for the GPPs and qualitatively comparing the results to a generic tight binding model for Kekulé-textured graphene-like structures. While second-order topological effects have been demonstrated in macroscopic (photonic and phononic)^39–42 and waveguide-based^34,43 systems, where the lattice spacing is much greater than the photon/phonon wavelengths, our work presents the first implementation of subwavelength polaritonic SOTI phases hosted by a two-dimensional (2-D) material. We also estimate how plasmonic band structures are affected by quantum nonlocal effects of Dirac electrons.⁴

Chemical potential (Fermi level) landscapes from field-effect gating using a metagate were obtained by electrostatic simulations in COMSOL. In order to incorporate the quantum capacitance effect of graphene, we considered the electrochemical equilibrium condition, which was done by setting a graphene sheet as a surface charge density element.³¹ The static surface charge density on graphene was given as a function of the local potential [n=(eV/ℏvF)2/π]; therefore, an iterative solver was used. Photonic bands of the GPPs were calculated using both the in-house (Fourier-based quasielectrostatic) code³¹ and the finite-elements (COMSOL) simulations. The two codes agree well with each other, and we used them according to their computational efficiency: the in-house code for obtaining the results plotted in Figs. 1 and 4(c) and Fig. 6, COMSOL for the large-domain simulation shown in Figs. 4(d) and 5.

The dispersion relation q(ω) of the GPPs propagating in homogeneous graphene with a uniform chemical potential EF surrounded by an embedding medium with a dielectric constant ε is given by |q|=2πϵℏ2e2|EF|ω(ω+iγ),(1)where the semiclassical Drude conductivity,⁴⁴σDrude=i(e2/πℏ2)|EF|/(ω+iγ), of graphene is assumed. The finite temperature effect would be negligible as long as kBT≪|EF|. Here q and ω are the GPP’s 2-D wave vector and frequency, respectively, γ is the phenomenological Drude loss corresponding to free-carrier scattering, e is the electron charge, and ℏ is Planck’s constant. Graphene’s chemical potential (corresponding to the Fermi level of its free carriers) can be varied within ±1eV by tuning the carrier density n through field-effect gating: EF=sgn(n)ℏvFπ|n|.⁴⁵ Throughout this paper, we use the n-dependent renormalized expression for the Fermi velocity vF.⁴⁶

Equation (1) implies that a periodic modulation in EF 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 EF 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 r1 at its center and two identical holes of radius r2 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 EF produced by metagate-induced doping is shown in Fig. 1(b) (a specific case of r1=r2 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 r1=r2 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 (k and k′) 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 |q|≪a−1 (i.e., near the γ-point of the BZ). Also, it is previously reported that the screening through this periodically structured gate, along with EF 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=600nm is much smaller than the free-space wavelength of the THz radiation. For example, at the k(k′) point, we obtain an ultrahigh confinement ratio λvacuum/λGPP∼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 (r1≠r2) cases. The double degeneracy of the Dirac cones at Γ 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 Γ shown here are the fundamental building blocks^26,27 for the topological crystalline phases discussed below.

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 C6 point group symmetry. If r1>r2 (r1<r2), 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 tin is weaker (stronger) than that across adjacent unit cells tout in the case of r1>r2 (r1<r2).

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 (|p1⟩=[1;1;0;−1;−1;0], |p2⟩=[1;−1;−2;−1;1;2]), two quadrupole-like (|d1⟩=[1;−1;0;1;−1;0], |d2⟩=[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 tin and tout). 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 tin=tout, |p1,2⟩ and |d1,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 tin≠tout, as shown in Fig. 2(b). Specifically, the dipole (quadrupole) modes can be shown to reside below (above) the bandgap for tout<tin. Band inversion occurs between these four orbitals around the Γ-point for the tout>tin case, thus creating an insulating phase that is topologically distinct from the “topologically trivial” insulating phase realized in the case of tout<tin.

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 Γ point is characterized by a Z2 topological invariant^16,26,27 (Z2=0 for tout<tin and Z2=1 for tout>tin) 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.

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: r1=90nm and r2=150nm (topologically trivial case), and the topologically nontrivial case of r1=160nm and r2=110nm. Here r1 and r2 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.

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 EF(x,y) in graphene. As depicted in Fig. 3(b), the plasmonic modes near the Γ point display dipole (quadrupole) charge distributions below the bandgap when r1 is smaller (greater) than r2. 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.

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 r0 of the peripheral metagate holes (i.e., the holes at the domain wall) to be the average of the two r2 values from each domain: r¯=(r2e+r2s)/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 C6 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 kx), 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.²⁷

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): J±=x^±iy^]. 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.

In addition to the Z2 invariant introduced above for characterizing the quantum orbital-Hall effect, the C6 crystalline insulator carries another binary invariant associated with topological corner charges:³⁵qc≡0 (mod 1) for tout<tin and 1/2 for tout>tin (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 C6 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 r¯=(r2s+r2e)/2] to large. The latter objective is accomplished by staggering the radii of the peripheral holes between r2s and r2e, 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.

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.34 The density of states (DOS) was defined as Im[∑n1ωn−ω], where ω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:¹⁰γ/2π=[1.5×(EFineV)×(TinK)]GHz (see the Supplementary Material). Thus, for the DOS calculation in Fig. 5(b) we used γ/2π=20GHz from EF=0.15eV and T=90K. 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.

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σRPA(ω,q)=σDrude×f[(vF|q|ω)2],(2)where f(z) is a dimensionless correction factor that describes the nonlocal effect, f(z)=2z(11−z−1)=1+34z+58z2+O(z3).(3)

Equation (2) is valid only when vF|q|ω<1 and ℏω≪|EF|. 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 vF, and plasmons, traveling at vp=ω/|q|; as z=(vF/vp)2→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σpheno[n](ω,q,q′)=ie2πℏ2E˜F(q−q′)ω+iγf[n](vF2q·q′ω2),(4)where f[n](z) is the Taylor expansion of f(z) in Eq. (3) up to the n’th order term, and E˜F(q) is the Fourier transform of EF(r). Equation (4) reduces back to Eq. (2) under a constant Fermi level: E˜F(q)=EFδq,0→σpheno(ω,q,q′)=δq,q′σRPA(ω,q). Also, the local Drude model is equivalent to the case n=0. The effectiveness of Eq. (4) with n=1 has been verified in a recent experimental work.³

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) V0=1V, as assumed so far, and (b) V0=0.05V, to demonstrate a case where nonlocal responses become significant. With V0=1→0.05V, GPP frequencies are roughly halved; ω∝V01/4 [Eq. (1)]. As stated earlier, Eqs. (2) and (4) are valid only when ρ≡vF|q|ω<1. Here the dimensionless quantity ρ 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 ρ.¹⁸ The primary wavevector for dipole and quadrupole states around the plasmonic bandgap shown in Fig. 3(a) is Γ with the first band folding; that is, its magnitude is |q|=4π3a. For case (a), ω/2π is around 4 THz, and thereby we get ρ≈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 ρ≈1.0 from ω/2π≈2THz. With ρ 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 ρ→1.

We emphasize that our metagate-graphene structure was designed for experimentally feasible low-loss propagation of topological GPP excitations. Thus our structure is designed to operate at a few THz frequency range, where the dielectric losses from hBN or SiO2 phonon absorptions are vanishing. Then the 1/e decay length of the topological chiral GPPs in Fig. 4 can be estimated to be ωγ1|q|≈30a if we assume γ/2π=20GHz. Considering that the lateral confinement of the edge states shown in Fig. 4(c) is 3 to 4 unit cells, it should be possible to observe the propagation along the sharp turns shown in Fig. 4(d). Also, we expect that the second-order topological corner state can be efficiently excited, featuring a high quality factor (Q∼400). Again, our metagate scheme plays an essential role in enabling this low-loss operation through avoiding direct patterning of graphene.

Our metagate-tuned graphene system serves as an interesting case where topological photonics becomes a useful tool for unraveling condensed matter properties. Figure 6 suggests that we can electrically tune the proposed system from a condition in which Eq. (4) is effective, to another in which a modification is required. In this sense, our platform can be used to study quantum nonlocal effects in periodically doped graphene. For example, it is known that when the nonlocal response strength parameter ρ becomes greater than 1, graphene plasmons experience severe Landau damping.^4,18 The Bloch states in our plasmonic crystal are described by multiple ρ values from multiple wavevector components. It is an intriguing problem to ask whether the reflection-free propagation of 1-D topological Bloch states would be disabled by this selective Landau damping acting on a particular set of wavevectors. Thus exotic optical properties of Dirac electrons under strongly modulated Fermi surfaces can be revealed by examining how our plasmonic crystal deviates from the predictions based on local descriptions or trivial nonlocal corrections.

In conclusion, we have designed a nanophotonic platform for achieving two types of topologically robust optical elements based on graphene plasmons: a topological crystalline insulator for making reflections-free waveguides and a second-order PTI for making robust nanocavities. The platform consists of a single-layer graphene that is electrically biased using a designer metagate: an array of holes forming a honeycomb lattice with a Kekulé distortion. We demonstrate that the change from the propagating GPPs guided in two dimensions to cavity-like plasmons localized in all three dimensions can be accomplished by engineering the interface between two plasmonic PTIs: trivial and topological. We envision that the metagate geometry introduced in this work can also be used for realizing fractionalized Kekulé-vortex states³³ with graphene plasmons as well as many other topological phenomena such as strain-induced pseudomagnetic fields.⁴⁸ Also, by adding the second gate (a top metagate or a simple backgate beneath the existing metagate), the electric control over topological plasmons can be greatly enhanced to enable switching between distinct low-dimensional topological phases.⁴⁹ This would constitute a major advance over the single-gate geometry analyzed here, e.g., switching between two domain wall configurations in Fig. 5(a). Our metagate approach, combined with engineered domain walls between topologically different PTIs, paves the way to novel methods of manipulating nanolight in graphene and for developing topologically robust devices.