# The dispersion and propagation of topological Langmuir-cyclotron waves in
cold magnetized plasmas

Yichen Fu<EMAIL_ADDRESS>Hong Qin<EMAIL_ADDRESS>Princeton
Plasma Physics Laboratory and Department of Astrophysical Sciences, Princeton
University, Princeton, NJ 08540

###### Abstract

Topological Langmuir-Cyclotron Wave (TLCW) is a recently identified
topological surface excitation in magnetized plasmas. We show that TLCW
originates from the topological phase transition at the Langmuir wave-
cyclotron wave resonance. By isofrequency surface analysis and 2D and 3D time-
dependent simulations, we demonstrate that the TLCW can propagate robustly
along complex phase transition interfaces in a unidirectional manner and
without scattering. Because of these desirable features, the TLCW could be
explored as an effective mechanism to drive current and flow in magnetized
plasmas. The analysis also establishes a close connection between the newly
instituted topological phase classification of plasmas and the classical CMA
diagram of plasma waves.

##  I Introduction

Surface waves in plasmas have been extensively studied [1, 2, 3] due to their
wide range of applications [4, 5, 6, 7, 8, 9]. Recently, a new class of
surface modes in plasma produced by topological phase transition have been
investigated [10, 11]. Here, we report studies on the dispersion and
propagation of the newly identified topological surface mode called
Topological Langmuir-Cyclotron Wave (TLCW) in cold magnetized plasmas.

Surface plasma waves are fluctuations propagating along the boundary between
two regions. The energy of surface waves is usually concentrated at the
boundary and attenuates rapidly in the normal direction. On the other hand,
bulk (or body) waves are not localized in the boundary and propagate through
plasmas. Bulk and surface waves were traditionally considered distinct waves
and studied separately. It was realized in recent decades that there exists a
connection between some surface waves and the topological properties of bulk
waves. In the study of the integer quantum Hall effect [12], it is found that
the bulk states in the quantum Hall systems can be characterized by a
topological invariant $n\in\mathbb{Z}$ called the Chern number. At the
interface between a quantum Hall state and vacuum, chiral surface (edge)
states exist [13] whose quantity is equal to the Chern number in the bulk
state. Such a correspondence between bulk topological invariants and the
number of chiral edge states can be generally summarized as the bulk-edge
correspondence [14]. The bulk topology and its relation with surface states
have been widely studied not only in quantum systems [14, 15, 16, 17], but
also in classical systems, such as those in photonics [18, 19, 20], acoustics
[21, 22, 23], mechanics [24, 25, 26] and fluid dynamics [27, 28, 29, 30, 31,
32, 33].

Magnetized plasmas have been studied from the topological point of view as
well. Cold plasmas in a constant magnetic field $B_{0}\hat{z}$ were found to
have non-trivial topology. When variation in the $z$-direction can be ignored,
i.e., $k_{z}=0$, Chern numbers for the X waves in 2D plasmas were calculated
[34, 35] and related to surface waves via the bulk-edge correspondence [36].
When $k_{z}\neq 0$, the non-trivial topology of the plasma as a 2D system were
found [10, 37] and its topological phase diagram was established [11]. In
particular, under-dense plasmas can be regarded as a Weyl semimetal [38, 39]
and host a surface state with the Fermi-arc structure [17]. The topological
matter properties of ideal MHD plasmas with magnetic shear [40], kinetic
plasmas [41], and surface waves within the continuous spectrum [42] were also
investigated.

The topological surface wave in the present study is called Topological
Langmuir-Cyclotron Wave (TLCW), because it is localized at the interface
between two plasmas in different topological phases separated by the phase
transition at the resonance between the Langmuir wave and the cyclotron wave
[11]. Based on a one-dimensional model, we study the dispersion of the TLCW
using its isofrequency contours, from which the main physical properties of
the TLCW, such as the unidirectional (chiral) propagation and the immunity of
scattering, can be directly observed. The isofrequency contours also
illustrate how the topology of the index-of-refraction surface of bulk plasma
waves evolves when the TLCW exists. This establishes an interesting
relationship between the topological classification by the Chern number and
the well-known CMA diagram [43, 44, 45] of plasma waves. Two- and three-
dimensional time-dependent simulations have been performed for the linearized
fluid equations. The TLCW are excited by a Gaussian source at a given
frequency and propagates in various configurations, confirming the physical
properties found in the analysis using the isofrequency contours. The momentum
and angular momentum carried by the TLCW are also studied. Because the
excitation and propagation of the TLCW are topologically protected, it could
explored as a robust mechanism to inject energy and momentum into magnetized
plasmas.

The paper is organized as follows. Section II briefly introduces the
theoretical model and the topological matter properties of cold magnetized
plasmas. The frequency range and isofrequency contours of the TLCW are given
in Sec. III. The connection between the isofrequency contours and the CMA
diagram is addressed. Section IV describes the numerical algorithms and 2D and
3D simulation results of the TLCW.

##  II Theoretical model

This section describes the governing equations for dynamics in cold magnetized
plasmas and briefly discusses the origin of the TLCW. We consider a cold
stationary plasma with immobile ions, and the background magnetic field is
constant, i.e., $\bm{B}_{0}=B_{0}\hat{z}$. Since there is no pressure in cold
plasmas, any given plasma density profile $n(\bm{r})$ is an equilibrium. With
proper renormalization, the linearized fluid equations can be written as [46,
11]:

$\displaystyle\partial_{t}\bm{v}=-\omega_{\mathrm{p}}\bm{E}-\Omega\bm{v}\times\hat{z},$
(1)
$\displaystyle\partial_{t}\bm{E}=\nabla\times\bm{B}+\omega_{\mathrm{p}}\bm{v},$
(2) $\displaystyle\partial_{t}\bm{B}=-\nabla\times\bm{E},$ (3)

where $\bm{v,}\bm{E},\bm{B}$ are perturbed velocity, electric and magnetic
fields,
$\omega_{\mathrm{p}}(\bm{r})=\sqrt{n(\bm{r})e^{2}/\epsilon_{0}m_{\mathrm{e}}}$
is the local plasma frequency, $\Omega=eB_{0}/m_{e}$ is the cyclotron
frequency, $e>0$ is the elementary charge, $m_{\mathrm{e}}$ is the electron
mass, and $\epsilon_{0}$ is the vacuum permittivity.

In a homogeneous bulk plasma, for each eigenmode with frequency $\omega$ and
wavenumber $\bm{k}$, Eqs. (1)-(3) reduce to

$H|\psi\rangle=\omega|\psi\rangle,$

where $|\psi\rangle=(\bm{v},\bm{E},\bm{B})^{\mathrm{T}}$ and

$\displaystyle H(\omega_{\text{p}},\Omega,\bm{k})$
$\displaystyle=\begin{pmatrix}-\mathrm{i}\Omega\hat{z}\times&-\mathrm{i}\omega_{\text{p}}&0\\\
\mathrm{i}\omega_{\text{p}}&0&-\bm{k}\times\\\ 0&\bm{k}\times&0\end{pmatrix}$
(4)

is a $9\times 9$ Hermitian matrix. For each $\bm{k}$, the system has 9
eigenmodes

$\omega_{n},|\psi_{n}\rangle,n=-4,-3,\cdots,3,4\thinspace.$

ordered by its value, i.e., $\omega_{i}\leq\omega_{j}$ for $i<j$. The spectrum
is symmetric with respect to the real axis, i.e., $\omega_{-n}=-\omega_{n}$
and $\omega_{0}=0$. The dispersion relations of $\omega_{n}$ $(n=1,2,3,4)$ for
over-dense and under-dense plasmas are plotted in Fig. 1. Because the spectrum
is also symmetric with respect to the rotation of $\bm{k}$ in the plane
perpendicular to the magnetic field, $\omega_{n}$ is plotted only as functions
of $k_{z}$ and $k_{y}$.

Figure 1: The dispersion relation $\omega_{n}(k_{z},k_{y})$ $(n=1,2,3,4)$ of
cold plasma with immobile ions for (a) overdense and (b) underdense plasmas.
Different colors represent different values of $k_{y}$. Because the spectrum
is also symmetric with respect to the rotation of $\bm{k}$ in the plane
perpendicular to the magnetic field, $k_{x}$ is assumed to be zero. The
dispersion relation $\omega(k_{y},k_{z})$ can be visualized by fixing $k_{z}$
and scanning all $k_{y}$ values.

Topological surface modes might exist at the boundary between two bulk plasmas
that share a common eigen-frequency gap and are in different topological
phases. Consider two adjacent regions of plasmas with density $n_{1}$ and
$n_{2}$, $n_{1}>n_{2}>0$, referred to as regions one and two, and assume that
the interface between two regions lies in the $y$-$z$ plane. It has been
demonstrated that [11] the TLCW exists if and only if the plasma parameters
satisfy the inequality,

$\displaystyle\frac{\omega_{\mathrm{p},1}}{|\Omega|}+\frac{\omega_{\mathrm{p},1}^{2}}{c^{2}k_{z}^{2}}>1>\frac{\omega_{\mathrm{p},2}}{|\Omega|}+\frac{\omega_{\mathrm{p},2}^{2}}{c^{2}k_{z}^{2}},$
(5)

where $\omega_{\mathrm{p},i}$ are the plasma frequencies in each region. Two
regions of plasmas satisfying the inequality above are topologically
different, or equivalently, there is a topological phase transition from one
region to the other.

The condition was established by invoking the principle of bulk-edge
correspondence in condensed matter physics for two bulk materials in different
topological phases. And the topological phases are characterized by the
numerically calculated Chern numbers over the 2D space of perpendicular
wavenumbers [10, 11]. Because the phase transition responsible for this
particular surface mode is due to the resonance between the Langmuir wave and
cyclotron wave, a proper name for it is Topological Langmuir-Cyclotron Waves
(TLCW).

However, unlike the scenarios in condensed matters, the wavenumber space for
plasma waves and many other waves in classical media is topologically
contractible. And the topology of vector bundles over a topologically
contractible base manifold is trivial. In addition, the Atiyah-Patodi-Singer
index theorem [47] of spectral flow, which is the rigorous mathematical
statement for the bulk-edge correspondence, was only proved for the periodic
wavenumber space commonly found in condensed matters. But for plasmas, the
wavenumber space is not periodic in general. Recently, a rigorous analysis
[48] has been developed using the tools of algebraic topology and a spectral
flow index theorem formulated by Faure [49] on non-periodic wavenumber space.
The analysis confirms the existence of TLCW under condition (5).

##  III The isofrequency contours of the TLCW

To quantitatively study the properties of the TLCW, we assume in this section
that the background plasma density $n(\bm{r})$ is non-uniform only in the
$x$-direction. The dispersion relation $\omega(k_{y},k_{z})$ of a homogeneous
bulk can be displayed by two different methods. The first method is to plot
$\omega(k_{y},k_{z})$ as a function of $k_{y}$ and $k_{z}$ as in Fig. 1, which
shows the range of band gaps and the location of the gap-closing. The other
method to exam the level sets of the eigen-frequencies in the space of
wavenumber, i.e., the isofrequency contours, which are more suitable for
visualizing the group velocity. The familiar CMA diagram [43, 44, 45] of bulk
plasma waves can also be established as isofrequency contours, indicating an
interesting connection between the CMA classification and the recently
instituted topological phase classification of plasmas. In this section, we
derive the frequency range where the TLCW exists and analyze the unique
properties of the TLCW by drawing the isofrequency contours for both bulk and
surface plasma waves in a one-dimensional model within the frequency range.

###  III.1 The frequency range of the TLCW

When the TLCW exists at the interface between two regions of different
topological phases, it will fall in the common frequency gap of the bulk waves
shared by the two regions. We denote the frequency of the first and second
branches of bulk waves by $\omega_{1}(k_{y},k_{z})$ and
$\omega_{2}(k_{y},k_{z})$, respectively. In each region, the frequency gap is
between the top of the first branch $\max_{k_{y}}(\omega_{1})$ and the bottom
of the second branch $\min_{k_{y}}(\omega_{2})$ [11]. Using this property, we
obtain the possible frequency gap in each region, which can be written as

$\displaystyle\min_{k_{z}}\left(\max_{k_{y}}(\omega_{1})\right)<\omega<\max_{k_{z}}\left(\min_{k_{y}}(\omega_{2})\right).$
(6)

Therefore, the common range given by the inequality above from the two regions
determines the frequency range of possible TLCW.

In Fig. 1, the resonance between $\omega_{1}$ and $\omega_{2}$ happens at

$\displaystyle
k_{i}^{\pm}:=\frac{\omega_{\mathrm{p},i}/c}{\sqrt{1\pm\omega_{\mathrm{p},i}/|\Omega|}}.$
(7)

Firstly, we consider region one with a higher density $n_{1}$. If region one
is overdense, $|\omega_{\mathrm{p,1}}/\Omega|>1$, the first inequality in
condition (5) holds for all $k_{z}$. From Fig. 1(a),
$\min_{k_{z}}\left(\max_{k_{y}}(\omega_{1})\right)$ is $0$ when $k_{z}=0$ and
$k_{y}=0,$ and $\max_{k_{z}}\left(\min_{k_{y}}(\omega_{2})\right)$ is
$\omega_{\text{p},1}$ at $k_{z}=k_{i}^{+}$ and $k_{y}=0.$ If region one is
underdense, $|\omega_{\mathrm{p,1}}/\Omega|<1$, the first inequality in
condition (5) gives $k_{z}<k_{1}^{-}$. From Fig. 1(b), we can see that the
frequency range is the same as in the overdense case. Thus, for region one,
gap range is $0<\omega<\omega_{\text{p},1}$. Next, for region two with a lower
density $n_{2}$, the second inequality in condition (5) shows that the plasma
has to be underdense and $k_{z}>k_{2}^{-}$. From Fig. 1(b),
$\min_{k_{z}}\left(\max_{k_{y}}(\omega_{1})\right)$ is $\omega_{\mathrm{p,2}}$
at $k_{z}=k_{2}^{-}$ and $k_{y}=0,$ and
$\max_{k_{z}}\left(\min_{k_{y}}(\omega_{2})\right)$ is $|\Omega|$ at
$k_{z}\to+\infty$ and $k_{y}=0.$ Thus, the gap range for region two is
$\omega_{\mathrm{p,2}}<\omega<\Omega$. Combining these two ranges, we infer
that the frequency range of possible TLCW is

$\displaystyle\omega_{\mathrm{p,2}}<\omega<\min\big{(}|\Omega|,\omega_{\mathrm{p,1}}\big{)}\,.$
(8)

###  III.2 The isofrequency contours from 1D model

After determining the frequency range of possible TLCW, we numerically
calculate the isofrequency contours for both the bulk and the surface waves
using a one-dimensional model. We choose the density profile to be continuous,
given by
$n(x)=\frac{1}{2}(n_{1}-n_{2})\\{\tanh[-(x-L)/\delta]+\tanh[(x+L)/\delta]\\}+n_{2}$,
where $n_{1}$ and $n_{2}$ are the plasma densities of region one and two,
respectively, and $L$ and $\delta$ are the location and width of the
interface. The periodic condition is used at $x=\pm 2L$. The equilibrium
configureation is sketched in Fig. 2. To draw the isofrequency contours of the
system, the dispersion relation $\omega(k_{y},k_{z})$ and the eigenmodes are
calculated using the method introduced in Ref. [11].

Figure 2: The equilibrium structure. The dotted red and blue lines are the
boundaries between two regions.

Figure 3 displays the isofrequency contours of the system when the TLCW
exists. We choose cyclotron frequency to be $|\Omega|=1$, and both regions are
underdense with $\omega_{\mathrm{p,1}}=0.8$ and $\omega_{\mathrm{p,2}}=0.3$.
The frequency range of TLCW according to Eq. (8) is $0.3<\omega<0.8$. The
isofrequency contours at frequency $\omega=0.5$ are shown in Fig. 3a for
$k_{z}>0$, and the $k_{z}<0$ part can be obtained from the symmetry condition
$\omega(k_{y},k_{z})=\omega(k_{y},-k_{z})$. The dispersion relations
$\omega(k_{y};k_{z})$ at $k_{z}=0.2,0.8$ and $1.0$ are shown in Fig. 3b-d for
$\omega>0$, and the $\omega<0$ part can be obtained from symmetry condition
$\omega(-k_{y},k_{z})=-\omega(k_{y},k_{z})$. In these plots, colored lines
represent the surface modes or the bulk modes with $k_{x}=0$, while the gray
lines/regions represent other bulk modes with $k_{x}=n\pi/2L$, where $n=\pm
1,\pm 2,\cdots$.

Figure 3: The dispersion relations for the system with
$\omega_{\mathrm{p,1}}=0.8$ and $\omega_{\mathrm{p,2}}=0.3$. The magenta and
green lines represent the isofrequency contours of bulk waves in regions one
and two, assuming $k_{x}=0$. The blue and red lines represent the surface
waves on the left and right boundaries. The gray areas and lines represent the
bulk waves with $k_{x}\neq 0$. (a) The isofrequency contours at frequency
$\omega=0.5$. The dashed lines represent the specific values of $k_{z}$.
(b)-(d) The $\omega(k_{y})$ at $k_{z}=0.3,0.8$ and $1.0$, respectively. The
dashed lines represent the location of $\omega=0.5$. The orange area are the
frequency gap of bulk waves.

For the isofrequency contours in Fig. 3a, the contours (green) for bulk waves
in region two are both ellipsoids, while one of the contours (magenta) for
region one is a hyperboloid. Well-defined surface waves at the left (red) and
right (blue) boundaries can be found approximately around $0.5\leq k_{z}\leq
1.1$, connecting the ellipsoidal and hyperboloidal contours. Outside this
range of $k_{z}$, since the density profile is continuous, the eigenmodes of
the surface waves decay slowly away from the boundaries into the bulk modes. A
closer look can be taken at specific values of $k_{z}$. At $k_{z}=0.3$ from
Figs. 3a and 3b, only bulk modes from region two exist. Near $k_{z}=0.8$,
there is a gap for bulk modes in both Fig. 3a (around $0.58\leq k_{z}\leq
0.95$) and Fig. 3c (around $0.42\leq\omega\leq 0.67$). The surface modes that
exist and fill-up such frequency gap are what we called TLCWs, predicted by
the topological phase transition and the bulk-edge correspondence. At
$k_{z}=1.0$, there exist both surface waves and bulk waves in region one.
These surface waves are not topological surface waves because they are not in
the frequency gap, though they are continuations of the topological surface
waves into the frequency ragne of the bulk waves.

Some physical properties of the TLCW are directly observable from the
isofrequency contours. Firstly, the TLCW is unidirectional. Since the group
velocity $\bm{v}_{\mathrm{g}}=\partial\omega/\partial\bm{k}$ is perpendicular
to the isofrequency contours, the TLCW always propagates along the
$(\hat{z},\hat{y})$ or $(\hat{z},-\hat{y})$ direction on the left or right
boundary, respectively, for $k_{z}>0$, and along the $(-\hat{z},\hat{y})$ or
$(-\hat{z},-\hat{y})$ direction for $k_{z}<0$. Secondly, since the TLCW exists
in a frequency gap, it is immune to back-scattering when the surface is
perturbed, at least when the scale-length of the perturbation is much larger
than the wavelength of the waves. Intuitively, this is because no other bulk
wave exists for the surface waves to couple with in such a frequency gap.
Although two surface waves exist in the frequency gap, they are spatially
separated, i.e., one on the left and the other on the right boundary, and
cannot interact. In section IV, these properties will be further verified by
numerical simulations in both 2D and 3D.

Figure 4: The isofrequency contours and dispersion relation
$\omega(k_{y};k_{z})$ for (a) and (b) $\omega_{\mathrm{p,1}}=0.8$ and
$\omega_{\mathrm{p,2}}=0.3$ , and for (c) and (d) $\omega_{\mathrm{p,1}}=1.1$
and $\omega_{\mathrm{p,2}}=0.9$.

For comparison, two isofrequency contours at two frequencies where no TLCW
exists are shown in Fig. 4. Again, we choose $\omega_{\mathrm{p,1}}=0.8$ and
$\omega_{\mathrm{p,2}}=0.3$, and the frequency range for TLCW is
$0.3<\omega<0.8$, if exists. The isofrequency contours at $\omega=0.85$ is
plotted in Fig. 4a. Here, the isofrequency contours of the bulk waves for both
regions are ellipsoid. There are multiple surface waves when $k_{z}<2.1$, some
of which are the continuation of the TLCW, while others are just surface waves
without topological origin. The dispersion relation $\omega(k_{y})$ at
$k_{z}=1$ is shown in Fig. 4b, which shows that $\omega=0.85$ is not in the
frequency gap of bulk modes. For the case in Fig. 4c,
$\omega_{\mathrm{p,1}}=1.1$ and $\omega_{\mathrm{p,2}}=0.9$, and the frequency
range for possible TLCW is $0.9<\omega<1$. At frequency $\omega=0.75$, the
isofrequency contours in Fig. 4c show that both regions have an ellipsoidal
and a hyperboloidal contours for bulk waves, and there is a bulk gap around
$0.7<k_{z}<2.3$. As expected, there is no surface wave in the gap. At
$k_{z}=1.0$, both Figs. 4c and 4d show that no surface mode exists at this
$k_{z}$.

###  III.3 Relation with the CMA diagram

In the previous numerical calculation, we see that when the TLCW exists (Fig.
3a), the shapes of the isofrequency contours of bulk waves in regions one and
two are different. On the contrary, when no TLCW exists (Figs. 4a and 4c), the
shapes of contours of bulk waves in two regions are the same. Here, we
demonstrate that this phenomenon originates from the relationship between the
topological classification by Chern invariants and the well-known CMA
classification of the plasma waves [43, 44, 45].

In the study of homogeneous cold plasmas, waves at given frequencies can be
classified based on the shapes of the wave normal surfaces
($\bm{v}_{\mathrm{p}}=\omega/c\bm{k}$), or equivalently the shape of the
index-of-refraction surface ($\bm{n}=c\bm{k}/\omega$). We use the second
method here since at a given frequency $\omega$ the shape of the index-of-
refraction surface is the same as the isofrequency contours. Such surfaces
have three possible shapes: ellipsoid, hyperboloid of two sheets, and
hyperboloid of one sheet. The CMA diagram plots the different shapes of the
index-of-refraction surface in the parameter space of the plasma density and
magnetic field strength. A CMA diagram for cold plasmas with immobile ions is
shown in Fig. 5. In different regions of the CMA diagram, the shapes of the
index-of-refraction surface may change. This is similar to the variation of
Fermi surfaces, known as the Lifshitz phase transitions [50] in condensed
matter physics. This transition closely connects with the topological phase
transition discussed in Sec. II. The frequency range of the TLCW in Eq. (8)
can be equivalently written as

$\displaystyle|{\Omega}/{\omega}|>1,\quad{\omega_{\mathrm{p,1}}^{2}}/{\omega^{2}}>1,\quad{\omega_{\mathrm{p,2}}^{2}}/{\omega^{2}}<1.$
(9)

It turns out that on the CMA diagram, the inequalities in (9) also determine
the possible shapes of bulk waves. Region one belongs to the two parts at the
top right, where one bulk wave is hyperboloid, while region two belongs to the
top left, where all bulk waves are ellipsoid. Therefore, we find that the
topological phase transition that generates the TLCW occurs simultaneously
with the transition of the shapes of the index-of-refraction surfaces. In
other words, whenever the TLCW exists, the shapes of the bulk waves on the
isofrequency contours in regions one and two must be different.

Figure 5: CMA diagram for cold plasmas with immobile ions. The shapes of the
index-of-refraction surface are sketched in each region.

##  IV Numerical simulations in 2D and 3D

In this section, we present numerical simulation results of time-dependent
propagation of the TLCW in two- and three-dimensions. We first introduce the
numerical method adopted for time-dependent simulation and then discuss the
simulation results and verify the key physics properties of the TLCW,
including the unidirectional propagation and the immunity to scattering. The
momentum and angular momentum of the waves will also be addressed

###  IV.1 Numerical algorithms

To simulate the dynamics governed by Eqs. (1)-(3), we solve the momentum
equation (1) using the Caylay transformation [51], and solve the Maxwell
equations (2) and (3) on the Yee grid [52]. The 3D space is discretized by a
$N_{x}\times N_{y}\times N_{z}$ grid, where the grid size in three directions
are $\Delta x,\Delta y$ and $\Delta z$. In each direction, we define integer
grid points and half-integer grid points, e.g., $x_{i}=i\Delta x$ and
$x_{i+\frac{1}{2}}=(i+\frac{1}{2})\Delta x$. The time grid has a similar
structure, i.e., $t^{n}=n\Delta t$ and
$t^{n+\frac{1}{2}}=(n+\frac{1}{2})\Delta t$, where $\Delta t$ is the time-step
size. The nine independent variables in Eqs. (1)-(3) are discretized on the
grid as

$\displaystyle\begin{pmatrix}v_{x,i,j,k}^{n+\frac{1}{2}}\\\\[10.0pt]
v_{y,i,j,k}^{n+\frac{1}{2}}\\\\[10.0pt]
v_{z,i,j,k}^{n+\frac{1}{2}}\end{pmatrix},\,\begin{pmatrix}E_{x,i+\frac{1}{2},j,k}^{n}\\\\[10.0pt]
E_{y,i,j+\frac{1}{2},k}^{n}\\\\[10.0pt]
E_{z,i,j,k+\frac{1}{2}}^{n}\end{pmatrix},\,\begin{pmatrix}B_{x,i,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}\\\\[8.0pt]
B_{y,i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}\\\\[8.0pt]
B_{z,i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}\end{pmatrix},$

where superscribe indices are for the time grid, and subscript indices are for
the spatial grid. For variables defined on integer grids point $f_{i}$, denote
finite difference and average between two adjacent grid points as

$\displaystyle\Delta_{x}[f_{i}]:=\dfrac{f_{i+1}-f_{i}}{\Delta x},\quad
f_{\,\overline{i+\frac{1}{2}}}:=\dfrac{f_{i+1}+f_{i}}{2}.$

For variables defined on half-integer grid points $g_{i+\frac{1}{2}}$,
similarly,

$\displaystyle\Delta_{x}\left[g_{i+\frac{1}{2}}\right]:=\dfrac{g_{i+\frac{1}{2}}-g_{i-\frac{1}{2}}}{\Delta
x},\quad g_{\,\bar{i}}:=\dfrac{g_{i+\frac{1}{2}}+g_{i-\frac{1}{2}}}{2}.$

Using these notations, the discretization of Eqs. (1)-(3) is

$\displaystyle\Delta_{t}\left[v_{x,i,j,k}^{n+\frac{1}{2}}\right]$
$\displaystyle=\omega_{\mathrm{p},i,j,k}\,E_{x,\bar{i},j,k}^{n}-\Omega\,v_{y,i,j,k}^{\bar{n}},$
(10) $\displaystyle\Delta_{t}\left[v_{y,i,j,k}^{n+\frac{1}{2}}\right]$
$\displaystyle=\omega_{\mathrm{p},i,j,k}\,E_{y,i,\bar{j},k}^{n}+\Omega\,v_{x,i,j,k}^{\bar{n}},$
$\displaystyle\Delta_{t}\left[v_{z,i,j,k}^{n+\frac{1}{2}}\right]$
$\displaystyle=\omega_{\mathrm{p},i,j,k}\,E_{z,i,j,\bar{k}}^{n}.$
$\displaystyle\Delta_{t}\left[E_{x,i+\frac{1}{2},j,k}^{n}\right]$
$\displaystyle=\Delta_{y}\left[B_{z,i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}\right]-\Delta_{z}\left[B_{y,i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}\right]-\left(\omega_{\mathrm{p}}v_{x}^{n+\frac{1}{2}}\right)_{\overline{i+\frac{1}{2}},j,k},$
(11) $\displaystyle\Delta_{t}\left[E_{y,i,j+\frac{1}{2},k}^{n}\right]$
$\displaystyle=\Delta_{z}\left[B_{x,i,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}\right]-\Delta_{x}\left[B_{z,i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}\right]-\left(\omega_{\mathrm{p}}v_{y}^{n+\frac{1}{2}}\right)_{i,\overline{j+\frac{1}{2}},k},$
$\displaystyle\Delta_{t}\left[E_{z,i,j,k+\frac{1}{2}}^{n}\right]$
$\displaystyle=\Delta_{x}\left[B_{y,i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}\right]-\Delta_{y}\left[B_{x,i,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}\right]-\left(\omega_{\mathrm{p}}v_{z}^{n+\frac{1}{2}}\right)_{i,j,\overline{k+\frac{1}{2}}}.$
$\displaystyle\Delta_{t}\left[B_{x,i,j+\frac{1}{2},k+\frac{1}{2}}^{n+\frac{1}{2}}\right]$
$\displaystyle=\Delta_{z}\left[E_{y,i,j+\frac{1}{2},k}^{n}\right]-\Delta_{y}\left[E_{z,i,j,k+\frac{1}{2}}^{n}\right],$
(12)
$\displaystyle\Delta_{t}\left[B_{y,i+\frac{1}{2},j,k+\frac{1}{2}}^{n+\frac{1}{2}}\right]$
$\displaystyle=\Delta_{x}\left[E_{z,i,j,k+\frac{1}{2}}^{n}\right]-\Delta_{z}\left[E_{x,i+\frac{1}{2},j,k}^{n}\right],$
$\displaystyle\Delta_{t}\left[B_{z,i+\frac{1}{2},j+\frac{1}{2},k}^{n+\frac{1}{2}}\right]$
$\displaystyle=\Delta_{y}\left[E_{x,i+\frac{1}{2},j,k}^{n}\right]-\Delta_{x}\left[E_{y,i,j+\frac{1}{2},k}^{n}\right].$

Here, $\omega_{\mathrm{p}}(\bm{r})$ is a prescribed function,
$\omega_{\mathrm{p},i,j,k}:=\omega_{\mathrm{p}}(x_{i},y_{j},z_{k})$,
${\displaystyle\left(\omega_{\mathrm{p}}v_{x}\right)_{\overline{i+\frac{1}{2}}}}:=(\omega_{\mathrm{p},i}v_{x,i}+\omega_{\mathrm{p},i+1}v_{x,i+1})/2$.
Notice that although its right-hand side depends on $v^{n+\frac{1}{2}}$
through $v^{\bar{n}}$, Eq. (10) is explicitly solvable by the Cayley
transformation that appears in the Boris algorithm [51] and other structure-
preserving algorithm for charged particle dynamics [53, 54, 55, 56, 57].

The algorithm in 2D is a special case of the 3D algorithm described above.
When the plasma density is invariant along the background magnetic field,
$\omega_{\mathrm{p}}(\bm{r})=\omega_{\mathrm{p}}(x,y)$, we can replace the $z$
dependence of all variables by $\exp(\mathrm{i}k_{z}z)$. Then, all variables
in Eqs. (10)-(12) no longer depend on the index $k$, and a constant
$\mathrm{i}k_{z}$ replaces the difference operators $\Delta_{z}[\cdot]$.

###  IV.2 Propagation of the TLCW in 2D and 3D

We now apply the algorithm to simulate the TLCW excited by a localized time-
dependent source in two- and three-dimensions. An external force that
represents the source is added to the RHS of the momentum equation (10). In
the 3D simulation, the force we used is

$\displaystyle\bm{F}(\bm{r},t)\sim\exp\left(-\dfrac{|\bm{r}-\bm{r}_{s}|^{2}}{\delta^{2}}\right)e^{\mathrm{i}(k_{s}z-\omega_{s}t)},$
(13)

where $\bm{r}_{s}$ and $\delta$ are the spatial center and width of the
source, $k_{s}$ and $\omega_{s}$ are the wavenumber and frequency of the
source. The direction of the force does not affect the excitation of the
surface waves. In the 2D simulation, the force does not depend on $z$, i.e.,
$k_{s}=0,$ and the factor $\mathrm{i}k_{z}$ replaces the difference operators
$\Delta_{z}[\cdot]$ in Eqs. (11) and (12).

Figure 6: Propagation of TLCW in (a) 2D and (b) 3D, where
$\omega_{\mathrm{p,1}}=0.8$, $\omega_{\mathrm{p,2}}=0.3$, $k_{s}=0.8$, and
$\omega_{s}=0.5$. In (a), the color map indicates the strength of real part of
$E_{z}$, the yellow star marks the location of the source, and the dashed line
is the boundary between two regions. In (b), the red and blue contours
represent the locations of $\mathrm{Re}[E_{z}]=\pm
0.15\,\mathrm{max}[\mathrm{Re}(E_{z})]$, the yellow sphere marks the source,
and the gray surface is the boundary between two regions.

The parameters for the first case simulated is $\omega_{\mathrm{p,1}}=0.8$,
$\omega_{\mathrm{p,2}}=0.3$, $k_{s}=0.8$, and $\omega_{s}=0.5$. The dispersion
relation calculated in 1D was shown in Figs. 3a and 3c. The center of the
source is on the boundary. The time-dependent propagation of the TLCW in 2D
and 3D are shown in Fig. 6. Notice that the configuration is equivalent to the
right boundary shown in Fig. 2, where the blue lines in Fig. 3 represent the
surface waves in Fig. 3. In the 2D simulation displayed in Fig. 6a, the TLCW
is excited, and the excitation of bulk waves is negligible. Furthermore, the
surface wave only travels to the left side of the source, confirming its
unidirectional propagation. The same phenomenons can also be observed in 3D
simulation shown in Fig. 6b, where the TLCW is excited and propagates towards
the $(-\hat{y},+\hat{z})$ direction.

###  IV.3 The momentum of the surface waves

In this subsection, we discuss the momentum and angular momentum carried by
the TLCW. There exist two different definitions of momentum for plasma waves
[58, 59], i.e., the (Minkowski) canonical momentum and the (Abraham) kinetic
momentum. The canonical momentum of the TLCW is proportional to the wavenumber
$\bm{k}$, and it can vary significantly, as illustrated in Fig. 3. However,
the kinetic momentum, which is proportional to it group velocity
$\bm{v}_{\mathrm{g}}$, is unidirectional. For example, the surface waves on
the right boundary can only have kinetic momentum in the
$(-\hat{y},\pm\hat{z})$ direction.

Figure 7: (a) The Poynting vectors in the 3D simulation. The yellow sphere
marks the source, the arrows indicate the Poynting vectors, the gray surface
is the boundary. (b) The Poynting vectors at the boundary between two regions
(at $x=30$). The yellow star marks the source, the red arrows indicate the
Poynting vectors in the $y$-$z$ plane, whose magnitudes are given by the color
of the background.

For the cold plasma model given by Eqs. (1)-(3), the dielectric tensor is
Hermitian and independent of $\bm{k}$. Therefore, the kinetic momentum, the
group velocity, and the flux of canonical energy are all proportional to the
Poynting vectors [46] $\bm{S}\sim\mathrm{Re}[\bm{E}\times\bm{B}^{*}]$. The
Poynting vectors of the TLCW excited on a planar boundary is shown in Fig. 7.
As expected, the kinetic momentum is along the direction of wave propagation

###  IV.4 Propagation of the TLCW on non-planar boundaries

The main properties of the TLCW is robust in more complex configurations. We
now show some simulation results of the TLCW on non-planar boundaries in both
two- and three-dimensions. Shown in Fig. 8 is the TLCW propagation on a
zigzagged boundary, which clearly demonstrates immunity of back-scattering and
unidirectional propagation. When the TLCW meets sharp turns on the boundary,
it propagates along the boundary without reflection or transition to bulk
waves.

Figure 8: (a) 2D and (b) 3D simulations of the TLCW excited at a zigzagged
boundary.

Additional examples are shown in Figs. 9 and 10, where the boundary between
two regions is a square or a circle. In these configurations, the left and
right boundaries are connected and indistinguishable. As anticipated, the TLCW
propagates along the boundaries, regardless their shapes, in a unidirectional
manner and without any scattering. When viewed against the direction of the
magnetic field, the TLCW propagates clockwise if
$\omega_{\mathrm{p,1}}>\omega_{\mathrm{p,2}}$, and therefore carries a non-
zero kinetic angular momentum proportional to
$\int\bm{r}\times\bm{S}\,\mathrm{d}V$. Notice that the source exciting the
TLCW does not carry any angular momentum. The generation of an angular-
momentum-carrying surface wave reflects the topological property of the bulk
material and the bulk-edge correspondence. Because of these desirable
properties, the TLCW could be explored as an effective mechanism for driving
current and flow in magnetized plasmas.

Figure 9: (a) 2D and (b) 3D simulations of the TLCW excited on a square
boundary.

Figure 10: (a) 2D and (b) 3D simulations of the TLCW excited on a circular
boundary.

##  V Conclusions and discussion

We have studied the dispersion relation and propagation of the recently
identified Topological Langmuir-Cyclotron Wave (TLCW). From the excitation
condition (5), we derived the frequency range in Eq. (8) where the TLCW can be
observed. It is shown that the topological phase transition responsible for
the TLCW excitation occurs at the resonance between the Langmuir wave and
cyclotron wave, and concurrently with the transition of the shapes of the
index-of-refraction surfaces in the familiar CMA diagram.

Based on the isofrequency contours, the group velocity of the TLCW is
unidirectional, and so is its kinetic momentum. Furthermore, the TLCW is
immune to back-scattering because they lie in the frequency gap of bulk waves.
The unidirectional propagation and immunity of back-scattering are verified
using time-dependent simulations in two- and three-dimensions. The excitation
of the TLCW could be an effective mechanism for driving current and flow in
magnetized plasmas.

In the present study, we theoretically investigated the dispersion and
propagation of the TLCW using a linearized cold plasma model. For experimental
observation, realistic effects need to be carefully considered to better
understand the properties and applications of the TLCW. For example, the
kinetic effects of wave-particle interaction will play an essential role in
the deposition of momentum carried by the TLCW. Other non-Hermitian [60] and
non-linear [61, 62, 60] effects may also be important for topological waves in
more complex fluid and plasma models [63, 64, 65, 33].

###### Acknowledgements.

This research was supported by the U.S. Department of Energy (DE-
AC02-09CH11466).

## References

  * Trivelpiece and Gould [1959] A. Trivelpiece and R. Gould, Space charge waves in cylindrical plasma columns, Journal of Applied Physics 30, 1784 (1959).
  * Gradov and Stenflo [1983] O. Gradov and L. Stenflo, Linear theory of a cold bounded plasma, Physics Reports 94, 111 (1983).
  * Vladimirov  _et al._ [1994] S. V. Vladimirov, M. Yu, and V. Tsytovich, Recent advances in the theory of nonlinear surface waves, Physics Reports 241, 1 (1994).
  * Moisan  _et al._ [1979] M. Moisan, Z. Zakrzewski, and R. Pantel, The theory and characteristics of an efficient surface wave launcher (surfatron) producing long plasma columns, Journal of Physics D: Applied Physics 12, 219 (1979).
  * Ida and Hayashi [1971] Y. Ida and K. Hayashi, Analytical and experimental study on landau damping of surface wave along a plasma column, Journal of Applied Physics 42, 2423 (1971).
  * Lopez-Rios and Vuye [1979] T. Lopez-Rios and G. Vuye, Use of surface plasmon excitation for determination of the thickness and optical constants of very thin surface layers, Surface Science 81, 529 (1979).
  * Agranovich [2012] V. M. Agranovich,  _Surface polaritons_ (Elsevier, 2012).
  * Kupersztych  _et al._ [2004] J. Kupersztych, M. Raynaud, and C. Riconda, Electron acceleration by surface plasma waves in the interaction between femtosecond laser pulses and sharp-edged overdense plasmas, Physics of Plasmas 11, 1669 (2004).
  * Ionson [1978] J. A. Ionson, Resonant absorption of alfvénic surface waves and the heating of solar coronal loops, The Astrophysical Journal 226, 650 (1978).
  * Parker  _et al._ [2020a] J. B. Parker, J. Marston, S. M. Tobias, and Z. Zhu, Topological gaseous plasmon polariton in realistic plasma, Physical Review Letters 124, 195001 (2020a).
  * Fu and Qin [2021] Y. Fu and H. Qin, Topological phases and bulk-edge correspondence of magnetized cold plasmas, Nature Communications 12, 3924 (2021).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized hall conductance in a two-dimensional periodic potential, Physical Review Letters 49, 405 (1982).
  * Halperin [1982] B. I. Halperin, Quantized hall conductance, current-carrying edge states, and the existence of extended states in a two-dimensional disordered potential, Physical Review B 25, 2185 (1982).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: topological insulators, Reviews of Modern Physics 82, 3045 (2010).
  * Bernevig [2013] B. A. Bernevig,  _Topological insulators and topological superconductors_ (Princeton university press, 2013).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Reviews of Modern Physics 83, 1057 (2011).
  * Armitage  _et al._ [2018] N. Armitage, E. Mele, and A. Vishwanath, Weyl and dirac semimetals in three-dimensional solids, Reviews of Modern Physics 90, 015001 (2018).
  * Haldane and Raghu [2008] F. Haldane and S. Raghu, Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry, Physical Review Letters 100, 013904 (2008).
  * Raghu and Haldane [2008] S. Raghu and F. D. M. Haldane, Analogs of quantum-hall-effect edge states in photonic crystals, Physical Review A 78, 033834 (2008).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, _et al._ , Topological photonics, Reviews of Modern Physics 91, 015006 (2019).
  * Yang  _et al._ [2015] Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, Topological acoustics, Physical Review Letters 114, 114301 (2015).
  * Wang  _et al._ [2015] P. Wang, L. Lu, and K. Bertoldi, Topological phononic crystals with one-way elastic edge waves, Physical Review Letters 115, 104302 (2015).
  * He  _et al._ [2016a] C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, Acoustic topological insulator and robust one-way sound transport, Nature Physics 12, 1124 (2016a).
  * Kane and Lubensky [2014] C. Kane and T. Lubensky, Topological boundary modes in isostatic lattices, Nature Physics 10, 39 (2014).
  * Nash  _et al._ [2015] L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner, and W. T. Irvine, Topological mechanics of gyroscopic metamaterials, Proceedings of the National Academy of Sciences 112, 14495 (2015).
  * Huber [2016] S. D. Huber, Topological mechanics, Nature Physics 12, 621 (2016).
  * Delplace  _et al._ [2017] P. Delplace, J. Marston, and A. Venaille, Topological origin of equatorial waves, Science 358, 1075 (2017).
  * Perrot  _et al._ [2019] M. Perrot, P. Delplace, and A. Venaille, Topological transition in stratified fluids, Nature Physics 15, 781 (2019).
  * Tauber  _et al._ [2019] C. Tauber, P. Delplace, and A. Venaille, A bulk-interface correspondence for equatorial waves, Journal of Fluid Mechanics 868 (2019).
  * Souslov  _et al._ [2019] A. Souslov, K. Dasbiswas, M. Fruchart, S. Vaikuntanathan, and V. Vitelli, Topological waves in fluids with odd viscosity, Physical Review Letters 122, 128001 (2019).
  * Marciani and Delplace [2020] M. Marciani and P. Delplace, Chiral maxwell waves in continuous media from berry monopoles, Physical Review A 101, 023827 (2020).
  * Venaille and Delplace [2021] A. Venaille and P. Delplace, Wave topology brought to the coast, Physical Review Research 3, 043002 (2021).
  * Zhu  _et al._ [2021] Z. Zhu, C. Li, and J. Marston, Topology of rotating stratified fluids with and without background shear flow, arXiv preprint arXiv:2112.04691 (2021).
  * Silveirinha [2015] M. G. Silveirinha, Chern invariants for continuous media, Physical Review B 92, 125153 (2015).
  * Gangaraj  _et al._ [2017] S. A. H. Gangaraj, M. G. Silveirinha, and G. W. Hanson, Berry phase, berry connection, and chern number for a continuum bianisotropic material from a classical electromagnetics perspective, IEEE journal on multiscale and multiphysics computational techniques 2, 3 (2017).
  * Silveirinha [2016] M. G. Silveirinha, Bulk-edge correspondence for topological photonic continua, Physical Review B 94, 205105 (2016).
  * Parker [2021a] J. B. Parker, Topological phase in plasma physics, Journal of Plasma Physics 87 (2021a).
  * Gao  _et al._ [2016] W. Gao, B. Yang, M. Lawrence, F. Fang, B. Béri, and S. Zhang, Photonic weyl degeneracies in magnetized plasma, Nature Communications 7, 12435 (2016).
  * Yang  _et al._ [2016] B. Yang, M. Lawrence, W. Gao, Q. Guo, and S. Zhang, One-way helical electromagnetic wave propagation supported by magnetized plasma, Scientific Reports 6, 1 (2016).
  * Parker  _et al._ [2020b] J. B. Parker, J. Burby, J. Marston, and S. M. Tobias, Nontrivial topology in the continuous spectrum of a magnetized plasma, Physical Review Research 2, 033425 (2020b).
  * Parker [2021b] J. Parker, Berry curvature and topological phase: From cold to hot plasmas, Bulletin of the American Physical Society (2021b).
  * Shvets  _et al._ [2021] G. Shvets, R. S. Rajawat, and V. Khudik, Continuum damping of topologically-protected edge modes at the boundary of magnetized plasma, Bulletin of the American Physical Society 66 (2021).
  * Clemmow and Mullaly [1955] P. C. Clemmow and R. F. Mullaly, The dependence of the refractive index in magneto-ionic theory on the direction of the wave normal, in  _Physics of the Ionosphere: Rep. Phys. Soc. Conf._ (London Physical Society, 1955) p. 340.
  * Allis [1959] W. P. Allis, Waves in a plasma, in  _Sherwood Conf. Contr. Fusion, Gatlinburg, TN, Apr. 27-28_ (1959) p. 32.
  * Allis  _et al._ [1963] W. P. Allis, S. J. Buchsbaum, and A. Bers,  _Waves in anisotropic plasmas_ (MIT Press, Cambridge, 1963).
  * Stix [1992] T. H. Stix,  _Waves in plasmas_ (Springer Science & Business Media, 1992).
  * Atiyah  _et al._ [1976] M. F. Atiyah, V. K. Patodi, and I. M. Singer, Spectral asymmetry and riemannian geometry. III, Mathematical Proceedings of the Cambridge Philosophical Society 79, 71 (1976).
  * [48] H. Qin and Y. Fu, Topological Langmuir-cyclotron waves, to be published.
  * Faure [2019] F. Faure, Manifestation of the topological index formula in quantum waves and geophysical waves (2019), arXiv:1901.10592 [math-ph] .
  * Volovik [2017] G. Volovik, Topological Lifshitz transitions, Low Temperature Physics 43, 47 (2017).
  * Qin  _et al._ [2013] H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M. Tang, Why is boris algorithm so good?, Physics of Plasmas 20, 084503 (2013).
  * Yee [1966] K. Yee, Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media, IEEE Transactions on antennas and propagation 14, 302 (1966).
  * He  _et al._ [2015] Y. He, Y. Sun, J. Liu, and H. Qin, Volume-preserving algorithms for charged particle dynamics, Journal of Computational Physics 281, 135 (2015).
  * He  _et al._ [2016b] Y. He, Y. Sun, J. Liu, and H. Qin, Higher order volume-preserving schemes for charged particle dynamics, Journal of Computational Physics 305, 172 (2016b).
  * He  _et al._ [2016c] Y. He, Y. Sun, R. Zhang, Y. Wang, J. Liu, and H. Qin, High order volume-preserving algorithms for relativistic charged particles in general electromagnetic fields, Physics of Plasmas 23, 092109 (2016c).
  * Zhang  _et al._ [2015] R. Zhang, J. Liu, H. Qin, Y. Wang, Y. He, and Y. Sun, Volume-preserving algorithm for secular relativistic dynamics of charged particles, Physics of Plasmas 22, 044501 (2015).
  * Fu  _et al._ [2022] Y. Fu, X. Zhang, and H. Qin, An explicitly solvable energy-conserving algorithm for pitch-angle scattering in magnetized plasmas, Journal of Computational Physics 449, 110767 (2022).
  * Barnett [2010] S. M. Barnett, Resolution of the abraham-minkowski dilemma, Physical Review Letters 104, 070401 (2010).
  * Dodin and Fisch [2012] I. Dodin and N. Fisch, Axiomatic geometrical optics, abraham-minkowski controversy, and photon properties derived classically, Physical Review A 86, 053834 (2012).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-hermitian systems, Reviews of Modern Physics 93, 015005 (2021).
  * Smirnova  _et al._ [2020] D. Smirnova, D. Leykam, Y. Chong, and Y. Kivshar, Nonlinear topological photonics, Applied Physics Reviews 7, 021306 (2020).
  * Zhou  _et al._ [2020] D. Zhou, D. Rocklin, M. Leamy, and Y. Yao, Topological invariant and anomalous edge states of fully nonlinear systems, arXiv preprint arXiv:2012.14103 (2020).
  * Qin  _et al._ [2019] H. Qin, R. Zhang, A. S. Glasser, and J. Xiao, Kelvin-helmholtz instability is the result of parity-time symmetry breaking, Physics of Plasmas 26, 032102 (2019).
  * Fu and Qin [2020] Y. Fu and H. Qin, The physics of spontaneous parity-time symmetry breaking in the kelvin–helmholtz instability, New Journal of Physics 22, 083040 (2020).
  * Qin  _et al._ [2021] H. Qin, Y. Fu, A. S. Glasser, and A. Yahalom, Spontaneous and explicit parity-time-symmetry breaking in drift-wave instabilities, Physical Review E 104, 015215 (2021).