# Nonreciprocal surface plasmons in angularly varying, magnetized, metasurface
tubes

Yarden Mazor  School Of Electrical Engineering, Tel Aviv University, Tel Aviv,
69978, Israel<EMAIL_ADDRESS>
###### Abstract

We analytically and numerically study the nonreciprocal surface waves guided
by a magnetized metasurface tube. When applying the magnetic bias
perpendicularly to the cylinder axis, the conductivity profile cross-section
is nonuniform, which enables us to obtain pronounced nonreciprocal modes with
different field distributions for propagation in opposite directions due to
the coupling between different values of orbital angular momentum. We show
that this property can be leveraged for directional power delivery when
changing the source location, and we study the isolation ratio as a function
of different parameters. The simple, continuous dependence of the isolation on
the direction of magnetization allows simple yet robust control over the wave
propagation properties.

##  I Introduction

Cylindrical structures are among the most common waveguide geometries being
used. Prominent examples include coaxial waveguides and hollow metallic tubes
[1]. The dispersion of guided waves can be controlled by varying the cylinder
size (potentially also the internal cylinder in a coaxial geometry) and the
distribution and properties of the filling material. These waveguides have
many uses, such as communications [1], open-ended aperture antennas [2, 3],
material sensing probes, and more.

More recently, incorporating metamaterials and metasurfaces within a
cylindrical geometry was examined as a platform for various applications.
Scattering from cylindrical metasurface structures was studied in [4] to
achieve cloaking, and in [5], the authors show how to systematically analyze
and synthesize a layered cylindrical structure to achieve different scattering
functionalities. Another venue is tailoring guided wave properties in
cylindrical waveguides using metasurfaces coating the waveguide walls or
serving as the waveguide themselves. This is a particularly interesting venue
since cylindrical guided modes exhibit complex polarization states, helicity,
and spin-OAM coupling [6]. This was pioneered in [7], where a helically
slotted PEC cylinder was studied. A modal theory for such structures with
general electric impedance walls was formulated in [8]. In [9], a metamaterial
lining was used to control the cut-off and modal properties, and in [10] it
was shown how general bianisotropic walls can be used to tailor the modal
dispersion, field distribution, and spin. A step forward was taken in [11],
where the authors discuss an actual implementation and consider realistic
constraints and the arising dispersion effects. Such structures also appear in
more natural systems, where the most prominent example is carbon nanotubes.
The electrodynamics of cylindrical graphene-based structures was first studied
in [12], where the appropriate boundary conditions were formulated, and the
dispersion of guided surface waves was carefully mapped. In the last decade,
several groups have fabricated and measured artificial cylindrical graphene
structures, for enhancing light-matter interaction [13, 14], four-wave mixing
[15] and all-optical modulation [16]. Multilayered folding is shown in [17],
and additional applications are reviewed in [18].

An equally important ingredient of this work is the nonreciprocal propagation
of waves. Nonreciprocal waveguides have many applications, such as routing and
sorting of signals, and serving as building blocks in the design of
nonreciprocal components such as isolators and circulators [19, 20]. To obtain
significant nonreciprocal behavior, such as different field distributions for
oppositely propgating fields which can be leveraged for isolation [1] (e.g.,
asymetrically ferrite loaded waveguides), or one-way guiding, one has to
incorporate nonreciprocity with some form of structural asymmetry [21, 22]. In
the context of tube-like structures, [23] studied guided waves in several
cylindrical geometries incorporating graphene, where electrical bias
(controlling the chemical potential $\mu_{c}$) and magnetic bias were used as
control parameters to manipulate the wave characteristics. Importantly, this
work also incorporated nonreciprocity since the magnetic field bias gives rise
to an asymmetric conductivity tensor of the graphene layer. Nevertheless, this
nonreciprocity did not visibly manifest in the modal dispersion or field
distribution due to the highly ”symmetric” magnetization scheme considered - a
radial magnetic field bias.

In this work, we would like to examine the nonreciprocal wave propagation on a
magnetized cylindrical impedance surface, such as graphene. Unlike other
examples, where the structural asymmetry requires ingomogeneous cross section
materials or scatterer parameters, here it stems from the magnetization scheme
itself. By applying a simple magnetic bias perpendicular to the waveguiding
surface (represented by a thin impedance sheet), an azimuthally dependant
conductivity tensor arises. Therefore, the magnetization itself provides both
the nonreciprocity and the spatial inhomogeneity. This, in turn, creates
coupling between different values of orbital angular momentum (OAM) of the
guided wave components. The nonreciprocity makes this coupling asymmetric,
giving rise to complex nonreciprocal waves composed of a superposition of
different cylindrical eigenfunctions (each with its characteristic OAM),
resulting in different field distributions modes propagating in opposite
directions. This will be shown to enable directional excitation using simple
dipolar sources as a function of the source location and spin-momentum
coupling.

##  II Formulation

Let us start from the basic formulation of the problem of waveguiding by thin
cylindrical shells. Since the waves propagate essentially in free space (with
the addition of the impedance/surface susceptibility boundary condition), we
can represent them as a superposition of the ”standard” TE and TM modes

$E_{z}=\begin{Bmatrix}A_{n}^{i}I_{n}(\alpha r)&;r<a\\\ A_{n}^{o}K_{n}(\alpha
r)&;r>a\end{Bmatrix}e^{-j\beta z}e^{-jn\varphi},$ (1a)
$H_{z}=\begin{Bmatrix}B_{n}^{i}I_{n}(\alpha r)&;r<a\\\ B_{n}^{o}K_{n}(\alpha
r)&;r>a\end{Bmatrix}e^{-j\beta z}e^{-jn\varphi},$ (1b)

with $\beta^{2}=\alpha^{2}+k_{0}^{2}$, and $k_{0}=\omega\sqrt{\mu\epsilon}$
and $I_{n},K_{n}$ are the modified Bessel functions of the first and second
kind. The other field components are given can be directly derived from
equations (1a),(1b), and are given in appendix A for completeness.

In the present problem, the waves are guided by a tube of radius $R$, which
satisfies the surface impedance boundary condition

$\hat{\boldsymbol{r}}\times\left(\boldsymbol{H}_{out}-\boldsymbol{H}_{in}\right)_{r=R}=\boldsymbol{J}_{s}=\underline{\underline{\boldsymbol{\sigma}}}\boldsymbol{E}_{tan,r=R},$
(2)

where $\underline{\underline{\boldsymbol{\sigma}}}$ is a $2\times 2$ matrix,
that can be written explicitly as

$\underline{\underline{\boldsymbol{\sigma}}}=\begin{bmatrix}\sigma_{\varphi\varphi}&\sigma_{\varphi
z}\\\ \sigma_{z\varphi}&\sigma_{zz},\end{bmatrix}$ (3)

and $\boldsymbol{E}_{tan}$ is the electric field components tangent to the
tube. In addition to this, since the cylindrical surface only hosts electric
surface currents, the tangential electric fields need to be continuous

$E_{z}(R^{+})=E_{z}(R^{-})\;,\;E_{\varphi}(R^{+})=E_{\varphi}(R^{-}),$ (4)

where $()^{+}and()^{-}$ indicate sampling of the fields on the surface of the
tube, on the outer and inner side, respectively. Using this condition, the
field representation on the surface of the tube can be reduced to rely on two
coefficients $A_{n}^{i},\bar{B}_{n}^{i}=\eta_{0}B_{n}^{i}$, (see appendix A).
If we now define

$\boldsymbol{D}_{n}=\begin{pmatrix}A_{n}^{i}\\\ \bar{B}_{n}^{i}\end{pmatrix},$
(5)

we can represent the operations required for substitution into equation (2) as

$\boldsymbol{E}^{o}_{n,tan}(r=R)=\boldsymbol{E}^{i}_{n,tan}(r=R)=\underline{\underline{\boldsymbol{M}}}_{E,n}\boldsymbol{D}_{n},$
(6a)
$\hat{\boldsymbol{r}}\times(\boldsymbol{H}^{o}_{n}-\boldsymbol{H}^{i}_{n})_{(r=R)}=\underline{\underline{\boldsymbol{M}}}_{H,n}\boldsymbol{D}_{n},$
(6b)

with

$\underline{\underline{\boldsymbol{M}}}_{E,n}=\begin{bmatrix}\frac{n\bar{\beta}}{\bar{\alpha}^{2}}I_{n}&-\frac{j\bar{k}_{0}}{\bar{\alpha}}I_{n}^{\prime}\\\
I_{n}&0\end{bmatrix}\;,\underline{\underline{\boldsymbol{M}}}_{H,n}=\begin{bmatrix}0&-\frac{1}{\bar{\alpha}K_{n}^{\prime}}\\\
-\frac{j\bar{k}_{0}}{\bar{\alpha}K_{n}}&\frac{n\bar{\beta}}{\bar{\alpha}^{3}K_{n}^{\prime}}\end{bmatrix}.$
(7)

To shorten the notation, we will omit the argument of the modified Bessel
function $I_{n},K_{n}$ and their derivatives $I_{n}^{\prime},K_{n}^{\prime}$
when we refer to sampling these on the surface of the tube. In addition, we
have defined $\bar{\alpha}=\alpha R\;,\;\bar{\beta}=\beta
R\;,\;\bar{k}_{0}=k_{0}R$. Now, if we substitute this representation into the
boundary condition given in equation (2) and assume the conductivity
$\underline{\underline{\boldsymbol{\sigma}}}$ is independant of $\varphi$, we
obtain the eigenmode equation

$\left(\eta_{0}\underline{\underline{\boldsymbol{\sigma}}}\,\underline{\underline{\boldsymbol{M}}}_{E,n}-\underline{\underline{\boldsymbol{M}}}_{H,n}\right)\boldsymbol{D}_{n}=0.$
(8)

For this equation to have a non-trivial solution for the coefficient vector
$\boldsymbol{D}_{n}$, the determinant of the bracketed matrix must vanish,
which yields the dispersion relation as a determinant of a $2\times 2$ matrix.
For any value of $n$ (the OAM of the mode) and angular frequency $\omega$ we
substitute, we will obtain $\beta(\omega)$ \- the dispersion relation of the
$n$’th mode. Substituting this value back will yield the modal coefficients
$A_{n}^{i},\bar{B}_{n}^{i}$, forming a mixed TE/TM mode, with the amplitude of
each constituent depending on the conductivity.

Based on this, our goal now would be to describe the waves propagating in a
system where $\underline{\underline{\boldsymbol{\sigma}}}$ is a function of
the azimuthal angle $\varphi$, generally expressed as
$\underline{\underline{\boldsymbol{\sigma}}}(\varphi)$. We start by describing
all the relevant quantities
$\boldsymbol{E},\boldsymbol{H},\underline{\underline{\boldsymbol{\sigma}}}$
using their Fourier series in $\varphi$. For each specific mode $n$ we use the
matrices defined in equations (6a),(6b),(7) to represent the fields using the
corresponding coefficients $\boldsymbol{D}_{n}$, yieding

$\begin{pmatrix}E_{\varphi}^{i}\\\
E_{z}^{i}\end{pmatrix}=\sum_{n=-\infty}^{\infty}\underline{\underline{\boldsymbol{M}}}_{E,n}\boldsymbol{D}_{n}e^{-jn\varphi},$
(9a)
$\hat{\boldsymbol{r}}\times\begin{pmatrix}H_{\varphi}^{o}-H_{\varphi}^{i}\\\
H_{z}^{o}-H_{z}^{i}\end{pmatrix}=\sum_{n=-\infty}^{\infty}\underline{\underline{\boldsymbol{M}}}_{H,n}\boldsymbol{D}_{n}e^{-jn\varphi},$
(9b)

and for the condictivity

$\eta\underline{\underline{\boldsymbol{\sigma}}}=\sum_{m=-\infty}^{\infty}\underline{\underline{\boldsymbol{Y}}}_{m}e^{-jm\varphi}.$
(10)

When
$\underline{\underline{\boldsymbol{\boldsymbol{\sigma}}}}=\underline{\underline{\boldsymbol{\boldsymbol{\sigma}}}}(\varphi)$,
there is coupling between different components with different $n$, through the
additional OAM provided by the geometrical structure of
$\underline{\underline{\boldsymbol{\sigma}}}$, and we can derive the modal
equation by substituting these expansions into equation (2),

$\sum_{n=-\infty}^{\infty}\underline{\underline{\boldsymbol{Y}}}_{\ell-n}\underline{\underline{\boldsymbol{M}}}_{E,n}\boldsymbol{D}_{n}-\underline{\underline{\boldsymbol{M}}}_{H,\ell}\boldsymbol{D}_{\ell}=0\;\;,\;\;\forall\ell,$
(11)

where each value of $\ell$ represents a row in a $2\times 2$ infinite block
matrix. This can be reformulated into a matrix equation, as shown in appendix
B.

For this system to have a nontrivial solution, the determinant of the matrix
(29) must vanish, giving the dispersion equation for the propagating modes.
Naturally, since different values of $n$ are coupled through the angle-
dependent conductivity of the tube, the modes will be composed of many
azimuthal harmonics excited and propagating together. This is in addition to
the natural TE/TM coupling in higher order modes propagating on impedance
cylinders, giving rise to hybrid-type modes that can yield interesting
patterns of spin and helicity [10].

For the sake of simplicity, We will limit our discussion here to conductivity
matrices of the form

$\underline{\underline{\boldsymbol{\sigma}}}(\varphi)=\underline{\underline{\boldsymbol{\sigma}}}_{0}+\underline{\underline{\boldsymbol{\sigma}}}_{1}\cos\varphi=\underline{\underline{\boldsymbol{\sigma}}}_{0}+\frac{1}{2}\left(\underline{\underline{\boldsymbol{\sigma}}}_{1}e^{j\varphi}+\underline{\underline{\boldsymbol{\sigma}}}_{1}e^{-j\varphi}\right).$
(12)

This would allow us to simplify equation (11) and its matrix representation
(29) to a tridiagonal block form. In this case, one could formally express the
dispersion relation in a matrix continued fraction form, in a process similar
to the one described in [24]

$\underline{\underline{\boldsymbol{Y}}}_{1}^{-1}\left[-\underline{\underline{\boldsymbol{Y}}}_{0}+\underline{\underline{\boldsymbol{G}}}_{n}+\left[-\underline{\underline{\boldsymbol{Y}}}_{0}+\underline{\underline{\boldsymbol{G}}}_{n-1}+\left[\underline{\underline{\boldsymbol{Y}}}_{0}+\underline{\underline{\boldsymbol{G}}}_{n-2}+[...]^{-1}\right]^{-1}\right]^{-1}\right]=\\\
\left[\underline{\underline{\boldsymbol{Y}}}_{1}^{-1}\underline{\underline{\boldsymbol{Y}}}_{0}-\underline{\underline{\boldsymbol{Y}}}_{1}^{-1}\underline{\underline{\boldsymbol{G}}}_{n+1}+\left[\underline{\underline{\boldsymbol{Y}}}_{1}^{-1}\underline{\underline{\boldsymbol{Y}}}_{0}-\underline{\underline{\boldsymbol{Y}}}_{1}^{-1}\underline{\underline{\boldsymbol{G}}}_{n+2}+[...]^{-1}\right]^{-1}\right]^{-1}$
(13)

With
$\underline{\underline{\boldsymbol{G}}}_{n}=\underline{\underline{\boldsymbol{M}}}_{H,n}\underline{\underline{\boldsymbol{M}}}_{E,n}^{-1}$,
and we have assumed
$\underline{\underline{\boldsymbol{Y}}}_{1}=\underline{\underline{\boldsymbol{Y}}}_{-1}$
as in our simplified case.

Since we would like to study nonreciprocal propagation, it is worth noticing
that when the guiding structure is reciprocal (satisfying
$\underline{\underline{\boldsymbol{\sigma}}}=\underline{\underline{\boldsymbol{\sigma}}}^{T}$),
inverting the propagation direction would result in time-reversed modal
fields. For a specific mode, the modal fields will be transformed as
$\boldsymbol{E}(-\beta)=\boldsymbol{E}(\beta)^{*}$,$\boldsymbol{H}(-\beta)=-\boldsymbol{H}^{*}(\beta)$.
However, this is no longer the case when the conductivity tensor does not
satisfy reciprocity. In many cases, merely violating the reciprocity of the
guiding structure is not enough for visible effects of non-reciprocity to
manifest in the guided modes (such as a non-symmetric dispersion curve and
different field profiles for oppositely propagating waves). This is one of the
reasons that [23] did not observe any significant nonreciprocal propagation
and excitation effects, although studying a nonreciprocal structure. However,
we will see that more significant nonreciprocal effects are possible in our
case due to the explicit geometrical asymmetry arising from the $\varphi$
dependence of the conductivity matrix.

###  II.1 Transversely magnetized gyrotropic tube

The case we would like to focus on is where the conductivity of the cylinder
can be expressed as

$\underline{\underline{\boldsymbol{\boldsymbol{\sigma}}}}(\varphi)=\frac{1}{\eta_{0}}\begin{pmatrix}jX&X_{o}\cos\varphi\\\
-X_{o}\cos\varphi&jX\end{pmatrix}.$ (14)

Although this choice seems quite specific, such a system naturally arises if
we envision a Graphene tube, and magnetize it perpendicularly to the cylinder,
as shown in figure 1(a) (in a setup similar to [13, 25] for instance). Since
the magnitude of the off-diagonal conductivity term is proportional to the
magnetic field component perpendicular to the surface, and assuming the
external magnetic field is uniform, say $\boldsymbol{B}=B_{0}\hat{x}$, we will
obtain such a conductivity matrix [26]. Due to the magnetization, the diagonal
terms would also depend slightly on $\varphi$. However, for magnetic fields up
to $~{}1.5T$ this dependence is relatively weak ($<10\%$ deviation from the
value without magnetic field). Therefore, we will neglect this dependence,
taking $X$ in equation (14) independent of $\varphi$. This variation can be
accounted for straightforwardly using the same formulation. Cylindrical
arrangements of plasmonic patches or other particles are also expected to
yield a similar conductivity structure (also shown in figure 1(a)). Extreme
wave propagation effects using plasmonic particles and patches was
demonstrated in [27], and nonreciprocity was shown in a planar setup in [28].
Here, we do not require the particle anisotropy (implemented through the usage
of elongated ellipsoidal pathces) since the spatial symmetry breaking required
for significant nonreciprocity occurs due to the cross-section inhomogeneity
in our case.

Figure 1: (a) Possible geometries for magnetizable cylindrical surface
waveguides. (b) Graphene conductivity under a magnetization of $B_{0}=1.5T$.
(c) Surface wave dispersion. Blue - no magnetization, red - with
magnetization, green - accounting for the dispersion of graphene conductivity.

To analyze this system using the developed formulation, we represent the
conductivity tensor using its $\varphi$ Fourier transform. The matrix
coefficients are

$\underline{\underline{\boldsymbol{Y}}}_{0}=\begin{pmatrix}jX&0\\\
0&jX\end{pmatrix},\;\;\underline{\underline{\boldsymbol{Y}}}_{-1}=\underline{\underline{\boldsymbol{Y}}}_{1}=\frac{1}{2}\begin{pmatrix}0&-X_{o}\\\
X_{o}&0\end{pmatrix}.$ (15)

Once substituted into equation (11) (or 29), we obtain a tridiagonal block
system of equations. Although the magnetic bias makes this system
nonreciprocal, we know from previous works that nonreciprocity can manifest in
several ways, which are often quite subtle.

The conductivity matrix elements for Graphene are shown in figure 1(b). The
parameters we take are the Fermi velocity $v_{F}=10^{6}m/s$, the chemical
potential $\mu_{c}=0.35V$, and a biasing magnetic field of $B_{0}=1.5T$. Since
the response of Graphene can be quite lossy in the frequency regions where the
gyrotropic effect is significant (yielding meaningful off-diagonal
conductivity terms), we want to choose our frequency such that the losses are
not too significant as to mask the guided wave propagation. Here, the
operation frequency we choose is $4THz$, corresponding to
$X\approx-0.6j+0.135,X_{o}\approx-0.093-0.043j$.

##  III Dispersion

###  III.1 Unmagnetized tube

In the unmagnetized case, the dispersion for each value of OAM, $n$, can be
solved separately (corresponding to the matrix in equation 29 being block
diagonal). Following equation (8), the dispersion equation can be written as

$det\left[\eta\underline{\underline{\boldsymbol{\sigma}}}\;\underline{\underline{\boldsymbol{M}}}_{E,n}-\underline{\underline{\boldsymbol{M}}}_{H,n}\right]=0.$
(16)

The dispersion curves for this case are shown using thick dotted blue lines in
figure 1(c). There exists a degeneracy for positive and negative values of
$n$, and both $n=1$ and $n=-1$ have the same dispersion. In our case, the
unmagnetized tube conductivity is scalar, and therefore for $n=0$ (and only
for this value), equation (16) can be completely separated to TE/TM, where
only one of them will constitute the propagating mode according to the sign of
the imaginary value of $\sigma$ (the sign of $X$). Here, we use $Im[\sigma]<0$
(an inductive surface), and therefore the $n=0$ mode will be a TM mode, with
$B_{0}^{i}=0$. For $n\neq 0$, there is always some ”mixing” between the TE and
TM constituents, although a TM component will still dominate the modes due to
the inductive response.

###  III.2 Magnetized tube

When applying the magnetic bias, several effects come into play. First,
$\underline{\underline{\boldsymbol{\sigma}}}$ is no longer uniform but depends
on $\varphi$. This creates coupling between different values of OAM, so we can
no longer separate the dispersion equations for distinct values of $n$.
Consequently, each modal field distribution, derived from the eigenvector of a
certain solution of the dispersion equation, will be composed of contributions
from several $n$’s. Usually, significant contributions will only come from one
or a few $n$, but at least theoretically, all of them are a part of the modal
field. Due to this, the magnitude of the modal field will no longer be uniform
across $\varphi$ (as expected in the unmagnetized case, for a simple
$e^{jn\varphi}$ dependence). It is worth mentioning that this is not related
to the nonreciprocal response of the magnetized tube but to the cross-
sectional inhomogeneity of the conductivity matrix. However, due to the
nonreciprocity, this is expected to happen in a non-symmetric way -
propagating modes with $\beta>0$ and $\beta<0$ are expected to possess
differently distributed fields. This opens the possibility of directional
excitation of waves by arbitrary sources, as will be shown later.

The dispersion curve for a magnetized cylinder, with $B_{0}=1.5T$, is shown in
solid magenta lines in figure 1(c), where we have neglected the associated
losses to obtain a clear modal dispersion. Due to the coupling between
different OAM values and the nonreciprocity, the degeneracy is now broken, and
the $n=\pm 1$ branch splits into two separate branches. The $n=0$ branch,
which was not degenerate in the first place, experiences some frequency shift.
For this calculation, the infinite matrix dispersion equation was truncated at
$n=\pm 7$, which is enough for the specific modes in our examples. Higher-
order modes might require a different truncation.

The calculation of this dispersion assumes that the conductivity values did
not change as a function of $k_{0}R$. Therefore, the vertical axis can be seen
as a variation of the cylinder radius, for instance. A variation in the
frequency would require that we also vary the conductivity matrix (both the
diagonal and off-diagonal terms depend on frequency). An example of the
surface wave dispersion when considering this dependence
$\underline{\underline{\boldsymbol{\sigma}}}(\omega)$ is shown in green in
figure 1(c). We see that qualitatively, the curves are similar, yet there is a
difference in the obtained wavenumbers due to the varying conductivity values
for different frequencies. The radius of the tube is taken as $a=6\mu m$ here.
The dashed black line indicates the operation frequency for which we later
calculate the fields in different scenarios.

Figure 2: Mode field profiles for a magnetized graphene cylinder. (a) and (b)
show the modes that split from the degenerate $n=\pm 1$ mode in the
unmagnetized case. (c) The perturbed $n=0$ mode. Thin lines show the field
profiles for $\beta<0$, displaying significant asymmetry. (d) The modal
coefficients $A_{n},B_{n}$ for the mode shown in (c).

To illustrate the modal fields, we calculate the modal coefficients. The
operation point was chosen as $k_{0}R\approx 0.5$, corresponding to a cylinder
with radius $R=6[\mu m]$ in $f=4THz$. Figure 2 shows the longitudinal fields
$|E_{z}|,|H_{z}|$ for the three modes that solve the presented dispersion.
Panels (a) and (b) show the modes that split from the $n=\pm 1$ branch,
presenting a mixed TE/TM nature. The mode close to the $n=0$ branch, shown in
panel (c), is TM dominated (notice the relatively weaker $H_{z}$ content)
since around this regime, the TE/TM coupling is weaker, a remnant of the pure
TM $n=0$ mode in the unmagnetized case. The thin lines in panel (c) show the
modal fields for $\beta<0$, and we see a significant difference in the field
distribution compared to $\beta>0$, a clear consequence of the nonreciprocity
(in the reciprocal case, we would expect a simple conjugation operation, which
would preserve the profile of $|E_{z}|,|H_{z}|$). The electric field
oscillates ”slowly” since the $n=0$ coefficient still dominates it, but the
magnetic field oscillates more rapidly since the magnetic fields are dominated
by the $n=\pm 1$ coefficients, coupled through the nonuniform conductivity.
This is confirmed in panel (d), where we see the coefficients composing the
field corresponding to the mode shown in panel (c).

##  IV Excitation of waves

We performed full-wave simulations of different excitation schemes to
illustrate the various aspects of this system using COMSOL multiphysics [29].
Throughout this section, we use the operation point previously introduced -
$f=4[THz],R=6[\mu m],B_{0}=1.5[T]$, unless explicitly stated otherwise.

###  IV.1 Excitation centered around specific OAM

We start by inspecting the fields excited by a source centered on a specific
value of $n$. To this end, we used a ring surrounding the magnetized tube,
with a radius of $1.1R$, as shown in figure 3(a). On this ring, we apply a
surface current $\boldsymbol{J}_{s}=\hat{\boldsymbol{z}}e^{-jn\varphi}$ and
examine the excited fields. Figure 3(b) shows the profile of $E_{z},H_{z}$,
similar to Fig. 2. The dashed line shows the distribution obtained from the
full-wave simulation, and the solid line shows the analytical result for
comparison. We see a good match between the field profiles since the
analytically obtained electric field mainly comprises the $n=0$ component. It
is worth mentioning that while the field distributions along $\varphi$ are
very similar, they are not identical (in terms of the ratio between the field
magnitudes, for instance). This is because we cannot excite a single mode this
way, and the propagating fields contain at least some contribution from other
modes.

Figure 3: (a) Surface current source is placed on the yellow ring around the
cylinder. (b) In dashed, the simulated $\boldsymbol{z}$ component of the
propagating fields excited by the source with $n=0$, as a function of
$\varphi$, measured on a circle in the XY plane with radius $1.1R$, placed
$L/3$ from the source along the $+\hat{\boldsymbol{z}}$ axis. Solid lines show
the analytical field distribution of the ’$n=0$’ mode. (c) Distribution of the
excited electric field longitudinal component $E_{z}$. The guiding cylinder is
between the dashed lines. (d) Same as (c) but for $|\boldsymbol{E}|$.

In figure 3(c), we see how propagation to the $+\hat{\boldsymbol{z}}$
direction has fields that are most dominant on the upper side (corresponding
to $\varphi=\pi/2$) whereas in propagation to the $-\hat{\boldsymbol{z}}$ the
lower side is more dominant. This is, again, a direct consequence of the
system nonreciprocity and wouldn’t occur if the system response was reciprocal
(even if it is strongly dependant on $\varphi$, or even anisotropic). From
figure 3(c) we are also able to approximate the wavelength of the propagating
wave, and it corresponds to $~{}\beta R\approx 1.6$ which is close to the
value we extract from the analytical formulation, $\beta R\approx 1.7$. The
inaccuracy is caused by the coupling of various modes in the numerical
simulation.

###  IV.2 Excitation by a point source

Next, we would like to understand how the modal properties manifest when
exciting the system with a point source. When using a point source, in
principle, a superposition of all possible wave species is excited, each with
a different amplitude. With such excitation, we can feel the system properties
as a whole, essentially ”probing” the Green’s function. One of the
consequences of nonreciprocity is that even with simple sources, we get
pronounced asymmetry in the propagation of power since the field distributions
of the modes that propagate in opposite directions are different. We simulated
the fields excited by a dipole source $\boldsymbol{p}=p\hat{\boldsymbol{z}}$
in different locations $[x_{s},y_{s}]$ around the center of the cylindrical
waveguide, $x_{s}=1.1R\cos\varphi_{s},y_{s}=1.1R\sin\varphi_{s}$, with $z_{s}$
being in the central plane, $z_{s}=L/2$. We then calculated the total power
flowing through two cross-sections at $z=L/4$ and $z=3L/4$.

Figure 4: (a) Ratio between the total power passing through the tube cross-
section to +z and -z directions, as a function of the location of the dipole
source around the cylinder, $[x_{s},y_{s}]$. (b) $E_{z}$ for
$\phi_{s}=\pi/2,3\pi/2$, source dipole is shown in magenta. The asymmetry is
visible. (c) $E_{z}$ for a smaller cylinder, with $k_{0}R\approx 0.25$.

Figure 4(a) shows the ratio between the power that propagates to the $+z$ and
$-z$ directions - the isolation ratio - as a function of the source location
angle $\varphi_{s}$, for different values of the magnetic bias. Significant
nonreciprocity is present, with $\phi_{s}=\pi/2,3\pi/2$ having the strongest
asymmetry in power propagation. In 4(b) we see the propagating fields for
$\phi_{s}=\pi/2$ ($\phi_{s}=3\pi/2$) in the top (bottom) panel. The asymmetry
is visible, and it is worth mentioning that this propagating power is not a
result of a single propagating mode but all of the possible ones (3 in this
case). Evidence of this can be seen from the ”inconsistent” wavefronts in
4(b), since the excited modes also have different attenuation and propagation
constants. Examining the dispersion, we see that for $k_{0}R<~{}0.25$ we can
obtain a ”single-mode” operation regime. Figure 4(c) shows the propagating
fields for this case when the source is placed at $\phi_{s}=\pi/2$, and we can
see the much more uniform wavefronts. Although the asymmetry in the field
picture might not seem as pronounced, calculating the $+z$ to $-z$ propagating
power reveals a similar picture, as shown in the dashed green line in figure
4(a). This demonstrates that isolation is a fundamental phenomenon here and is
not dependent the ”mixed” excitation of modes.

Figure 5: (a) Surface wave dispersion when the Graphene sheet is wrapped
around a dielectric core for $\epsilon_{r,core}=[2,3,10]$. (b) Modal TE/TM
composition for three modes (out of 5), when $\epsilon_{r}=3$. (c) The
isolation ratio for different values of $\epsilon_{r}$

##  V A magnetized tube wrapped around a dielectric core

When the conductive sheet is wrapped around a dielectric core, the value of
$\alpha$ in equations (1a) and (1b) differs between the inner and outer parts,
which we will term $\alpha_{i,o}$. If the core has a dielectric constant of
$\epsilon_{r}$ we will have $\beta^{2}=\alpha_{i}^{2}+\epsilon_{r}k_{0}^{2}$
and $\beta^{2}=\alpha_{o}^{2}+k_{0}^{2}$. Substituting these into the same
equations and performing similar manipulations to the ones shown in appendix
A, yields the same infinite matrix relation shown in appendix B, with the
reformulated building blocks

$\underline{\underline{\boldsymbol{M}}}_{E,n}=\begin{bmatrix}\frac{n\bar{\beta}}{\bar{\alpha}_{i}^{2}}I_{n}&-\frac{j\bar{k}_{0}\sqrt{\epsilon_{r}}}{\bar{\alpha}_{i}}I_{n}^{\prime}\\\
I_{n}&0\end{bmatrix}\;\;,\;\;\underline{\underline{\boldsymbol{M}}}_{H,n}=\underline{\underline{\boldsymbol{M}}}_{H,n}^{o}-\underline{\underline{\boldsymbol{M}}}_{H,n}^{i},$
(17)

where

$\underline{\underline{\boldsymbol{M}}}_{H,n}^{i}=\frac{1}{\eta_{0}}\begin{bmatrix}0&-I_{n}\sqrt{\epsilon_{r}}\\\
\frac{j\bar{k}_{0}\sqrt{\epsilon_{r}}}{\bar{\alpha}_{i}}I_{n}^{\prime}&\frac{n\bar{\beta}\sqrt{\epsilon_{r}}}{\bar{\alpha}^{2}}I_{n}\end{bmatrix},$
(18)

and

$\underline{\underline{\boldsymbol{M}}}_{H,n}^{o}=\frac{1}{\eta_{0}}\begin{bmatrix}\frac{K_{n}}{K_{n}^{\prime}}\frac{n\bar{\beta}\bar{\alpha}_{o}}{j\bar{k}_{0}}I_{n}\Delta&-\sqrt{\epsilon_{r}}\frac{\alpha_{o}}{\alpha_{i}}\frac{I_{n}^{\prime}K_{n}}{K_{n}^{\prime}}\\\
\frac{j\bar{k}_{0}}{\bar{\alpha}_{o}}\frac{I_{n}K_{n}^{\prime}}{K_{n}}-\frac{(n\bar{\beta})^{2}}{j\bar{k}_{0}\eta_{0}\bar{\alpha}_{o}}\frac{K_{n}I_{n}}{K_{n}^{\prime}}\Delta&\frac{n\bar{\beta}\sqrt{\epsilon_{r}}}{\bar{\alpha}_{o}\bar{\alpha}_{i}}\frac{I_{n}^{\prime}K_{n}}{K_{n}^{\prime}}\end{bmatrix}.$
(19)

Here
$\Delta=\left(\frac{1}{\bar{\alpha}_{i}^{2}}-\frac{1}{\bar{\alpha}_{o}^{2}}\right)$
which is indicative of the material contrast between the core and the
surroundings, $I_{n}=I_{n}(\alpha_{i}R),K_{n}=K_{n}(\alpha_{o}R)$ and
$I_{n}^{\prime}=dI_{n}(x)/dx|_{x=\alpha_{i}R},K_{n}^{\prime}=dK_{n}(x)/dx|_{x=\alpha_{o}R}$.
The parameters $\eta,k$ are the intrinsic impedance and wavenumber of the
dielectric medium, and $k_{0},\eta_{0}$ are those of vacuum. One can see that
when substituting $\epsilon_{r}=1$ we get
$\alpha_{i}=\alpha_{o},\eta=\eta_{0}$, and we revert back to equations (7).
The solutions are obtained in the same way as in the previous case, by
substituting the Graphene conductivity and
$\underline{\underline{\boldsymbol{M}}}_{E,n},\underline{\underline{\boldsymbol{M}}}_{H,n}$
into equation (29).

Figure 5(a) shows the dispersion relation for a dielectric core with
$\epsilon_{r}=[2,3,10]$. The first difference we notice with respect to the
all-vacuum case is that new modes emerge that did not exist previously. These
modes are much less confined to the tube located in the shaded region in the
top-left corner of the dispersion plot. We see almost no ”split” in their
wavenumber since the split is induced by the properties of the angularly
varying conductivity, and being less confined, their interaction with the
Graphene shell is much weaker. As expected, increasing the dielectric constant
of the core gives larger values of $\alpha$ and shorter wavelengths of the
guided waves. In figure 5(b), we see the contribution of the different angular
components. On the top, we see the new, low-confinement mode. Due to its weak
interaction with the Graphene layer, there is much less coupling between
different values of $n$ (we see mainly the $\pm 1$ components), and the mode
is much more balanced in the TE/TM sense. In the more confined modes (bottom
two panels in (b)), the interaction with the Graphene is stronger, which leads
to a more pronounced TM nature and stronger coupling between different $n$.
Figure 5 shows the confinement as a function of the dielectric core
permittivity. We see that increasing $\epsilon_{r}$ may result in stronger
isolation, but there is a tradeoff since increasing $\epsilon_{r}$ increases
the mode confinement, which makes the propagating waves more prone to losses
induced by the Graphene layer.

##  VI Conclusions

In this work, we have analyzed the wave propagation and the nonreciprocal
characteristics of surface waves guided on a perpendicularly magnetized
Graphene tube. This magnetization scheme is more realistic than the previously
studied, but it results in a system with a nonuniform cross-section. The
resulting propagating modes form a mixture of TE/TM components and several
dominant values of OAM. The geometrical asymmetry results in more pronounced
nonreciprocal effects, such as a different cross-sectional field distribution
for oppositely propagating modes. These were leveraged for a directional
excitation of waves using any source as a function of the source placement
around the tube. Generalizing the analysis for a Graphene sheet wrapped on a
dielectric core, we see new modes emerge, poorly confined and weakly
interacting with the Graphene sheet, in addition to the highly confined
surface waves. When exciting the waves with a point source, the dielectric
core can yield more pronounced isolation, with demonstrated ratios of around
$\approx 3$. This can be further improved by adding additional elements, such
as an inhomogeneous core or structuring of the graphene later, which is left
for future work. The continuous dependence of the isolation ratio on the
direction of the magnetization lets us tailor the power propagation properties
using the magnetization direction as an easy control knob without needing to
solve for a different magnetization intensity or redesign the geometry and
materials. This component can have several applications, such as a building
block for miniature isolators and circulators in the THz range, control of
signal flow, and enhance light-matter interaction in graphene-wrapped
microfiber systems [18], and as a nonreciprocal antenna [23].

##  Appendix A Modal field derivation

Given the $\hat{z}$ components of the electric and magnetic fields in the TM
and TE components respectively, we can use the following equations to derive
the rest of the field components in cylindrical coordinates [1]

$E_{r}=\frac{j}{\alpha^{2}}\left(\beta\frac{\partial E_{z}}{\partial
r}+\frac{\omega\mu}{r}\frac{\partial H_{z}}{\partial\varphi}\right)$ (20)
$E_{\varphi}=\frac{j}{\alpha^{2}}\left(\frac{\beta}{r}\frac{\partial
E_{z}}{\partial\varphi}-\omega\mu\frac{\partial H_{z}}{\partial r}\right)$
(21) $H_{r}=-\frac{j}{\alpha^{2}}\left(\frac{\omega\epsilon}{r}\frac{\partial
E_{z}}{\partial\varphi}-\beta\frac{\partial H_{z}}{\partial r}\right)$ (22)
$H_{\varphi}=\frac{j}{\alpha^{2}}\left(\omega\epsilon\frac{\partial
E_{z}}{\partial r}+\frac{\beta}{r}\frac{\partial
H_{z}}{\partial\varphi}\right)$ (23)

Substituting equations (1a) and (1b) into these relations and applying the
continuity boundary condition in equation (4) yields

$I_{n}A_{n}^{i}=K_{n}A_{n}^{o}.$ (24)

The continuity of $E_{\varphi}$ gives

$\frac{n\beta}{\alpha^{2}R}\left[I_{n}A_{n}^{i}-K_{n}A_{n}^{o}\right]-\frac{jk_{0}}{\alpha}\left[I_{n}^{\prime}B_{n}^{i}-K_{n}^{\prime}B_{n}^{o}\right]=0$
(25)

from which we get

$I_{n}^{\prime}B_{n}^{i}=K_{n}^{\prime}B_{n}^{o}$ (26)

In Eq. (5), we have defined a truncated coefficient column vector
$\boldsymbol{D}_{n}$. Using these relations, the representation of the guided
wave fields can be written as

$E_{z}^{i}(r=R)=E_{z}^{o}(r=R)=\left[I_{n},0\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27a)
$E_{\varphi}^{i}(r=R)=E_{\varphi}^{o}(r=R)=\left[\frac{n\beta}{\alpha^{2}R},-\frac{jk_{0}}{\alpha}I_{n}^{\prime}\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27b)
$H_{z}^{i}(r=R)=\left[0,\frac{1}{\eta_{0}}I_{n}\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27c)
$H_{z}^{o}(r=R)=\left[0,\frac{1}{\eta_{0}}\frac{I_{n}^{\prime}}{K_{n}^{\prime}}K_{n}\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27d)
$H_{\varphi}^{i}(r=R)=\left[\frac{jk_{0}}{\alpha\eta_{0}}I_{n}^{\prime},\frac{\beta
n}{\alpha^{2}R\eta_{0}}I_{n}\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27e)
$H_{\varphi}^{o}(r=R)=\left[\frac{jk_{0}}{\alpha\eta_{0}}\frac{I_{n}}{K_{n}}K_{n}^{\prime},\frac{\beta
n}{\alpha^{2}R\eta_{0}}\frac{I_{n}^{\prime}}{K_{n}^{\prime}}K_{n}\right]\boldsymbol{D}_{n}e^{-jn\varphi}e^{-j\beta
z}$ (27f)

These relations can be transformed into a matrix representation which will be
instrumental to the convenient analysis of the guided modes, and using the
wronskian identity

$I_{n}^{\prime}(x)K_{n}(x)-K_{n}^{\prime}(x)I_{n}(x)=\frac{1}{x}$ (28)

one can arrive at the definitions given in equation (7).

##  Appendix B Infinite matrix eigenmode equation

Equation 11 can be written in an infinite matrix form

$\begin{pmatrix}\ddots&\vdots&\vdots&\vdots&\ddots\\\
\ldots&\underline{\underline{\boldsymbol{Y}}}_{0}\underline{\underline{\boldsymbol{M}}}_{E,-1}-\underline{\underline{\boldsymbol{M}}}_{H,-1}&\underline{\underline{\boldsymbol{Y}}}_{-1}\underline{\underline{\boldsymbol{M}}}_{E,0}&\underline{\underline{\boldsymbol{Y}}}_{-2}\underline{\underline{\boldsymbol{M}}}_{E,1}&\ldots\\\
\ldots&\underline{\underline{\boldsymbol{Y}}}_{1}\underline{\underline{\boldsymbol{M}}}_{E,-1}&\underline{\underline{\boldsymbol{Y}}}_{0}\underline{\underline{\boldsymbol{M}}}_{E,0}-\underline{\underline{\boldsymbol{M}}}_{H,0}&\underline{\underline{\boldsymbol{Y}}}_{-1}\underline{\underline{\boldsymbol{M}}}_{E,1}&\ldots\\\
\ldots&\underline{\underline{\boldsymbol{Y}}}_{2}\underline{\underline{\boldsymbol{M}}}_{E,-1}&\underline{\underline{\boldsymbol{Y}}}_{1}\underline{\underline{\boldsymbol{M}}}_{E,0}&\underline{\underline{\boldsymbol{Y}}}_{0}\underline{\underline{\boldsymbol{M}}}_{E,1}-\underline{\underline{\boldsymbol{M}}}_{H,1}&\ldots\\\
\ddots&\vdots&\vdots&\vdots&\ddots\end{pmatrix}\begin{pmatrix}\vdots\\\
\boldsymbol{D}_{-2}\\\ \boldsymbol{D}_{-1}\\\ \boldsymbol{D}_{0}\\\
\boldsymbol{D}_{1}\\\ \boldsymbol{D}_{2}\\\ \vdots\end{pmatrix}=0.$ (29)

When a conductivity of the form given in equation (12) is used, this matrix
reduces to a tridiagonal block matrix.

## References

  * [1] R. Collin, Foundations for microwave engineering.  1992\.
  * [2] J. D. Kraus and R. J. Marhefka, Antennas for all applications.  McGraw-Hill series in electrical engineering, New York: McGraw-Hill, 3rd ed ed., 2002.
  * [3] H. D. Levine and C. H. Papas, “Theory of the circular diffraction antenna,” Journal of Applied Physics, vol. 19, pp. 29–43, 1951.
  * [4] Y. R. Padooru, A. B. Yakovlev, P.-Y. Chen, and A. Alù, “Analytical modeling of conformal mantle cloaks for cylindrical objects using sub-wavelength printed and slotted arrays,” Journal of Applied Physics, vol. 112, 08 2012.  034907\.
  * [5] C.-W. Lin and A. Grbic, “Analysis and synthesis of cascaded cylindrical metasurfaces using a wave matrix approach,” IEEE Transactions on Antennas and Propagation, vol. 69, no. 10, pp. 6546–6559, 2021.
  * [6] F. Alpeggiani, K. Y. Bliokh, F. Nori, and L. Kuipers, “Electromagnetic helicity in complex media,” Phys. Rev. Lett., vol. 120, p. 243605, Jun 2018.
  * [7] J. R. Pierce and P. K. Tien, “Coupling of modes in helixes,” Proceedings of the IRE, vol. 42, no. 9, pp. 1389–1396, 1954.
  * [8] N. Raveu, B. Byrne, L. Claudepierre, and N. Capet, “Modal theory for waveguides with anisotropic surface impedance boundaries,” IEEE Transactions on Microwave Theory and Techniques, vol. 64, no. 4, pp. 1153–1162, 2016.
  * [9] J. G. Pollock, A. K. Iyer, D. Pratap, and S. Anantha Ramakrishna, “A class of circular waveguiding structures containing cylindrically anisotropic metamaterials: Applications from radio frequency/microwave to optical frequencies,” Journal of Applied Physics, vol. 119, 02 2016.  083103\.
  * [10] Y. Mazor and A. Alù, “Angular-momentum selectivity and asymmetry in highly confined wave propagation along sheath-helical metasurface tubes,” Physical Review B, vol. 99, no. 15, p. 155425, 2019.
  * [11] C. J. M. Barker, N. D. Zanche, and A. K. Iyer, “Dispersion and polarization control in below-cutoff circular waveguides using anisotropic metasurface liners,” IEEE Transactions on Microwave Theory and Techniques, pp. 1–12, 2023.
  * [12] G. Y. Slepyan, S. A. Maksimenko, A. Lakhtakia, O. Yevtushenko, and A. V. Gusakov, “Electrodynamics of carbon nanotubes: Dynamic conductivity, impedance boundary conditions, and surface wave propagation,” Phys. Rev. B, vol. 60, pp. 17136–17149, Dec 1999.
  * [13] J. long Kou, J. hui Chen, Y. Chen, F. Xu, and Y. qing Lu, “Platform for enhanced light–graphene interaction length and miniaturizing fiber stereo devices,” Optica, vol. 1, pp. 307–310, Nov 2014.
  * [14] K. Chen, X. Zhou, X. Cheng, R. Qiao, Y. Cheng, C. Liu, Y. Xie, W. Yu, F. Yao, Z. Sun, et al., “Graphene photonic crystal fibre with strong and tunable light–matter interaction,” Nature Photonics, vol. 13, no. 11, pp. 754–759, 2019.
  * [15] Y. Wu, B. C. Yao, Q. Y. Feng, X. L. Cao, X. Y. Zhou, Y. J. Rao, Y. Gong, W. L. Zhang, Z. G. Wang, Y. F. Chen, and K. S. Chiang, “Generation of cascaded four-wave-mixing with graphene-coated microfiber,” Photon. Res., vol. 3, pp. A64–A68, Apr 2015.
  * [16] J.-H. Chen, B.-C. Zheng, G.-H. Shao, S.-J. Ge, F. Xu, and Y.-Q. Lu, “An all-optical modulator based on a stereo graphene–microfiber structure,” Light: Science & Applications, vol. 4, no. 12, pp. e360–e360, 2015.
  * [17] I. D. Barcelos, L. G. Moura, R. G. Lacerda, and A. Malachias, “Observation of strain-free rolled-up cvd graphene single layers: Toward unstrained heterostructures,” Nano letters, vol. 14, no. 7, pp. 3919–3924, 2014.
  * [18] J.-h. Chen, Y.-f. Xiong, F. Xu, and Y.-q. Lu, “Silica optical fiber integrated with two-dimensional materials: towards opto-electro-mechanical technology,” Light: Science & Applications, vol. 10, no. 1, p. 78, 2021.
  * [19] C. Caloz, A. Alù, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Léger, “Electromagnetic nonreciprocity,” Phys. Rev. Appl., vol. 10, p. 047001, Oct 2018.
  * [20] A. Kord, D. L. Sounas, and A. Alù, “Microwave nonreciprocity,” Proceedings of the IEEE, vol. 108, no. 10, pp. 1728–1758, 2020.
  * [21] Z. Yu, G. Veronis, Z. Wang, and S. Fan, “One-way electromagnetic waveguide formed at the interface between a plasmonic metal under a static magnetic field and a photonic crystal,” Phys. Rev. Lett., vol. 100, p. 023902, Jan 2008.
  * [22] Y. Mazor, Y. Hadad, and B. Z. Steinberg, “Planar one-way guiding in periodic particle arrays with asymmetric unit cell and general group-symmetry considerations,” Phys. Rev. B, vol. 92, p. 125129, Sep 2015.
  * [23] D. Correas-Serrano, J. S. Gomez-Diaz, A. Alù, and A. Á. Melcón, “Electrically and magnetically biased graphene-based cylindrical waveguides: analysis and applications as reconfigurable antennas,” IEEE Transactions on Terahertz Science and Technology, vol. 5, no. 6, pp. 951–960, 2015.
  * [24] A. Oliner and A. Hessel, “Guided waves on sinusoidally-modulated reactance surfaces,” IRE Transactions on Antennas and Propagation, vol. 7, no. 5, pp. 201–208, 1959.
  * [25] C. Li, J.-h. Chen, S.-c. Yan, F. Xu, and Y.-q. Lu, “A fiber laser using graphene-integrated 3-d microfiber coil,” IEEE Photonics Journal, vol. 8, no. 1, pp. 1–7, 2016.
  * [26] D. L. Sounas and C. Caloz, “Gyrotropy and nonreciprocity of graphene for microwave applications,” IEEE Transactions on Microwave Theory and Techniques, vol. 60, no. 4, pp. 901–914, 2012.
  * [27] O. Y. Yermakov, A. I. Ovcharenko, M. Song, A. A. Bogdanov, I. V. Iorsh, and Y. S. Kivshar, “Hybrid waves localized at hyperbolic metasurfaces,” Phys. Rev. B, vol. 91, p. 235423, Jun 2015.
  * [28] Y. Mazor and B. Z. Steinberg, “Metaweaves: Sector-way nonreciprocal metasurfaces,” Phys. Rev. Lett., vol. 112, p. 153901, Apr 2014.
  * [29] COMSOL Multiphysics® v. 6.2, COMSOL AB, Stockholm, Sweden.


# Quench dynamics and bulk-edge correspondence in nonlinear mechanical systems

Motohiko Ezawa  Department of Applied Physics, University of Tokyo, Hongo
7-3-1, 113-8656, Japan

###### Abstract

We study a topological physics in a one-dimensional nonlinear system by taking
an instance of a mechanical rotator model with alternating spring constants.
This nonlinear model is smoothly connected to an acoustic model described by
the Su-Schrieffer-Heeger model in the linear limit. We numerically show that
quench dynamics of the kinetic and potential energies for the nonlinear model
is well understood in terms of the topological and trivial phases defined in
the associated linearized model. It indicates phenomenologically the emergence
of the edge state in the topological phase even for the nonlinear system,
which may be the bulk-edge correspondence in nonlinear system.

##  I Introduction

Topological insulatorsHasan ; Qi  were originally discovered in in-organic
solid state materials. However, it is now well recognized that topological
physics is ubiquitous in various systems such as acousticProdan ; TopoAco ;
Berto ; Xiao ; He ; Abba ; Xue ; Ni ; Wei ; Xue2 , mechanicalLubensky ; Chen ;
Nash ; Paul ; Sus ; Sss ; Huber ; Mee ; Kariyado ; Hannay ; Po ; Rock ;
Takahashi ; Mat ; Taka ; Ghatak ; Wakao , photonicKhaniPhoto ; Hafe2 ; Hafezi
; WuHu ; TopoPhoto ; Ozawa ; Hassan ; Li  and electric circuitTECNature ;
ComPhys ; Hel ; Lu ; YLi ; EzawaTEC ; Research ; Zhao ; EzawaLCR ; EzawaSkin ;
Garcia ; Hofmann ; EzawaMajo ; Tjunc ; Lee ; Kot  systems. They are called
artificial topological systems. The merit of them is that it is possible to
fabricate ideal systems comparing to natural solid state materials. Another
merit of artificial topological systems is that nonlinearity is naturally
introduced into them.

Topological physics is mostly studied in linear systems. There are only a few
studies on it in nonlinear systemsKot ; Smi ; Sone ; TopoToda  because the
study of the topological properties is not straightforward. One of the reasons
is that it is a formidable problem to obtain a band structure and a
topological number in a generic nonlinear system. Recently, it is proposed to
account for the topological properties in nonlinear systems phenomenologically
based on the bulk-edge correspondence well established in the linear
theoryTopoToda . It seems to work provided that the nonlinear theory is
continuously connected to a linear theory where the topological number is well
defined. One may say that the topological properties are inherited from a
linear theory to a nonlinear theory. It is an interesting problem to explore
other systems which share similar properties.

A good signal to detect the topological phase transition is to study quench
dynamics starting from one of the edges in the case of one dimension based on
the bulk-edge correspondenceQWalk . There is almost no time evolution and the
state remains almost localized at the edge for a topological phase. It is
because the initial state is almost given by the topological localized edge
state, which has no dynamics. On the other hand, the state rapidly spreads
into the bulk for a trivial phase because the initial state is composed of
bulk eigen functions. Although the usefulness of quench dynamics has been
established in a linear system, it is also useful in a nonlinear
systemTopoToda .

In this paper, we investigate a one-dimensional mechanical rotator model with
alternating spring constants: See Fig.1. It is a nonlinear generalization of
the Su-Schrieffer-Heeger (SSH) model. We solve the quench dynamics of a
mechanical rotator model as an initial condition problem, where only one
rotator at the left-most edge is excited initially. The potential and kinetic
energies exhibit distinct behaviors depending on whether the system is in the
topological or trivial phase defined in the linearized model. Namely, after
enough time, they are well localized in the topological phase, while they are
spread over the bulk in the trivial phase. These phenomena are understood in
terms of the emergence of the topological edge state in the topological phase.
It would represent the bulk-edge correspondence in nonlinear system.

Figure 1: Illustration of the mechanical rotator model. Two adjacent rotators
are connected by a spring with the spring constant $\kappa_{j}$. A rotator
colored in green rotates around a center colored in orange with the radius
$\ell_{0}$. All centers are on a line, and the distance between two adjacent
centers is $L_{0}$.

##  II Mechanical rotator model

We consider a mechanical rotator modelLubensky  as illustrated in Fig.1. We
place rotators indexed by the number $j$, $j=1,2,\cdots,N$. Each rotator has
the radius $\ell_{0}$, and its center is fixed at the position
$\left(jL_{0},0\right)$ on the $x$ axis, with $L_{0}$ the distance between two
adjacent centers. The dynamical variables are angles $\theta_{j}$ for rotators
$j$.

The Hamiltonian of the system consists of the kinetic energy $U^{\text{K}}$,
the potential energy of the springs $U^{\text{spring}}$ and the gravitational
energy of the $U^{\text{g}}$,

$H=U^{\text{K}}+U^{\text{spring}}+U^{\text{g}}.$ (1)

They are given by

$U^{\text{K}}=\sum_{j}U_{j}^{\text{K}},\quad
U_{j}^{\text{K}}=\frac{m}{2}\dot{\theta}_{j}^{2}$ (2)

with the mass $m$,

$U^{\text{g}}=-g\sum_{j}\cos\theta_{j}$ (3)

with the gravitational constant $g$, and

$\displaystyle U^{\text{spring}}$ $\displaystyle=$
$\displaystyle\sum_{j}U_{j}^{\text{spring}},$ (4) $\displaystyle
U_{j}^{\text{spring}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\kappa_{j}\left(L_{j}\left(\theta_{j},\theta_{j+1}\right)-L\right)^{2},$
(5)

where $L_{j}\left(\theta_{j},\theta_{j+1}\right)$ is the length of the spring
between $j$ and $j+1$ nodes,

$L_{j}\left(\theta_{j},\theta_{j+1}\right)=\sqrt{\left[L_{j}^{x}\left(\theta_{j},\theta_{j+1}\right)\right]^{2}+\left[L_{j}^{y}\left(\theta_{j},\theta_{j+1}\right)\right]^{2}}$
(6)

with

Figure 2: Illustration of the trivial and topological phases in the SSH model.
There are two isolated edge states in the topological phase, while all of the
sites are dimerized in the trivial phase. Magenta and cyan lines represent
strong and weak bondings $\kappa_{j}$ in Eq.(10), respectively. This
illustration is argued to hold even in nonlinear system described by Eq.(9).

$\displaystyle L_{j}^{x}\left(\theta_{j},\theta_{j+1}\right)$ $\displaystyle=$
$\displaystyle\ell_{0}\sin\theta_{j+1}-\ell_{0}\sin\theta_{j}+L_{0},$ (7)
$\displaystyle L_{j}^{y}\left(\theta_{j},\theta_{j+1}\right)$ $\displaystyle=$
$\displaystyle\ell_{0}\cos\theta_{j+1}-\ell_{0}\cos\theta_{j}.$ (8)

The equation of motion is given by

$m\ddot{\theta}_{j}=-g\sin\theta_{j}-\frac{\partial
U^{\text{spring}}}{\partial\theta_{j}}.$ (9)

The spring constant is assumed to be alternating,

$\kappa_{j}=\kappa\left(1+\lambda\left(-1\right)^{j}\right),$ (10)

where the dimerization is controlled by $\lambda$ with $|\lambda|\leq 1$. For
$\lambda<0$, the spring constant $\kappa_{j}$ with odd (even) $j$ is strong
(weak). On the other hand, for $\lambda>0$, the spring constant $\kappa_{j}$
with odd (even) $j$ is weak (strong): See Fig.2.

Figure 3: Time evolution of the angle $\sin\theta_{j}$ for (a1)$\sim$(a8) and
(d1)$\sim$(d8), the potential energy $U_{j}^{\text{spring}}$ for
(b1)$\sim$(b8) and (e1)$\sim$(e8), and the kinetic energy $U_{j}^{\text{K}}$
for (c1)$\sim$(c8) and (f1)$\sim$(f8). (a1)$\sim$(c8) topological phase with
$\lambda=-0.5$, and (d1)$\sim$(f8) trivial phase with $\lambda=0.5$.
(a1)$\sim$(f2) $\theta_{0}=\pi/10$, (a3)$\sim$(f4) $\theta_{0}=\pi/4$,
(a5)$\sim$(f6) $\theta_{0}=\pi/2$, and (a7)$\sim$(f8) $\theta_{0}=3\pi/4$.
(a1)$\sim$(f1), (a3)$\sim$(f3), (a5)$\sim$(f5) and (a7)$\sim$(f7) Bird’s eye’s
view, and (a2)$\sim$(f2), (a4)$\sim$(f4), (a6)$\sim$(f6) and (a8)$\sim$(f8)
Color plot of the time evolution. We have used a chain containing 100 sites.
We have set $m=0.01$, $g=0.1$, $L_{0}=10$ and $\ell_{0}=1$. Color indicates
the site index $j$, which is manifest in bird’s eye’s views.

##  III Linearized model

Provided the angle $\theta_{j}$ is small enough, the potential energy is well
approximated by the harmonic potential,

$U^{\text{spring}}=\frac{\ell_{0}^{2}}{2}\sum_{j}\kappa_{j}\left(\theta_{j+1}-\theta_{j}\right)^{2},$
(11)

and the equation of motion is obtained as

$m\ddot{\theta}_{j}=-g\theta_{j}-\ell_{0}^{2}\sum_{j}\kappa_{j}\left(\theta_{j}-\theta_{j+1}\right)+\kappa_{j-1}\left(\theta_{j}-\theta_{j-1}\right),$
(12)

which is a linearized model.

Figure 4: (a1)$\sim$(a4) Kinetic energy at the edge $U_{1}^{\text{K}}$, and
(b1)$\sim$(b4) the potential energy at the edge $U_{1}^{\text{spring}}$ as a
function of the initial position $\theta_{0}$ after enough time. (a1) and (b1)
$\theta_{0}=\pi/10$, (a2) and (b2) $\theta_{0}=\pi/4$, (a3) and (b3)
$\theta_{0}=\pi/2$, (a4) and (b4) $\theta_{0}=3\pi/4$. See the caption of
Fig.3 for the parameters.

Eq.(12) is rewritten in the form of

$m\ddot{\theta}_{j}=\sum_{j}K_{ij}\theta_{j},$ (13)

where

$K_{ij}=-\delta_{ij}\left(g+\ell^{2}\right)+\ell_{0}^{2}\left(\kappa_{j}\delta_{i,j+1}+\kappa_{j-1}\delta_{i,j-1}\right)$
(14)

is identical to the SSH Hamiltonian. After a Fourier transformation, we have

$K\left(k\right)=-\left(g+\ell_{0}^{2}\right)I_{2}+\kappa\left(\begin{array}[]{cc}0&1-\lambda
e^{-ik}\\\ 1-\lambda e^{ik}&0\end{array}\right).$ (15)

It is knownLubensky ; Chen  that the system is topological for $\lambda<0$ and
trivial for $\lambda>0$. There are two isolated edge states in the limit
$\lambda\simeq-1$, while all of the states are dimerized in the limit
$\lambda\simeq 1$: See Fig.2.

The topological number associated with the SSH Hamiltonian is the chiral index
defined by

$\Gamma=\int_{0}^{2\pi}\text{Tr}\left[\sigma_{z}K^{-1}\partial_{k}K\right]dk,$
(16)

where $K$ is given by Eq.(15). We obtain $\Gamma=1$ for $\lambda<0$ and
$\Gamma=0$ for $\lambda>0$.

Our interest is the effect of the nonlinearity based on Eq.(9) instead of the
linear equation of motion (12).

##  IV Nonlinear quench dynamics and bulk-edge correspondence

The quench dynamics starting from a state localized at one edge well captures
the topological phase transition in linear systemsQWalk , where the
eigenfunctions are easily obtained and a topological number is well defined.
In the topological phase, there are topological edge states well localized at
edges, which is known as the bulk-edge correspondence. If we excite only one
edge site, the most component is dominated by a topological edge state. The
topological edge state remains as it is after time evolution. On the other
hand, there is no such localized edge state in the trivial phase. Hence, all
of the components of the initial state are bulk eigenfunctions. They are
spread into the bulk after the time evolution. Thus, it is possible to
distinguish topological and trivial phases by checking whether the state is
localized or spread into the bulk.

The above observation is also applicable even for nonlinear systemsTopoToda .
The existence of the localized topological edge state is obscure in nonlinear
systems because it is not possible to diagonalize the Hamiltonian and obtain
eigenfunctions. Nevertheless, the quench dynamics shows distinct behaviors
between the topological and trivial phases as in the case of the linear
system. Thus, the quench dynamics starting from one of the edges is a strong
signal to detect a topological phase transition.

We numerically solve the equation of motion (9) with the initial condition
$\theta_{j}\left(t=0\right)=\theta_{0}\delta_{j,1}$ and
$\dot{\theta}_{j}\left(t=0\right)=0$. If $\left|\theta_{0}\right|\ll 1$, the
system is well described by the linear equation (12). Otherwise, the system is
nonlinear. We study a case $\theta_{0}=\pi/10,\pi/4,\pi/2$ and $3\pi/4$ as
typical examples.

We show the time evolution of $\sin\theta_{j}$, the potential energy
$U_{j}^{\text{spring}}$ and the kinetic energy $U_{j}^{\text{K}}$ in Fig.3.
The quench dynamics of $\sin\theta_{j}$ is significantly different between the
topological and trivial phases. The amplitude of $\sin\theta_{1}$ at the left
edge site remains finite in the topological phase. On the other hand, the
amplitude of $\sin\theta_{1}$ decreases in the trivial phase. The potential
energy $U_{j}^{\text{spring}}$ remains well localized at the left edge in the
topological phase although the absolute value of $\sin\theta_{j}$ spreads into
the bulk as shown in Fig.3(b). On the other hand, it moves as if it were a
soliton in the trivial phase as shown in Fig.3(e). If $\theta_{0}=\pi/2$,
there remains a finite component in $U_{j}^{\text{spring}}$ as well as a
soliton-like propagation. This will be because that $\sin\theta_{j}$ takes a
maximum for $\theta_{0}=\pi/2$, where the dynamics is most enhanced and it
takes longer time to reach a steady state.

We also show the time evolution of the kinetic energy $U_{j}^{\text{K}}$ in
Fig.3(c) and (f). The behavior of $U_{j}^{\text{K}}$ is quite similar to that
of $U_{j}^{\text{spring}}$.

The dynamics between $\theta_{0}=\pi/4$ and $3\pi/4$ are very similar as shown
in Fig.3. It will be due to the fact that $\sin\pi/4=\sin 3\pi/4$ although we
have $\cos\pi/4\neq\cos 3\pi/4$. These results indicate a symmetry between the
angle $\theta_{1}$ and $\theta_{2}$ satisfying the condition
$\sin\theta_{1}=\sin\theta_{2}$, where the dynamics is similar.

The potential and kinetic energies after enough time present a good signal for
the topological phase transition comparing to the dynamics of
$\sin\theta_{j}$. We show the kinetic energy $U_{j}^{\text{K}}$ and the
potential energy $U_{j}^{\text{spring}}$ after enough time as a function of
$\lambda$ in Fig.4. It is finite for the topological phase with $\lambda<0$,
while it is almost zero for the trivial phase with $\lambda>0$. These features
are common for all initial conditions $\theta_{0}=\pi/10,\pi/4,\pi/2$ and
$3\pi/4$. It indicates the bulk-edge correspondence in nonlinear systems. It
suggests the validity of the topological number (16) even in strong nonlinear
regime.

A comment is in order with respect to finite components around $\lambda\simeq
1$ in Fig.4. They are interpreted as follows. In the limit $\lambda\simeq 1$,
the system is almost dimerized as shown in Fig.2, where there are no isolated
sites at the both edges. In this limit, the energy is not well transferred to
the bulk because the energy is localized in the dimer located at the edge to
which the energy is injected, resulting in finite components around
$\lambda\simeq 1$.

##  V Discussion

We have shown that the topological and trivial phases are well differentiated
in the mechanical rotator model even in a strong nonlinear regime based on the
bulk-edge correspondence. Our finding is that the topological properties are
inherited to the nonlinear model from the associated linearized model provided
they are smoothly connected. Then, it would be possible to use a topological
number defined in the linearized model. This phenomenon is quite similar to
the one in the Toda latticeTopoToda , which is a typical exactly solvable
model containing a soliton.

The author is very much grateful to N. Nagaosa for helpful discussions on the
subject. This work is supported by the Grants-in-Aid for Scientific Research
from MEXT KAKENHI (Grants No. JP17K05490 and No. JP18H03676). This work is
also supported by CREST, JST (JPMJCR16F1 and JPMJCR20T2).

## References

  * (1) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * (4) Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong and B. Zhang, Phys. Rev. Lett. 114, 114301 (2015).
  * (5) P. Wang, L. Lu and K. Bertoldi, Phys. Rev. Lett. 115, 104302 (2015).
  * (6) M. Xiao, G. Ma, Z. Yang, P. Sheng, Z. Q. Zhang and C. T. Chan, Nat. Phys. 11, 240 (2015).
  * (7) C. He, X. Ni, H. Ge, X.-C. Sun,Y.-B. Chen1 M.-H. Lu, X.-P. Liu, L. Feng and Y.-F. Chen, Nature Physics 12, 1124 (2016).
  * (8) H. Abbaszadeh, A. Souslov, J. Paulose, H. Schomerus and V. Vitelli, Phys. Rev. Lett. 119, 195502 (2017).
  * (9) H. Xue, Y. Yang, F. Gao, Y. Chong and B.Zhang, Nature Materials 18, 108 (2019).
  * (10) X. Ni, M. Weiner, A. Alu and A. B. Khanikaev, Nature Materials 18, 113 (2019).
  * (11) M. Weiner, X. Ni, M. Li, A. Alu, A. B. Khanikaev, Science Advances 6, eaay4166 (2020)
  * (12) H. Xue, Y. Yang, G. Liu, F. Gao, Y. Chong and B. Zhang, Phys. Rev. Lett. 122, 244301 (2019).
  * (13) C. L. Kane and T. C. Lubensky, Nature Phys. 10, 39 (2014).
  * (14) B. Gin-ge Chen, N. Upadhyaya and V. Vitelli, PNAS 111, 13004 (2014)
  * (15) L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner and W. T. M. Irvine, PNAS 112, 14495 (2015).
  * (16) J. Paulose, A. S. Meeussen and V. Vitelli, PNAS 112, 7639 (2015)
  * (17) R. Susstrunk, S. D. Huber, Science 349, 47 (2015).
  * (18) R. Susstrunk and S. D. Huber, Proc. Natl. Acad. Sci. USA 113, E4767 (2016).
  * (19) S. D. Huber, Nature Physics 12, 621 (2016).
  * (20) A. S. Meeussen, J. Paulose and V. Vitelli, Phys. Rev. X 6, 041029 (2016).
  * (21) T. Kariyado and Y. Hatsugai, Sci. Rep. 5, 18107 (2016).
  * (22) T. Kariyado and Y. Hatsugai, J. Phys. Soc. Jpn. 85, 043001 (2016).
  * (23) H. C. Po, Y. Bahri and A. Vishwanath, Phys. Rev. B 93, 205158 (2016).
  * (24) D. Zeb Rocklin, Bryan Gin–ge Chen, Martin Falk, Vincenzo Vitelli, and T. C. Lubensky, Phys. Rev. Lett. 116, 135503 (2016)
  * (25) Y. Takahashi, T. Kariyado and Y. Hatsugai, New J. Phys. 19, 035003 (2017).
  * (26) K. H. Matlack, M. Serra-Garcia, A. Palermo, S. D. Huber and C. Daraio, Nature Mat. 17, 323 (2018).
  * (27) Y. Takahashi, T/ Kariyado and Y. Hatsugai, Phys. Rev. B 99, 024102 (2019).
  * (28) A. Ghatak, M. Brandenbourger, J. van Wezel and C. Coulais, Proc. Natl. Ac. Sc. U.S.A. 117, 29561 (2020)
  * (29) H. Wakao, T. Yoshida, H. Araki, T. Mizoguchi and Y. Hatsugai, Phys. Rev. B 101, 094107 (2020).
  * (30) A. B. Khanikaev, S. H. Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, G. Shvets, Nature Materials 12, 233 (2013).
  * (31) M. Hafezi, E. Demler, M. Lukin, J. Taylor, Nature Physics 7, 907 (2011).
  * (32) M. Hafezi, S. Mittal, J. Fan, A. Migdall, J. Taylor, Nature Photonics 7, 1001 (2013).
  * (33) L.H. Wu and X. Hu, Phys. Rev. Lett. 114, 223901 (2015).
  * (34) L. Lu. J. D. Joannopoulos and M. Soljacic, Nature Photonics 8, 821 (2014).
  * (35) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg and L. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * (36) A. E. Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, M. Bourennane, Nature Photonics 13, 697 (2019)
  * (37) M. Li, D. Zhirihin, D. Filonov, X. Ni, A. Slobozhanyuk, A. Alu and A. B. Khanikaev, Nature Photonics 14, 89 (2020)
  * (38) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, R. Thomale, Nat. Phys. 14, 925 (2018).
  * (39) C. H. Lee , S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling and R. Thomale, Communications Physics, 1, 39 (2018).
  * (40) T. Helbig, T. Hofmann, C. H. Lee, R. Thomale, S. Imhof, L. W. Molenkamp and T. Kiessling, Phys. Rev. B 99, 161114 (2019).
  * (41) Y. Lu, N. Jia, L. Su, C. Owens, G. Juzeliunas, D. I. Schuster and J. Simon, Phys. Rev. B 99, 020302 (2019).
  * (42) Y. Li, Y. Sun, W. Zhu, Z. Guo, J. Jiang, T. Kariyado, H. Chen and X. Hu, Nat. Com. 9, 4598 (2018)
  * (43) M. Ezawa, Phys. Rev. B 98, 201402(R) (2018).
  * (44) K. Luo, R. Yu and H. Weng, Research, ID 6793752. (2018)
  * (45) E. Zhao, Ann. Phys. 399, 289 (2018).
  * (46) M. Ezawa, Phys. Rev. B 99, 201411(R) (2019).
  * (47) M. Ezawa, Phys. Rev. B 99, 121411(R) (2019).
  * (48) M. Serra-Garcia, R. Susstrunk and S. D. Huber, Phys. Rev. B 99, 020304 (2019).
  * (49) T. Hofmann, T. Helbig, C. H. Lee, M. Greiter, R. Thomale, Phys. Rev. Lett. 122, 247702 (2019).
  * (50) M. Ezawa, Phys. Rev. B 100, 045407 (2019).
  * (51) M. Ezawa, Phys. Rev. B 102, 075424 (2020).
  * (52) C. H. Lee, T. Hofmann, T. Helbig, Y. Liu, X. Zhang, M. Greiter and R. Thomale, Nature Communications 11, 4385 (2020).
  * (53) T. Kotwal, H. Ronellenfitsch, F. Moseley, A. Stegmaier, R. Thomale, J. Dunkel, PNAS 118, e2106411118 (2021)
  * (54) D. Smirnova, D. Leykam, Y. Chong and Y. Kivshar, Applied Physics Reviews 7, 021306 (2020)
  * (55) K. Sone, Y. Ashida, T. Sagawa, arXiv:2012.09479
  * (56) M. Ezawa, cond-mat/arXiv:2105.10851
  * (57) M. Ezawa, Phys. Rev. B 100, 165419 (2019).