Topological electromagnetic waves in dispersive and lossy plasma crystals

Topological photonic crystals, which offer topologically protected and back-scattering-immune transport channels, have recently gained significant attention for both scientific and practical reasons. Although most current studies focus on dielectric materials with weak dispersions, this study focuses on topological phases in dispersive materials and presents a numerical study of Chern insulators in gaseous-phase plasma cylinder cells. We develop a numerical framework to address the complex material dispersion arising from the plasma medium and external magnetic fields and identify Chern insulator phases that are experimentally achievable. Using this numerical tool, we also explain the flat bands commonly observed in periodic plasmonic structures, via local resonances, and how edge states change as the edge termination is periodically modified. This work opens up opportunities for exploring band topology in new materials with non-trivial dispersions and has potential radio frequency (RF) applications, ranging from plasma-based lighting to plasma propulsion engines.


Introduction
The discovery that topological phases can exist beyond electronics [1][2][3] has led to significant interest in other wave systems such as photonics [4][5][6], plasmonics [7][8][9], polaritonics [10], acoustics [11,12], and even water waves [13,14].Of particular interest is the quantum anomalous Hall effect, also known as Chern insulators, which can protect unidirectional and back-scatteringimmune transport channels at the interface with normal insulators.Such transport channels have both fundamental and practical applications in fields such as optical communication, low-loss waveguides, and circulators.
To create photonic Chern insulators, time-reversal symmetry must be broken while retaining near-Hermiticity.This requires a permittivity or permeability tensor that breaks reciprocity, such that  T ≠  or  T ≠ .The properties of the underlying materials, along with their geometric aspects, determine the topological invariants (Chern numbers) of the electromagnetic band gaps and their transport properties.
Most current studies of photonic Chern insulators have relied on gyromagnetic materials (e.g.yttrium iron garnet [2,[15][16][17]) and external magnetic fields.These materials are mostly dielectric in nature, meaning their permeability remains positive in the frequency range of interest, resulting in Bloch modes that are delocalized in the photonic crystals.There has been recent interest in gyroelectric materials, which have the potential for large Faraday effects in magnetized plasmas in metals.This focus has centered on continua [18] and plasmonic crystals [7][8][9], where the Chern insulator phase is mostly composed of coupled plasmonic resonances localized at individual sites with Drude-like material responses and negative permittivity.
Here, we present a numerical study of photonic Chern insulators in plasma crystals with gaseous-phase plasma, which simultaneously exhibit both extended photonic bands and localized plasmonic modes in the RF regime.Our design is based on a 2D crystal of plasma cells placed in an external magnetic field.Without the magnetic field, the responses of the plasma elements are Drude-like, and the associated structure is known to support coexisting de-localized and localized modes in suitable polarizations [19][20][21].We explore the time-reversal broken generalization, which exhibits an interesting interplay between the Drude and Lorentzian dispersion due to the applied magnetic field causing cyclotron motions in the plasma.We propose a plasma crystal design that features a Chern insulator gap between de-localized photonic bands, coexisting with nearby dense groups of flat bands associated with localized plasmons of fixed-handedness.On termination of this crystal, we observe a rich interplay between localized plasmonic bands and de-localized chiral edge states inside the Chern insulating gap.Finally, we explore how local and de-localized edge modes evolve under continuous deformations of the interface between the Chern insulator and perfect magnetic conductors.The mode evolution can then be interpreted as a manifestation of the filling anomalies.
Our work is organized as follows: in Sections 2 and 3, we review the plasma dispersion without and with an external magnetic field.In Sections 4 and 5, we present the band structure of plasma crystals without and with an external magnetic field and their associated topological invariants.In Section 6, we explain the origin of the observed flat bands in calculations as localized surface plasmon polariton resonances.In Section 7, we explore how the chiral edge state dispersion evolves when the edge termination changes.Finally, we discuss the limitations existing in our calculations and practical aspects related to the experimental verification of our proposal.

Dispersion of plasma without magnetic field: Drude model
We start by reviewing the dispersion of the plasma permittivity , following the standard Drude model [22].Assuming that the positive ions are too heavy to move, the volume current density is solely contributed by electrons: J = −v.Here  is the volume density of electrons, and − is the electron charge.The equation of motion for electrons in the plasma crystals reads: Here,  e is the electron mass,  is the damping rate,  is the number density of electrons, and E is the electric field.For harmonic solutions at a fixed angular frequency , all temporal derivatives can be substituted via   → −i.Accordingly, the frequency-dependent plasma conductivity  can be written as Here, e  0 is the plasma frequency.Noting that the volume current density is also related to the electric polarization: J =   P = −iP, the Drude permittivity of plasma  can be defined as:

Dispersion of magnetized plasma
Following similar steps, we next show the permittivity tensor describing gaseous phase plasma placed in an external magnetic field (Fig. 1a).Following Eq. (3), we need to re-write the conductivity tensor ( σ), based on the updated equation of motion for electrons.
Taking into account the Lorentz force, the equation of motion becomes: Re-writing the equation using circular bases in the  plane: ,   ) and ( + ,  − ,   ) = ( ,   ), both matrices σ and ε become diagonal.Specifically, the conductivity tensor σ can be expressed as: Here  c =   e is the cyclotron resonance.Accordingly, the permittivity tensor of the magnetized plasma ε has the following form in a Cartesian basis: where the column vectors of the matrix  label the directions of the optical principle axes: Meanwhile, along the principal axes, the material dispersion can be expressed as: This result can be intuitively understood as follows: under the influence of the external magnetic field along ẑ, electrons undergo cyclotron motion and orbit at a fixed angular frequency of  c in the  plane.In the frame co-rotating (counter-rotating) with the electrons, the driving electric field of the same (opposite) circular polarization is thus offset in frequency by  c (− c ), giving rise to  + ( − ).Meanwhile, the applied electric field in the  direction is not affected by the cyclotron motion, so   remains the standard Drude dispersion.Next, we numerically compute the material dispersion as functions of frequency  = /2 using the typical values in gaseous phase plasma.The results are shown in Fig. 1, where the external magnetic field is set at 0.054 T and the corresponding cyclotron resonance is at  c =  c /2 = 1.5 GHz.The plasma frequency is controlled by the number density of electrons, which is set to be a practical value of  = 2.8×10 11 cm −3 throughout the calculations.Accordingly, the plasma frequency is at  p =  p /2 = 5 GHz.The real and imaginary parts of the permittivity are shown in blue and red, respectively.The parameter  refers to the damping rate, which mostly originates from the electron-ion collisions at room temperature.For helium plasma gas, the damping rate is roughly linearly proportional to the gas pressure: Using a practical pressure of  = 0.05 Torr, the damping is set to be  = 0.1 GHz throughout our calculations.
As expected from Eq. ( 8), both  + and  − are affected by the cyclotron resonance and deviate from the standard Drude model (  ).For example, at low frequencies (  ≈ 0),  + is positive and diverges as 1/  ;  − is negative and diverges as 1/  ; meanwhile, the Drude model   is negative and diverges as 1/  2 .These fundamental differences in scaling lead to challenges when fitting the dispersion to standard formalism in commercial software, as described later.We note that our description of the plasma medium is limited by a few approximations.Overcoming these approximations would lead to modifications of our results and will be discussed elsewhere.First, only electrons are assumed to move under the influence of external electromagnetic fields, while ions are assumed to be always stationary.Second, our permittivity ignores non-local effects, which leads to a frequency gap between surface plasmon polaritons traveling in opposite directions [23][24][25].Finally, we neglect the inhomogeneous broadening effect, and thus the bandwidth of our permittivity is defined solely by electron damping.

Band structures of plasma crystals without external magnetic field
Based on the dispersion in Eq. ( 8) [26], we compute the band structure of a square lattice of plasma cylinders placed in the air.Specifically, we focus on quadratic point degeneracy, which is protected by spatial symmetry and time-reversal symmetry.
The plasma photonic crystal unit cell is shown in Fig. 2a, where the lattice constant  is 6 cm, and the radius of the cylinder  is 1.5 cm.The external magnetic field along the  direction preserves the mirror symmetry in  (  ) and separates the electromagnetic modes into two mode types: to avoid possible terminology confusion, we follow the definition of Sakoda [27] to distinguish the E-polarization case (with the E-field parallel to the -axis) from the H-polarization case (with the H-field parallel to the -axis).We also use the E-case and H-case for short.Given that the gyroelectric response has response components only in the -plane but not along the -direction [Eq.( 5)], we focus solely on the H-case and disregard the associated E-case.
Without breaking time-reversal symmetry (), i.e.,  = 0, the plasma dispersion follows the standard Drude model, and the H-case can be calculated using a standard finite element method (FEM).The results are plotted along high-symmetry lines in the Brillouin zone (Fig. 2b).As shown, a pair of quadratic degeneracies are found at the  point in the Brillouin zone around 2.4 GHz, which is protected by the 90-degree rotation symmetry  4 and .Specifically, the two modes, marked as '+' and '−', have  4 indices of ±,.They are connected to each other by .The phases of the corresponding mode profiles, arg(  ), confirm the  4 indices of the two modes.We note that a set of flat bands are observed in the calculation (blue ribbon) with an upper-frequency bound of  p / √ 2 = 3.5 GHz, which is further discussed in the next section.

Band structure of the magnetized-plasma photonic crystals and Chern insulators
When an external magnetic field is applied, the plasma dispersion ( + and  − ) deviates from the Drude model, and the band structure proves difficult to be faithfully captured by using the standard material dispersion in various commercial software.Hence, we modify a standard FEM method to trace the desired number of bands, reliably resolving the band crossings.Similar to conventional approaches (see e.g., [28][29][30][31][32]), we employ the Floquet periodic boundary conditions to solve a quadratic eigenvalue problem in a square unit cell depicted in Fig. 2a; we also discretize the cell to get a weak FEM formulation.We use commercial FEM software (COMSOL M ™ ) -a flexible platform for constructing a system of customized coupled equations that employs the Portable Large Scale Eigenvalue Package (P_ARPACK).[33] P_ARPACK, a parallel version of the ARPACK software, [34] implements the Implicitly Restarted Arnoldi Method (IRAM) capable of solving large sparse eigenvalue problems for a select number of additionally constrained eigenvalues.COMSOL M ™ also provides a convenient user interface for controlling the key tunable parameters of its IRAM-based scalable eigensolver.Advantages of solving customized quadratic eigenproblems with P_ARPACK in COMSOL have been proven by numerous studies, see e.g., [31,35,36].
In contrast to the known methods largely employing the auxiliary equations for the polarization vector, we couple the weak-form equations for the current density and the E-field, achieving a stable performance with flexible tracking of a desired number of bands and unambiguous resolution of the band crossings.Section 3 provides a general equation for current density under an external magnetic field.In the H-case, we directly use the External Current Density Interface with an auxiliary algebraic equation (AE) for in-plane vectors, We then employ a normalized current density j   = J   /  0  2 p in a normalized Drude model The weak form is obtained by integrating the dot product of Eq. ( 10) with the test function j   = test(j   ) over the Drude material domain, An auxiliary equation for the current density is introduced to the COMSOL framework through the Weak Contribution interface, where is the E-field test function and  p =  p / is the plasma wave number.The approach exhibits good error convergence with accuracy controlled through the meshing density and the FE order.In contrast with [31,[35][36][37][38] we completely exclude the polarization vector, reducing the order of the auxiliary equations and improving the numerical accuracy at  → 0. The full details on the numerical implementation and verification of the IRAM-based eigensolver are not the main focus of the paper and these details will be published elsewhere.An example of the calculated band structures is shown in Fig. 3a when the external magnetic field is set to be  = 0.054 T. As time-reversal symmetry  is broken, the -point degeneracy is lifted, opening an 8% full energy gap, from 2.18 to 2.37 GHz (green ribbon).As the structure still maintains  4 symmetry, the Chern number  of the first band is necessarily non-trivial [39], where  2 (Γ) = 1 since the phase of an electromagnetic wave is locked at zero frequency, and  2 () = −1 because it originates from the time-reversal symmetry breaking of the degenerate modes shown in Fig. 2b.As a result, the first gap highlighted in green corresponds to a Chern insulator [40] and supports unidirectional transport channels, as shown next.Due to the magnetic field, the flat bands split into two regions (two blue ribbons), as explained in detail in the next section.

Flat bands from localized surface plasmon polariton resonances
In this section, we elucidate the origin of the flat bands observed above, from the viewpoint of localized surface plasmon polaritonic (SPP) resonances.We note that such flat band features are also commonly observed elsewhere, such as in metallic photonic crystals.The one unusual feature is related to the splitting of flat band regions under an external magnetic field (Fig. 3a).
It is more straightforward to understand the flat bands if we consider the local SPP resonances supported by a single plasma cylinder [7,8].As the cylinder has full rotation symmetry, the SPP Such a trend can be well observed in the comparison between the mode profile of  = −5 (more localized) versus the model profile of  = 1 (more extended, inset of Fig. 3b).Taken together, at large ||, the SPP resonances are well-localized near individual cylinders.Their modal overlaps are small, giving rise to little dispersion and thus the flat bands.Furthermore, the frequency of the flat bands (blue ribbons) also splits into two regions, approaching the two frequency bounds mentioned above.
Here we note that the number of calculated flat bands using our numerical method increases with increased mesh density -a common feature also observed in the literature [28] -although the frequency of the flat bands is always confined to the blue regions.Furthermore, the reliable azimuthal number of calculated modes is always limited by the numerical resolution.In our specific setting, modes with || > 10 are no longer reliable.See Supplementary Information for more details.
Finally, we verify the local SPP resonance frequencies in each plasma cylinder using the SPP dispersion along the interface between air and magnetized plasma, where a good agreement is found (Fig. 3c).The momentum of each local SPP resonance is determined as /, where  is the azimuthal number and  is the radius of the cylinder.

Evolution of the chiral edge state dispersion with changing edge termination
While the existence of chiral edge states (CES) is guaranteed at the interface between a Chern insulator and a trivial insulator, their exact dispersion depends on the details of the interface.In this section, we continuously change the interface configuration and study how the CES dispersion evolves accordingly.Our finding suggests that CES dispersion essentially reflects that a localized plasma resonance emerges at the interface when the plasma is cut through, and the frequency of the antenna state decreases with the shrinking of the plasma region.
The CES dispersions are calculated from the interfaces between a Chern insulator super-cell and a pair of perfect magnetic conductors (PMC), as shown in Fig. 4.Each unit cell in the Chern insulator has the same design and parameters in Fig. 3.The lower interface (blue) is fixed, while the top interface (red) is continuously modified as the distance  increases from 0 to .When / = 0, a pair of CES are found in the super-cell dispersion, including one band at the top interface (red) and another at the bottom interface (blue).As / increases, the CES on top (red) moves up in frequency until the PMC gets in contact with the plasma cell at / = 0.25.As / keeps increasing, only the low-frequency portion of the CES continues to move up, while the high-frequency part remains roughly unchanged until / = 0.5 when the PMC cuts the plasma cell in half.For even larger /, the CES dispersion becomes flatter and flatter, until when / reaches 0.7 and the CES splits into two bands: a CES band that traverses the topological gap, and a trivial band that is outside the topological gap.A second splitting occurs when / = 0.735, and a second trivial band emerges.As / keeps increasing, both trivial bands keep going down in frequency, traveling through the first bulk continuum and eventually disappearing at zero frequency  = 0 when / = 0.75.See Supplementary Information for the case when / = 0.74996.Finally, when / increases to 1, the interface configuration becomes the same as / = 0, and the CES dispersion goes back to the initial setting.
Taken together, as / increases from 0 to 1, both the first and the second band continua lose 1 mode, through the CES and the trivial edge state evolution.This can intuitively be understood as: when the super-cell loses a unit cell, each bulk continuum should also lose a mode, accordingly.

Discussion
Our simulation has a few limitations that should be noted.Firstly, the ions are assumed to be stationary and not moving, i.e., only electrons are allowed to move under the external fields.Introducing ions' motion would lead to an effective mass and lower plasma frequency.Secondly, our permittivity model ignores non-local effects, which results in a frequency gap between surface plasmon polaritons traveling in opposite directions [23,24].Thirdly, we have neglected the inhomogeneous broadening effect in the plasma, and the bandwidth of our permittivity is contributed solely by electron damping.Introducing the inhomogeneous broadening effect would further broaden the energy bands and reduce the effective size of the band gaps.
Despite these limitations, we believe that our proposal is feasible to demonstrate in an experiment.The required magnetic field of 0.05 T can be achieved, even over large areas, via commercial electromagnets or permanent magnets.Meanwhile, the required carrier number density of 2.8 × 10 11 cm −3 and the pressure of 0.05 Torr are also within the typical range in experiments.

Fig. 1 .
Fig.1.Permittivity dispersion of plasma medium in an external magnetic field. + ( − ) is the permittivity of right-handed (left-) circular polarization under an external magnetic field of 0.054 T.   is for polarization along the  direction.The real and imaginary parts of the permittivity are shown in blue and red, respectively.Both  + and  − are affected by the cyclotron resonance at  c = 1.5 GHz, while   is unaffected and remains to be given by the standard Drude model.

Fig. 2 .
Fig. 2. Band structure of a plasma photonic crystal without external magnetic field.(a) Schematic drawing of a photonic crystal made of gaseous plasma cylinders placed in the air.(b) Calculated H-case without magnetic field ( = 0), where a quadratic degeneracy, modes '+' and '−', with  4 = ± respectively, is found at the Brillouin zone corner.A set of flat bands (blue ribbon) are also observed in the calculation.

Fig. 3 .
Fig. 3. Band structure of magnetized-plasma photonic crystal featuring a Chern insulator gap.(a) Under an external magnetic field of  = 0.054 T, a full energy gap is opened (green ribbon), featuring a non-zero Chern number.Meanwhile, the flat bands split into two groups (blue ribbons).(b) Local resonances, labeled by different azimuthal numbers , are responsible for observed flat bands in a.Two example mode profiles ( = 1 and  = −5) are shown, both in amplitude (hot color map) and phase (gray-scale).(c) The frequency of local resonances agrees well with the surface plasmon polariton (SPP) dispersion at the interface between air and magnetized plasma.

Fig. 4 .
Fig. 4. Evolution of chiral edge states as the interface changes.A super-cell of electromagnetic Chern insulator (same design as in Fig. 3) is terminated by PMC.The bottom interface (blue) is fixed, while the top interface (red) changes from / = 0 to / = 1.The corresponding chiral edge states (CES) dispersion under different configurations of / is shown on the right.Red arrows indicate the evolution direction of the CES.