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 Reference Faure2019; Delplace Reference Delplace2022; Qin & Fu Reference Qin and Fu2022), 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 Reference Fu and Qin2021). 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 Reference Qin and Fu2022) 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 (Reference Fu and Qin2021). 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 (Reference Fu and Qin2021).

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 Reference Clemmow and Mullaly1955; Allis Reference Allis1959; Allis et al. Reference Allis, Buchsbaum and Bers1963).

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 Reference Volovik2017) 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 Reference Clemmow and Mullaly1955; Allis Reference Allis1959; Allis et al. Reference Allis, Buchsbaum and Bers1963). 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.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. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013), and solve the Maxwell equations (2.2) and (2.3) on the Yee grid (Yee Reference Yee1966). 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. Reference Qin, Zhang, Xiao, Liu, Sun and Tang2013) and other structure-preserving algorithms for charged particle dynamics (He et al. Reference He, Sun, Liu and Qin2015; Zhang et al. Reference Zhang, Liu, Qin, Wang, He and Sun2015; He et al. Reference He, Sun, Liu and Qin2016b, ; Fu, Zhang & Qin Reference Fu, Zhang and Qin2022).

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 Reference Barnett2010; Dodin & Fisch Reference Dodin and Fisch2012), 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 Reference Stix1992) $\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$.