Topologically protected Dirac plasmons in a graphene superlattice

Topological optical states exhibit unique immunity to defects, rendering them ideal for photonic applications. A powerful class of such states is based on time-reversal symmetry breaking of the optical response. However, existing proposals either involve sophisticated and bulky structural designs or can only operate in the microwave regime. Here we show a theoretical demonstration for highly confined topologically protected optical states to be realized at infrared frequencies in a simple two-dimensional (2D) material structure—a periodically patterned graphene monolayer—subject to a magnetic field of only 2 tesla. In our graphene honeycomb superlattice structures, plasmons exhibit substantial nonreciprocal behavior at the superlattice junctions under moderate static magnetic fields, leading to the emergence of topologically protected edge states and localized bulk modes. This approach is simple and robust for realizing topologically nontrivial optical states in 2D atomic layers, and could pave the way for building fast, nanoscale, defect-immune photonic devices.

T opologically protected photonic states 1,2 exhibit remarkable robustness against disorder and backscattering, which make them extremely useful for the realization of defect-immune photonic devices. So far, such states have been explored based on reciprocal metamaterials 3,4 and photonic crystals 5,6 . But their realization has relied on sophisticated structural designs, which greatly complicate their fabrication and limit their miniaturization for nanoscale optical integration. Additionally, those topologically protected edge modes could couple to backscattering channels along opposite directions compatible with time-reversal symmetry (T-symmetry).
A more robust method to realize topologically protected photonic states consists in introducing T-symmetry breaking in periodic structures [7][8][9][10][11] , which leads to the opening of topologically nontrivial bandgaps at degeneracy points (e.g., Dirac points). T-symmetry breaking can be introduced via dynamical modulation of the refractive index 10 , although this approach is extremely challenging from the experimental viewpoint. A simpler alternative consists in using the magnetization of the material by exposing it to a static magnetic field. While this approach has been demonstrated in the microwave regime 9 , the weak magnetooptical (MO) response of most materials at visible and infrared frequencies renders it difficult to achieve substantial T-symmetry breaking in these technologically important spectral ranges.
Recent studies have shown that Dirac fermion (DF) systems possess a giant MO response in the infrared regime, such as the conducting surface of topological insulators 12 and a monolayer graphene [13][14][15] . Moreover, the plasmons supported by DF systems [16][17][18][19][20][21][22][23] exhibit deep-subwavelength confinement. This confinement makes the Dirac plamons extremely sensitive to external modulations 24 , and possibly more susceptible to the MO response than using light plane waves incident on structureless graphene, therefore providing substantial T-symmetry breaking in the infrared regime. Although edge magnetoplasmons already show obvious nonreciprocal behavior [13][14][15] and even form topologically protected states 25 at low frequencies, versatile applications are enabled by constructing topologically nontrivial bandgaps in the high-frequency regime and additionally steering plasmon propagation by using periodically patterned graphene.
Here we theoretically demonstrate that topologically protected plasmonic states can be robustly realized in DF superlattices constructed from single-layer graphene, thanks to the giant MO response of this material under exposure to static magnetic fields of only a few tesla. The superlattice is a honeycomb network constructed by graphene nanoribbons. The applied magnetic field induces asymmetry in the guided ribbon plasmon modes, thus resulting in directional coupling at the junctions of the structure. We show that, as a direct consequence of this directional coupling, localized modes are formed inside the superlattice, as well as topologically protected edge states at the boundary.

Results
Structure and working principle. We focus on graphene 26 because it sustains ultra-confined infrared plasmons that have been already observed in experiments 21,22 . Topological protection of plasmons in our proposed graphene superlattice can be intuitively understood as schematically illustrated in Fig. 1a. In contrast to previous designs based on photonic crystals, our structure consists of a network of waveguides (ribbons) in which topological protection emerges by analogy to the classical phenomenological picture of the quantum Hall effect. In this picture, free electrons in a Fermi gas follow circular cyclotron orbits under the influence of a magnetic field, and consequently they form a bulk insulating state and a topologically protected conducting state at the edge boundary. Although photons are charge-free, and consequently, a magnetic field cannot change their direction of motion, the splitting ratio of the plasmons at the junctions of the superlattice can be controlled by the MO effect due to the breaking of T-symmetry, leading to nonreciprocal directional coupling (Fig. 1b-f).
In the limit of total directional coupling (e.g., nearly 100% of the power exiting through the right side at each junction, as  Fig. 1 Topologically protected plasmons in a graphene superlattice. a Illustration of a honeycomb graphene superlattice with 100% right-coupling efficiency at each junction. Plasmons form localized vortex modes (1 and 2) and a unidirectional topologically protected edge state (3). b Dispersion relations of the two lowest-order modes sustained by a graphene nanoribbon (width W = 200 nm, Fermi energy E F = 0.2 eV) with and without a normal static magnetic field B (inset: schematic and electric-field distributions of guided modes for plasmon energy E p = 0.1 eV). c Schematic of a graphene nanoribbon junction with plasmons incident from branch 1 (unit incident amplitude), and scattered toward the three branches 1-3 with coefficients S 1,2−3 . d-f Electricfield distributions (50 nm above graphene for E p = 0.06 eV) produced by plasmon scattering at the nanoribbon junction with and without magnetic field, assuming the same ribbon parameters as in b and neglecting inelastic losses. Dashed lines delineate the graphene edges shown in Fig. 1f), the guided waves circulate around single lattice hexagons and become localized vortex states in the interior of the superlattice, mimicking the cyclotron motion of free electrons (routes 1 and 2 in Fig. 1a). Additionally, a unidirectional and topologically protected guided mode is produced at the edge boundary (route 3 in Fig. 1a). In a general situation characterized by nonzero backward and leftward scattering at each junction (e.g., Fig. 1e), localized bulk modes and a topologically protected edge wave are formed at a frequency within the bandgap opened by the T-symmetry breaking under magnetic field exposure 7 at the Dirac point of the honeycomb superlattice-these localized bulk modes and the topologically protected edge state exhibit analogous behaviors to the phenomenological picture in Fig. 1a, and further elaborated in our calculations below. We note that other periodic patterns different from the honeycomb structure, such as triangular or square lattices with degenerate band points, are equally applicable. The above design for realizing topologically protected plasmonic states is universal for waveguide networks incorporating a MO material, but the graphene nanoribbons are ideally suited because radiative losses can be neglected due to the large mismatch between plasmon and photon wavelengths. Additionally, experiments have demonstrated rather low propagation losses in this material 27,28 . More importantly, graphene exhibits a giant MO response in the infrared regime, so the supported plasmons are highly susceptible to the external magnetic field and allow us to realize the required directional coupling (Fig. 1b-f).
We use a finite-element method (FEM) to numerically solve Maxwell's equations and calculate the plasmon dispersion relations and associated field distributions (see "Methods" section). The response of graphene under a static magnetic field B is described by the Drude conductivity 13 as where E F is the Fermi energy, γ is the plasmon damping rate, ω c ¼ e B Á z ð Þv 2 F =E F is the cyclotron frequency 29 and v F ≈ 10 6 ms −1 is the graphene Fermi velocity. For simplicity, we ignore losses (γ = 0), which is reasonable in view of the long plasmon lifetimes (>500 fs) observed in graphene 27,28 . The local dielectric formalism is a reasonable approximation because the nanoribbons are hundreds of nanometers wide, and consequently, the Fermi energy exceeds the plasmon energy in all of our calculations, thus eliminating any dependence on the crystallographic orientation of atomic edges 30 .
The two lowest-order plasmonic modes of the nanoribbon are bonding and antibonding combinations of induced charge pileup at the ribbon edges (inset to Fig. 1b) 31 , with the second order mode showing a wavelength cutoff. Under an externally applied static magnetic field (Fig. 1b, dashed curves), the optical field distributions of the two modes become asymmetric, with field piling up toward different sides of the ribbon (inset to Fig. 1b). This nonreciprocal behavior is caused by the MO response, which is captured by the off-diagonal elements of the conductivity (Eq. (1)), similar to what happens for edge magnetoplasmons in a two-dimensional electron gas 32,33 and graphene 34 (Supplementary Note 1 and Supplementary Fig. 1). Below the cutoff frequency of the second mode, the ribbon only holds a single mode, so symmetry breaking can lead to directional coupling at the ribbon junctions, which is required for the realization of topological protection in the superlattice (Fig. 1c-f). The unit cell of the superlattice is formed by two nanoribbon junctions with lengths of the branches equal to L/2. For each junction, we calculate the scattering coefficients S 1,1-3 , starting from a plasmon wave of unit amplitude launched on the terminal of branch 1 and going to ports at the branches 1-3 (Fig. 1c). The results confirm that without magnetic field the powers at the outputs 2 and 3 are equal due to geometrical symmetry, while an applied magnetic field of 5 T leads to clear directional coupling (Fig. 1e, for a photon energy of 0.06 eV) and the right-coupling efficiency exceeds 96% for a magnetic field of 15 T.
Band diagrams of the graphene supperlattice. The magneticfield-induced directional coupling at the junctions leads to the opening of a Dirac point (Fig. 2). Separate bands at this point have nonzero Chern numbers 7 , which yield topologically protected edge modes within the bandgap, with the number of edge modes equal to the Chern numbers of the bands below the gap. To obtain the band diagram of the superlattice, we first calculate the scattering coefficients of the junction structure illustrated in Fig. 1c for different frequencies and magnetic field strengths via FEM simulations (Fig. 2a). Energy conservation leads to the relation S 11 j j 2 þ S 12 j j 2 þ S 13 j j 2 ¼ 1 between the scattering coefficients, which is satisfied by our numerical simulations.
Using Bloch's theorem, we calculate the band diagram of the superlattice (see "Methods" section)- Fig. 2b shows the results for an excursion along symmetry points within the first Brillouin zone. The nanoribbons in the unit cell are not identical because they have different orientations, so the band diagram is exactly the same as in a Kagome lattice, which contains a flat band and two Dirac points 35,36 . Our results confirm that, when a magnetic field of 4 T is applied, a large bandgap opens at the Dirac points K and K′ (Fig. 2b, d, f), which results in the formation of localized bulk modes within the bandgap region-in contrast, no bandgap is present in the absence of magnetic field (Fig. 2b, c, e). The bandgap opened by the magnetic field is topologically nontrivial, as demonstrated by the nonzero Chern number (C = 1) of the first band below the gap. This band has a nonzero Chern number by exchanging topological charge with upper bands separated by the Dirac point (see Supplementary Note 2 and Supplementary   Fig. 2 for a proof of this). The lower bandgap (in between orange and blue bands in Fig. 2b) is topologically trivial, corresponding to a zero sum of Chern numbers over all bands below 44 meV; the trivial character of this bandgap is evidenced by the fact that it does not close even when increasing the magnetic field up to 4 T. According to the bulk-edge correspondence principle 37 , the Chern number C = 1 of the first band below E 0 implies a unidirectional anticlockwise edge mode on the external boundary of the superlattice in the opened bandgap. The dispersion relations of the edge modes are revealed by the projected band diagrams, which are calculated along the zig-zag (Fig. 2c, d) and armchair (Fig. 2e, f) directions for a superlattice of finite period (N = 20) (see "Methods" section). The dispersion curves (red curve in Fig. 2d) in the bandgap correspond to edge modes propagating on the upper and lower boundaries of the superlattice-these have unidirectional group velocities, which imply that these edge modes are topologically protected. The directions of the group velocities also indicate that the edge mode on the external boundary of a finite surperlattice is anticlockwise, in agreement with that predicted from the Chern number.
Topologically protected localized and edge modes. To reveal the localized mode and prove the topological protection of the edge modes in the proposed graphene superlattice, we construct a numerical network with honeycomb topology and simulate the field evolution using the scattering coefficients represented in Fig. 2a (see "Methods" section). The field distributions for a directional point excitation (indicated by the red arrows) at the center of the network with magnetic field of 4 T is shown in Fig. 3a. For a photon energy of E 0 = 56.37 meV, which lies in the bandgap opened by the magnetic field, the directional coupling at the junctions causes the excitation power to circulate around a single lattice hexagon with enhanced intensity-this confirms the localization of the mode and corroborates the phenomenological picture we introduced in Fig. 1a. As this picture suggested, the localized vortex mode is the foundation for the topologically protected states, which can be clearly demonstrated by Fig. 3b. We now introduce a vacant defect by removing several nanoribbons from the center of the superlattice. Then, a directional point excitation adjacent to the defect generates a localized wave circulating around the defect unidirectionally without backscattering, which reveals the topological protection of the edge mode on the interior boundary of the superlattice.
To demonstrate that the edge mode on exterior boundary of the superlattice is topologically protected, we further cut the honeycomb network into a complex shape (Fig. 4a) and simulate the field evolution in successive time steps. In the presence of a B = 4 T field, a point source of energy E 0 placed at the boundary of the network can generate the sought-after one-way boundary mode, which efficiently propagate energy over sharp corners without back reflection (Fig. 4b). When the magnetic field is reversed in sign, the direction of the one-way edge state is also reversed (Fig. 4c). Finally, when the magnetic field decreases to 2 T, the edge mode is still preserved, although it is less localized at the boundary (Fig. 4d).
The effect of inelastic optical losses. We further show that inelastic losses do not deteriorate topological protection of the plasmons in the structure of Fig. 4, and just induce an attenuation of propagation, as shown in Fig. 5a. More precisely, we introduce losses into the Drude model through a damping rate γ = μE F /ev F 2 , where μ is mobility of graphene, for which we adopt the experimentally measured value in high-quality suspended graphene 38 . Although this value is measured at a temperature of 5 K, it can be maintained at room temperature by encapsulating the graphene in hexagonal boron nitride 27 . Figure 5a clearly shows that the excitation point source still generates unidirectional edge modes with good contrast between the two opposite directions, and that the edge modes retain a good capacity of topological protection. Inelastic losses only result in the attenuation of the modes along their propagation. In Fig. 5a, we show that edge plasmons can propagate over~30 periods. Since the individual branch structures under consideration can serve as classical and quantum interferometers [39][40][41] or as logic gates 42,43 , the proposed graphene superlattice provides a versatile platform to explore various complex nonreciprocal optical computing functions. Importantly, similar graphene superlattices of different sizes are equivalent according to basic electrostatic scaling laws of plasmons in two-dimensional (2D) materials 44,45 (Supplementary Note 3 and Supplementary Fig. 3), which reveal the possibility to further miniaturize the graphene superlattice in order to increase the propagation distance in units of plasmon wavelengths. In particular, when the geometry of the graphene superlattice in Fig. 4d shrinks by a factor of 2 and the magnetic field increases by a factor of ffiffi ffi 2 p , for a photon energy of ffiffi ffi 2 p E 0 , the response of the superlattice and the field distribution is equivalent to that Fig. 4d. However, when considering inelastic losses (Fig. 5b), the propagating distance increases considerably in units of the plasmon wavelength 31 compared with Fig. 5a, which means that the plasmon wave can propagate over many more periods in the network. Therefore, the noted scaling leads to a dramatic reduction in the effect of inelastic losses that should facilitate the design of practical devices.

Discussion
The realization of the topologically protected plasmonic states in the proposed graphene superlattice is experimentally feasible using state-of-the-art fabrication and measurement techniques, similar to those recently used to observe edge plasmon modes on graphene nanoribbons 46 (i.e., the elementary unit based on which we obtain topologically protected plasmons). To reduce inelastic losses, the sample can be fabricated by dielectric structuring rather than direct lithographic patterning of the graphene, thus relying on extended graphene encapsulated in hexagonal boron nitride, for which the measured mobility exceeds the values assumed here 27 . Additionally, our results can be readily extrapolated to include the effect of the substrate, the size of the structure and the level of graphene doping via an electrostatic scaling law formulated for 2D structures 44,45 (Supplementary Note 3 and Supplementary Fig. 3). An extension of this law to include the MO response reveals that the magnetic field strength required to achieve topological protection is below 1 T-such field strengths can be achieved using commercially available permanent magnets.
In conclusion, we have proposed and theoretically demonstrated how topologically protected plasmon modes can be realized robustly in single-layer graphene honeycomb surperlattice structures using experimentally attainable static magnetic fields. The key ingredients for this realization are the strong MO response of graphene, which induces T-symmetry breaking, and the ensuing directional coupling at geometrically symmetric x (μ m ) x (μ m ) B = 2√2 T Fig. 5 The effect of inelastic optical losses. a Optical electric-field distributions in the graphene superlattice in Fig. 4 incorporating loss, where the magnetic field is 2 T and the photon energy is E 0 . b Edge mode propagating in a superlattice with length size decreased by a factor of 2, assuming a magnetic field of 2 ffiffiffi 2 p T, for a photon energy of 2 ffiffiffi 2 p E 0 . The graphene mobility is 100,000 cm 2 V 1 s −1 in both a, b NATURE COMMUNICATIONS | DOI: 10.1038/s41467-017-01205-z ARTICLE NATURE COMMUNICATIONS | 8: 1243 | DOI: 10.1038/s41467-017-01205-z | www.nature.com/naturecommunications junctions between nanoribbons in the lattice. This nonreciprocal propagation leads to the formation of localized bulk modes and a topologically protected edge mode. Although we have focused on graphene, we expect other kinds of 2D materials such as topological insulators, semiconductor junctions, and conducting 2D materials, to exhibit similar phenomena, for example, by exploiting their long-lived phonon-polaritons. Our design provides an a simple, yet robust platform for fast speed, ultracompact, nonreciprocal optical computing networks, thus paving the way toward realistic applications of topological photonics.

Methods
Numerical simulations. We adopt FEM (COMSOL) to numerically calculate the electromagnetic field distribution and plasmon dispersion relations of uniform graphene nanoribbons, as well as the reflection and left/right-transmission coefficients of plasmons at the nanoribbon junction structures. We model graphene as a thin layer of thickness d and permittivity ε ω ð Þ ¼ 1 þ iσ ω ð Þ= ε 0 dω ð Þ, where σ(ω) is the frequency-dependent 2D conductivity in the Drude model (Eq. (1)). In the simulation, we take d = 0.5 nm as a reasonable value close to the d ! 0 limit. The scattering coefficients are calculated as the ratios between the complex amplitudes of the absorbed and input fields at ports on each of the three branches of the junction, as shown in Fig. 1c.
Calculation of band diagrams. In general, when plasmon waves are simultaneously launched at all three branches of a junction, the amplitude of input and output fields at the three branch terminals are linearly connected through a scat- where S is a scattering matrix, which contains three independent coefficients: S 11 = S 22 = S 33 , S 12 = S 23 = S 31 and S 13 = S 21 = S 32 , considering the symmetry of the structure. Due to T-symmetry breaking, S is asymmetric for B ≠ 0. Incidentally, out-coupling at the junction to free space is rather inefficient due to impedance mismatch (e.g., out-coupling losses are found to be below 0.01%). Neglecting out-coupling, the scattering matrix should be unitary, and therefore, energy conservation leads to the condition S 11 j j 2 þ S 12 j j 2 þ S 13 j j 2 ¼ 1. Using the simulation results for S 11 , S 22 , and S 13 shown in Fig. 2a, we can calculate the band diagrams of the graphene superlattice. We label the five pieces of nanoribbons within the unit cell as i = 1, …, 5. Each of them supports two counter-propagating waves, with field amplitudes A j i , where j = 1,2 denotes the directions of the waves corresponding to the red and black arrows in Fig. 6a, respectively. The amplitudes of these waves are connected by the scattering matrix of the single junction discussed above as ½A 2 T . Now, Bloch's theorem on the boundaries leads to ½A 1 4 A 2 4 T ¼ expðiK Á a 2 Þ½A 1 1 A 2 1 T and ½A 1 5 A 2 5 T ¼ expðiK Á a 1 Þ½A 1 2 A 2 2 T , where K = k x x + k y y is the Bloch wave vector. The band diagram of the superlattice in Fig. 2b is obtained from the condition of vanishing determinant of the above linear equations, while the wave functions of the nth band n K ð Þ j iare given by the corresponding eigenvectors.
Using a similar method, we also calculate the projected band diagram for a superlattice of finite width, as shown in Fig. 6b, c. Here, the unit cell is indicated by the two parallel dashed lines, which contain a larger number of nanoribbons. Again, each nanoribbon supports two counter-propagating waves. The scattering between these waves is determined by similar relations as above. Here, apart from the coefficients S 11 , S 12 , and S 13 of Fig. 2a, we also need to obtain the scattering coefficients for V-shape junctions appearing at the boundary of the superlattice. We use the same simulation method as for the three-lobbed junction. The resulting linear set of equations, combined with Bloch's theorem along the direction of periodicity, allows us to generate the projected band diagram shown in Fig. 2c-f.
Calculation of Chern numbers. With the whole set of wave eigenfunctions n K ð Þ j i obtained for all wave vectors K in the first Brillouin zone for any given nth band, the Chern numbers of this band can be easily calculated from the integral of the Berry connections along the close path around the boundary of the first Brillouin zone 8 i. According to the Stokes theorem, the Chern numbers can be equivalently calculated as a surface integral of the Berry curvature Ω n ¼ ∇ K n K ð Þ h ji∇ K n K ð Þ j iover the first Brillouin zone 8 : C n ¼ 2π ð Þ À1 R SdS Á Ω n . The Chern numbers of the bands in our study as calculated with these two well-established methods are in excellent mutual numerical agreement.
Calculation of field distributions in a superlattice. We focus on a superlattice of specific shape constructed out of N nanoribbons. We thus use a 2 N element vector A(t) to describe the distribution of field amplitude on the superlattice at time t, with A 2i (t) and A 2i−1 (t) (i = 1,…, N) representing the amplitude of two counterpropagating waves at the ith nanoribbon. Therefore, A 2i t ð Þ þ A 2iÀ1 t ð Þ j j gives the modulus of the field amplitude at the ribbons.
For the simulations, we introduce a point source at the initial time and calculate the field distribution in sequential time steps upon iteration. The point source is introduced at the ith nanoribbon through a prescribed pair of values A 2i (0) and A 2i−1 (0), by setting one of them to 1 (unidirectional point source for Fig. 3) or both of them to 1 (omnidirectional point source for Figs. 4 and 5), with all other elements of A(0) set to 0. Then, the field distribution of the waves scattered by the junctions of the superlattice after a time step (arbitrarily denoted 1) necessary for propagation along a ribbon length are approximated by A(1) = TA(0), where T is a 2 × 2 N matrix constructed from the S 11 , S 12 and S 13 coefficients (for the bulk), as well as the scattering coefficients for V-shape junctions, which describe scattering and interference between waves in the superlattice and are determined by the linear connecting relations of the nanoribbons (see above). Subsequent field distributions are obtained by iteration A(t + 1) = T[A(t) + A(0)], where the initial distribution A(0) is added at every step, so that the point source produces a persistent excitation.
Data availability. The data that support the findings of this study are available from the corresponding author upon request.  Fig. 6 Superlattice used for the calculations of band diagrams. a We consider an infinite superlattice constructed by repetitive translation of the unit cell along the principal axes a 1 and a 2 shown in Fig. 1a. b, c Superlattices of finite width with zig-zag (b) and armchair (c) boundaries, where the two vertical dashed lines indicate the unit cell