Topological magnetoplasmon

Classical wave fields are real-valued, ensuring the wave states at opposite frequencies and momenta to be inherently identical. Such a particle–hole symmetry can open up new possibilities for topological phenomena in classical systems. Here we show that the historically studied two-dimensional (2D) magnetoplasmon, which bears gapped bulk states and gapless one-way edge states near-zero frequency, is topologically analogous to the 2D topological p+ip superconductor with chiral Majorana edge states and zero modes. We further predict a new type of one-way edge magnetoplasmon at the interface of opposite magnetic domains, and demonstrate the existence of zero-frequency modes bounded at the peripheries of a hollow disk. These findings can be readily verified in experiment, and can greatly enrich the topological phases in bosonic and classical systems.

Plasmon is a unique type of bosonic excitations. Microscopically, it consists of collective motion of electron-hole pairs in a Coulomb interaction electron gas, whereas macroscopically it appears as coherent electron-density oscillations. Under most circumstances, it can be well described by classical density (and velocity) fields on the hydrodynamic level. These fields, like all other classical fields, are intrinsically realvalued and respect an unbreakable particle-hole symmetry. For a particle-hole symmetric system, its Hamiltonian H transforms under an antiunitary particle-hole conjugate operation C via C À 1 HC ¼ À H. (Here C 2 ¼ þ 1 for bosons.) This property ensures a symmetric spectrum o q ð Þ with respect to zero frequency, o q ð Þ ¼ À o À q ð Þ, in which q is the wavevector. The associated wave fields are thus superpositions of complex-conjugate pairs, F r; The spin-0 real-scalar field governed by the Klein-Gordon equation and the spin-1 real-vector field governed by Maxwell's equations also share the same feature.
Particle-hole symmetry greatly expands the classification of topological phases, according to the results of tenfold classification 28,29 . For example, the 2D quantum Hall phase belongs to the class-A in 2D with broken time-reversal symmetry T . The Su-Schrieffer-Heeger model and the recently studied phonon zero modes 12,15 belong to the class-BDI in one-dimension (1D) with particle-hole symmetry C 2 ¼ þ 1 and time-reversal symmetry T 2 ¼ þ 1. However, the class-D 2D topological phase with C 2 ¼ þ 1 and broken T has so far only been proposed in p þ ip superconductors 30 , which carry chiral Majorana edge states and Majorana zero modes.
In this work, we show that the historically studied 2D magnetoplasmon (MP) [31][32][33][34][35][36][37][38][39][40][41] belongs to the class-D 2D topological phase with unbreakable C and broken T . It is governed by three-component linear equations carrying a similar structure as that of the two-band Bogoliubov-de Gennes (BdG) equations of the p þ ip topological superconductor 30,42 . It contains a gapped bulk spectrum around zero frequency, and possesses gapless topological edge states and zero-frequency bound states (zero modes). Many properties of p þ ip superconductor can find their analogy in 2D MP. Therefore, 2D MP provides the first realized example of class-D in the 10-class table 28,29 . The experimentally observed edge MP states dated back to 1985 (ref. 37) are in fact topologically protected, analogous to the 1D chiral Majorana edge states. Some simulations of Majorana-like states in photonics have been reported [43][44][45][46][47] ; however, they only appear at finite frequencies around which the particle-hole symmetry does not rigorously hold, unless using high nonlinearity 48 . In addition, we are able to derive non-zero Chern numbers adhering to the band topology of 2D MP. On the basis of this, we are able to predict a new type of one-way edge MP states on the boundary between opposite magnetic domains, and the existence of MP zero modes on a hollow disk. Our prediction can be experimentally verified in any 2D electron gas (2DEG) systems, such as charged liquid-helium surface 33 , semiconductor junctions 34,38,49 and graphene 50,51 .

Results
Governing equations. We consider a 2DEG confined in the z ¼ 0 and r ¼ xe x þ ye y plane. The dynamics of MPs are governed by the linearized charge-continuity equation and the Lorentz force equation 52,53 . In the frequency domain, they are where r r; o ð Þ is the variation of 2D electron-density off equilibrium and j r; o ð Þ is the induced 2D current density. a r ð Þ is a space-dependent coefficient that gives a 2D local longitudinal conductivity, s r; o ð Þ¼ia r ð Þ=o. o c ðrÞ is the cyclotron frequency from a perpendicularly applied static magnetic field, B 0 (r) ¼ B 0 (r)e z . For the massive electrons chosen for our demonstration in this paper, a r ð Þ ¼ e 2 n0 r ð Þ mÃ and o c r ð Þ ¼ eB 0 r ð Þ m Ã c , in which n 0 (r) is the equilibrium electron-density distribution, m * is the effective mass and c is the speed of light 32,54 Here q ¼ |q|, and x q ð Þ ¼ 1 g is a q-dependent screening function 54 , in which d A and d B are the thicknesses, and E A and E B are the permittivities of the dielectrics on the two sides of 2DEG. A generalization of equation 3 can be made to include photon retardation (plasmon becomes plasmon polariton) 32,55,56 (see Supplementary Note 2). The Drude loss by electron-phonon collision can be included in this formulation by introducing a finite lifetime in equation 2 (ref. 54).
Bulk Hamiltonian and bulk states. We analytically solve for the case of homogeneous bulk plasmons, where n 0 (r), B 0 (r), a r ð Þ and o c r ð Þ are all constants. B 0 and o c can be positive or negative, depending on the direction of the magnetic field to be parallel or antiparallel to e z . We find that equations 1-3 can be casted into a Hermitian eigenvalue problem. The Hamiltonian H is blockdiagonalized in the momentum space, acting on a generalized current density vector J, Þg are the left-and right-handed chiral components of current density, and is the dispersion relation of conventional 2D bulk plasmon without a magnetic field. u p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the effective plasmon velocity originated from the screening of Coulomb interaction at long wavelengths by surrounding metals (see Fig. 1a). j D q; o ð Þ asymptotically equals u p r q; o ð Þ at long wavelengths.
The bulk Hamiltonian H q ð Þ and its long-wavelength limit read Strikingly, this long-wavelength three-band H q ð Þ has a very similar structure to that of the Bogoliubov-de Gennes hamiltonian of a 2D p þ ip topological superconductor 42 , which has two bulk bands. Here the odd number of bands is crucial for the topological consequences shown below. For an even number of bosonic bulk bands, it has been proved that the summed Chern number is always zero for the bands below or above the zero frequency 16 .
We can define the antiunitary particle-hole operation C and the antiunitary time-reversal operation T here, where K is the complex-conjugate operator and j q! À q stands for a momentum flipping. As can be checked, C acts as the complex conjugation for all the field components j x , j y and j D , while T additionally reverses the sign for j x and j y (refer to Methods: Cartesian representation). Our Hamiltonian has C symmetry, but broken T symmetry by the non-zero o c . The calculated homogeneous bulk spectra, plotted in Fig. 1b,c, are The corresponding (unnormalized) wavefunctions are The positive-and negative-frequency bands are gapped by o c . The zero-frequency band represents purely rotational currents (r r Á j ¼ 0, r r Â ja0) balanced by static charge density distribution.
It is worthwhile to point out that the 2D bulk plasmon discussed here is distinctively different from the surface plasmon (or surface plasmon polariton) in a 3D boundary (see Supplementary Note 3).
Chern number on an infinite momentum plane. Unlike the regular lattice geometries whose Brillouin zone is a torus, our system is invariant under continuous translation, where the unbounded wavevector plane can be mapped on a Riemann sphere 30,57 . As long as the Berry curvature decays faster than q À 2 as q-N, the Berry phase around the north pole (q ¼ N) of Riemann sphere is zero and the Chern number C is quantized 57 . For our 2D MP, we can verify this analytically, where dS q is the momentum-space surface element, and the term in the curly bracket is the Berry curvature. Since the Berry curvature changes its sign under the particle-hole symmetry, the positive and negative-frequency bands must have opposite Chern numbers, and the zero-frequency band, with odd Berry curvature, must have zero Chern number.
Majorana-type one-way edge states. We are able to calculate the analytical edge solutions in the long-wavelength limit q-0, and numerical edge solutions for unrestricted wavelengths, as plotted in Fig. 2. When q-0, the Hamiltonian in equation 6 becomes local, allowing us to replace q x and q y with the operators À i@ x and À i@ y and to solve for the edge states by matching boundary conditions 25,53 . We consider a 1D edge situated at x ¼ 0 caused by a discontinuity either in n 0 (x) (Configuration-I) or in B 0 (x) (Configuration-II). In Configuration-I as shown in Fig. 2a, we let the x40 region be filled with 2DEG, u p x40 ð Þ¼u p À Á , while the xo0 region be vacuum, u p xo0 ð Þ¼0 À Á and the magnetic field be uniformly Þ . Since we know that the change of Chern number across the edge is DC ¼ ± 1, there must be a single topological edge state present in this configuration. We look for solutions that behave like $ e À kx þ iqyy À iot , q x ¼ ik ð Þ , in the x40 region, where k40 is an evanescent wavenumber. The boundary condition is j x j x¼0 þ ¼ 0, meaning that the component of the current normal to the edge must vanish. Therefore, the edge solutions have to satisfy q y o c ¼ k ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi o 2 c þ u 2 p q 2 y À k 2 r , whose right-hand side is always positive. Since we have chosen o c ¼ o c j j40, q y must be positive too. The edge state is hence one-way propagating. (We discard the unphysical solution k ¼ q y that yields a null wavefunction.) The physical solution is k ¼ oc up , which gives an edge state spectrum, This is a gapless state running across the zero-frequency band as shown in Fig. 2c. It is a bosonic analogue to the chiral Majorana edge states in the topological superconductor. Its wavefunction is The zero-frequency edge state at q y ¼ 0 has the same finite decay length k ¼ o c u p and is well localized at the edge, despite being degenerate with the flat middle band.
In Configuration-II as shown in Fig. 2b, we let the 2DEG have a constant electron density throughout the whole space, u p x ð Þ ¼ u p , but the magnetic field have opposite signs in the two regions (o c . This is a novel configuration with DC ¼ ± 2 across the edge, and permits two topological edge states. The boundary conditions are j x j x¼0 À ¼ j x j x¼0 þ and j D j x¼0 À ¼ j D j x¼0 þ , meaning that both the normal current and density must be continuous across the edge. In our long-wavelength approximation, the first edge solution has the same spectrum and wavefunction as those of Configuration-I, except for a symmetric extension of the wavefunction into the xo0 region (compare the plots of r and j y in red in Fig. 2a,b). The second edge solution satisfies k ¼ À q y , and so q y must be negative to ensure k positive. This edge state is also one-way with a spectrum, as shown in Fig. 2d. Its wavefunction is It is antisymmetric for j y as plotted in green in Fig. 2b. The antisymmetry plus the continuity condition for j D here make the density oscillations vanish identically; therefore, the second edge state carries only current oscillations about the magnetic domain wall. Figure 2c,d gives our numerically calculated bulk and edge spectra that are not restricted by the small-q local approximation (refer to Methods: Numerical scheme). When q is small, the numerical results accurately agree with our analytical derivation above. When q is large, the first edge dispersion asymptotically approaches the bulk bands and eventually connects to them at q-N. The second edge dispersion drops gradually towards the zero-frequency bulk band and connects to it at q-N. It does show a (tiny) positive group velocity indicating the correct chiral direction, in spite of its negative phase velocity (refer to Fig. 2b,d). Physically though, the plasmon picture breaks down in the large-q regime, where electron-hole pair production takes place.
Besides, we should note that there exists a non-topological gapped edge state in Configuration-I (see the blue curve in Fig. 2c), which is consistent with the literature 54 . At long wavelengths q-0, this state resembles a bulk state bearing an excitation gap of o c . Its edge confinement shows up only at short wavelengths. As elaborated below, under a parameter evolution from Configuration-I to II, this edge state evolves, remarkably, into the second gapless topological edge state in Configuration-II.
One-way propagation of edge states. Owing to the topological protection, one-way edge states are immune to backscattering from a random defect. In Fig. 3, we display two snapshots from our real-time simulation. Edge states are excited by a point source oscillating at a frequency of 1 2 o c within the bulk gap. The generated waves propagate only towards the right. A sharp zigzag defect has been purposely inserted into the route of propagation. The waves then exactly follow the edge profile and insist on propagating forward without undergoing any backscattering. The slightly visible fluctuation to the left of the point source is completely local (non-propagating). It is caused by the unique long-range Coulomb interaction in this system, different from the more familiar photonic and acoustic systems. Figure 3a corresponds to the traditional Configuration-I, where the physical edge is formed by the density termination of 2DEG at vacuum. Although the one-way nature in this configuration is known historically, its absolute robustness due to the topological protection is less known and is manifested here. Figure 3b corresponds to our new Configuration-II. Only the first edge state classified in Fig. 2 is excited here; the second edge state has an energy too close to the band edge o c and hence is not excited here. On the basis of our argument on topology above, Configuration-II permits protected one-way edge states on a magnetic domain boundary even if the 2DEG may be homogeneous. Figure 3b vividly demonstrates this scenario. Furthermore, it shows the perfect immunization to backscattering when the two magnetic domains drastically penetrate each other.
Evolution of edge states. We find that the adiabatic mode evolution of the MP edge states described by a three-band model here is rather different from that in 2D topological insulators or superconductors, which can be commonly described by a two-band or four-band model. The zero-frequency bulk band here plays a critical role for the appearance and disappearance of the topological edge states. We have investigated the evolution corresponding to the aforementioned two configurations shown in Fig. 2.
We first study the interaction between two topological edge states propagating along the opposite directions on two 2DEG edges in the same magnetic field. The result is shown in Fig. 4a. When we gradually narrow the spacer, we let the equilibrium density in the spacer to gradually change from 0 to the bulk density n 0 . This allows the two topological edge states to couple and end up with a bulk configuration. Instead of opening a gap and entering the top band (which would happen in the conventional two-band system), the two topological edge states here (marked in red) fall down into the flat middle band. In addition, there are two non-topological edge states (marked in blue, and refer to Fig. 2c) in this configuration. They move completely into the bulk during this process.
We then study the interaction between two topological edge states propagating along the same direction on the two 2DEG edges under opposite magnetic fields. The result is shown in Fig. 4b. When we gradually narrow the spacer, we find that one of the edge states (marked in red) goes slightly upwards, while the other falls downwards into the flat middle band. In the meanwhile, one of the non-topological edge states (marked in blue) detaches from the top band, continuously bends downwards MP on a hollow disk. For infinitely long edges as discussed above, the first topological edge state when k y -0 merges into the zero-frequency middle band, as can be seen in Fig. 2c,d. Importantly, we want to quest whether such zero-frequency modes preserve (not being gapped) when we wrap a long edge into a small circle. We want to investigate the behaviours of 2D MP on, for instance, a hollow disk geometry depicted in Fig. 5a.
We write the low-energy long-wavelength Hamiltonian in realspace polar coordinates and still use the chiral representation The azimuthal symmetry ensures the eigensolutions of j r , j f and j D to pick up a common factor e inf , where n is the integer of discretized angular momentum. This low-energy long-wavelength problem can be solved semi-analytically. In all situations, j r and j f are derived quantities from j D , Depending on the energy to be above or below the bandgap, we have j D r; n; o ð Þ¼ A n J n k r r ð ÞþB n Y n k r r ð Þ; o 2 oo 2 c À Á ; C n K n k r r ð Þ; 0oo 2 oo 2 c À Á : where k r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi =u p , J n k r r ð Þ and Y n k r r ð Þ are the n-th order Bessel and Neumann functions, and K n k r r ð Þ is the modified Bessel function. A n , B n and C n are coefficients. The no-normal-current boundary condition j r r; n; o ð Þ j r¼a ¼ 0 determines the allowed eigenfrequencies o for every given n. In this approximate model, the solutions above the bandgap o 2 oo 2 c are all continuous bulk-like and within the bandgap 0oo 2 oo 2 c are all discrete edge-like, as shown in Fig. 5b.
Majorana-type bound state: zero mode. We are particularly interested in the limit o ! 0. According to equations 19 and 21, the boundary condition can be satisfied if and only if n ¼ 0, which completely annihilates the radial current, j r (r, 0, 0) ¼ 0. The zero-mode profile is given by ; ð22Þ j D describes a static charge distribution; j f describes a d.c. circulating current. The electric field generated by j D balances out the Lorentz force due to j f and B 0 . The overall amplitude and phase in A 0 can be arbitrary, but j D and j f are phase-locked to each other by a ratio of sgn(o c ), that is, the sign of magnetic field. In the chiral representation, with the arbitrary phase ignored, the zero-mode wavefunction behaves in the below manner: We see that the Majorana-type topological edge state within the bandgap indeed preserve as o ! 0 and n ¼ 0, on the finite-sized inner edge of a hollow disk. It is easy to verify that shrinking the hole radius will only reduce the number of discrete edge states at a finite n, but not kill the zero mode. These observations are consistent with the k y -0 limit of the infinite straight edge case (see Fig. 2c,d). Nevertheless, we should cautiously note that the zero mode here is not isolated. It degenerates with a large number of other zero-frequency modes in the middle band (see Supplementary Note 4).
The guaranteed existence of zero mode in a hollow disk geometry for 2D MP is reminiscent of the Majorana zero modes in the two-band p þ ip topological superconductor 30 . In the latter case, a p flux inside the vortex core is needed to balance the antiperiodic boundary condition in the f direction. By comparison, our three-band Hamiltonian here hosts zero mode with periodic boundary condition, that is the wavefunction in equation 23 returns to itself when f rotates 2p. More generally, similar kind of discretized edge states and zero modes must exist on a noncircular geometry as well, and can be on both the inner and outer edges, if the geometry is finite. For a finite geometry, the inner and outer edge zero modes must come out in pairs. If one closes the inner boundary and turns the hollow disk into a solid disk, then both modes must vanish. Otherwise, the vortex currents associated with both zero modes are singular at the centre, and are unphysical.

Discussion
Although the configuration of 2D MP looks almost identical to that of the quantum Hall effect (QHE), there are fundamental differences 39 . QHE deals with electron (fermionic) transport, whereas 2D MP deals with collective electron-density oscillations (bosonic excitations). In QHE, the bulk is electronically insulating, whereas in 2D MP, the bulk is electronically conducting. QHE experiments usually require a high magnetic field (\5 T) and a high electron density (\10 11 cm À 2 ). By comparison, MP experiments normally does not require a high magnetic field or a high electron density. Moreover, the bandgap in QHE is for electron transport, while the bandgap for 2D MP is for electron-density oscillations.
2D MP belongs to the topological class-D, which is different from the quantum Hall class-A. The topological edge states in MP are not governed by the traditional bulk-edge correspondence of the quantum Hall states. The traditional rule states that the number of gapless edge states inside a gap is equal to the sum of all the bulk Chern numbers from the energy zero up to the gap. In our topological MP, the lowest bulk band (zero-frequency flat band) contributes zero Chern number. A generalized bulk-edge correspondence must be established by considering the particlehole symmetry that extends the spectrum into the 'redundant' negative-frequency regime. This accounts for the Berry flux exchanged across the zero frequency. Only by doing so, the sum of all the bulk Chern numbers, spanning both the positive and negative-frequency regimes, can be correctly kept zero.
Experiments for 2D plasmon were performed by pioneering researchers dated back to 1970s on the charged liquid helium surface 33  Our prediction of new one-way edge MP in Configuration-II has not been reported, despite there have been similar studies in the QHE systems 63,64 . Our prediction can, in principle, be verified in any 2DEG system as long as a pair of opposite ARTICLE magnetic fields is applied across. The magnetic domain boundary can be experimentally created by, for example, two concentric solenoids of opposite currents or ferromagnetic materials. The gapless one-way edge modes can be mapped out by microwave bulk transmission and by Fourier transforming near-field scans along the edges, as demonstrated in ref. 65 for a topological photonic system. In Fig. 6, we provide a schematic experimental design. A semiconductor heterojunction serves as the platform for 2DEG. A ferromagnetic film is grown on the top and is polarized into up-and down domains. For the circular domain boundary as sketched, discretized one-way edge states can travel around. A meandering coplanar waveguide can be fabricated on the back of the substrate. One can then send in and receive microwave signals at a frequency below o c , and in the meanwhile monitor the characteristic absorption, which signifies the excitation of one-way edge plasmons.
Our proposed zero-frequency bound modes feature localized charge accumulation (inducing a static electric field) and circulating current (inducing a static magnetic field). These effects are in principle measurable when the zero mode is actuated by charge or current injection.
In summary, we have revealed the salient topological nature of 2D MP as the first realization of the class-D topological phase. The predicted new Majorana-type one-way edge states and zero modes can be verified experimentally. Our work opens a new dimension for topological bosons via introducing the bosonic phase protected by the particle-hole symmetry. We anticipate similar realizations for other bosonic particles 66 and discoveries of new topological phases, after combining the particle-hole symmetry with other symmetries 67 , for example, the timereversal and many kinds of spatial symmetries 10,68 .

Methods
Theoretical model. Our theoretical model of 2D MP is based on a hydrodynamic formalism. In equilibrium, the 2D electron number density is n 0 (r), which is determined by the device structure, doping concentration and gating condition 38 . Off equilibrium, the density changes to n(r, t). For 2D MP, we are only concerned with the small deviation of 2D charge density, r(r, t) ¼ À e{n(r,t) À n 0 (r)}, and the induced 2D current density up to the linear order, j(r, t) ¼ À en 0 (r)v(r,t), where v(r, t) is the local velocity of electron gas restricted to move in the z ¼ 0 plane only. For massive electrons, the conductivity coefficient is a r ð Þ ¼  69,70 . Hence, our theory works for both massive and massless electrons. We consider a model system shown in Fig. 1a. With different choices of the material constants and structural parameters, it can lead to different practical 2DEG systems, such as charged liquid-helium surface, metal-insulatorsemiconductor junction, top-gated graphene transistor, and so on. In Fig. 1a, the 2DEG lies in the z ¼ 0 plane. The 0ozod A region is filled with an insulator (or vacuum) of the dielectric constant E A . The À d B ozo0 region is filled with another insulator (or semiconductor) of the dielectric constant E B . The z4d A and zo À d B regions are filled with perfect metals whose dielectric constant is E M ¼ À N in the (radio to microwave) frequency range of MP. Encapsulating the system by two metals mimics the experimental configurations with top and bottom electrodes, and in the meanwhile, cutoffs the theoretically infinitely long-ranged Coulomb interaction.
For the results shown in the main article, we simply choose Cartesian representation. In the main article we have used the chiral representation for the bulk Hamiltonian and generalized current vector, in order to demonstrate the similarity between our 2D MP and the 2D p þ ip topological superconductor. However, it is sometimes more convenient to adopt the traditional Cartesian representation, especially when we solve for the edge modes and apply boundary conditions on the edge. In the Cartesian representation, the bulk 0 u p q x À io c u p q x 0 u p q y þ io c u p q y 0 0 @ 1 A : The particle-hole and time-reversal operators in this Cartesian representation are Clearly, C is just a complex conjugation on the matrix or vector elements (combined with a momentum reversal q-À q here, when written on the singleparticle momentum bases).
Numerical scheme. To numerically solve the nonuniform edge state problem, which is unrestricted by the long-wavelength approximation, we expand the governing equations 1 and 2 in the main article into plane waves and use the form of interaction in the momentum space equation 3. We can derive a matrix equation in the Cartesian representation, The submatrices U x q; q 0 ð Þ, U y q; q 0 ð Þ, Q x q; q 0 ð Þ, Q y q; q 0 ð Þ, W q; q 0 ð Þ, and the subcolumn vectors J x q 0 ð Þ, J D q 0 ð Þ, J y q 0 ð Þ all have the dimension of the number of plane waves used for expansion.
The elements of submatrices are U x q; q 0 ð Þ¼2pã q À q 0 ð Þ q 0 x q 0 x q 0 ð Þ ; ð27Þ U y q; q 0 ð Þ¼2pã q À q 0 ð Þ q 0 y q 0 x q 0 ð Þ ; ð28Þ Q x q; q 0 ð Þ¼d q; q 0 q 0 x ; ð29Þ Q y q; q 0 ð Þ¼d q; q 0 q 0 y ; ð30Þ and W q; q 0 ð Þ¼o c q À q 0 ð Þ : whereã Dq ð Þ and o c Dq ð Þ with Dq ¼ q À q 0 are the Fourier transform of a r ð Þ and o c r ð Þ. For the system with edges along the y axis, the expansion only needs to be done in the x direction while keeping q y ¼ q 0 y as good quantum numbers. Typically, we use 128 plane waves to do the expansion, and the standard eigensolver to obtain all the eigenfrequencies and eigenvectors.
The 2D real-space calculation for Fig. 3 is performed by a homemade finitedifference time-domain code. The field quantities are defined on a square grid. Equations 1 and 2 are adopted for the real-time evolution after replacing À io with q t . A point-source term of a given driving frequency is added to the right-hand side of equation 2.
Data availability. All relevant data are available from the authors on request.