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

Abstract

The 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 two- and three-dimensional 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 Clemmow-Mullaly-Allis (CMA) diagram of plasma waves.

1. Introduction

Surface waves in plasmas have been extensively studied (Trivelpiece & Gould 1959; Gradov & Stenflo 1983; Vladimirov, Yu & Tsytovich 1994) due to their wide range of applications (Ida & Hayashi 1971; Ionson 1978; Lopez-Rios & Vuye 1979; Moisan, Zakrzewski & Pantel 1979; Kupersztych, Raynaud & Riconda 2004; Agranovich 2012). Recently, a new class of surface modes in plasma produced by topological phase transition has been investigated (Parker et al. 2020b; Fu & Qin 2021). Here, we report on studies of the dispersion and propagation of the newly identified topological surface mode called the 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 (Thouless et al. 1982), it is found that the bulk states in 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 (Halperin 1982) 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 (Hasan & Kane 2010). The bulk topology and its relation to surface states have been widely studied not only in quantum systems (Hasan & Kane 2010; Qi & Zhang 2011; Bernevig 2013; Armitage, Mele & Vishwanath 2018), but also in classical systems, such as those in photonics (Haldane & Raghu 2008; Raghu & Haldane 2008; Ozawa et al. 2019; Marciani & Delplace 2020), acoustics (Wang, Lu & Bertoldi 2015; Yang et al. 2015; He et al. 2016a), mechanics (Kane & Lubensky 2014; Nash et al. 2015; Huber 2016) and fluid dynamics (Delplace, Marston & Venaille 2017; Perrot, Delplace & Venaille 2019; Souslov et al. 2019; Tauber, Delplace & Venaille 2019; Venaille & Delplace 2021; Zhu, Li & Marston 2021).

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 two-dimensional (2-D) plasmas were calculated (Silveirinha 2015; Gangaraj, Silveirinha & Hanson 2017) and related to surface waves via the bulk–edge correspondence (Silveirinha 2016). When $k_{z}\neq 0$, the non-trivial topology of the plasma as a 2-D system was found (Parker et al. 2020b; Parker 2021b) and its topological phase diagram was established (Fu & Qin 2021). In particular, under-dense plasmas can be regarded as a Weyl semimetal (Gao et al. 2016; Yang et al. 2016) and host a surface states with the Fermi-arc structure (Armitage et al. 2018). The topological matter properties of ideal magnetohydrodynamics (MHD) plasmas with magnetic shear (Parker et al. 2020a), kinetic plasmas (Parker 2021a) and surface waves within the continuous spectrum (Rajawat, Khudik & Shvets 2022) were also investigated.

The topological surface wave in the present study is called the 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 (Fu & Qin 2021). 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 Clemmow-Mullaly-Allis (CMA) diagram (Clemmow & Mullaly 1955; Allis 1959; Allis, Buchsbaum & Bers 1963) of plasma waves. Two- and three-dimensional time-dependent simulations have been performed for the linearized fluid equations. The TLCW is 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 be explored as a robust mechanism to inject energy and momentum into magnetized plasmas.

The paper is organized as follows. Section 2 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 § 3. The connection between the isofrequency contours and the CMA diagram is addressed. Section 4 describes the numerical algorithms and 2-D and 3-D simulation results of the TLCW.

2. Theoretical model

This section describes the governing equations for the 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. $\boldsymbol {B}_{0}=B_{0}\hat {z}$. Since there is no pressure in cold plasmas, any given plasma density profile $n(\boldsymbol {r})$ is an equilibrium. With proper renormalization, the linearized fluid equations can be written as (Stix 1992; Fu & Qin 2021)

(2.1)\begin{gather} \partial_{t}\boldsymbol{v}={-}\omega_{\mathrm{p}}\boldsymbol{E}-\varOmega\boldsymbol{v}\times\hat{z}, \end{gather}

(2.2)\begin{gather}\partial_{t}\boldsymbol{E}=\boldsymbol{\nabla}\times\boldsymbol{B}+\omega_{\mathrm{p}}\boldsymbol{v}, \end{gather}

(2.3)\begin{gather}\partial_{t}\boldsymbol{B}={-}\boldsymbol{\nabla}\times\boldsymbol{E}, \end{gather}

where $\boldsymbol {v},\boldsymbol {E},\boldsymbol {B}$ are perturbed velocity, electric and magnetic fields, $\omega _{\mathrm {p}}(\boldsymbol {r})=\sqrt {n(\boldsymbol {r})\,{\rm e}^{2}/\epsilon _{0}m_{\mathrm {e}}}$ is the local plasma frequency, $\varOmega =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 $\boldsymbol {k}$, (2.1)–(2.3) reduce to

(2.4)\begin{equation} H|\psi\rangle=\omega|\psi\rangle, \end{equation}

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

(2.5)\begin{equation} H(\omega_{\text{p}},\varOmega,\boldsymbol{k}) =\begin{pmatrix} -\mathrm{i}\varOmega\hat{z}\times & -\mathrm{i}\omega_{\text{p}} & 0\\ \mathrm{i}\omega_{\text{p}} & 0 & -\boldsymbol{k}\times\\ 0 & \boldsymbol{k}\times & 0 \end{pmatrix}\end{equation}

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

(2.6)\begin{equation} \omega_{n},|\psi_{n}\rangle,\quad n={-}4,-3,\ldots,3,4. \end{equation}

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 figure 1. Because the spectrum is also symmetric with respect to the rotation of $\boldsymbol {k}$ in the plane perpendicular to the magnetic field, $\omega _{n}$ is plotted only as functions of $k_{z}$ and $k_{y}$. The resonance between $\omega _{1}$ and $\omega _{2}$ happens at

(2.7)\begin{equation} k^{{\pm}}(\omega_{\mathrm{p}},\varOmega):= \frac{\omega_{\mathrm{p}}/c}{\sqrt{1\pm\omega_{\mathrm{p}}/|\varOmega|}}.\end{equation}

Figure 1. The dispersion relation $\omega _{n}(k_{z},k_{y})$ $(n=1,2,3,4)$ of cold plasma with immobile ions for (a) an over-dense plasma, i.e. $|\omega _{\mathrm {p}}/\varOmega \!>\!\!1|$, and (b) an under-dense plasma, i.e. $|\omega _{\mathrm {p}}/\varOmega \!<\!\!1|$. Here, $k^{\pm }$ are the crossing points of $\omega _{1}$ and $\omega _{2}$, defined in (2.7). Different colours represent different values of $k_{y}$. Because the spectrum is also symmetric with respect to the rotation of $\boldsymbol {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 eigenfrequency gap and are in different topological phases. Consider two adjacent regions of plasma 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 the two regions lies in the $y$–$z$ plane. It has been demonstrated that (Fu & Qin 2021) the TLCW exists if and only if the plasma parameters satisfy the inequality

(2.8)\begin{equation} \frac{\omega_{\mathrm{p},1}}{|\varOmega|}+\frac{\omega_{\mathrm{p},1}^{2}}{c^{2}k_{z}^{2}}>1> \frac{\omega_{\mathrm{p},2}}{|\varOmega|}+\frac{\omega_{\mathrm{p},2}^{2}}{c^{2}k_{z}^{2}},\end{equation}

where $\omega _{\mathrm {p},i}$ are the plasma frequencies in each region. Two regions of plasma 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. The topological phases are characterized by the numerically calculated Chern numbers over the 2-D space of perpendicular wavenumbers (Parker et al. 2020b; Fu & Qin 2021). 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 the TLCW.

However, unlike the scenarios in condensed matters, where the wavenumber space is periodic, the wavenumber space for plasma waves and many other waves in classical media is topologically contractible. It is known that the topology of vector bundles over a topologically contractible base manifold is trivial. In addition, the Atiyah–Patodi–Singer index theorem (Atiyah, Patodi & Singer 1976) of spectral flow, which is a rigorous mathematical statement for the bulk–edge correspondence, was only proved for the periodic wavenumber space. But for plasmas, the wavenumber space is not periodic in general. Recently, a rigorous analysis (Qin & Fu 2022) has been developed using the tools of algebraic topology and a spectral flow index theorem formulated by Faure (Faure 2019) over an $\mathbb {R}$-valued wavenumber space. The analysis confirms the existence of TLCW under condition (2.8).

3. 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(\boldsymbol {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 figure 1, which shows the range of band gaps and the location of the gap closing. The other method is to examine the level sets of the eigenfrequencies in wavenumber space, i.e. the isofrequency contours, which are more suitable for visualizing the group velocity. The familiar CMA diagram (Clemmow & Mullaly 1955; Allis 1959; Allis et al. 1963) 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 analyse the unique properties of the TLCW by drawing the isofrequency contours for both bulk and surface plasma waves in a 1-D model within the frequency range.

3.1. The frequency range of the TLCW

When the TLCW exists at the interface between two regions of different topological phases, it will fall into the common frequency gap of the bulk waves shared by the two regions. At the interface between two regions of different topological phases, a topological edge mode, also known as spectral flow (Faure 2019; Delplace 2022; Qin & Fu 2022), exists and it transits between the spectrum corresponding to the two regions. Specifically, in this paper, we define TLCW to be the topological edge mode in the frequency gap between the bulk spectrum. As shown later, mode within the bulk frequency gap is immune to backscattering. Strictly speaking, frequencies of topological edge modes can overlap with the bulk spectrum, which is evident from the numerical results shown in figure 4. This is expected since the topological mode connects the gaped bands of the bulk spectrum. But to simplify the discussion with respect to the topological robustness, in the present context we choose to define the TLCW in a narrow sense as the topological mode whose frequency is in the gap of the bulk bands. This restriction simplifies the discussion with respect to the topological robustness. We denote the frequencies 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})$ (Fu & Qin 2021). Using this property, we obtain the possible frequency gap in each region, which can be written as

(3.1)\begin{equation} \min_{k_{z}}\left(\max_{k_{y}}(\omega_{1})\right)<\omega<\max_{k_{z}} \left(\min_{k_{y}}(\omega_{2})\right).\end{equation}

Therefore, the common range given by the inequality above from the two regions determines the frequency range of possible TLCW. In other words, we need to solve inequality (3.1) when the TLCW exists, i.e. under the constraint of condition (2.8). For simplicity, we define $k_{i}^{\pm }\equiv k^{\pm }(\omega _{\mathrm {p},i},\varOmega )$ for $i=1,2$.

Firstly, consider region one with a higher density $n_{1}$. The wavenumber $k_{z}$ should satisfy the first inequality in condition (2.8). If region one is over-dense, i.e. $|\omega _{\mathrm {p,1}}/\varOmega |>1$, the first inequality in condition (2.8) holds for all $k_{z}$. From figure 1(a), we can see that

(3.2a,b)\begin{equation} \min_{k_{z}\in\mathbb{R}}\left(\max_{k_{y}\in\mathbb{R}}(\omega_{1})\right)= \omega_{1}(0,0)=0,\quad\max_{k_{z}\in\mathbb{R}}\left(\min_{k_{y}\in\mathbb{R}}(\omega_{2})\right)= \omega_{2}(0,\infty)=\omega_{\mathrm{p},1}.\end{equation}

If region one is under-dense, i.e. $|\omega _{\mathrm {p,1}}/\varOmega |<1$, the first inequality in condition (2.8) gives $|k_{z}|< k_{1}^{-}$. Then from figure 1(b) it is clear that

(3.3a,b)\begin{equation} \min_{|k_{z}|< k_{1}^{-}}\left(\max_{k_{y}\in\mathbb{R}}(\omega_{1})\right)= \omega_{1}(0,0)=0,\quad\max_{|k_{z}|< k_{1}^{-}}\left(\min_{k_{y}\in\mathbb{R}}(\omega_{2})\right)= \omega_{2}(0,k_{1}^{-})=\omega_{\text{p},1}.\end{equation}

Thus, for region one, the gap range is $0<\omega <\omega _{\text {p},1}$.

Next, for region two with a lower density $n_{2}$, the second inequality in condition (2.8) shows that the plasma has to be under-dense, i.e. $|\omega _{\mathrm {p,2}}/\varOmega |<1$, and $|k_{z}|>k_{2}^{-}$. From figure 1(b), we have

(3.4a,b)\begin{equation} \min_{|k_{z}|>k_{2}^{-}}\left(\max_{k_{y}\in\mathbb{R}}(\omega_{1})\right)= \omega_{1}(0,k_{2}^{-})=\omega_{\mathrm{p,2}},\quad\max_{|k_{z}|>k_{2}^{-}} \left(\min_{k_{y}\in\mathbb{R}}(\omega_{2})\right)=\omega_{2}(0,\infty)=|\varOmega|.\end{equation}

Thus, the gap range for region two is $\omega _{\mathrm {p,2}}<\omega <|\varOmega |$. Combining these two ranges, we infer that the frequency range of possible TLCW is

(3.5)\begin{equation} \omega_{\mathrm{p,2}}<\omega<\min(|\varOmega|,\omega_{\mathrm{p,1}}).\end{equation}

3.2. Numerical demonstration of the frequency range

After deriving the frequency range of TLCW in (3.5), we now numerically demonstrate it using a 1-D 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 regions one and two, respectively, and $L$ and $\delta$ are the location and width of the interface. The periodic boundary condition is used at $x=\pm 2L$. The equilibrium configuration is sketched in figure 2. The width $\delta$ will only change the spectrum of bulk modes, but not affect the existence of TLCW, see (Qin & Fu 2022) for a detailed discussion. In the present study, we used $L=40$ and $\delta =3$.

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

Since all the extreme values of $\omega _{1,2}$ are reached at $k_{y}=0$ in (3.2a,b)–(3.4a,b), we show the numerically calculated dispersion relation of the non-uniform system $\omega (k_{z})$ at $k_{y}=0$. Two different cases are shown in figure 3, which are adapted from Fu & Qin (2021). The TLCWs are shown by the red curves, while the grey curves show all other modes. The green and magenta curves represent the bulk modes when $k_{x}=k_{y}=0$. In figure 3(a), both regions 1 and 2 are under-dense. We see that TLCW only exists when $k_{2}^{-}< k_{z}< k_{1}^{-}$, and its frequency range is $\omega _{\mathrm {p,2}}<\omega <\omega _{\mathrm {p,1}}$. In figure 3(b), region 1 is over-dense while region 2 is under-dense. Now TLCW exists when $k_{z}>k_{2}^{-}$ within the frequency range $\omega _{\mathrm {p,2}}<\omega <|\varOmega |$. Notice that the frequency of TLCW converges to $|\varOmega |$ when $k_{z}\to \infty$. Figure 3 confirms the frequency range derived in (3.5).

Figure 3. The dispersion relation $\omega (k_{y}=0,k_{z})$ in the non-uniform system. The green and magenta lines represent the bulk modes in regions 1 and 2 at $k_{x}=k_{y}=0$; the red lines represent the TLCW; the grey curves represent the other modes. The frequency ranges where TLCW exists are highlighted on the right of each panel.

3.3. The isofrequency contours from 1-D model

After determining the frequency range of possible TLCW, we numerically calculate the isofrequency contours for both the bulk and the surface waves. 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 Fu & Qin (2021).

Figure 4 displays the isofrequency contours of the system when the TLCW exists. We choose the cyclotron frequency to be $|\varOmega |=1$, and both regions are under-dense with $\omega _{\mathrm {p,1}}=0.8$ and $\omega _{\mathrm {p,2}}=0.3$. The frequency range of TLCW according to (3.5) is $0.3<\omega <0.8$. The isofrequency contours at frequency $\omega =0.5$ are shown in figure 4(a) 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 figure 4(b–d) for $\omega >0$, and the $\omega <0$ part can be obtained from the symmetry condition $\omega (-k_{y},k_{z})=-\omega (k_{y},k_{z})$. In these plots, coloured lines represent the surface modes or the bulk modes with $k_{x}=0$, while the grey lines/regions represent other bulk modes with $k_{x}=n{\rm \pi} /2L$, where $n=\pm 1,\pm 2,\ldots$.

Figure 4. 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 grey 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 value of $\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 areas are the frequency gaps of bulk waves.

For the isofrequency contours in figure 4(a), 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 figures 4(a) and 4(b), only bulk modes from region two exist. Near $k_{z}=0.8$, there is a gap for bulk modes in both figure 4(a) (around $0.58\leq k_{z}\leq 0.95$) and figure 4(c) (around $0.42\leq \omega \leq 0.67$). The surface modes that exist and fill up such a frequency gap are what we call 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, although they are continuations of the topological surface waves into the frequency range 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 $\boldsymbol {v}_{\mathrm {g}}=\partial \omega /\partial \boldsymbol {k}$ is perpendicular to the isofrequency contours, we can see that $v_{\mathrm {g},z}>0$ ($v_{\mathrm {g},z}<0$) for $k_{z}>0$ ($k_{z}<0$) and $v_{\mathrm {g},y}>0$ ($v_{\mathrm {g},y}<0$) on the left (right) boundary. Secondly, since by definition the TLCW refers to the topological modes in a frequency gap of bulk modes, it is immune to backscattering when the surface is perturbed, at least when the density perturbation is a function of the $x$ and $y$ coordinates only, and 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 § 4, these properties will be further verified by numerical simulations in both two and three dimensions.

For comparison, two isofrequency contours at two frequencies where no TLCW exists are shown in figure 5. 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 it exists. The isofrequency contour at $\omega =0.85$ is plotted in figure 5(a). 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 figure 5(b), which shows that $\omega =0.85$ is not in the frequency gap of bulk modes. For the case in figure 5(c), $\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 figure 5(c) show that both regions have an ellipsoidal and a hyperboloidal contour 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 figures 5(c) and 5(d) show that no surface mode exists at this $k_{z}$.

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

3.4. Relation to the CMA diagram

In the previous numerical calculation, we see that when the TLCW exists (figure 4a), the shapes of the isofrequency contours of bulk waves in regions one and two are different. On the contrary, when no TLCW exists (figure 5a,c), 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 (Clemmow & Mullaly 1955; Allis 1959; Allis et al. 1963).

In the study of homogeneous cold plasmas, waves at given frequencies can be classified based on the shapes of the wave normal surfaces ($\boldsymbol {v}_{\mathrm {p}}=\omega /c\boldsymbol {k}$), or equivalently the shape of the index-of-refraction surface ($\boldsymbol {n}=c\boldsymbol {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 figure 6. 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 (Volovik 2017) in condensed matter physics. This transition closely connects with the topological phase transition discussed in § 2. The frequency range of the TLCW in (3.5) can be equivalently written as

(3.6a–c)\begin{equation} |{\varOmega}/{\omega}|>1,\quad {\omega_{\mathrm{p,1}}^{2}}/{\omega^{2}}>1,\quad {\omega_{\mathrm{p,2}}^{2}}/{\omega^{2}}<1.\end{equation}

It turns out that, in the CMA diagram, the inequalities in (3.6a–c) also determine the possible shapes of bulk waves. Region one belongs to the two parts at the top right filled with magenta colour, where one bulk wave is a hyperboloid. Region two belongs to the top left filled with green colour, 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, the existence of TLCM implies different shapes of the bulk waves on the two sides of the interface and vice versa.

Figure 6. CMA diagram for cold plasmas with immobile ions (Clemmow & Mullaly 1955; Allis 1959; Allis et al. 1963). The shapes of the index-of-refraction surface are sketched in each region. When TLCW exists, region one must belong to the top right part filled by magenta, while region two must belong to the top left part filled by green.

The CMA diagram includes all transitions of wave shapes, one of which corresponds to the TLCW. However, it is not known whether every transition of wave shape will lead to a topological edge mode. For magnetized cold plasmas, TLCW is the only known topological mode identified so far. Nevertheless, we suspect that similar mechanism might exist for other possible topological modes and the CMA diagram could be used as a guide for searching for these modes.

4. Numerical simulations in two and three dimensions

This section presents numerical simulation results of time-dependent propagation of the TLCW in two and three dimensions, based on the linearized incompressible fluid system specified by (2.1)–(2.3). We first introduce the numerical method adopted for time-dependent simulation and then discuss the simulation results and verify the key physical 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. The time-dependent wave propagation in this section can be found in the supplementary movies available at https://doi.org/10.1017/S0022377822000629.

4.1. Numerical algorithms

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

(4.1a–c)\begin{equation} \begin{pmatrix} v_{x,i,j,k}^{n+{1}/{2}}\\ v_{y,i,j,k}^{n+{1}/{2}}\\ v_{z,i,j,k}^{n+{1}/{2}} \end{pmatrix},\quad \begin{pmatrix} E_{x,i+{1}/{2},j,k}^{n}\\ E_{y,i,j+{1}/{2},k}^{n}\\ E_{z,i,j,k+{1}/{2}}^{n} \end{pmatrix},\quad \begin{pmatrix} B_{x,i,j+{1}/{2},k+{1}/{2}}^{n+{1}/{2}}\\ B_{y,i+{1}/{2},j,k+{1}/{2}}^{n+{1}/{2}}\\ B_{z,i+{1}/{2},j+{1}/{2},k}^{n+{1}/{2}} \end{pmatrix}, \end{equation}

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

(4.2a,b)\begin{equation} \varDelta_{x}[f_{i}]:=\frac{f_{i+1}-f_{i}}{{\rm \Delta} x},\quad f_{\overline{i+{1}/{2}}}:=\frac{f_{i+1}+f_{i}}{2}. \end{equation}

For variables defined on half-integer grid point $g_{i+{1}/{2}}$, similarly,

(4.3a,b)\begin{equation} \varDelta_{x}\left[g_{i+{1}/{2}}\right]:= \frac{g_{i+{1}/{2}}-g_{i-{1}/{2}}}{{\rm \Delta} x},\quad g_{\bar{i}}:=\frac{g_{i+{1}/{2}}+g_{i-{1}/{2}}}{2}. \end{equation}

Using these notations, the discretization of (2.1)–(2.3) is

(4.4)\begin{gather} \left.\begin{gathered} \varDelta_{t}\left[v_{x,i,j,k}^{n+{1}/{2}}\right] = \omega_{\mathrm{p},i,j,k}E_{x,\bar{i},j,k}^{n}-\varOmega v_{y,i,j,k}^{\bar{n}},\\ \varDelta_{t}\left[v_{y,i,j,k}^{n+{1}/{2}}\right] = \omega_{\mathrm{p},i,j,k}E_{y,i,\bar{j},k}^{n}+\varOmega v_{x,i,j,k}^{\bar{n}},\\ \varDelta_{t}\left[v_{z,i,j,k}^{n+{1}/{2}}\right] = \omega_{\mathrm{p},i,j,k}E_{z,i,j,\bar{k}}^{n}. \end{gathered}\right\} \end{gather}

(4.5)\begin{gather} \left.\begin{gathered} \varDelta_{t}\left[E_{x,i+{1}/{2},j,k}^{n}\right] =\varDelta_{y} \left[B_{z,i+{1}/{2},j+{1}/{2},k}^{n+{1}/{2}}\right]-\varDelta_{z} \left[B_{y,i+{1}/{2},j,k+{1}/{2}}^{n+{1}/{2}}\right]- \left(\omega_{\mathrm{p}}v_{x}^{n+{1}/{2}}\right)_{\overline{i+{1}/{2}},j,k},\\ \varDelta_{t}\left[E_{y,i,j+{1}/{2},k}^{n}\right] = \varDelta_{z}\left[B_{x,i,j+{1}/{2},k+{1}/{2}}^{n+{1}/{2}}\right]-\varDelta_{x} \left[B_{z,i+{1}/{2},j+{1}/{2},k}^{n+{1}/{2}}\right]- \left(\omega_{\mathrm{p}}v_{y}^{n+{1}/{2}}\right)_{i,\overline{j+{1}/{2}},k},\\ \varDelta_{t}\left[E_{z,i,j,k+{1}/{2}}^{n}\right] = \varDelta_{x}\left[B_{y,i+{1}/{2},j,k+{1}/{2}}^{n+{1}/{2}}\right]-\varDelta_{y} \left[B_{x,i,j+{1}/{2},k+{1}/{2}}^{n+{1}/{2}}\right]- \left(\omega_{\mathrm{p}}v_{z}^{n+{1}/{2}}\right)_{i,j,\overline{k+{1}/{2}}}. \end{gathered}\right\} \end{gather}

(4.6)\begin{gather} \left.\begin{gathered} \varDelta_{t}\left[B_{x,i,j+{1}/{2},k+{1}/{2}}^{n+{1}/{2}}\right] = \varDelta_{z}\left[E_{y,i,j+{1}/{2},k}^{n}\right]-\varDelta_{y}\left[E_{z,i,j,k+{1}/{2}}^{n}\right],\\ \varDelta_{t}\left[B_{y,i+{1}/{2},j,k+{1}/{2}}^{n+{1}/{2}}\right] = \varDelta_{x}\left[E_{z,i,j,k+{1}/{2}}^{n}\right]-\varDelta_{z}\left[E_{x,i+{1}/{2},j,k}^{n}\right],\\ \varDelta_{t}\left[B_{z,i+{1}/{2},j+{1}/{2},k}^{n+{1}/{2}}\right] = \varDelta_{y}\left[E_{x,i+{1}/{2},j,k}^{n}\right]-\varDelta_{x}\left[E_{y,i,j+{1}/{2},k}^{n}\right]. \end{gathered}\right\} \end{gather}

Here, $\omega _{\mathrm {p}}(\boldsymbol {r})$ is a prescribed function, $\omega _{\mathrm {p},i,j,k}:=\omega _{\mathrm {p}}(x_{i},y_{j},z_{k})$, ${(\omega _{\mathrm {p}}v_{x})_{\overline {i+{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+{1}/{2}}$ through $v^{\bar {n}}$, (4.4) is explicitly solvable by the Cayley transformation that appears in the Boris algorithm (Qin et al. 2013) and other structure-preserving algorithms for charged particle dynamics (He et al. 2015; Zhang et al. 2015; He et al. 2016b, ; Fu, Zhang & Qin 2022).

The algorithm in two dimensions is a special case of the 3-D algorithm described above. When the plasma density is invariant along the background magnetic field, $\omega _{\mathrm {p}}(\boldsymbol {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 (4.4)–(4.6) no longer depend on the index $k$, and a constant $\mathrm {i}k_{z}$ replaces the difference operators $\varDelta _{z}[\cdot ]$.

4.2. Propagation of the TLCW in two and three dimensions

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 right-hand side of the momentum equation (4.4). In the 3-D simulation, the force we used is

(4.7)\begin{equation} \boldsymbol{F}(\boldsymbol{r},t)\sim\exp \left(-\frac{|\boldsymbol{r}-\boldsymbol{r}_{s}|^{2}}{\delta^{2}}\right) \exp({\mathrm{i}(k_{s}z-\omega_{s}t)}), \end{equation}

where $\boldsymbol {r}_{s}$ and $\delta$ are the spatial centre and width of the source and $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 2-D simulation, the force does not depend on $z$, i.e. $k_{s}=0,$ and the factor $\mathrm {i}k_{z}$ replaces the difference operators $\varDelta _{z}[\cdot ]$ in (4.5) and (4.6).

The parameters for the first case simulated are $\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 one dimension was shown in figures 4(a) and 4(c). The centre of the source is on the interface. The time-dependent propagation of the TLCW in two and three dimensions after a finite time are shown in figure 7. Notice that the configuration is equivalent to the right interface shown in figure 2, where the blue lines in figure 4 represent the surface waves in figure 4. In the 2-D simulation displayed in figure 7(a), 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 phenomena can also be observed in 3-D simulation shown in figure 7(b), where the TLCW is excited and propagates towards the upper-left direction in the $y$–$z$ plane.

Figure 7. Propagation of TLCW in (a) two dimensions and (b) three dimensions at $|t\varOmega |=250$, 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 colour map indicates the strength of the real part of $E_{z}$, the yellow star marks the location of the source and the dashed line is the interface between two regions. In (b), the red and blue contours represent the locations of $\mathrm {Re}[E_{z}]=\pm 0.15\max [\mathrm {Re}(E_{z})]$, the yellow sphere marks the source and the grey surface is the interface between two regions.

4.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 (Barnett 2010; Dodin & Fisch 2012), i.e. the (Minkowski) canonical momentum and the (Abraham) kinetic momentum. The canonical momentum of the TLCW is proportional to the wavenumber $\boldsymbol {k}$, and it can vary significantly, as illustrated in figure 4. However, the kinetic momentum, which is proportional to it group velocity $\boldsymbol {v}_{\mathrm {g}}$, is unidirectional. For example, the surface waves on the right interface can only have kinetic momentum in the upper-left or lower-left direction in the $y$–$z$ plane.

For the cold plasma model given by (2.1)–(2.3), the dielectric tensor is Hermitian and independent of $\boldsymbol {k}$. Therefore, the kinetic momentum, the group velocity and the flux of canonical energy are all proportional to the Poynting vectors (Stix 1992) $\boldsymbol {S}\sim \mathrm {Re}[\boldsymbol {E}\times \boldsymbol {B}^{*}]$. The Poynting vectors of the TLCW excited on a planar interface are shown in figure 8. As expected, the kinetic momentum is along the direction of wave propagation

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

4.4. Propagation of the TLCW on non-planar interface

The main properties of the TLCW are 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 figure 9 is the TLCW propagation on a zig-zagged interface after a finite time, which clearly demonstrates immunity of backscattering and unidirectional propagation. When the TLCW meets sharp turns on the interface, it propagates along the interface without reflection or transition to bulk waves.

Figure 9. (a) Two- and (b) three-dimensional simulations of the TLCW excited at a zig-zagged interface The source was turned on at $t=0$ and the field strength is plotted at $|t\varOmega |=600$.

Additional examples are shown in figures 10 and 11, where the interface 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 \boldsymbol {r}\times \boldsymbol {S}\,\mathrm {d}V$. Notice that the source exciting the TLCW does not carry any angular momentum. The existence of an angular-momentum-carrying surface wave reflects the topological property of the bulk material and is predicted by 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 10. (a) Two- and (b) three-dimensional simulations of the TLCW excited on a square interface. The source was turned on at $t=0$ and the field strength is plotted at $|t\varOmega |=600$.

Figure 11. (a) Two- and (b) three-dimensional simulations of the TLCW excited on a circular interface. The source was turned on at $t=0$ and the field strength is plotted at $|t\varOmega |=600$.

5. Conclusions and discussion

We have studied the dispersion relation and propagation of the recently identified TLCW. From the excitation condition (2.8), we derived the frequency range in (3.5) 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 backscattering because it lies in the frequency gap of bulk waves. The unidirectional propagation and immunity of backscattering 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 (Bergholtz, Budich & Kunst 2021) and nonlinear (Smirnova et al. 2020; Zhou et al. 2020; Bergholtz et al. 2021) effects may also be important for topological waves in more complex fluid and plasma models (Qin et al. 2019; Fu & Qin 2020; Qin et al. 2021; Zhu et al. 2021).