# Twisted topology and Bipolar Non-Hermitian Skin Effect induced by long-range
asymmetric coupling

S. M. Rafi-Ul-Islam<EMAIL_ADDRESS>Department of Electrical and
Computer Engineering, National University of Singapore, Singapore 117583,
Republic of Singapore  Zhuo Bin Siu<EMAIL_ADDRESS>Department of
Electrical and Computer Engineering, National University of Singapore,
Singapore 117583, Republic of Singapore  Haydar Sahin<EMAIL_ADDRESS>Department of Electrical and Computer Engineering, National University of
Singapore, Singapore 117583, Republic of Singapore  Institute of High
Performance Computing, A*STAR, Singapore 138632, Republic of Singapore  Md.
Saddam Hossain Razo<EMAIL_ADDRESS>Department of Electrical and Computer
Engineering, National University of Singapore, Singapore 117583, Republic of
Singapore  Mansoor B. A. Jalil<EMAIL_ADDRESS>Department of Electrical
and Computer Engineering, National University of Singapore, Singapore 117583,
Republic of Singapore

###### Abstract

We investigate the twisted topology of the complex eigenspectrum of a one-
dimensional non-Hermitian system under the influence of long-range
unidirectional coupling. Unlike the complex energy spectrum of the
conventional Hatano-Nelson chain, which takes the form of a single loop with a
topological winding index of a definite sign, the introduction of long-range
unidirectional hopping results in the creation of multiple twisted loops.
These twisted loops exhibit opposite signs of the topological winding index,
which correlate to alternating clockwise and anticlockwise energy windings.
The simultaneous presence of both signs of the winding index translates into a
bipolar non-Hermitian skin effect (NHSE), which challenges the conventional
wisdom that the NHSE localization is dependent on the direction of the
dominant nearest-neighbor interactions. In this bipolar NHSE, the exponents of
the complex energy eigenvectors corresponding to clockwise and anti-clockwise
windings, lie inside and outside of the complex unit circle, respectively.
Interestingly, at the intersections of oppositely oriented energy loops where
the sign of the topological winding index flips, the energy becomes real-
valued, leading to a suppression of the NHSE. This marks the emergence of
Bloch-like contact points, where both the bipolar NHSE and the traditional
NHSE vanish. Based on the non-Hermitian model we provide analytical insights
into the effects of long-range unidirectional coupling to the winding topology
of its complex energy spectra and their broader implications for the field of
condensed matter physics.

##  I Introduction

In recent years, topological phases [1, 2, 3, 4, 5, 6, 7, 8, 9] and non-
Hermitian physics [10, 11, 12, 13, 14, 15] have become rapidly growing fields
in condensed matter physics [16, 17, 18] and attracted the attention of
researchers owing to their unique properties [19, 20] and potential
applications in various areas [21, 22, 23, 24].

Fascinatingly, the investigation of non-Hermitian systems, rooted in the early
days of quantum mechanics, has flourished into a dynamic research field,
marked by notable advancements over the past few decades [25, 26, 27, 28, 29,
30]. These non-Hermitian systems manifest novel and exotic physical phenomena
absent in their Hermitian counterparts. In essence, non-Hermiticity in these
systems arises from coupling asymmetry or the presence of gain or loss terms
at the onsite level [31, 32, 33], contributing to a rich tapestry of diverse
and intriguing behaviors. The coupling asymmetry in non-Hermitian lattice
models lies at the origin of many non-Hermitian phenomena ranging from the
non-hermitian skin effect (NHSE)[31, 34, 35, 36, 37, 38], and exceptional
points [39, 40], which may be utilized in ultra-sensitive sensing [24, 41,
42], exponential signal enhancement [22, 43, 44], and unidirectional transport
[45, 46, 47].

However, the effect of long-range coupling and its asymmetry on non-Hermitian
systems have not been fully analyzed. Incorporating unidirectional long-range
coupling into these platforms could open up a new realm of opportunities to
investigate various intriguing aspects of non-Hermitian band topology.
Recently, it was found that such long-range coupling asymmetry introduces
peculiar phenomena associated with non-Hermiticity such as Type-II corner
modes [48, 49, 50], enhancement of topological boundary modes [51, 52] and
complex energy braiding [53, 54]. Furthermore, the inclusion of asymmetric
long-range coupling can lead to dramatic changes in the complex energy
eigenspectra and their corresponding topology. The interplay between
unidirectional long-range coupling and non-Hermiticity has also been studied
in various platforms such as photonics [55], metamaterial [56, 57], optics
[58, 59], condensed matter [38, 60], and topolectrical circuit systems [37,
61, 62, 63, 64, 65, 66, 67, 68]. In photonic systems, asymmetric long-range
coupling can be achieved by using waveguides with asymmetric coupling
coefficients [69], while in metamaterials, it is implemented using
asymmetrical split-ring resonators [70, 71]. Likewise, in condensed matter
systems, such coupling can be modelled via asymmetric long-range electronic
hopping between atoms [72, 73], while in topolectrical circuits, it can be
realized through asymmetric circuit components such as operational amplifiers
[14, 63, 74].

In this work, we investigate the creation of complex energy spectra with
multiple and arbitrary number of twisted loops in the presence of long-range
unidirectional hopping in a Hatano-Nelson (HN) [35, 75, 76] chain. The complex
energy spectrum of a finite HN chain has a non-trivial topology [77, 78] that
can be characterized by a topological invariant called the winding number [13,
79, 80]. In the absence of long-range coupling, the energy spectrum of the HN
chain takes the form of a single closed loop and the winding number assumes a
constant value with a definite sign throughout the complex energy plane
enclosed by the loop. This leads to the accumulation of eigenstates at only
one end of the chain depending on the direction of the dominant asymmetric
coupling, which is reflected in the sign of the winding number.

This simple scenario can be dramatically modified by incorporating
unidirectional long-range coupling, which modifies the trajectory of the
complex eigenvalues and results in multiple twisted loops in the complex
energy plane. The twisted topology translates into an unconventional NHSE
configuration in which eigenstates associated with the clockwise and anti-
clockwise complex energy loops are localized at opposite ends of the chain.
The resulting bipolar NHSE localization overturns the conventional expectation
that NHSE only occurs at one end of the NH chain corresponding to the dominant
coupling direction of the asymmetric nearest-neighbor interactions.
Furthermore, novel behavior was observed at the contact points between
oppositely oriented energy loops, where the energy becomes real-valued. Here,
the topological invariant changes sign leading to the emergence of Bloch-like
contact points at which the NHSE vanishes. Our analysis of a NH model with
long-range unidirectional coupling offers fresh insights into the interplay
between non-Hermiticity and long-range coupling, and the resulting topology
described by its complex energy eigenspectra.

##  II Effects of long-range coupling on complex energy distribution

We consider a periodic one-band system with non-reciprocal nearest-neighbor
couplings $t_{1}\neq t_{-1}$, i.e., the Hatano-Nelson (HN) model [26, 28, 32],
whose eigenenergy is given by

$\displaystyle E_{\mathrm{HN}}$ $\displaystyle=t_{1}\exp(ik)+t_{-1}\exp(-ik)$
$\displaystyle=(t_{1}+t_{-1})\cos(k)+i(t_{1}-t_{-1})\sin(k).$ (1)

The locus of $E_{\mathrm{HN}}$ takes the form of an ellipse in the complex
energy plane with the axis lengths of $|t_{1}+t_{-1}|$ and $|t_{1}-t_{-1}|$
(Fig. 1a). The ellipse cuts across the real energy axis twice under periodic
boundary condition (PBC).

The energy eigenspectrum exhibits distinct distributions under open boundary
conditions (OBC) from the PBC eigenspectrum owing to the asymmetrical coupling
between lattice sites. The tight-binding Hamiltonian corresponding to Eq. (1)
for a finite HN chain that extends from $x=1$ to $x=N$ with OBC is given by

$\displaystyle H_{\mathrm{HN;OBC}}=\sum_{x=1}^{N-1}|x\rangle t_{1}\langle
x+1|+|x+1\rangle t_{-1}\langle x|$ (2)

where $|x\rangle$ and $\langle x|$ are ket and bra vectors representing the
basis states at site $x$.

We derive the eigenenergy spectrum of Eq. (2) using the imaginary gauge
approach in Appendix A and the generalized Brillouin zone (GBZ) approach in
Appendix B, and show that in the eigenenergies $E$ lies on the real energy
axis with $|E|<2\sqrt{t_{1}t_{-1}}$ (Fig. 1a). The marked disparity between
the OBC and PBC eigenspectra (line vs. ellipse on the complex energy plane)
heralds the breakdown of the conventional bulk boundary correspondence (BBC)
in a non-Hermitian system. Specifically, in non-Hermitian systems with
coupling asymmetry, the eigenstates under OBC become localized near a single
edge of the system in the non-Hermitian skin effect (NHSE). As discussed in
detail in Appendix A, the eigenstate localization direction depends on the
relative magnitudes of $t_{1}$ and $t_{-1}$ \- the NHSE localization occurs at
the left (right) edge when $\ln|t_{1}/t_{-1}|$ is positive (negative).

We now consider the introduction of a long-range unidirectional coupling
$t_{-n}\exp(-ink)$ along the left direction a distance of $n$ nodes away. The
eigenenergy $E$ for a periodic system now takes the form of

$\displaystyle E=$ $\displaystyle
t_{1}\exp(ik)+t_{-1}\exp(-ik)+t_{-n}\exp(-ink)$ $\displaystyle=$
$\displaystyle(t_{1}+t_{-1})\cos(k)+t_{-n}\cos(nk)$
$\displaystyle+i\left((t_{1}-t_{-1})\sin(k)-t_{-n}\sin(nk)\right).$ (3)

The locus of the eigenenergy now intersects the real energy axis more than the
two times it does in the conventional HN model. The intersection between the
eigenenergy locus and the real energy axis is governed by the following
equation:

$\displaystyle(t_{1}-t_{-1})\sin(k)-t_{-n}\sin(nk)=0$ (4)
$\displaystyle\Rightarrow$
$\displaystyle\sin(nk)=\frac{t_{1}-t_{-1}}{t_{-n}}\sin(k)$ (5)

Figure 1:  Complex eigenenergy spectra induced by long-range unidirectional
coupling (a) PBC spectra (thick red line), and OBC spectra and spatial
probability density distribution of a conventional Hatano-Nelson chain, 80
nodes in length and without unidirectional long-range coupling
($t_{-1}=0.4,t_{-n}=0$).(b), (c) Plots of $t^{\prime}\sin x$ (green line) and
$\sin(nx)$ (blue line) for (b) $n=5$, $t^{\prime}=C_{5,1}=1$ (top) and
$t^{\prime}=C_{5,2}=1.250$ (bottom), and (c) $n=4$, $t^{\prime}=C_{4,1}=1.089$
and $t^{\prime}=C_{4,2}=4$. (d) Plot of the number of times that the PBC
eigenenergy spectrum curve intersects the real energy axis as a function of
$t_{-1}$ and $t_{-n}$ for $n=5$. The boundaries between regions with different
number of intersection points are demarcated by the white dotted lines and the
corresponding values of $t^{\prime}$ indicated. (e) and (f) Exemplary PBC
eigenenergy spectra and OBC eigenenergy spectra and spatial density
distribution for $n=5$ with (e) $t_{-1}=1.2$, $t_{-5}=-1$, for which the
eigenstates are all localized towards the right, and (f) $t_{-1}=0.4$,
$t_{-5}=0.5$, which exhibits bipolar NHSE. The $(t_{-1},t_{-5})$ values for
(e) and (f) are marked by dots on (d). (g)-(i). Corresponding plots for $n=4$
with (h) $t_{-1}=0.5$, $t_{-4}=-2$ and (i) $t_{-1}=0.2$, $t_{-4}=-0.5$. The
darker dots in the spatial probability density distributions represent high
density values. The states that are localized towards the left (smaller node
numbers) are denoted by green dots and those localized towards the right
(larger node numbers) by blue dots. $t_{1}=1$ for all the plots in the figure.

Introducing $t^{\prime}\equiv\frac{t_{1}-t_{-1}}{t_{-n}}$, the number of real
solutions for of $k$ that satisfies Eq. (5) in the range of $-\pi<k<\pi$,
depends on the value of $t^{\prime}$. For illustration, $\sin(k)$ is plotted
together with $t^{\prime}\sin(nk)$ for an odd value of $n=5$ and even value of
$n=4$ (see Fig. 1b and 1c, respectively). It can be seen that the maximum
number of times the $t^{\prime}\sin(nk)$ curve can intersect the $\sin(k)$
curve is $2n$ times regardless of whether $n$ is even or odd. Moreover, the
$\sin(k)$ curve always touches the $t^{\prime}\sin(nk)$ curve at $k=0,\pi$
regardless of the value of $t^{\prime}$ or $n$.

For odd $n$, the $t^{\prime}\sin(nk)$ curve intersects the $\sin(k)$ curve the
maximum number of $2n$ times when $0<t^{\prime}<1$. However, when $t^{\prime}$
exceeds 1 slightly, the maximum value of $|t^{\prime}\sin(k)|$ at $k=\pm\pi/2$
now exceeds the maximum value of $\sin(k)$ at $\sin(\pm\pi/2)=1$. This results
in the number of intersections between the two curves decreasing by four since
the $t^{\prime}\sin(nk)$ curve now no longer intersects the $\sin(k)$ curve at
$k=\pi/2\pm\delta k$ and at $k=-\pi/2\pm\delta k$. We denote this first
critical value of $t^{\prime}$ which results in the lowering of the number of
intersection points as $C_{n,1}=1$. When $t^{\prime}$ is increased further to
beyond $t^{\prime}=n$, the $t^{\prime}\sin(nk)$ curve no longer intersects the
$\sin(k)$ curve near $k=\pm\pi$ and both sides of $k=0$. This results in a
further decrease by four in the number of intersection points. Moreover, when
$t^{\prime}$ is negative, and its magnitude is increased, we approach a second
critical value of $t^{\prime}=C_{n,2}$ beyond which the $t^{\prime}\sin(k)$
curve no longer intersects the $\sin(nk)$ curve near
$k=\pm\pi(\frac{1}{2}\pm\frac{1}{2n})$ (see lower plot of Fig. 1b), and the
number of intersection points between the $t^{\prime}\sin(k)$ and $\sin(nk)$
curves drops by further a step of eight. For the particular case of $n=5$
shown in Fig. 1b, $C_{5,2}$ was numerically determined to be -1.25.

The abovementioned trends for the number of intersection points between the
$\sin(k)$ and $t^{\prime}\sin(nk)$ curves result in the phase diagram shown in
Fig. 1. Notice that the number of intersection points drop in steps of 4 or 8
as explained above, and that the different phases are not symmetrically
distributed about $t_{-n}=0$. In general, for larger odd values of $n$ there
will be further critical values of $C_{n,m}$ at which the number of
intersection points between the $t^{\prime}\sin(k)$ and $\sin(nk)$ curves
drops in steps of 4 or 8 from the maximum value of $2n$ down to the minimum
value of 2.

Conversely, for even $n$, the $t^{\prime}\sin(k)$ curve intersects across the
$\sin(nk)$ curve for $2n$ times for small $|t^{\prime}|$. As $|t^{\prime}|$
increases, it reaches a critical value of $|t^{\prime}|=C_{n,1}$ beyond which
the $t^{\prime}\sin(x)$ curve no longer intersects the $\sin(k)$ curve near
$k=\pm\pi(\frac{1}{2}+\frac{1}{2n})$. For even $n\geq 4$, as $|t^{\prime}|$
increases further, there are further values of $|t^{\prime}|=C_{n,m}$ beyond
which the number of times the $t^{\prime}\sin(k)$ curve intercepts the
$\sin(k)$ curve decreases further because the $|t^{\prime}\sin(k)|>|\sin(k)|$
at increasing values of $|k-(\pm\pi/2)|$. In particular, at $|t^{\prime}|$
values slightly larger $C_{n,n/2}=n$ the $t^{\prime}\sin(k)$ curve no longer
intersects the $\sin(k)$ curve near $k=0$ for all even $n$ (lower plot of Fig.
1c ). In contrast to the case of odd $n$, the critical values of $t^{\prime}$
exist in $\pm|t^{\prime}|$ pairs, which gives rise to the phase diagram shown
in Fig. 1d, which is symmetrical about $t^{\prime}_{-n}=0$.

With these considerations on the number of intersection points and hence real
solutions to the eigenenergy for the even ($n=4$) and odd case ($n=5$), we
turn our attention to the phase diagrams shown in Fig. 1d and 1g. Figs. 1e, f
and 1h, i show the eigenenergy curves for representative examples of systems
with intermediate (Fig. 1e and h) and maximal (Fig. 1f and i) number of
intersections with the real axis, respectively. Note that these eigenenergy
curves are symmetric about the real axis. This symmetry may be readily
explained by the fact that since the coefficients $t_{i}$s are real,
$E(k)=E^{*}(-k)$ from Eq. (3). Thus, for each value of energy $E$
corresponding to crystal momentum $k$ on the eigenenergy curve, its reflection
about the real axis $E^{*}$ corresponding to $-k$ would also be on the
eigenenergy curve. In particular, the time-reversal invariant momenta $k=0$
and $k=\pi$ lie on the real axis.

For odd $n$, the eigenenergy curves are also symmetric about the imaginary
axis, as shown in Figs. 1e and f. This symmetry is to due to the fact that
$\sin(n(x\pm\pi))=-\sin(nx)$, $\cos(n(x\pm\pi))=-\cos(nx)$ for odd $n$, which
implies that $E(k)=-E(k+\pi)$, i.e., for each energy $E$ corresponding to $k$
that exists on the eigenenergy curve, its reflection about the origin $-E$
corresponding to $k+\pi$ is also on the eigenenergy curve. This results in the
inversion symmetry of the eigenenergy curve about $E=0$ on the complex energy
plane. Correspondingly, the real part of the eigenenergy curve spans between
$-(t_{1}+t_{-1}+t_{-n})$ at $k=\pi$ to $(t_{1}+t_{-1}+t_{-n})$ at $k=-\pi$. In
contrast, for even $n$, the eigenenergy curves are not symmetric about the
imaginary axis (Fig. 1i, j). This is because of the fact that
$\sin(n(x\pm\pi))=\sin(nx)$ and $\cos(n(x\pm\pi))=\cos(nx)$ which implies that
the correspondence $E(k)=-E(k+\pi)$ no longer holds for odd $n$. We find
instead that the real values of $E$ now span from $-(t_{1}+t_{-1}-t_{-n})$ at
$k=\pi$ to $t_{1}+t_{-1}+t_{-n}$ at $k=0$.

We introduce OBC to a modified HN chain of length $N$ with lattice sites
located at $x=1$ to $x=N$ by setting the couplings that extend outside the
extent of the chain to 0. The real-space tight-binding Hamiltonian of the
chain is then given by

$H_{\mathrm{HNL;OBC}}=\sum_{x=1}^{N-1}\left(t_{1}|x+1\rangle\langle
x|+t_{-1}|x\rangle\langle
x+1|\right)+\sum_{x=1}^{N-n}\left(t_{-n}|x\rangle\langle x+n|\right).$ (6)

Remarkably, the long-range coupling can cause the OBC eigenspectrum to become
complex, as depicted in Fig. 1f, h, and i. This differs from the OBC
eigenspectrum of a HN chain without long-range coupling, which lies completely
on the real energy axis. We illustrate how a complex eigenspectrum can emerge
in Appendix B. Notably, the NHSE persists as long as the eigenenergy spectra
under OBC and PBC remain dissimilar. The persistence of the NHSE localization
when the OBC and PBC spectra differ from each other can be intuitively
understood through the following argument: As described in detail in Appendix
B, the wavefunction of an OBC eigenstate $\psi(x)$ at an eigenenergy $E$ has
the general form of $\sum_{j}c_{j}\beta_{j}^{x}$ where
$\beta_{j}\equiv\exp(ik_{j})$ and the $(n+1)$ $k_{j}$s, which are generally
complex, are related to $E$ via Eq. (3). The PBC eigenenergy spectrum is
essentially the loci of $E$ values at which at least one of the $\beta_{j}$
values has a modulus of 1, which in turn corresponds to a real value of $k$.
This implies that none of the $\beta$ values of the OBC eigenenergies has a
modulus of 1 when the OBC eigenspectrum differs from the PBC eigenspectrum,
and the wavefunction grows exponentially towards the left or the right
depending on the sign of the dominant $\mathrm{ln}|\beta_{j}|$ component. We
show in Appendix B that the OBC eigenspectrum in a sufficiently large (on the
order of 5 sites for the parameter ranges here) system is in turn given by the
locus of energy values at which the two largest $|\beta_{j}|$s have the same
moduli.

Unlike the conventional HN system in which the eigenstates are localized near
one edge of the system (Figs. 1a), we see in Figs. 1f and i that for some
parameter ranges of $t_{-1}$ and $t_{-n}$, a peculiar bipolar NHSE appears in
which the OBC eigenstates are localized around both edges of the system. (The
states localized nearer the left edge are denoted by yellowish dots and
located near $E=0$, while the states localized nearer the right edge are
denoted by bluish circles and located further away from $E=0$. ) This unique
localization of eigenstates overturns the conventional expectation that the
edge at which the NHSE localization occurs is determined by the dominant
nearest-neighbor asymmetric coupling direction [31]. In the next section, we
will analyze this bipolar NHSE in more detail by considering the specific
example of the $n=5$ long-range non-Hermitian coupling system, whose eigenmode
localization is illustrated in Figs. 1e and f.

Figure 2:  Winding number and NHSE localization : (a) PBC energy spectrum and
winding number distribution of $n=5$ system without bipolar NHSE corresponding
to Fig. 1e with $t_{-1}=1.2$ and $t_{-5}=-1$. A lighter color on the PBC
energy curve corresponds to smaller values of $k$, which range from $-\pi$ to
$\pi$. The arrow on the PBC curve indicates the direction of increasing $k$.
(b) The $\beta$ distribution on the complex $\beta$ plane at the complex
energies of (top) $E=0.9+0.4i$ and (bottom) $E=2+0.4i$, which are indicated by
the labelled circles in (a). (c) A schematic illustration of the winding
number distribution around an internal loop, as exemplified by the area within
the white dotted rectangle in Fig. 2a. (d) PBC energy spectrum and winding
number distribution of $n=5$ system exhibiting bipolar NHSE with $t_{-1}=0.4$
and $t_{-5}=0.5$, and (e) the $\beta$ distribution at $E=0.1$ (top) and
$E=1.5$ (bottom). (f) A schematic illustration of the winding number around an
external loop, as exemplified by the dotted rectangle in Fig. 2d.

###  II.1  Topological invariant in presence of long range coupling

In the preceding section, we explored how the introduction of long-range
coupling gives rise to skin modes at one or both boundaries. In this section,
we introduce a non-Hermitian topological invariant to analyze the conventional
and bipolar NHSE. The NHSE localization direction at an arbitrary reference
energy $E_{\mathrm{ref}}$ can be determined from the non-Bloch winding number
$\eta$, which is as

$\eta=\frac{1}{2\pi}\oint_{|\beta|=1}\mathrm{d}\beta\
\mathrm{Arg}(E(\beta)-E_{\mathrm{ref}}).$ (7)

A positive (negative) finite value of $\eta$ indicates that a semi-infinite
system extending from $-\infty$ to 0 (0 to $\infty$) would host bulk states
that are NHSE-localized at the right (left) edge of the system. A positive
(negative) winding number also indicates that if bulk OBC states exist within
the energy region around $E_{\mathrm{ref}}$ is bounded by the PBC curve, then
these OBC states will be localized at the right (left) boundaries of the
system. We explain why this is so below using the concept of the GBZ. Figs. 2a
and c show the winding numbers at different energies on the complex energy
planes for the $n=5$ systems shown in Figs. 1e and f, respectively. A
comparison between Fig. 2a and Fig. 1e, and between Fig. 2d and Fig. 1f shows
that this correspondence between the sign of the winding numbers and the
localization direction of the OBC eigenstates indeed holds.

Furthermore, the winding number can be visually determined from the PBC
eigenenergy curve on the complex energy plane if the direction of increasing
$k$ on the eigenenergy curve is known. From Eq. (7), the winding number at the
reference energy $E_{\mathrm{ref}}$ is simply the number of times the
eigenenergy curve winds around $E_{\mathrm{ref}}$ as $k$ is increased from
$-\pi$ to $\pi$ with positive (negative) values of the winding number
corresponding to counter-clockwise (clockwise) windings. As an illustration,
let us consider the reference energy $E_{\mathrm{ref}}=0.9+0.4i$, indicated by
the open circle labelled “A”, in Fig. 2a. Its associated winding number is -1
because the vertical line $x=0.9$ passing through the centre of this circle
cuts across the eigenenergy curve twice, with the eigenenergy curve
progressing from left to right with increasing $k$ above the circle and the
curve progressing from right to left below the circle, indicating a clockwise
winding. Similarly, the winding number around $E_{\mathrm{ref}}=2+0.4i$,
indicated by the open circle labelled “B”, is -2 because the vertical line
$x=2$ passing through the circle cuts across the eigenenergy curve twice above
the circle and twice below it.

Besides this visual interpretation of the winding number, the winding number
can also be expressed in an alternative form that gives more insight into the
connection between the winding number and the NHSE localization direction. Eq.
(7) can be recast into

$\eta=\frac{-i}{2\pi}\oint_{|\beta|=1}\mathrm{d}\beta\
\frac{\partial_{\beta}\left(E(\beta)-E_{\mathrm{ref}}\right)}{E(\beta)-E_{\mathrm{ref}}}$
(8)

where $\beta\equiv\exp(ik)$. Applying the argument principle, $\eta$ is then
given by $\eta(E_{\mathrm{ref}})=Z(E_{\mathrm{ref}})-P(E_{\mathrm{ref}})$
where $P(E_{\mathrm{ref}})$ is the number of poles of
$E(\beta)-E_{\mathrm{ref}}$ while $Z(E_{\mathrm{ref}})$ is the corresponding
number of zeros lying within the complex unit circle. From Eq. (3),
$E(\beta)-E_{\mathrm{ref}}=t_{1}\beta+t_{-1}\beta^{-1}+t_{-n}\beta^{-n}-E_{\mathrm{ref}}$.
Therefore $P=n$, while $Z(E_{\mathrm{ref}})$ is the number of $\beta$s that
satisfy $E(\beta)=E_{\mathrm{ref}}$ and fall within the unit circle on the
complex $\beta$ plane. Fig. 2b shows the distribution of $\beta$s satisfying
$E(\beta)=E_{\mathrm{ref}}$ for $E_{\mathrm{ref}}=0.9+0.4i$ (top) and
$E_{\mathrm{ref}}=2+0.4i$ (bottom). From these plots, it is evident that
$Z(E_{\mathrm{ref}})$ are 4 and 3, respectively, which correspond to the
respective winding numbers of $\eta=4-5=-1$ and $\eta=3-5=-2$, respectively,
for the two values of $E_{\mathrm{ref}}$, noting that n = 5. The winding
numbers obtained here by considering the number of poles and zeros are in
agreement with those obtained by visually counting the number of times the
eigenenergy curve winds around $E_{\mathrm{ref}}$ following Eq. 7 and depicted
in Fig. 2a. Fig. 2e shows corresponding examples for the two other exemplary
$E_{\mathrm{ref}}$ values of 0.1 and 1.5 for the $n=4$ case.

From the above analysis, we can correlate the winding number $\eta$ to the
NHSE localization direction under OBC: The fact that $\eta=Z-P=Z-n$ implies
that the number of $\beta$ values that satisfy $E(\beta)=E_{\mathrm{ref}}$ and
lie within the unit circle in the complex $\beta$ plane is given by
$Z=n+\eta$. We show in Appendix B that for a system with $n$th-order long-
range coupling to the left, the condition for $E_{\mathrm{ref}}$ to lie on the
GBZ and thus be an OBC eigenenergy in the thermodynamic limit, is
$|\beta_{n+1}|=|\beta_{n}|$ where
$|\beta_{1}|\leq|\beta_{2}|...\leq|\beta_{n+1}|$. An $\eta$ value equal to $1$
indicates that all $(n+1)$ $\beta$ values lie within the unit circle at energy
values ‘near’ $E_{\mathrm{ref}}$, which implies that $\beta_{n}$ and
$\beta_{n+1}$ would both be within the unit circle if the GBZ condition
$|\beta_{n}|=|\beta_{n+1}|$ holds. Here, we define an arbitrary energy $E$ to
be ‘near’ $E_{\mathrm{ref}}$ if the PBC eigenenergy curve does not fall
between $E$ and $E_{\mathrm{ref}}$ on the complex energy plane. We explain the
reason for this definition in the next paragraph. Thus, if an OBC eigenstate
exists near $E_{\mathrm{ref}}$, then $|\beta_{n}|=|\beta_{n+1}|<1$ for that
eigenstate, and the eigenstate would be localized at the left edge.
Conversely, if $\eta<0$ at $E_{\mathrm{ref}}$, then both $|\beta_{n}|$ and
$|\beta_{n+1}|$ have magnitudes larger than 1, which implies that any OBC
eigenstate that exists near $E_{\mathrm{ref}}$ would be localized at the right
edge. Finally, an $\eta$ value of exactly 0 would indicate that exactly $n$ of
the smaller $\beta$ values near $E_{\mathrm{ref}}$, including $\beta_{n}$, lie
strictly inside the unit circle while $\beta_{n+1}$ lies outside. This implies
that the condition $|\beta_{n}|=|\beta_{n+1}|$ for an OBC eigenstate to exist
cannot be satisfied, so there are no OBC eigenstates near $E_{\mathrm{ref}}$.

One corollary of $\eta=Z-P$, i.e., the difference between the number of poles
and zeros of $E(\beta)-E_{\mathrm{ref}}$, is that changes in the winding
number as $E_{\mathrm{ref}}$ is varied across the complex plane, would always
entail a crossing of the PBC eigenenergy curve. This is because any change in
the winding numbers will involve the transit of at least one of the $\beta$
roots of the $E(\beta)=E_{\mathrm{ref}}$ equation across the unit circle
boundary in the complex $\beta$ plane. Thus, during this transition, at least
one of the $\beta$ values must lie exactly on the complex unit circle at some
point. When $\beta=\beta^{\prime}$ lies exactly on the complex unit circle,
its corresponding $k=-i\ln(\beta^{\prime})$ value is real, which in turn
implies that $E(\beta^{\prime})$ lies on the PBC energy spectrum.
Consequently, the winding number changes by 1 across a simple linear section
of the PBC spectrum, such as depicted by the blue line segment in Fig. 2c. We
can thus conclude that the NHSE must exist within any non-Hermitian system in
which the eigenenergy curve encloses a finite area . This is because the
winding number within the interior regions bounded by the eigenenergy curve
has a finite value, since it differs by 1 from the winding number (which is 0)
in the region outside the curve that extends to infinity. Conversely, there
will be no NHSE when the eigenenergy curve takes the form of an open curve
which does not enclose any finite area in the complex energy plane.

###  II.2 Mechanism of complex energy loop crossing and its correspondence
with NHSE

In the preceding subsection, we have examined the impact of the interplay
between complex energy and long-range coupling on the sign and magnitude of
the winding number, which in turn affects the localization of the NHSE.
Intriguingly, the presence of long-range coupling can alter the complex energy
spectra and lead to intersecting eigenenergy curves and various loop
configurations. In this subsection, we discuss how long-range coupling also
induces and modifies the complex eigenenergy loop crossings and discuss their
correspondence with the winding number and NHSE.

As an illustration, Fig. 2c shows a schematic representation of an internal
loop demarcated by the white dotted rectangle in Fig. 2a in which the
eigenenergy curve intersects itself in a crossing that resembles an “X” shape
to form an inner loop enclosed within a larger outer loop in the complex
energy plane area bounded by the curve. We shall refer to such self-crossings
as X crossings for conciseness henceforth. The intersection point in an X
crossing corresponds to having two values of $\beta$ lying on the unit circle
on the complex plane. Consequently, when one traverses across the two opposite
sides of the X crossing, the winding number must change by either $\pm 2$
corresponding to having both $\beta$ values moving into or out of the unit
circle, or remain unchanged which corresponds to one $\beta$ value moving into
the unit circle while the other $\beta$ moves out of it.

By considering the direction of increasing $k$ on the eigenenergy curve
denoted by the arrowheads in Fig. 2c, it can be seen that the latter (i.e. $0$
change in the winding number) corresponds to the two quadrants of the
X-crossing bounded by arrowheads pointing both towards or away from the
crossing point (i.e., the upper and lower quadrants of the X-crossing), and
the former (i.e. $\pm 2$ change in the winding number) involves the quadrants
bounded by one arrowhead pointing towards and the other pointing away from the
crossing point (the upper and lower quadrants of the X-crossing). Additonally,
as shown in Fig. 2c, the winding number within the internal loop has the same
sign as that the larger outer loop bounded by the eigenenergy curve.

We now consider the self-intersection of the PBC eigenenergy curve to form an
external loop that protrudes into the region exterior to the eigenenergy
curve. (Note that “intersection” here refers to the eigenenergy curve
intersecting itself, while the related term “crossing” refers to the
eigenenergy curve crossing the real axis.). One example is demarcated by the
dotted white rectangle in Fig. 2d, whose schematic representation is shown in
Fig. 2f. Because an external loop is directly adjacent to the region exterior
to the PBC curve, its winding number is necessarily $\pm 1$ (since the winding
number associated with the exterior region is $0$). Furthermore, by
considering the direction of increasing $k$ on the eigenenergy curve, it can
be seen that the winding number changes by $\pm 2$ as we traverse between the
two quadrants of the X-crossing that link the external loop to the interior of
the region bounded by the eigenenergy curve. This results in the former having
a winding number of opposite sign to that in the latter. The different signs
of the winding number translates to the phenomenon of bipolar NHSE in which
the OBC eigenstates are localized along both edges of the finite system.

Figure 3:  Complex energy loops with alternating signs of winding numbers (a)
Schematic representation of system with $n$ complex energy loops with
alternating signs of winding numbers and $(n-1)$ X crossings denoted as BP1 to
BPn-1. (b) Number of self-intersections in complex eigenenergy spectrum at
$n=5$ as function of $t_{-1}$ and $t_{-5}$. The area demarcated by the solid
red lines represents the region in which the eigenenergy spectrum takes the
form of alternating loops shown in Fig. 3a. The dotted white line demarcates
the region on the $t_{-1}$–$t_{-5}$ plane in which the eigenenergy
distribution crosses the real energy axis 10 times. (c) - (e) (Top) PBC
eigenenergy spectrum, and OBC eigenenergy spectrum and spatial density
distribution, and (bottom) $\beta$-plane GBZ, of (c) $n=5$, (d) $n=6$, and (e)
$n=7$ systems with $t_{1}=1$, $t_{-1}=1.05$, and $t_{-n}=0.05$.

Having established the correlation between the external eigenenergy loops and
the emergence of the bipolar NHSE, we will proceed to show that the number of
switches between the NHSE localization at the two ends of the chain can be
engineered by tuning the system parameters so to increase the number of
X-crossings. We have earlier explained how the eigenenergy of an $n$-th order
long-range coupled system described by Eq. 3 can cross the imaginary axis up
to a maximum of $2n$ times. Hence, a large number of alternations between NHSE
localizations can be achieved by tuning the system parameters such that all
the $2(n-1)$ real-axis eigenenergy crossings (disregarding the two crossings
at the extreme ends of the eigenenergy curve) take the form of X crossings. In
other words, we have $(n-1)$ X-crossings (n.b. each X-crossing corresponds to
$2$ crossings of the real axis by the two degenerate states at each X-point)
resulting in $n$ loops, as shown schematically in Fig. 3a. Such a
configuration corresponds to the maximum number of times the complex energy
loop crosses the real energy axis and the minimal number that the loop
intersects itself where all the self-intersections occur at the crossings.
Fig. 3b plots the number of self-intersections for a $n=5$ long-range coupled
system as a function of the coupling parameters $t_{-n}$ and $t_{-1}$. The
intersection between the region with the maximal number of real energy axis
crossings demarcated by the dotted white lines (refer also to Fig. 1d) and the
region with the minimal number of self-intersections is denoted by the solid
white lines marks. Within this intersection region, we have the desired
eigenenergy spectrum which takes the form of consecutive loops along the real
axis (schematically drawn in Fig. 3a). From Fig. 3b, we see that such a
configuration exists around the vicinity of $t_{1}=1$, $t_{-n}=0$ (although
not exactly at $t_{-n}=0$, which would correspond to a conventional HN
system.)

###  II.3 Formation of series of twisted complex energy loops with opposite
windings and Bloch points

To illustrate the concept of series of twisted complex energy loops with
opposite winding index (shown schematically in Fig. 3a), we consider various
order of unidirectional long range coupling. With the optimal choice of the
coupling parameters (e.g. corresponding to the triangular areas bounded by the
white dotted lines in Fig. 3b for the case $n=5$), we depict the respective
examples for $n=5$ to $7$ for which the eigenenergy spectrum takes the form of
$n$ consecutive loops. As we traverse each X-crossing, the winding direction
(clockwise / anticlockwise) of the eigenenergy curve flips. This leads to a
change in the sign of the winding number and the flipping of the NHSE
localization direction of the OBC states between successive loops. The
existence of bipolar NHSE localization implies that the GBZ would contain
segments that lie both within and outside the unit circle on the complex
$\beta$ plane, as verified in the bottom plots of Figs. 3b to 3e.

Let us consider the Bloch-like points that lie exactly on the complex unit
circle in Figs. 3b to 3e, which correspond to the $k$ values at which the
eigenenergy curve intersects the real axis and forms the X-crossings.
Interestingly, these OBC eigenstates that lie exactly at each X crossing are
not localized at either of the edges. This absence of localization can be
explained in terms of the following: Each internal intersection point at a X
crossing on the complex energy plane represents a transition between a region
with $\eta=1$ which possesses $(n+1)$ $\beta$ values lying within the unit
circle, and a region with $\eta=-1$ which possesses $(n-1)$ $\beta$ values
within the unit circle (and hence two $\beta$ values outside it). Thus, each
$X$ crossing corresponds to an energy value where two $\beta$ states lie
exactly on the complex unit circle as one state moves from outside the unit
circle to inside, and the other state moves in the opposite direction. At this
intersection point, there are $(n-2)$ $\beta$ values inside the complex unit
circle and one $\beta$ value outside. This means that the $n$ and $(n+1)$th
largest $|\beta|$ values are both on the complex unit circle. In other words,
the $\beta$-plane GBZ on the X-crossing falls on the complex unit circle.
Thus, the NHSE localization length is 0, which results in the absence of NHSE
localization.

##  III Conclusion

In conclusion, we have analyzed the generalized Hatano-Nelson chain which
incorporates the effects of asymmetric long-range coupling of arbitrary
orders. The introduction of unidirectional long-range coupling results in
multiple twisting topology of its PBC eigenspectra in complex energy space.
The number of eigenenergy loops equates to the order of coupling. Furthermore,
we showed the close correlation between the loop topology and the NHSE
localization of the eigenstates of the system. The loop topology can be
engineered by tuning the long-range coupling parameters to exhibit maximum
number of external loops of alternating winding numbers. This leads to the
phenomenon of bipolar NHSE, where the position of the NHSE localization is
determined by the sign of the winding index rather than the direction of
coupling asymmetry according to the conventional understanding. Interestingly,
the crossing points of the eigenenergy loops are characterized by real energy
Bloch-like points, at which the NHSE localization vanishes. We provide the
underlying theoretical basis for the observed NHSE phenomena induced by the
eigenenergy loop topology, as well as numerical verification of the theory.
Finally, our results provide a flexible and accessible method to modulate the
complex energy spectra of a long-range coupled non-Hermitian system, thereby
realizing multiple twisted eigenenergy loop topology and inducing the novel
NHSE localization as described above.

#### Acknowledgements

This work is supported by the Ministry of Education (MOE) Tier-II Grant
MOE-T2EP50121-0014 (NUS Grant No. A-8000086-01-00), and MOE Tier-I FRC Grant
(NUS Grant No. A-8000195-01-00).

## Appendix A: Localization of eigenstates in Hatano-Nelson model: imaginary
gauge transformation approach

The imaginary gauge transformation (similarity transformation) is a
historically important approach in understanding the exponential localization
and eigenenergy spectrum of the HN chain without long-range coupling. In this
section, we apply this approach on the finite chain represented by the
Hamiltonian Eq. (2). An imaginary gauge transformation

$\displaystyle|x\rangle$ $\displaystyle\rightarrow|x\rangle\rangle e^{-gx},$
(9) $\displaystyle\langle x|$ $\displaystyle\rightarrow e^{gx}\langle\langle
x|,$ (10) $\displaystyle g$ $\displaystyle=\ln\sqrt{\frac{t_{1}}{t_{-1}}}$
(11)

can be applied to Eq. (2) to convert it into the equivalent Hermitian
Hamiltonian

$H_{\mathrm{OBC}}=\sqrt{t_{1}t_{-1}}\sum_{x=1}^{N-1}(|x\rangle\rangle\langle\langle
x+1|+|x\rangle\rangle\langle\langle x+1|),$ (12)

where the double left / right angular brackets denote the transformed basis
states.

Importantly, this transformation preserves the energy spectrum because it is a
similarity transformation. Consequently, the eigenvalue spectrum of Eq. (2) is
identical to that of Eq. (12). This implies that the eigenspectrum remains
real for any degree of non-Hermiticity as long as $\sqrt{t_{1}t_{-1}}$ is
real.

At the same time, the introduction of the imaginary gauge potential induces a
position-dependent scaling of the eigenfunctions. From Eq. (10), we read off
that the spatial wavefunction of the eigenstate $|\psi\rangle$ in the basis of
Eq. (2), $\langle x|\psi(x)\rangle$, is related to that in the basis of Eq.
(12), $\langle\langle x|\psi(x)\rangle$ by

$\langle x|\psi\rangle=e^{gx}\langle\langle x|\psi\rangle.$ (13)

The $e^{gx}$ term on the right of Eq. (13) results in an exponential growth or
decrease of the right wavefunction depending on the sign of $g$. All the
eigenstates therefore become localized at one of the chain edges under OBC
when the coupling is asymmetrical ($t_{1}\neq t_{-1}$). Eq. (12) is
essentially the lattice version of Hamiltonian one-dimensional free electron
gas in an infinite potential well, which has the well-known solution
$\langle\langle x|\psi\rangle=e^{ikx}-e^{-ikx}$ where $k=(2n\pi/(N+1))$,
$n=1,...,N$ for a chain that extends from $x=1$ to $x=N$. From the Bloch
counterpart of Eq. (12), $H(k)=2\sqrt{t_{1}t_{-1}}\cos(k)$, it can be deduced
that the OBC eigenenergies of Eq. (13) lie on the real line $|E|\leq
2\sqrt{t_{1}t_{-1}}$.

The relationship between $\langle x|$ and $\langle\langle x|$ then implies
that the wavefunction of an OBC eigenstate $\psi_{HN}(x)\equiv\langle
n|\psi_{\mathrm{HN}}$ in the basis of Eq. (2) is explicitly given by

$\psi_{\mathrm{HN}}(x)=\left(\frac{t_{-1}}{t_{1}}\right)^{x/2}(e^{ikx}-e^{-ikx}).$
(14)

A notable observation from Eq. (14) is that the magnitude of the wavefunction
is directly proportional to $\left(\frac{t_{-1}}{t_{1}}\right)^{N/2}$. This
dependence has a significant implication: when $|t_{-1}|>|t_{1}|$, the
eigenstates become exponentially localized near the right edge of the chain at
$x=N$. Conversely, if $|t_{-1}|<|t_{1}|$, the eigenstates exhibit exponential
localization near the left edge of the chain at $x=1$.

Although the imaginary gauge transformation described above sheds light on the
effects of asymmetrical coupling and the spatial distribution of the
wavefunction in simple non-Hermitian systems like the non-reciprocal Hatano-
Nelson system above or the Su-Schrieffer-Heeger chain, it is not universally
applicable to more complicated systems such as the ones with asymmetrical
long-range coupling studied here. This is because there is no similarity
transformation with a constant value of $g$ like that in Eq. (9) and (10) that
can be performed on a generic non-reciprocal long-range Hamiltonian to convert
it into a Hermitian Hamiltonian with purely real eigenvalues. The non-
existence of such similarity transformations is hinted at by the fact that
unlike the HN system where $|\beta|$ has the constant value of $g$ in Eq. (11)
throughout the entire GBZ, which then takes the form of a circle in the
complex $\beta$ plane, the $|\beta|$ values in general vary at different
points on the GBZ, as can be seen from the lower plots in Fig. 3c–e. Another
indication that such similarity transformations do not exist is the fact that
the OBC energy spectra for these systems are complex rather than real: if a
similarity transformation that converts the non-Hermitian Hamiltonian to a
Hermitian one and preserves the eigenvalues exists, it would not have been
possible to obtain complex eigenvalues from a Hermitian Hamiltonian. A more
modern and universally applicable approach that has been commonly adopted to
explain the non-Hermitian skin effect in more recent works over the past three
years is the GBZ, which we explain in more detail in the next section.

## Appendix B: Generalized Brillouin zone

As noted above, although the OBC eigenspectra of prototypical systems as the
Su-Schrieffer-Heeger (SSH) and Hatano-Nelson chains with nearest-neighbor
couplings are consistently real, this is not always the case for more complex
models involving long-range couplings or gain/loss terms [31, 73, 32, 81, 82,
35].

The emergence of complex eigenenergies is not surprising because the
eigenvalues of a non-Hermitian matrix, such as the Hamiltonian of a finite-
length chain with asymmetrical coupling in Eq. (6), are not restricted to real
values but can, in general, be complex. Fundamentally, the eigenvalue $E$ must
satisfy the requirement that the Schrödinger equation $\langle
x|H|\psi\rangle=\langle x|\psi\rangle E$ is satisfied by the eigenstate of a
Hamiltonian of a finite-length chain (i.e., under OBC) at all the lattice
sites lying within the extent of the chain, i.e., $x=1,...,N$.

Consider the Hamiltonian Eq. (6). For a lattice site $x$ that lies within the
interior in the chain for which all the sites it is coupled to by $t_{-1}$,
$t_{1}$, and $t_{-n}$ lie within the chain, i.e., $n<x<N-1$, the Schrödinger
equation at $x$ reads

$t_{1}\psi(x+1)+t_{-1}\psi(x-1)+t_{-n}\psi(x-n)=E\psi(x).$ (15)

This is the same equation that is obeyed at any lattice site inside an
infinitely long chain. In non-Hermitian systems, the Bloch theorem for
Hermitian system, which states that the wavefunction of a periodic system with
a unit cell containing a single lattice point has the form of $\exp(ikx)$, is
extended so that $k$ is no longer limited to real values but can, in general,
be complex (see, for example, [83]). It is conventional to introduce
$\beta\equiv\exp(ik)$. Writing $\psi(x)=\beta^{x}$ in Eq. (15) gives

$t_{1}\beta^{n+1}+t_{-1}\beta^{n-1}+t_{-n}-E\beta^{n}=0,$ (16)

which is an $(n+1)$th-order polynomial in $\beta$. For a given $E$, Eq. (16)
has $n+1$ solutions for $\beta$, which we label as
$\beta_{1},\beta_{2},...,\beta_{n+1}$ where
$|\beta_{1}|\leq|\beta_{2}|\leq...\leq|\beta_{n+1}|$. Eq. (15) is then
satisfied by any linear combination of the $n+1$ $\beta$ values

$\psi(x)=\sum^{n+1}_{j=1}\beta_{j}^{x}c_{j}$ (17)

where the $n+1$ $c_{j}$s are position-independent constant coefficients. In
particular, $\psi(x)$ in Eq. (17) is also an eigenstate of the OBC Hamiltonian
Eq. (6) when appropriate boundary conditions are applied as follows: We note
that Eq. (15) will also hold for an eigenstate of Eq. (6) at $1\leq x\leq n$
and $x=N$ if we introduce the $n+1$ constraints that
$\psi(-n+1)=...=\psi(0)=0$ and $\psi(N+1)=0$ [for example, $\langle
x=1|H\rangle=E\psi(1)=t_{1}\psi(2)$ in Eq. (6) is equal to
$t_{1}\psi(2)+t_{1}\psi(0)+t_{-n}\psi(1-n)$ if $\psi(0)=\psi(1-n)=0$].
Substituting the expression for $\psi(x)$ in Eq. (17) into these $n+1$
constraints results in a homogenous system of $n+1$ linear equations in the
$n+1$ unknown $c_{j}$s. The eigenenergies of a chain with any finite value of
$N$ can then be solved for exactly by finding the values of $E$ at which the
determinant of this system of linear equations is zero.

The GBZ approach provides a simpler approach for obtaining the loci of the
eigenenergies on the complex energy plane in the thermodynamic limit
$N\rightarrow\infty$ compared to computing the zeros of the determinant
explicitly. The key idea in this is that a certain pair of the $\beta$ values
is required to have the same moduli so that the boundary conditions can be
satisfied at both ends of the chain simultaneously, as explained in the
following. To facilitate the explanation, we shift the $x$ position labels of
the chain from $x=1$ – $x=N$ to $x=-(N-1)/2$ – $x=(N-1)/2$. The boundary
conditions then become
$\psi(-(N-1)/2-n)=\psi((N-1)/2-n+1)=...=\psi(-(N-1)/2-1)=0$ and
$\psi((N-1)/2+1)=0$.

Consider first the boundary condition at $x=(N-1)/2+1$. As
$N\rightarrow\infty$, the absolute value of $\beta_{j}^{(N-1)/2+1}$ for the
smaller $\beta_{j}$ values with $j=1,..,n-1$ become negligibly small compared
to those of $\beta_{n}$ and $\beta_{n+1}$. The corresponding
$c_{j}\beta_{j}^{(N-1)/2+1},j=1,...,n-1$ terms in Eq. (17) can then be
approximated to zero and we have

$\psi((N-1)/2+1)\approx
c_{n}\beta_{n}^{(N-1)/2+1}+c_{n+1}\beta_{n+1}^{(N-1)/2+1}.$ (18)

Note that we cannot approximate the $c_{n}\beta_{n}^{(N-1)/2+1}$ term in Eq.
(18) to 0 because otherwise, the $c_{n+1}\beta_{n+1}^{(N-1)/2+1}$ cannot be
cancelled off to make $\psi((N-1)/2+1)$ zero.

Consider next the $n$ boundary conditions at the left end of the chain
$\psi(x)=0,x=-(N-2)/2-n,...,-(N-1)/2-1$. To guarantee the existence of a
solution for these $n$ equations, we need all $n+1$ terms in Eq. (17) to be of
approximately the same order of magnitude at these values of $x$ so that there
are more non-negligible free variables (i.e., the $n+1$ $c_{j}$s) than
constraints (the $n$ boundary conditions). Now considering all the $n+1$
boundary conditions at the left and right ends collectively, we note that
there is a requirement for

$|c_{n}\beta_{n}^{x}|\simeq|c_{n+1}\beta_{n+1}^{x}|$ (19)

at both negative values of $x$ at $x=-(N-2)/2-n,...,-(N-1)/2-1$ and at a
positive value of $x$ at $x=(N-1)/2+1$. Eq. (19) can hold at both the negative
and positive values of $x$ as $N\rightarrow\infty$ only when
$|\beta_{n}|=|\beta_{n+1}|$; otherwise, if $|\beta_{n+1}|$ is slightly larger
than $|\beta_{n}|$, the left side of Eq. (19) will become exponentially
smaller than the right side at $x=(N-1)/2+1$, and exponentially larger than
the right at $x=-(N-2)/2-n,...,-(N-1)/2-1$ as $N\rightarrow\infty$ for any
finite values of $c_{n}$ and $c_{n+1}$. The loci of the OBC energy eigenvalues
at large values of $N$ therefore approaches the loci of $E$ at which
$|\beta_{n}|=|\beta_{n+1}|$, which gives the GBZ. (This condition differs from
the usual criteria that it is the moduli of the middle pair of $|\beta|$ that
needs to have the same value rather than the largest pair of $|\beta|$ values
here because the former applies only for systems at which the furthest
coupling to the left and right have the same distances, whereas the long-range
coupling here is unidirectional. )

We illustrate the application of the GBZ approach through the example of the
system with second-order unidirectional coupling in Eq. (6). The solutions for
$\beta$ of Eq. (16) at $n=2$ are then given by

$\displaystyle\beta_{1}=\frac{E}{3t_{1}}-\frac{2^{\frac{1}{2}}\lambda_{2}}{3t_{1}\lambda}+\frac{\lambda}{3\times
2^{\frac{1}{3}}t_{1}},$ (20)
$\displaystyle\beta_{2}=\frac{E}{3t_{1}}+\frac{(1+\sqrt{3}i)\lambda_{2}}{3\times
2^{\frac{2}{3}}t_{1}\lambda}-\frac{(1-\sqrt{3}i)\lambda}{6\times
2^{\frac{1}{3}}t_{1}},$ (21)
$\displaystyle\beta_{3}=\frac{E}{3t_{1}}+\frac{(1-\sqrt{3}i)\lambda_{2}}{3\times
2^{\frac{2}{3}}t_{1}\lambda}-\frac{(1+\sqrt{3}i)\lambda}{6\times
2^{\frac{1}{3}}t_{1}},$ (22)

where $\lambda_{1}=2E^{3}-27t_{1}^{2}t_{-2}^{2}-9Et_{1}t_{-1}$,
$\lambda_{2}=3t_{1}t_{-1}-E^{2}$ and
$\lambda=\left(\lambda_{1}+\sqrt{4\lambda_{2}^{3}+\lambda_{1}^{2}}\right)^{1/3}$.
Following the arguments above, the OBC eigenspectrum is given by the loci of
$E$ where $|\beta_{2}|=|\beta_{3}|$. Although the loci of $E$ that satisfies
this requirement is obviously too complicated to solve for analytically, it
can be appreciated from the presence of the complex coefficients in Eq. (21)
and (22) that the solutions for $E$ are, in general, complex and not purely
real. Moreover, the common value of $|\beta_{2}|$ and $|\beta_{3}|$ is, in
general, not necessarily 1. This results in an exponential localization of the
wavefunction, i.e., the NHSE via Eq. (17).

For comparison, we also derive the OBC eigenenergy spectrum of the HN chain
without long-range coupling (i.e., $t_{n}=0$) using the GBZ approach. In this
case, the two values of $\beta$ are given by

$\beta_{\pm}=\frac{E\pm\sqrt{E^{2}-4t_{1}t_{-1}}}{2t_{1}},$ (23)

and it is required that $|\beta_{+}|=|\beta_{-}|$ on the GBZ. A key difference
between Eqs. (23), for which there is no long-range coupling, and (21) and
(22), for which there is a second-order long-range coupling, is that there are
no complex coefficients in the former. This opens the possiblity for the
solutions of $E$ in $|\beta_{+}|=|\beta_{-}|$ to be purely real in the HN
chain rather than complex. Indeed, $|\beta_{+}|=|\beta_{-}|$ holds when the
$\beta_{\pm}$s form a complex conjugate pair. This occurs when the term in the
square root, viz. $E^{2}-4t_{1}t_{-1}$, is negative. The OBC spectrum of the
HN chain in the thermodynamic limit is hence given by
$|E|<2\sqrt{t_{1}t_{-1}}$, which matches the OBC spectrum obtained using the
imaginary gauge approach.

## References

  * Senthil [2015] T. Senthil, Annu. Rev. Condens. Matter Phys. 6, 299 (2015).
  * Rafi-Ul-Islam  _et al._ [2022a] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, and M. Jalil, Front. Phys. 10, 1202 (2022a).
  * Fujita  _et al._ [2011] T. Fujita, M. Jalil, S. Tan, and S. Murakami, J. Appl. Phys. 110 (2011).
  * Tan  _et al._ [2020] S. G. Tan, M. B. Jalil, C.-R. Chang, and S. Murakami, in  _Spin_ , Vol. 10 (World Scientific, 2020) p. 2050018.
  * Tan and Jalil [2012] S. G. Tan and M. B. Jalil,  _Introduction to the Physics of Nanoelectronics_ (Elsevier, 2012).
  * Ma  _et al._ [2019] G. Ma, M. Xiao, and C. T. Chan, Nat. Rev. Phys. 1, 281 (2019).
  * Rafi-Ul-Islam  _et al._ [2023a] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, and M. B. Jalil, J. Magn. Magn. Mater. , 170559 (2023a).
  * Hafezi  _et al._ [2013] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. Taylor, Nat. Photon. 7, 1001 (2013).
  * Rafi-Ul-Islam  _et al._ [2020a] S. Rafi-Ul-Islam, Z. Bin Siu, and M. B. Jalil, Commun. Phys. 3, 72 (2020a).
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020).
  * Rafi-Ul-Islam  _et al._ [2021a] S. Rafi-Ul-Islam, Z. B. Siu, and M. B. Jalil, Phys. Rev. B 103, 035420 (2021a).
  * Yao and Wang [2018] S. Yao and Z. Wang, Phys. Rev. Lett. 121, 086803 (2018).
  * Leykam  _et al._ [2017] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Phys. Rev. Lett. 118, 040401 (2017).
  * Rafi-Ul-Islam  _et al._ [2021b] S. Rafi-Ul-Islam, Z. B. Siu, and M. B. Jalil, New J. Phys. 23, 033014 (2021b).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Wang  _et al._ [2019] H. Wang, J. Ruan, and H. Zhang, Phys. Rev. B 99, 075130 (2019).
  * Martinez Alvarez  _et al._ [2018] V. Martinez Alvarez, J. Barrios Vargas, M. Berdakin, and L. Foa Torres, Eur. Phys. J.: Spec. Top. 227, 1295 (2018).
  * Okuma and Sato [2022] N. Okuma and M. Sato, Annu. Rev. Condens. Matter Phys. 14 (2022).
  * Ghorashi  _et al._ [2021] S. A. A. Ghorashi, T. Li, and M. Sato, Phys. Rev. B 104, L161117 (2021).
  * Bao  _et al._ [2021] L. Bao, B. Qi, D. Dong, and F. Nori, Phys. Rev. A 103, 042418 (2021).
  * McDonald and Clerk [2020] A. McDonald and A. A. Clerk, Nat. Commun. 11, 5382 (2020).
  * Chen  _et al._ [2021] W. Chen, M. Abbasi, Y. N. Joglekar, and K. W. Murch, Phys. Rev. Lett. 127, 140504 (2021).
  * Budich and Bergholtz [2020] J. C. Budich and E. J. Bergholtz, Phys. Rev. Lett. 125, 180403 (2020).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Phys. Rev. Lett. 77, 570 (1996).
  * Hatano and Nelson [1997] N. Hatano and D. R. Nelson, Phys. Rev. B 56, 8651 (1997).
  * Hatano and Nelson [1998] N. Hatano and D. R. Nelson, Phys. Rev. B 58, 8384 (1998).
  * Rudner and Levitov [2009] M. S. Rudner and L. Levitov, Phys. Rev. Lett. 102, 065703 (2009).
  * Jazaeri and Satija [2001] A. Jazaeri and I. I. Satija, Phys. Rev. E 63, 036222 (2001).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Phys. Rev. Lett. 124, 086801 (2020).
  * Lin  _et al._ [2023] R. Lin, T. Tai, L. Li, and C. H. Lee, Front. Phys. 18, 53605 (2023).
  * Li  _et al._ [2022a] Y. Li, C. Liang, C. Wang, C. Lu, and Y.-C. Liu, Phys. Rev. Lett. 128, 223903 (2022a).
  * Rafi-Ul-Islam  _et al._ [2022b] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, C. H. Lee, and M. B. Jalil, Phys. Rev. Res. 4, 013243 (2022b).
  * Li  _et al._ [2020] L. Li, C. H. Lee, S. Mu, and J. Gong, Nat. Commun. 11, 5491 (2020).
  * Song  _et al._ [2019] F. Song, S. Yao, and Z. Wang, Phys. Rev. Lett. 123, 170401 (2019).
  * Rafi-Ul-Islam  _et al._ [2022c] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, C. H. Lee, and M. B. Jalil, Phys. Rev. Res. 4, 043108 (2022c).
  * Rafi-Ul-Islam  _et al._ [2022d] S. Rafi-Ul-Islam, H. Sahin, Z. B. Siu, and M. B. Jalil, Phys. Rev. Res. 4, 043021 (2022d).
  * Kawabata  _et al._ [2019] K. Kawabata, T. Bessho, and M. Sato, Phys. Rev. Lett. 123, 066405 (2019).
  * Minganti  _et al._ [2019] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Phys. Rev. A 100, 062131 (2019).
  * Hokmabadi  _et al._ [2019] M. P. Hokmabadi, A. Schumer, D. N. Christodoulides, and M. Khajavikhan, Nature 576, 70 (2019).
  * De Carlo  _et al._ [2022] M. De Carlo, F. De Leonardis, R. A. Soref, L. Colatorti, and V. M. Passaro, Sensors 22, 3977 (2022).
  * Bao  _et al._ [2022] L. Bao, B. Qi, and D. Dong, Phys. Rev. Appl. 17, 014034 (2022).
  * Koch and Budich [2022] F. Koch and J. C. Budich, Phys. Rev. Res. 4, 013113 (2022).
  * Longhi  _et al._ [2015a] S. Longhi, D. Gatti, and G. Della Valle, Phys. Rev. B 92, 094204 (2015a).
  * Longhi  _et al._ [2015b] S. Longhi, D. Gatti, and G. D. Valle, Sci.Rep. 5, 1 (2015b).
  * Du  _et al._ [2020] L. Du, Y. Zhang, and J.-H. Wu, Sci. Rep. 10, 1113 (2020).
  * Rafi-Ul-Islam  _et al._ [2022e] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, and M. B. Jalil, Phys. Rev. B 106, 245128 (2022e).
  * Chen  _et al._ [2022] C. Chen, T. Chen, W. Ding, and J. Zhu, ArXiv preprint arXiv:2203.13149 (2022).
  * Yang  _et al._ [2022a] H. Yang, L. Song, Y. Cao, and P. Yan, Phys. Rev. B 106, 075427 (2022a).
  * Viyuela  _et al._ [2018] O. Viyuela, L. Fu, and M. A. Martin-Delgado, Phys. Rev. Lett. 120, 017001 (2018).
  * Wang  _et al._ [2022] Y.-C. Wang, H. Jen, and J.-S. You, ArXiv preprint arXiv:2211.16565 (2022).
  * Wang  _et al._ [2021] K. Wang, A. Dutt, C. C. Wojcik, and S. Fan, Nature 598, 59 (2021).
  * Li  _et al._ [2022b] Y. Li, X. Ji, Y. Chen, X. Yan, and X. Yang, Phys. Rev. B 106, 195425 (2022b).
  * Zhu  _et al._ [2020] X. Zhu, H. Wang, S. K. Gupta, H. Zhang, B. Xie, M. Lu, and Y. Chen, Phys. Rev. Res. 2, 013280 (2020).
  * Thevamaran  _et al._ [2019] R. Thevamaran, R. M. Branscomb, E. Makri, P. Anzel, D. Christodoulides, T. Kottos, and E. L. Thomas, J. Acoust. Soc. Am. 146, 863 (2019).
  * Wen  _et al._ [2020] X. Wen, X. Zhu, A. Fan, W. Y. Tam, J. Zhu, F. Lemoult, M. Fink, and J. Li, ArXiv preprint arXiv:2011.01006 (2020).
  * Kominis  _et al._ [2017] Y. Kominis, T. Bountis, and S. Flach, Phys. Rev. A 95, 063832 (2017).
  * Ivars  _et al._ [2022] S. B. Ivars, M. Botey, R. Herrero, and K. Staliunas, Phys. Rev. A 105, 033510 (2022).
  * Jin and Song [2019] L. Jin and Z. Song, Phys. Rev. B 99, 081103 (2019).
  * Sahin  _et al._ [2023] H. Sahin, Z. B. Siu, S. Rafi-Ul-Islam, J. F. Kong, M. B. Jalil, and C. H. Lee, Phys. Rev. B 107, 245114 (2023).
  * Rafi-Ul-Islam  _et al._ [2023b] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, and M. B. Jalil, Phys. Rev. Res. 5, 013107 (2023b).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Nat. Phys. 16, 747 (2020).
  * Rafi-Ul-Islam  _et al._ [2022f] S. Rafi-Ul-Islam, Z. B. Siu, H. Sahin, C. H. Lee, and M. B. Jalil, Phys. Rev. B 106, 075158 (2022f).
  * Rafi-Ul-Islam  _et al._ [2020b] S. Rafi-Ul-Islam, Z. B. Siu, C. Sun, and M. B. Jalil, New J. Phys. 22, 023025 (2020b).
  * Hofmann  _et al._ [2020] T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši, _et al._ , Phys. Rev. Res. 2, 023265 (2020).
  * Rafi-Ul-Islam  _et al._ [2020c] S. Rafi-Ul-Islam, Z. B. Siu, and M. B. Jalil, Appl. Phys. Lett. 116, 111904 (2020c).
  * Zhang  _et al._ [2023] X. Zhang, B. Zhang, H. Sahin, Z. B. Siu, S. Rafi-Ul-Islam, J. F. Kong, B. Shen, M. B. Jalil, R. Thomale, and C. H. Lee, Commun. Phys. 6, 151 (2023).
  * Bell  _et al._ [2017] B. A. Bell, K. Wang, A. S. Solntsev, D. N. Neshev, A. A. Sukhorukov, and B. J. Eggleton, ArXiv preprint arXiv:1709.01267 (2017).
  * Decker  _et al._ [2011] M. Decker, N. Feth, C. M. Soukoulis, S. Linden, and M. Wegener, Phys. Rev. B 84, 085416 (2011).
  * Roy Chowdhury  _et al._ [2012] D. Roy Chowdhury, R. Singh, A. J. Taylor, H.-T. Chen, W. Zhang, and A. K. Azad, Int. J. Opt. 2012 (2012).
  * Longhi [2021] S. Longhi, Phys. Rev. B 104, 125109 (2021).
  * Kawabata  _et al._ [2020] K. Kawabata, M. Sato, and K. Shiozaki, Phys. Rev. B 102, 205118 (2020).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nat. Phys. 14, 925 (2018).
  * Suthar  _et al._ [2022] K. Suthar, Y.-C. Wang, Y.-P. Huang, H.-H. Jen, and J.-S. You, ArXiv preprint arXiv:2202.12675 (2022).
  * Schindler and Prem [2021] F. Schindler and A. Prem, Phys. Rev. B 104, L161106 (2021).
  * Yang  _et al._ [2022b] F. Yang, C. S. Prasad, W. Li, R. Lach, H. O. Everitt, and G. V. Naik, Nanophoton. 11, 1159 (2022b).
  * Ghatak and Das [2019] A. Ghatak and T. Das, J. Phys.: Conden. Matter 31, 263001 (2019).
  * Yokomizo and Murakami [2019a] K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123, 066404 (2019a).
  * Zhang  _et al._ [2020] K. Zhang, Z. Yang, and C. Fang, Phys. Rev. Lett. 125, 126402 (2020).
  * Tai and Lee [2023] T. Tai and C. H. Lee, Phys. Rev. B 107, L220301 (2023).
  * Lee and Thomale [2019] C. H. Lee and R. Thomale, Phys. Rev. B 99, 201103 (2019).
  * Yokomizo and Murakami [2019b] K. Yokomizo and S. Murakami, Phys. Rev. Lett. 123, 066404 (2019b).