# Uncovering Exceptional Contours in non-Hermitian Hyperbolic Lattices

Nisarg Chadha  Undergraduate Programme, Indian Institute of Science, Bangalore
560012, India  Awadhesh Narayan<EMAIL_ADDRESS>Solid State and
Structural Chemistry Unit, Indian Institute of Science, Bangalore 560012,
India

###### Abstract

Hyperbolic lattices are starting to be explored in search of novel phases of
matter. At the same time, non-Hermitian physics has come to the forefront in
photonic, optical, phononic, and condensed matter systems. In this work, we
introduce non-Hermitian hyperbolic lattices and elucidate its exceptional
properties in depth. We use hyperbolic Bloch theory to investigate band
structures of hyperbolic lattices in the presence of non-Hermitian on-site
gain and loss as well as non-reciprocal hopping. Using various analytical and
numerical approaches we demonstrate widely accessible and tunable exceptional
points and contours in {10,5} tessellations, which we characterize using phase
rigidity, energy scaling, and vorticity. We further demonstrate the occurrence
of higher-order exceptional points and contours in the {8,4} tessellations
using the method of Newton polygons, supported by vorticity and phase rigidity
computations. Finally, we investigate the open boundary spectra and densities
of states to compare with results from band theory, along with a demonstration
of boundary localization. Our results unveil an abundance of exceptional
degeneracies in hyperbolic non-Hermitian matter.

##  I Introduction

Spaces with negative curvature emerge naturally in general relativity
Friedmann (1979), and find applications in graph theory Brooks (1996), random
walks Monthus and Texier (1996), complexity theory Krioukov  _et al._ (2010),
and quantum information theory Breuckmann and Terhal (2016); Breuckmann  _et
al._ (2017). The enhanced bulk connectivity of the system Mertens and Moore
(2017) makes hyperbolic surfaces an efficient candidate for data storage and
communication, making these geometries ubiquitous in data science and
electrical engineering. Recent works on networks of coupled microwave
resonators and superconducting qubits Kollár  _et al._ (2019); Chen  _et al._
(2023) provide more tangible insight into these ideas from the perspective of
band theory.

The theoretical extension of ideas from band theory and topology to hyperbolic
geometries Maciejko and Rayan (2021) and their experimental realisations
Kollár  _et al._ (2019); Chen  _et al._ (2023) on tabletop experimental
platforms have propelled hyperbolic lattices to the forefront in the search
for novel phases of matter. Progress in circuit quantum electrodynamics has
made it possible to simulate hyperbolic lattices and explore band theory on
these models, therefore being a direct probe of the influence of the geometry
and metric on the properties of the system. This has motivated the study of
hyperbolic lattices, which are regular tilings of hyperbolic space. Tilings of
regular $p$-gons with coordination number $q$ are denoted in the Schläfli
notation as {$p$,$q$}.

The negative curvature of hyperbolic surfaces endows interesting properties to
these systems, such as the finite ratio of boundary sites to the total sites
and the non-Abelian nature of the translation group.

Figure 1: Visualising the hyperbolic plane in two dimensions through the
Poincaré disk model. (a) The {8,3} tessellation with different coloured sites
denoting different epochs generated recursively. The area in green shows the
16-site unit cell, also called the Bolza lattice, which itself is an {8,8}
tessellation. (b) The {10,5} tessellation up to the first four epochs, along
with a zoomed-in image showing the sites near the boundary.

Recent works have investigated hyperbolic analogues of topological insulators
Urwyler  _et al._ (2022); Liu  _et al._ (2022). The use of real-space
topological invariants on flat projections of hyperbolic tessellations has
shown Chern insulator phases accompanied by topological edge modes and a
quantized conductance Liu  _et al._ (2022). Topological Anderson insulator
phases have also been shown to arise in the presence of symmetry-preserving
disorder. Furthermore, higher-order topological insulator phases have been
discovered in 2D hyperbolic tessellations Liu  _et al._ (2023); Tao and Xu
(2023). These phases show zero-dimensional corner edge states whose degeneracy
depends on the symmetry of the crystal. The realisation of hyperbolic lattices
with arbitrary rotational symmetries allows these corner edge states to have
higher degeneracies than possible for Euclidean systems. The effects of flux
threading through the hyperbolic lattice have also been studied Mosseri  _et
al._ (2022); Stegmaier  _et al._ (2022). Remarkably, a universality in the
shape of the Hofstadter butterfly is obtained when the unit cells in the
compactified manifold are threaded by a flux Stegmaier  _et al._ (2022). These
Hofstadter butterfly patterns turn out to be independent of the coordination
number $q$ for tessellations of the same $p$-gons.

Unlike Euclidean lattices, hyperbolic lattices show a finite fraction of
boundary sites in the thermodynamic limit. Consequently, an open boundary
analysis gets prominent contributions from the boundary sites. As a result,
open boundary spectra show considerable deviations from the bulk spectra even
in the thermodynamic limit. Refs. Maciejko and Rayan (2022); Stegmaier  _et
al._ (2022); Chen  _et al._ (2023) have implemented a compactification of the
hyperbolic manifold using regular maps to identify boundary sites to enforce
periodic boundaries on a finite to obtain the true bulk spectra.

In a parallel development, non-Hermitian systems have been in the limelight in
recent years due to their intriguing fundamental properties, with a potential
for interesting applications Ashida  _et al._ (2020); El-Ganainy  _et al._
(2018); Bergholtz  _et al._ (2021); Ghatak and Das (2019); Banerjee  _et al._
(2023). Beginning with the pioneering ideas of Bender and co-workers Bender
and Boettcher (1998); Bender (2007), the interplay with topology has
reinvigorated the field and has led to several seminal discoveries. The role
of topology in non-Hermitian systems is now at the forefront of research, with
implications for photonics, optics, metamaterials, topoelectric circuits and
many others. Exceptional points (EPs) are an important feature of non-
Hermitian systems. These are singularities where eigenvalues and eigenvectors
become identical Heiss (2012); Miri and Alù (2019); Özdemir  _et al._ (2019),
and have no counterpart in Hermitian systems. The order $n$ of the EP is the
number of eigenvectors coalescing together. In this regard, higher-order EPs
($n>$2) have attracted interest Hodaei  _et al._ (2017); Chen  _et al._
(2017); Wiersig (2020) as higher-order branch points in the energy spectrum.
Not only are they interesting from a fundamental point of view – with
remarkable associated properties such as Berry phases, Riemann sheet
structures, and bulk Fermi arcs – they are also beginning to be exploited in
applications. For instance, EPs have been used to design sensors with enhanced
sensitivity Wiersig (2014, 2020); De Carlo  _et al._ (2022). Another
intriguing phenomenon which has been recently identified in non-Hermitian
systems is the non-Hermitian skin effect (NHSE), where a macroscopic fraction
of states localize at the boundary Yao and Wang (2018); Alvarez  _et al._
(2018); Yokomizo and Murakami (2019); Borgnia  _et al._ (2020); Zhang  _et
al._ (2022). The NHSE is another topic of very active research in the past
couple of years, whose implications are only beginning to be understood. It
features intriguing connections to spectral topology, and a plethora of
experimental platforms have been used to implement its phenomenology Zhang
_et al._ (2022); Lin  _et al._ (2023).

In this work, we explore the interplay of non-Hermiticity and hyperbolic
geometry by using analytic band theory as well as numerical tight-binding
calculations. We find that the non-Abelian nature of the translational group
leads to a greater degree of freedom in parameter space, allowing the tuning
of parameters to readily obtain EPs as well as their higher dimensional
analogues, exceptional contours. We investigate the behaviour of EPs in the
presence of different kinds of non-Hermitian interactions and characterise
their properties.

The {10,5} tessellation lies in the same infinite crystalline family of
$\\{2(2g+1),2g+1\\}$ tessellations (with $g$ being the genus) as the honeycomb
lattice in Euclidean space, and has been termed ”hyperbolic graphene”
Boettcher  _et al._ (2022). Just as in graphene, the {10,5} tessellation has a
two-site unit cell, making it a convenient, analytically tractable two-band
model to study the behaviour of exceptional contours. However, the presence of
only two bands precludes higher-order EPs. To demonstrate the tunability of
higher-order EPs, we use the {8,4} lattice, which has a four-site unit cell.
We see that this four-band system gives us a simple geometry to realize second
as well as fourth order EPs.

The rest of this work is organised as follows: We briefly introduce the model
for projecting hyperbolic space onto Euclidean space in Sec. II. In Sec. III,
we introduce the tight-binding model on the {10,5} lattice and discuss the
spectra obtained for the Bloch Hamiltonian under different non-Hermitian
perturbations. Non-Hermiticity is introduced using gain and loss and non-
reciprocal hoppings in the unit cell. Sec. IV discusses the EPs obtained for
the non-Hermitian Hamiltonian. Phase rigidity and vorticity (spectral winding)
are used to characterize the behaviour of the exceptional contours by
quantifying the coalescing of the eigenstates. The {8,4} lattice with a four-
site unit cell provides a tractable model to show the onset of
higher(fourth)-order EPs in a hyperbolic lattice. We also utilize the recently
proposed method of Newton polygons Jaiswal  _et al._ (2023) to obtain the
conditions for higher-order EPs. In Sec. V, we construct the real-space
lattice and use exact diagonalisation to obtain the energy spectra and the
densities of states. Deviations in the results obtained using hyperbolic band
theory and exact diagonalisation are discussed. We summarise our work in Sec.
VI and motivate experimental realizations of hyperbolic non-Hermitian
lattices.

##  II Understanding Hyperbolic Geometry

We briefly summarize here essential concepts in hyperbolic geometry as a
foundation for the rest of our results. We will be interested in space-filling
tilings of regular $p$-gons, with each vertex having a coordination number of
$q$. This is denoted in the Schläfli notation as a {$p$,$q$} tessellation. In
this notation, the square lattice is a {4,4} tessellation, and similarly, the
honeycomb and triangular lattices are the {6,3} and {3,6} tessellations,
respectively. Due to the restriction imposed by the angle sum property, these
are the only permissible tessellations for Euclidean space. The absence of
this restriction allows hyperbolic systems to have an infinite number of
realisations of {$p$,$q$} tessellations where $(p-2)(q-2)>4$.

2D hyperbolic spaces are often represented using projections onto Euclidean 2D
spaces endowed with a non-Euclidean metric. One such model for projection is
the Poincaré disk model. This projects the hyperbolic space onto a unit disk
with distances measured in the Poincaré metric and is useful for compact
depictions of large lattice sizes as shown in Fig. 1 for the {8,3} and {10,5}
tessellations. We have discussed the procedure to generate the Poincaré model
for hyperbolic lattices in Appendix A.

The lattice translation operators for hyperbolic lattices are elements of the
discrete non-Abelian Fuchsian symmetry groupYuncken (2011). This marks a
striking deviation from Euclidean lattices, where the translation operators
commute, and leads to the presence of higher-dimensional irreducible
representations(irreps) for the symmetry groups of hyperbolic lattices.

Since Abelian groups can only have 1D irreps, the eigenstates of Euclidean
systems with translational symmetry transform with a U(1) phase under the
action of the translation operators. This is the statement of Bloch’s theorem.
Clearly, we cannot use Bloch’s theorem for an arbitrary hyperbolic lattice,
since there will be eigenvectors which transform as higher-dimensional irreps
of the symmetry group. Nevertheless, automorphic band theory has been used for
hyperbolic lattices to describe the eigenstates that transform as a U(1)
irrep.

Thus, hyperbolic Bloch theory is expected to be incomplete due to the presence
of eigenvectors which do not transform as 1D-irreps. We shall subsequently
visit the deviation between results from open boundary numerical
diagonalisation and automorphic Bloch theory for different non-Hermitian
parameters in Section V.

Figure 2: Visualising the energy spectra for Hermitian hyperbolic graphene.
(a) The energy spectrum with $M=0$ produces the characteristic Dirac cones
when $k_{3}=0$, $k_{4}=\pi$ in which case
$h(\mathbf{k})=1+e^{ik_{1}}+e^{ik_{2}}$ becomes the phase factor for Euclidean
graphene. (b) Nodal lines are obtained with $M=0$ for $k_{3}=2\pi/3$ and
$k_{4}=4\pi/3$, and the condition for band touching becomes
$k_{1}=k_{2}\pm\pi$. Modulating $k_{3}$, $k_{4}$ gives different shapes of
nodal surfaces. (c) The surface represents the allowed values of
$k_{1},k_{2},k_{3}$ at which there is a $k_{4}$ that produces a node. The
lines shown in red along the surface are the values of $k_{1},k_{2},k_{3}$
where $k_{4}=0$ gives a node. (d) Adding an on-site potential ($M\neq 0$)
opens a gap in the system. For $k_{3}=0$, $k_{4}=\pi$ a gapped spectrum
($\Delta E_{g}=2M$) is obtained. The linear dispersion near the extrema is
replaced by a quadratic dispersion ($|E|=\frac{(\Delta k)^{2}}{2M}$). (e) For
$k_{3}=2\pi/3$, $k_{4}=4\pi/3$, a gapped system is obtained with the linear
scaling being replaced by a quadratic scaling as in (d). (f) Energy spectrum
through the Brillouin zone. The points defined here are $\Gamma$(0,0,0,0),
A$(\pi,0,0,0)$, B$(\pi,0,\pi,-\pi)$, C$(2\pi/3,-2\pi/3,0,\pi)$, and
D$(\pi,\pi,\pi,\pi)$.

##  III The Hyperbolic Tight-Binding Model

To implement the Bloch ansatz, we need to obtain the unit cell that can
generate the entire lattice through the discrete translations. The imposition
of periodic boundary conditions leads to the identification of pairs of edges
such that the lattice becomes a closed surface with only one face(such as the
torus in 2D). The generalization of the notion of a unit cell to hyperbolic
lattices was carried out in Ref. Sausset and Tarjus (2007); Boettcher  _et
al._ (2022) to obtain the conditions on {$p$,$q$} such that we can obtain a
closed primitive unit cell by imposing periodic boundaries. This leads to the
condition that the identified unit cell must correspond to a surface with only
one face, giving rise to constraints on $p$ and $q$. Ref. Boettcher  _et al._
(2022) detail the families of lattices($\\{p,q\\}$) whose Bravais lattices
have a unit cell.

Of interest to us is the $\\{2g+1,2(2g+1)\\}$ family, where $g$ is the genus
of the Brillouin zone. $g=1$ gives us the honeycomb lattice, and $g=2$ gives
the {10,5} lattice. Like graphene({6,3}), we obtain two sublattices, each
forming a $\\{2g+1,2(2g+1)\\}$ tessellation of its own. In Sec. IV.2, we shall
also use the {8,4} lattice, which is part of the $\\{4g,4\\}$ family with a
$2g$ site unit cell.

To proceed with the hyperbolic Bloch ansatz in a tight binding model, we
choose a unit cell and identify the nearest neighbours in terms of the
translation operators. These have been worked out for the {10,5} tessellation
in Ref. Chen  _et al._ (2023) and have been discussed in Appendix B for a
self-consistent narrative.

We introduce a simple tight-binding Hamiltonian with uniform hoppings between
nearest neighbours and an on-site potential. This Hamiltonian can be written
in the second quantized form as:

$H=\sum_{i}(Mc_{i,A}^{\dagger}c_{i,A}-Mc_{i,B}^{\dagger}c_{i,B})-t(\sum_{\langle
i,j\rangle}c_{i,A}^{\dagger}c_{j,B}+c.c.)$ (1)

Here $c_{i,A(B)}$ represents the annihilation operator at the A(B) sub-lattice
of the i-th site in the lattice. Likewise, $c_{i,A(B)}^{\dagger}$ represents
the creation operator at the A(B) sub-lattice of the i-th site in the lattice.

$M$ defines an on-site potential with opposite signs on either sub-lattice,
and $t$ is the strength of the inter-site hopping. We restrict ourselves to
nearest neighbour hoppings to get an analytically tractable model which is
also motivated by the circuit QED setups implemented in experimental
realizations Chen  _et al._ (2023). We shall show that this model is
sufficient to provide an abundance of exceptional contours in the presence of
non-Hermitian perturbations. Just like in Euclidean space, the $k$-space
Hamiltonian can be obtained using a U(1) Bloch theory by attaching a phase of
$e^{ik_{j}}$ corresponding to every translation operator $\gamma_{j}$ applied
to perform the hopping. This procedure has been detailed in App B to get the
Bloch Hamiltonian in the sublattice basis as:

$\hat{H}_{0}=\begin{pmatrix}M&-t(1+e^{ik_{1}}+e^{ik_{2}}+e^{ik_{3}}+e^{ik_{4}})\\\
-t(1+e^{-ik_{1}}+e^{-ik_{2}}+e^{-ik_{3}}+e^{-ik_{4}})&-M\end{pmatrix},$ (2)

For the rest of this work, we set $t=1$. We can get the tight-binding graphene
Hamiltonian from this model by setting $k_{3}=0$ and $k_{4}=\pi$. We note that
tuning $k_{3}$ and $k_{4}$ provide additional degrees of freedom to the off-
diagonal term:

$h(\mathbf{k}):=1+e^{ik_{1}}+e^{ik_{2}}+e^{ik_{3}}+e^{ik_{4}}$ (3)

For ease of representation, we will treat the $\mathbf{k}$-vectors as
parameters to study the positions of EPs and nodal points in cross-sections of
the 4D Brillouin zone. The simplicity of the Hamiltonian will be helpful to
analytically study the onset and behaviour of exceptional contours upon
introducing non-Hermiticity.

Figure 3: Energy spectra for the non-Hermitian cases with gain and loss
($\delta$), and non-reciprocity ($\Gamma$). (a)-(c) The shape of the nodal
contour on the projected $k_{1}-k_{2}$ plane with varying $k_{3}$ and $k_{4}$
for the case of gain and loss. (a) For $\delta=1$, $k_{3}=k_{4}=0$ gives the
limiting case of a single band intersection point $\vec{k_{0}}=(\pi,\pi,0,0)$.
The dispersion near this point is quadratic in $\mathbf{k}-\mathbf{k_{0}}$.
(b) A single nodal loop is obtained for $\delta=1,k_{3}=k_{4}=\pi$. The
spectrum scales as $E\thicksim(|\mathbf{k}-\mathbf{k_{0}}|)^{\frac{1}{2}}$
near the nodal contour. (c) Nodal contours obtained for
$\delta=1,k_{3}=7\pi/6,k_{4}=\pi/4$ (d) The energy spectrum on the same path
(but with $\delta=1$) as in Fig. 2 (f). The points where both the energy
levels become degenerate are marked with hollow circles, and occur when
$|h(\mathbf{k})|=\delta$. (e)-(g) The energy bands and corresponding nodal
contours for different choices of parameters $(\Gamma,k_{3},k_{4})$, i.e.,
non-reciprocal hopping case. (e) Similar to (a), one obtains a single nodal
point where the bands touch. The choice of parameters is
$\Gamma=0.5,k_{3}=\arccos(3/4),k_{4}=2\pi-k_{3}$. The dispersion is quadratic
in $k$ near the nodal point. (f) Absence of nodal points, with a finite
complex gap in the spectrum for parameter values
$\Gamma=0.7,k_{3}=k_{4}=\pi/2.$ (g) For
$\Gamma=0.5,k_{3}=4\pi/3,k_{4}=4\pi/5$, open nodal contours are obtained,
where the spectrum scales as
$E\thicksim(|\mathbf{k}-\mathbf{k_{0}}|)^{\frac{1}{2}}$, where $k_{0}$ is a
point on the nodal contour. (h) The energy spectrum is shown on the above-
mentioned path, with $\Gamma=1$. There are two nodal points where both the
real and imaginary parts of the eigenvalues become equal.

###  III.1 Hermitian spectra

We briefly summarize the properties of the Hermitian model in Eq. (2). In the
Hermitian regime, the parameters at our disposal are the on-site potential
$\pm M$ and the hopping strength $t$ (assumed to be isotropic). The energy
spectrum is similar to Euclidean graphene and is given by Eq. (4), with the
additional freedom furnished by $k_{3}$ and $k_{4}$ providing additional
tuning parameters on the $k_{1}$-$k_{2}$ cross-section.

$E_{\pm}=\pm\sqrt{M^{2}+|h(\mathbf{k})|^{2}}$ (4)

Plotting cross-sections for the energy spectrum with $M=0$ for different
choices of $k_{3},k_{4}$ can give nodal points as well as nodal lines as shown
in Fig. 2(a)-(b). The nodal contours are surfaces in the Brillouin zone where
the bands touch. For a Hermitian system, the nodal surface corresponds to
points of degeneracies where both the eigenvalues are equal. Such nodal
surfaces have been shown for the ($k_{1}$,$k_{2}$,$k_{3}$) space in Ref. Chen
_et al._ (2023). For a non-Hermitian system, these correspond to exceptional
nodal contours (see Fig. 4), where the eigenvectors also coalesce. The locus
of the nodal surface is given by $h(\mathbf{k})=0$. The nodal surface is
visualised in Fig. 2(c), where the surface represents values of
$k_{1},k_{2},k_{3}$ for which there exists a node for some value of $k_{4}$.
Fig. 2(f) shows the variation of the energy spectrum through different points
in the four-dimensional Brillouin zone.

A sublattice potential $M$ introduces a gap (equal to 2$M$) in the system, and
the linear dispersions at the Dirac points are replaced by quadratic
dispersions at the band extrema, as shown in Fig. 2(d)-(e).

###  III.2 Introducing non-Hermiticity

We will introduce non-Hermiticity through on-site gain and loss (terms
proportional to $\sigma_{z}$) and non-reciprocal hopping (terms proportional
to $\sigma_{x,y}$), by adding terms of the form $i\delta\sigma_{z}$,
$i\Omega\sigma_{x}$, and $i\Gamma\sigma_{y}$, to the Hamiltonian $\hat{H}_{0}$
in Eq. (2). For simplicity, we will consider the non-Hermitian terms to be
independent of $k$. This implies that the non-reciprocity is contained within
the unit cell. We will see how this anisotropy in the non-Hermiticity gives a
rich complex energy spectrum with an abundance of EPs.

Further, the higher dimensionality of the $k$-space endows greater freedom to
tune the parameters of the system, allowing us to realise rich non-Hermitian
phenomena. This freedom will have an even greater significance in obtaining
higher-order EPs for tessellations with larger unit cells, such as the four
band {8,4} model, which will be discussed later.

Figure 4: Exceptional nodal surfaces for non-Hermitian hyperbolic graphene.
The figures display the nodal surface on the $k_{1}$,$k_{2}$,$k_{3}$ subspace
of the Brillouin zone where energy eigenvalues vanish. (a) The nodal surface
for the case of on-site gain and loss with $\delta=0.5$. In this case, the
condition for a nodal surface is $|h(\mathbf{k})|=|\delta|$. The pink region
shows the regions where there exists a node for any value of
$k_{4}\in[0,2\pi]$. The blue subregion is the nodal surface with
$k_{4}=\pi/2$. (b) The nodal surface for the case of non-reciprocal hoppings
with $\Gamma=0.5$. The condition for nodes is $h(\mathbf{k})=\pm\Gamma$. The
nodal surface becomes a two-dimensional surface compared to the three-
dimensional volume in (a). This is due to the additional constraint
$2i\Gamma\mathrm{Im}(h)=0$ in Eq. 9, thus reducing the degrees of freedom of
the spectrum.

####  III.2.1 On-Site Gain and Loss

First, we add an on-site gain and loss of strength $\delta$, by adding
$i\delta\sigma_{z}$ to $\hat{H}_{0}$ in Eq. (2). As a result of the non-
Hermiticity, the spectrum is no longer entirely real. In fact, the eigenvalues
$E_{\pm}=\pm\sqrt{|h(\mathbf{k})|^{2}-\delta^{2}}$ will be purely real or
purely imaginary.

For a non-zero $M$, the square root has a constant non-zero imaginary part
which prohibits any band touchings, even in the individual real and imaginary
spectra. This can be seen by the replacement $\delta\rightarrow\delta-iM$ in
the energy dispersion. The locus of the real and imaginary energies ($E_{r}$
and $E_{i}$) is the intersection of three hyperbolic surfaces, given by,

$\displaystyle E_{r}E_{i}$ $\displaystyle=M\delta,$ (5) $\displaystyle
E_{r}^{2}-E_{i}^{2}$ $\displaystyle\leq M^{2}+25-\delta^{2},$ (6)
$\displaystyle E_{r}^{2}-E_{i}^{2}$ $\displaystyle\geq M^{2}-\delta^{2}.$ (7)

Clearly, increasing the non-Hermiticity strength $\delta$ will lead to a more
prominent imaginary spectrum. For $M=0$, we find that either the real or the
imaginary part of the eigenvalues will vanish, and we can recover nodal
contours where the eigenvalues are zero, with the cross-sections shown in Fig.
3(a)-(d) for some choices of $k_{3},k_{4}$. Fig. 4(a) shows the surface in the
$k_{1}-k_{2}-k_{3}$ subspace with $\delta=0.5$ for which a node can be
obtained for some value of $k_{4}$ in the Brillouin zone.

The nodal contours mark the transition from purely real to purely imaginary
eigenvalues. In Fig. 3(d), a one-dimensional contour is parametrised in the
Brillouin zone, which shows the appearance of EPs (circled in black). Contours
connecting these EPs are non-Hermitian Fermi arcs characterised by purely real
or imaginary energy eigenvaluesBergholtz  _et al._ (2021). While these arcs
occur trivially due to the behaviour of the band structure for this case, we
will observe their existence even for non-reciprocal hopping in the next
section. In the limit $\delta>5t$, the spectrum becomes purely imaginary, and
the bands do not touch in the four-dimensional Brillouin zone.

####  III.2.2 Non-Reciprocal Hopping

Another way to introduce non-Hermiticity is through non-reciprocal hopping
within the unit cell. This can be done by adding terms of the form
$i\Omega\sigma_{x}$ or $i\Gamma\sigma_{y}$. These terms produce an imaginary
$k$-dependent term in the expression for the eigenvalues given by

$i\Omega\sigma_{x}:E_{\pm}=\pm\sqrt{|h|^{2}+M^{2}-\Omega^{2}+2i\Omega\mathrm{Re}(h)},\\\
$ (8)
$i\Gamma\sigma_{y}:E_{\pm}=\pm\sqrt{|h|^{2}+M^{2}-\Gamma^{2}-2i\Gamma\mathrm{Im}(h)}.$
(9)

To obtain a node, both the real and imaginary parts in the square root must go
to zero. The shape of the nodal surface is very sensitive to non-reciprocity
and changes drastically upon the addition of a small non-reciprocal hopping,
as shown in Fig. 4(b). Typical two-dimensional cross-sections of the
eigenspectra are shown in Fig. 3(e)-(h), with different nodal structures,
which could be a nodal point (e), a nodeless cross-section (f), or a nodal
contour (g). Fig. 3(h) shows the variation of the spectrum on a one-
dimensional contour in the Brillouin zone. The appearance of the nodal EPs is
shown in black, where both the real and imaginary parts of the eigenvalues are
equal. These nodes are connected through non-Hermitian Fermi arcs where the
real part of the energy goes to zero, as shown by our choice of contour. For
the case of the $i\Gamma\sigma_{y}$, the condition for the appearance of nodes
is obtained to be

$\displaystyle|h|^{2}=\Gamma^{2}-M^{2},$
$\displaystyle\Gamma\mathrm{Im}(h)=0.$

Therefore, for non-zero $\Gamma$, we require $h=\pm\sqrt{\Gamma^{2}-M^{2}}$.
The requirement for the imaginary part to be zero results in lower dimensional
nodal surfaces on the $k_{1},k_{2}$ cross-section as shown in Fig. 4(b).
Similarly, the nodal spectrum for non-reciprocity $i\Omega\sigma_{x}$ is
obtained by setting $h=\pm i\sqrt{\Omega^{2}+M^{2}}$.

##  IV Diagnostics of Exceptional Contours

Figure 5: Diagnosing exceptional points using phase rigidity. (a) The energy
spectrum as a function of $\delta$ for $\mathbf{k}=(\pi/2,\pi/3,3\pi/4,\pi)$.
The bands touch each other when $\delta=\delta_{EP}=|h(\mathbf{k})|\approx
2.58145.$ (b) The scaling of the energy very close to the EP
($\delta=\delta_{EP}+\epsilon$), where $\epsilon$ is a perturbation from the
EP). The inset shows the logarithmic scale plot, which has a slope of 0.5,
revealing a square root dependence on $\epsilon$ characteristic of second-
order EPs. (c) The magnitude of the phase rigidity, $r$, as a function of
$\delta$ for the same k. (d) For a system with both non-reciprocity ($\Gamma$)
and on-site gain and loss ($\delta$) parametrised by
$\delta=r_{c}\cos\theta,\Gamma=r_{c}\sin\theta$, the spectrum is independent
of $\theta$. For $k_{1}=4\pi/3,k_{2}=0,k_{3}=2\pi/3,k_{4}=0$, the condition
for obtaining a gapless point is $r_{c}=2$. The phase rigidity is shown to go
to zero as a function of $r_{c}$, with $\theta=\pi/4$. (e) The circular
exceptional contour in the $\delta-\Gamma$ parameter space for
$\mathbf{k}=(0,0,2\pi/3,4\pi/3)$. Figure 6: Inter-Band winding for the {10,5}
tessellation. (a) Illustration of the winding of energy levels for the {10,5}
tessellation with the on-site gain and loss. For the parameter choice
$k_{1}=2\pi/3,k_{2}=\pi/3,k_{3}=0,k_{4}=0$ and $\delta=2$, the Hamiltonian is
not defective, and we get trivial winding of the energy levels with respect to
a parametric perturbation $e^{i\theta_{pert}}\sigma_{x}$ added to
$\hat{H}_{0}$. (b) $\delta=2\sqrt{3}$ makes the Hamiltonian rank-deficit and
gives an EP at the chosen $\vec{k}$-point. In this case, we can see the inter-
band winding visualised through the exchange of the energy levels in one
complete cycle of the perturbation.

###  IV.1 Identifying Exceptional Points

For non-Hermitian Hamiltonians, the usual orthogonality of eigenvectors is
replaced by the weaker condition of _biorthogonality_ Brody (2013). This gives
us the right and left eigenkets, defined as eigenvectors of $H$ and
$H^{\dagger}$ respectively obeying the biorthogonality condition

$\braket{L_{i}}{R_{j}}=\delta_{ij},$ (10)

where $\delta_{ij}$ is the Kronecker delta function, and the eigenvectors are
ordered such that $H\ket{R_{i}}=E_{i}\ket{R_{i}}$ and
$H^{\dagger}\ket{L_{i}}=E_{i}^{*}\ket{L_{i}}$. EPs represent branch points in
the spectral manifold and are endowed with topological properties as a
consequence of the spectral winding around the branch point Bergholtz  _et
al._ (2021); Ding  _et al._ (2022). At an EP, two or more eigenvectors
coalesce. Hence, Eq. (10) falls through since if $\ket{R_{i}}=\ket{R_{j}}$
then their inner products with $\bra{L_{i}}$ (expected to be 1 and 0
respectively) must also be equal, contrary to the statement of the
biorthogonality. Thus, at an EP, the state $\ket{R_{i}}$ cannot be left
normalized to unity since $\braket{L_{i}}{R_{i}}=0$. The phase rigidity, $r$,
measures this biorthogonality to quantify the approach to an EP Bulgakov  _et
al._ (2006); Eleuch and Rotter (2016). The phase rigidity, $r_{j}$ for the
j-th eigenstate is defined as Banerjee  _et al._ (2023)

$r_{j}=\frac{\braket{L_{j}}{R_{j}}}{\braket{R_{j}}{R_{j}}}.$ (11)

For a Hermitian system, $r_{j}=1$, since $H=H^{\dagger}$, implying the left
and right eigenkets are equal. On the other hand, $r_{j}=0$ when $\ket{R_{j}}$
coalesces with another eigenstate at an EP. Therefore, as we increase the non-
Hermiticity, the left and right eigenkets become increasingly orthogonal,
leading to a reduction of $r$. Remarkably, phase rigidity is not just a
theoretically defined quantity, but has also been experimentally measured Tang
_et al._ (2020).

EPs lead to degeneracies in the spectrum, which in our case can only occur if
the eigenvalues are zero (since the spectrum is symmetric about $E=0$). Thus,
the nodal surfaces that we obtained earlier correspond to the intersections of
the Fermi arcs for the real and imaginary parts of the energy leading to
second-order EPs. We will subsequently investigate the occurrence of higher-
order EPs in a four-band model for the {8,4} tessellation where we shall see
EPs with all four eigenvectors coalescing.

For second-order EPs, the eigenspectrum shows a square-root scaling, which is
displayed in Fig. 5(a)-(b). The complex eigenspectrum of the non-Hermitian
system allows the definition of the winding of energy eigenvalues, called the
vorticity, as a topological invariant Bergholtz  _et al._ (2021). EPs are
characterized by fractional values for the vorticity $\nu_{n}$ of the $n$-th
band as a consequence of the coalescing of more than two eigenvectors Ding
_et al._ (2022).

$\nu_{n}=\frac{1}{2\pi}\oint_{\mathcal{C}}\nabla_{\textbf{k}}\mathrm{arg}(E_{n}(\mathbf{k}))\cdot
d\mathbf{k}.$ (12)

On a contour $\mathcal{C}$ encircling the EP, the branch cut causes a net
winding of the coalescing eigenvectors. As a consequence, the sum of the
spectral windings (called the eigenvalue winding number Ding  _et al._ (2022))
of the eigenvectors is unity ($\Sigma_{j}\nu_{j}=1$).

The phase rigidity and vorticity quantify the behaviour of eigenvectors and
eigenvalues, respectively, near an EP. An $n$-th order EP creates a branch
point associated with the $n$-th complex root in parameter space. Assuming the
Hamiltonian near an EP to be a small perturbation from the Hamiltonian at the
EP, the eigenvalues and eigenvectors can be calculated using a non-Hermitian
perturbation theory Budich and Bergholtz (2020); Wiersig (2023). This yields
the branch point in the eigenvalue dispersion which we quantify using
vorticity. Similarly, the phase rigidity (which vanishes at the EP), can be
calculated near the EP to obtain either a $1/n$ or $(n-1)/n$ scaling depending
on the order of the EP Jaiswal  _et al._ (2023).

For a contour $\mathcal{C}$ encircling the EP, $\nu_{\pm}=1/2$ for both bands
at the second order EP in our model, giving a net winding of 1. This is shown
in terms of a parametrisation of a circular path centred at an EP in terms of
the vorticity in Fig. 6, and leads to an exchange of the eigenvalue and
eigenvectors across the branch cut.

Figure 7: Inter-band winding and higher-order EPs in {8,4} tessellation. (a)
The eigenvalue winding for the second order EP. The bands in red show an
exchange, whereas the bands in blue return to their original values under a
cyclic perturbation with $\epsilon=0.07e^{i\theta_{pert}}$ added as shown in
Eq. (15). This gives a second-order EP. The parameters taken here are
$k_{1}=\pi/4,k_{2}=2\pi/3,k_{3}=3\pi/5,k_{4}=61\pi/60$ and $\delta\approx
2.07923$. (b) The fourth order EP, where
$k_{1}=k_{2}=0,k_{3}=k_{4}=\arccos{(2\sqrt{5}-5)}$ and $\delta\approx
3.14461$. Panel (c) shows the Newton polygon diagram for the {8,4} model
considered above. The hollow circles show the terms that are absent in the
characteristic polynomial. The blue circles are the terms that need to vanish
to obtain a second-order EP, whose bounding line is shown in blue. The red
point is the additional term that needs to vanish to reach a fourth-order EP,
whose bounding line is shown in red. Panels (d) and (e) show the energy
spectrum near the EP for the second and fourth-order EPs, respectively. The
scaling for the spectra is shown in (f) and (g), respectively. The lowest-
order perturbation terms in the spectra have scaling coefficients 0.5 and
0.25, respectively.

###  IV.2 Higher order Exceptional Points in {8,4} Tesselations

In recent years, higher-order EPs – where more than two eigenvalues and
eigenvectors merge – have been of significant interest owing to their
remarkable properties beyond second-order EPs. Here, we next show how
hyperbolic lattice models allow such higher-order EPs in the presence of non-
Hermiticity. As an example, we consider the {8,4} tessellation, which is part
of the {4$g$, 4} infinite family of tessellations with genus 2 Boettcher  _et
al._ (2022). The unit cell has four sites and the Bolza lattice is an {8,8}
tessellation, with the discrete Fuchsian symmetry group generated by four
elements. Therefore, this four-band system gives us a simple system to realize
second and fourth-order EPs. We implement a similar nearest neighbour tight-
binding model using the translation operators for the {8,8} lattice. Adding a
balanced complex gain and loss $\pm i\delta$ on two sites gives the $k$-space
Hamiltonian:

$H_{0}(\mathbf{k})=\begin{pmatrix}i\delta&1+e^{i(k_{1}-k_{2})}&0&e^{ik_{1}}+e^{-ik_{4}}\\\
1+e^{i(k_{2}-k_{1})}&0&1+e^{i(k_{2}-k_{3}})&0\\\
0&1+e^{i(k_{3}-k_{2})}&0&1+e^{i(k_{3}-k_{4})}\\\
e^{-ik_{1}}+e^{ik_{4}}&0&1+e^{i(k_{4}-k_{3})}&-i\delta\\\ \end{pmatrix}.$ (13)

To study the EPs generated in this model, we use the Newton polygon method
elaborated in Ref. Jaiswal  _et al._ (2023). For a contour enclosing an $n$-th
order EP, we obtain a mutual winding of $n$ eigenvalues, leading to a cyclic
permutation of the eigenvalues across the branch cut. At the $n$-th order EP,
the eigenvalues of the perturbed Hamiltonian show a $1/n$ scaling with the
perturbation $\epsilon$. The Newton polygon method uses this scaling relation
in the binomial expansion for the secular equation, given by

$P(\omega,\epsilon)=\det(H+\epsilon
H_{1}-\omega\mathcal{I})=\sum_{n(\epsilon),n(\omega)}p_{n(\epsilon),n(\omega)}\epsilon^{n(\epsilon)}\omega^{n(\omega)},$
(14)

where $\epsilon H_{1}$ represents a perturbation of strength $\epsilon$ added
to the Hamiltonian $H$ with parameters tuned to give an EP. The points
$(n(\epsilon),n(\omega))$ with $p_{n(\epsilon),n(\omega)}\neq 0$ are plotted
in $\mathbb{R}^{2}$, which gives the Newton polygon as the smallest convex
polygon containing all the plotted points. The negative inverse of the slope
of a line with all points on or above it gives the lowest order scaling of
$\omega$ with $\epsilon$ ($1/n$), and thus the order of the EP ($n$).

We implement this method to study the {8,4} lattice. We begin by introducing a
perturbation, $\epsilon$, to the hopping between the first two sites in the
unit cell such that the perturbed Hamiltonian becomes

$H(\mathbf{k})=\begin{pmatrix}i\delta&1+e^{i(k_{1}-k_{2})}+\epsilon&0&e^{ik_{1}}+e^{-ik_{4}}\\\
1+e^{i(k_{2}-k_{1})}+\epsilon&0&1+e^{i(k_{2}-k_{3}})&0\\\
0&1+e^{i(k_{3}-k_{2})}&0&1+e^{i(k_{3}-k_{4})}\\\
e^{-ik_{1}}+e^{ik_{4}}&0&1+e^{i(k_{4}-k_{3})}&-i\delta\\\ \end{pmatrix}.$ (15)

The Newton polygon diagram for our Hamiltonian is shown in Fig. 7 (c). Using
the analytic form for the characteristic polynomial, we can fine-tune
different parameter choices to eliminate the (0,0) and (1,0) points on the
diagram to obtain a second-order EP in the system. Subsequently, we can remove
the (2,0) point to obtain a fourth-order EP. The energy spectra for the four
bands are shown in Fig. 7(d)-(e) for the second and fourth-order EPs,
respectively. The $E_{3},E_{4}$ bands in Fig. 7(d) do not mix and remain
purely real. The scaling for $E_{1},E_{2}$ bands is shown in Fig. 7(f),
revealing a square root dependence on the parameters, as expected for a
second-order EP.

For the fourth-order EP, all four bands participate in the mixing process and
scale more steeply as a function of $\delta$. The logarithmic scaling shows a
$\delta^{1/4}$ dependence on the parameters, as expected for a fourth-order
EP. The obtained scaling is shown for the spectrum near the EPs in both cases
in Fig. 7(f) and (g), thereby confirming the presence of higher(fourth)-order
EPs in this model.

Figure 8: Comparing the density of states between hyperbolic band theory (HBT)
and exact diagonalization (ED). (a) The density of states in a cumulative
probability distribution (The cumulative probability $p(E)$ is defined as the
fraction of states with energy less than or equal to $E$. The cumulative
density of states is $p(E)$.) for the Hermitian $k$-space model with periodic
boundaries and the open boundary system. (b) The cumulative density of
states(p(E)) for the system with $\delta=0.2$. The density is measured with
respect to only the real part of the energy, leading to a spike in the
probability distribution at $E=0$ corresponding to the spectrum along the
imaginary axis as shown in Fig. 8(d). (c) The cumulative density of states for
the system with non-reciprocity $\Omega=0.2$. (d) The complex energy spectra
with on-site gain and loss ($\delta=0.2$). The open boundary spectrum is a
subset of the periodic boundary condition spectrum. . Figure 9: localization
effects with non-reciprocal hopping. The density of states summed over all the
eigenfunctions for the system with two different values of the non-reciprocity
parameter (a) $\Omega=0.2$ and (b) $\Omega=0.4$. Although localization can be
seen at the outer epochs, the profile is not tenfold rotationally symmetric
due to the asymmetry while defining unit cells for non-Hermitian hopping. The
extent and magnitude of boundary localization increase with $\Omega$. The
asymmetry in the localization is more evident for larger values of $\Omega$.
(c) The inverse participation ratio (IPR) for eigenvectors with different
values of non-reciprocity. The macroscopic fraction of boundary sites leads to
a high IPR even for states of the Hermitian system. The absence of an
extensive shift of boundary localised modes with non-reciprocity rules out the
possibility of a skin effect in this class of models.

##  V Implications in real space: Poincaré disk

Till now, we have focused on the reciprocal space picture of the exceptional
contours. Now, we construct a real space system for the {10,5} tessellation
using the process of circular inversion as outlined in Ref. Liu  _et al._
(2022). For the sake of self-consistency, we have included this procedure in
Appendix A. This method allows the construction of successive epochs of the
hyperbolic lattice recursively. The number of sites in each successive epoch
increases exponentially, and we limit the system size to four epochs (which
amounts to 7040 lattice sites).

Defining the sublattice and unit cells is a recursive variation of the dimer
covering problem on the hyperbolic lattice and results in a few dangling sites
at the boundary which are not paired into any unit cell. This manifests itself
as a weak perturbation that breaks the tenfold rotational symmetry of the
lattice.

We numerically compute the density of states with respect to the real part of
the energy. Fig. 8(a) shows the comparison of the density of states obtained
for the Hermitian system under open boundaries and in the $k$-space. The
deviation between Bloch theory and exact diagonalization is due to the
macroscopic fraction of boundary sites for hyperbolic tessellations as well as
the presence of higher-dimensional irreps for the translational symmetry
group. Such deviations between exact diagonalization on finite flakes with
open boundaries and band theory have been demonstrated in Ref. Urwyler  _et
al._ (2022); Chen  _et al._ (2023). Recent works have used continued fraction
expansions Mosseri and Vidal (2023), supercells Lenggenhager  _et al._ (2023),
and converging periodic boundary conditions Lux and Prodan (2023a, b) to
remove spurious contributions from dangling edges and get true bulk spectra
from a finite-sized lattice. Recently, the process of using regular maps to
obtain a boundary-removed periodic spectrum has been proposed Stegmaier  _et
al._ (2022); Chen  _et al._ (2023) to get better agreements with the $k$-space
density of states and remove the effect of the dominating boundary
contribution. However, we will use open boundary spectra to probe localization
effects in the non-Hermitian hyperbolic system, since we are also interested
in investigating the presence of a non-Hermitian skin effect.

The spectra and density of states for the system with gain and loss are shown
in Fig. 8(b)-(d). The open boundary spectrum obtained by exact
diagonalization(ED), is a subset of the hyperbolic band theory(HBT) spectrum
as seen in Fig. 8(b) – this is consistent with the theorem outlined in Ref.
Okuma  _et al._ (2020). Whereas the local density of states profile does not
show any localization features, the density of states shows a jump at zero
energy due to a large number of states with purely imaginary energy.

Upon adding non-reciprocal hopping, we notice a localization at the outer
epoch of the system, as shown in Fig. 9(a)-(b). The localization increases
with non-reciprocity and has an asymmetric profile due to the method of
defining the unit cells, as discussed above. We use the inverse participation
ratio (IPR) to study this localization by measuring the fraction of wave-
function probability at the boundary. The IPR for an eigenket
$\ket{\psi_{n}}$, is defined as

$IPR=\frac{\sum_{r\in\rm{boundary}}|\braket{\mathbf{r}}{\psi_{n}}|^{2}}{\sum_{r}|\braket{\mathbf{r}}{\psi_{n}}|^{2}}.$
(16)

The numerator calculates the total occupation probability of boundary sites,
normalised by the total probability density. Plotting the IPR for different
values of non-reciprocity in Fig. 9(c) reveals that the non-reciprocity does
not cause a skin effect in the Hamiltonian and is simply a perturbation to the
density profile, enhanced by the microscopic probability density in the wave-
function density profile. In the Hermitian limit under the Bloch ansatz
itself, having a macroscopic number of sites at the boundary leads to states
with a large IPR, as can be seen from the plot with $\Omega=0.0$ in Fig. 9(c).

##  VI Summary and outlook

In this work, we introduce non-Hermiticity to simple nearest neighbour tight-
binding models on hyperbolic lattices to obtain an abundance of exceptional
degeneracies. We assume an automorphic Bloch ansatz for the hyperbolic system
to obtain the k-space Hamiltonian in the 4-dimensional Brillouin zone. The
two-band Bloch Hamiltonian for the {10,5} serves as a convenient and
analytically tractable model to show the exceptional tuning upon the
introduction of on-site gain and loss and non-reciprocal hoppings. Using the
phase rigidity and vorticity, we study the behaviour of the complex
eigenspectrum near the EPs.

To demonstrate the tuning for higher-order EPs, we use the {8,4} lattice with
a four-band Hamiltonian. Using the recently proposed method of Newton’s
polygons, we can fine-tune the parameters to obtain both second and fourth-
order EPs on the {8,4} model.

We compare the results from Bloch theory with results from exact
diagonalisation for a finite-sized flake with OBC, which shows a deviation in
the density of states expected due to the contributions from dangling boundary
sites and higher-dimensional irreps of the symmetry group. A boundary
localization is observed for the flake, which is shown to be different from
the non-Hermitian skin effect and arises from the non-Hermitian perturbation.

Given the recent experimental realisation of hyperbolic lattices Kollár  _et
al._ (2019); Chen  _et al._ (2023) using circuit quantum electrodynamics and
the simulation of non-Hermitian phases using optical and electrical circuit
elements Hodaei  _et al._ (2017); Chen  _et al._ (2017); Wiersig (2020),
engineering non-Hermitian lattices in hyperbolic geometry is an experimentally
viable endeavour. The abundance of EPs can have potential applications in
designing non-Hermitian sensors operating on the scaling of eigenvalues near
EPs Wiersig (2014, 2020); De Carlo  _et al._ (2022). We hope our work
motivates further theoretical and experimental explorations of non-Hermitian
hyperbolic matter.

Note Added: During the final stages of this work, we came across a
complementary preprint, which studies the skin effect in hyperbolic
topological lattice models Sun  _et al._ (2023).

## Acknowledgments

We thank Adhip Agarwala and Vijay Shenoy for illuminating discussions. N.C.
acknowledges a fellowship from the Kishore Vaigyanik Protsahan Yojana (KVPY).
A.N. is supported by the Indian Institute of Science.

##  Appendix A Poincaré disk projection

Projection models such as the Poincaré disk are used to conveniently represent
hyperbolic geometries. Such models are constructed using distance-preserving
maps (isometries) from the hyperbolic space to the representation space, in
this case, $\mathbb{D}$. The Poincaré disk projection is frequently used to
describe hyperbolic lattices since it allows a compact description of the
infinite hyperbolic plane. Distances in the hyperbolic space are calculated
using the Riemannian metric

$ds^{2}=\frac{4|dz|^{2}}{(1-|z|^{2})^{2}}.$ (17)

where $z=x+iy$. The translation operators for the Poincaré disk are isometries
belonging to the Fuchsian group with their operation defined as

$\begin{pmatrix}\alpha&\beta\\\ \gamma&\delta\end{pmatrix}z:=\frac{\alpha
z+\beta}{\gamma z+\delta},$ (18)

where $\alpha,\beta,\gamma,\delta\in\mathbb{R}$ such that
$\alpha\delta-\beta\gamma>0$ Here, we discuss the procedure to construct the
Poincaré disk projection for a hyperbolic tessellation following the method
outlined in Ref. Liu  _et al._ (2022). This method uses recursive circular
inversion to generate outer epochs using inner epochs. The choice of {$p$,$q$}
uniquely fixes the positions of all the sites in the disk $\mathbb{D}$ up to a
global rotation of the system.

  1. 1.

Taking the origin as the centre of the tessellation and the first site to lie
along the positive $x$-axis, the construction of the first epoch can be
trivially carried out as

$z_{n}=d_{pq}e^{2\pi in/p},\quad
d_{pq}=\sqrt{\frac{\cot{\pi/q}-\tan{\pi/p}}{\cot{\pi/q}+\tan{\pi/p}}},$ (19)

for $n=1,2...p$ and the positions of the sites are stored as complex numbers.

  2. 2.

The circular inversion of point A with respect to a circle centred at O with
radius $r$ is defined as the point A’ such that OA and OA’ are collinear with
OA.OA’=$r^{2}$.

  3. 3.

We take pairs of neighbouring sites, say A and B and invert them with respect
to $\mathbb{D}$ to get A’ and B’.

  4. 4.

Inverting the remaining sites in the regular polygon containing A and B with
respect to the circle circumscribing ABB’A’ to get the sites corresponding to
the next epoch.

  5. 5.

This process is repeated for pairs of neighboring sites to construct
subsequent epochs.

##  Appendix B Algebra of the generators of hyperbolic translations

The group structure of translational symmetry operators plays a crucial role
in the band theoretic description of systems. Since translations commute in
Euclidean space, they form an Abelian group which only permits 1D irreps.
Bloch’s theorem states that systems (Hamiltonians) with translational
symmetries will have eigenvectors that transform as the 1D irreps, and
therefore pick up a U(1) phase under the action of the operator.

Since hyperbolic translation operators, represented in Eq. (18) do not
commute, the translational group will, in general, have higher dimensional
irreps. This implies that the U(1) Bloch theorem cannot be applied to
arbitrary eigenvectors for Hamiltonians with hyperbolic translational
symmetry. For Euclidean geometries, the group of lattice translations forms a
normal subgroup for the $d$-dimensional manifold ($\mathbb{R}^{d}$) whose
quotient group gives the Brillouin zone ($\mathbb{T}^{d}$) – a $d$-dimensional
torus.

On the other hand, for a hyperbolic surface, the quotient group of the
hyperbolic manifold under the non-Abelian Fuchsian translational symmetry
group gives a higher genus Riemann surface. This leads to a higher genus
representation for the hyperbolic Brillouin zone, which emerges as the
Jacobian of the Riemann surface Maciejko and Rayan (2021).

Translational symmetry operators for hyperbolic lattices are represented as:

$\begin{pmatrix}\alpha&\beta\\\
\beta^{*}&\alpha^{*}\end{pmatrix}z:=\frac{\alpha
z+\beta}{\beta^{*}z+\alpha^{*}},$ (20)

with $|\alpha|^{2}-|\beta|^{2}=1$.

The expressions for the nearest neighbour translational group operators have
been discussed in Ref. Chen  _et al._ (2023). For our choice of lattice
orientation, they take the form

$\gamma_{n}=\frac{1}{\sqrt{1-\sigma^{2}}}\begin{pmatrix}1&\sigma
e^{2\pi(n-1/2)i/p}\\\ \sigma e^{-2\pi(n-1/2)i/p}&1\end{pmatrix},$ (21)

where $\sigma=\sqrt{\frac{\cos{2\pi/p}+\cos{2\pi/q}}{1+\cos{2\pi/q}}}$ and
$n=1,2,...,5$. These are not all independent since
$\gamma_{5}=-\gamma_{1}^{-1}\gamma_{2}\gamma_{3}^{-1}\gamma_{4}$ Chen  _et
al._ (2023). Thus, we get four independent generators for the discrete
translational symmetry group of the {10,5} tessellation, with the action of
each translation $\gamma_{n}$ associated with a phase addition $e^{i\phi_{n}}$
to the wave function. Ref. Chen  _et al._ (2023) found that the five nearest
neighbouring sites are identified by the action of
$\mathbb{I},\kappa_{1}=\gamma_{1}\gamma_{2}^{-1},\kappa_{2}=\gamma_{2}\gamma_{3}^{-1},\kappa_{3}=\gamma_{1}\gamma_{4}\gamma_{3}^{-1},\kappa_{4}=\gamma_{2}\gamma_{3}^{-1}\gamma_{4}\gamma_{3}^{-1}$.
We can use the independent operators
$\kappa_{1},\kappa_{2},\kappa_{3},\kappa_{4}$ as generators of the translation
group, with the associated crystal momenta $k_{n}$ such that the tight-binding
Hamiltonian (Eq.[2] in the manuscript) gives the Bloch Hamiltonian

$H(\mathbf{k})=\begin{pmatrix}M&t(1+e^{ik_{1}}+e^{ik_{2}}+e^{ik_{3}}+e^{ik_{4}})\\\
t(1+e^{-ik_{1}}+e^{-ik_{2}}+e^{-ik_{3}}+e^{-ik_{4}})&-M\end{pmatrix}.$ (22)

It is important to remember that since we are applying automorphic Bloch
theory, we attach a phase of $k_{n}$ upon the action of $\kappa_{n}$.
Therefore, the overall phase obtained in a translation will be independent of
the order in which the translations have been applied – it acts as if the
translations commute, which is true for 1D irreps. Note that the $\\{k_{n}\\}$
are related to the original $\\{\phi_{n}\\}$ as Chen  _et al._ (2023)

$\mathbf{k}=\begin{pmatrix}1&-1&0&0\\\ 0&1&-1&0\\\ 1&0&-1&1\\\
0&1&-2&1\end{pmatrix}\mathbf{\phi}.$ (23)

## References

  * Friedmann (1979) A. Friedmann, in  _A Source Book in Astronomy and Astrophysics, 1900–1975_ (Harvard University Press, 1979) pp. 838–843.
  * Brooks (1996) R. Brooks, in  _Proc. Ashkelon Workshop on Complex Function Theory_ (Citeseer, 1996) pp. 61–73.
  * Monthus and Texier (1996) C. Monthus and C. Texier, Journal of Physics A: Mathematical and General 29, 2399 (1996).
  * Krioukov  _et al._ (2010) D. Krioukov, F. Papadopoulos, M. Kitsak, A. Vahdat, and M. Boguná, Physical Review E 82, 036106 (2010).
  * Breuckmann and Terhal (2016) N. P. Breuckmann and B. M. Terhal, IEEE transactions on Information Theory 62, 3731 (2016).
  * Breuckmann  _et al._ (2017) N. P. Breuckmann, C. Vuillot, E. Campbell, A. Krishna, and B. M. Terhal, Quantum Science and Technology 2, 035007 (2017).
  * Mertens and Moore (2017) S. Mertens and C. Moore, Physical Review E 96, 042116 (2017).
  * Kollár  _et al._ (2019) A. J. Kollár, M. Fitzpatrick, and A. A. Houck, Nature 571, 45 (2019).
  * Chen  _et al._ (2023) A. Chen, H. Brand, T. Helbig, T. Hofmann, S. Imhof, A. Fritzsche, T. Kießling, A. Stegmaier, L. K. Upreti, T. Neupert, _et al._ , Nature Communications 14, 622 (2023).
  * Maciejko and Rayan (2021) J. Maciejko and S. Rayan, Science advances 7, eabe9170 (2021).
  * Urwyler  _et al._ (2022) D. M. Urwyler, P. M. Lenggenhager, I. Boettcher, R. Thomale, T. Neupert, and T. Bzdušek, Physical Review Letters 129, 246402 (2022).
  * Liu  _et al._ (2022) Z.-R. Liu, C.-B. Hua, T. Peng, and B. Zhou, Physical Review B 105, 245301 (2022).
  * Liu  _et al._ (2023) Z.-R. Liu, C.-B. Hua, T. Peng, R. Chen, and B. Zhou, Physical Review B 107, 125302 (2023).
  * Tao and Xu (2023) Y.-L. Tao and Y. Xu, Physical Review B 107, 184201 (2023).
  * Mosseri  _et al._ (2022) R. Mosseri, R. Vogeler, and J. Vidal, Phys. Rev. B 106, 155120 (2022).
  * Stegmaier  _et al._ (2022) A. Stegmaier, L. K. Upreti, R. Thomale, and I. Boettcher, Physical Review Letters 128, 166402 (2022).
  * Maciejko and Rayan (2022) J. Maciejko and S. Rayan, Proceedings of the National Academy of Sciences 119, e2116869119 (2022).
  * Ashida  _et al._ (2020) Y. Ashida, Z. Gong, and M. Ueda, Advances in Physics 69, 249 (2020).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * Bergholtz  _et al._ (2021) E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Reviews of Modern Physics 93, 015005 (2021).
  * Ghatak and Das (2019) A. Ghatak and T. Das, Journal of Physics: Condensed Matter 31, 263001 (2019).
  * Banerjee  _et al._ (2023) A. Banerjee, R. Sarkar, S. Dey, and A. Narayan, Journal of Physics: Condensed Matter (2023).
  * Bender and Boettcher (1998) C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender (2007) C. M. Bender, Reports on Progress in Physics 70, 947 (2007).
  * Heiss (2012) W. Heiss, Journal of Physics A: Mathematical and Theoretical 45, 444016 (2012).
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, eaar7709 (2019).
  * Özdemir  _et al._ (2019) Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nature materials 18, 783 (2019).
  * Hodaei  _et al._ (2017) H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Nature 548, 187 (2017).
  * Chen  _et al._ (2017) W. Chen, Ş. Kaya Özdemir, G. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * Wiersig (2020) J. Wiersig, Photonics Research 8, 1457 (2020).
  * Wiersig (2014) J. Wiersig, Phys. Rev. Lett. 112, 203901 (2014).
  * De Carlo  _et al._ (2022) M. De Carlo, F. De Leonardis, R. A. Soref, L. Colatorti, and V. M. Passaro, Sensors 22, 3977 (2022).
  * Yao and Wang (2018) S. Yao and Z. Wang, Physical review letters 121, 086803 (2018).
  * Alvarez  _et al._ (2018) V. M. Alvarez, J. B. Vargas, and L. F. Torres, Physical Review B 97, 121401 (2018).
  * Yokomizo and Murakami (2019) K. Yokomizo and S. Murakami, Physical review letters 123, 066404 (2019).
  * Borgnia  _et al._ (2020) D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Physical review letters 124, 056802 (2020).
  * Zhang  _et al._ (2022) X. Zhang, T. Zhang, M.-H. Lu, and Y.-F. Chen, Advances in Physics: X 7, 2109431 (2022).
  * Lin  _et al._ (2023) R. Lin, T. Tai, L. Li, and C. H. Lee, arXiv preprint arXiv:2302.03057 (2023).
  * Boettcher  _et al._ (2022) I. Boettcher, A. V. Gorshkov, A. J. Kollár, J. Maciejko, S. Rayan, and R. Thomale, Physical Review B 105, 125118 (2022).
  * Jaiswal  _et al._ (2023) R. Jaiswal, A. Banerjee, and A. Narayan, New Journal of Physics 25, 033014 (2023).
  * Yuncken (2011) R. Yuncken, arXiv preprint arXiv:1103.2051 (2011).
  * Sausset and Tarjus (2007) F. Sausset and G. Tarjus, J. Phys. A: Math. Theor. 40, 12873 (2007).
  * Brody (2013) D. C. Brody, Journal of Physics A: Mathematical and Theoretical 47, 035305 (2013).
  * Ding  _et al._ (2022) K. Ding, C. Fang, and G. Ma, Nature Reviews Physics 4, 745 (2022).
  * Bulgakov  _et al._ (2006) E. N. Bulgakov, I. Rotter, and A. F. Sadreev, Phys. Rev. E 74, 056204 (2006).
  * Eleuch and Rotter (2016) H. Eleuch and I. Rotter, Physical Review A 93, 042116 (2016).
  * Tang  _et al._ (2020) W. Tang, X. Jiang, K. Ding, Y.-X. Xiao, Z.-Q. Zhang, C. T. Chan, and G. Ma, Science 370, 1077 (2020).
  * Budich and Bergholtz (2020) J. C. Budich and E. J. Bergholtz, Phys. Rev. Lett. 125, 180403 (2020).
  * Wiersig (2023) J. Wiersig, Phys. Rev. Res. 5, 033042 (2023).
  * Mosseri and Vidal (2023) R. Mosseri and J. Vidal, arXiv preprint arXiv:2304.02382 (2023).
  * Lenggenhager  _et al._ (2023) P. M. Lenggenhager, J. Maciejko, and T. Bzdušek, arXiv preprint arXiv:2305.04945 (2023).
  * Lux and Prodan (2023a) F. R. Lux and E. Prodan, in  _Annales Henri Poincaré_ (Springer, 2023) pp. 1–40.
  * Lux and Prodan (2023b) F. R. Lux and E. Prodan, Phys. Rev. Lett. 131, 176603 (2023b).
  * Okuma  _et al._ (2020) N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Physical review letters 124, 086801 (2020).
  * Sun  _et al._ (2023) J. Sun, C.-A. Li, S. Feng, and H. Guo, arXiv preprint arXiv:2305.19810 (2023).