# Topological Effects in Two-Dimensional Quantum Emitter Systems

Miguel Bello<EMAIL_ADDRESS>Max-Planck-Institut für Quantenoptik,
Hans-Kopfermann-Straße 1, D-85748 Garching, Germany  J. I. Cirac  Max-Planck-
Institut für Quantenoptik, Hans-Kopfermann-Straße 1, D-85748 Garching, Germany
Munich Center for Quantum Science and Technology, Schellingstraße 4, 80799
München, Germany

###### Abstract

In this work we discuss particular effects that take place in systems of
quantum emitters coupled to two-dimensional bosonic topological insulators.
For a single emitter coupled to the Haldane model, we find a “fragile”
quasibound state that makes the emitter dynamics very sensitive to the model’s
parameters, and gives rise to effective long-range interactions that break
time-reversal symmetry. We then discuss one-dimensional arrangements of
emitters, emitter line defects, and how the topology of the bath affects the
effective polariton models that appear in the weak-coupling regime when the
emitters are spectrally tuned to a bandgap. In the Harper-Hofstadter model we
link the non-monotonic character of the effective interactions to the Chern
numbers of the surrounding energy bands, while in the Haldane model we show
that the effective models are either gapless or not depending on the topology
of the bath. Last, we discuss how the presence of emitters forming an ordered
array, an emitter superlattice, can produce polariton models with non-trivial
Chern numbers, and also modify the topology of the photonic states in the
bath.

##  I Introduction

Quantum emitters coupled to structured baths display a number of phenomena
that are not present in more conventional baths: bound states, non-Markovian
dynamics, etc. [1, 2, 3] In particular, topological baths offer the
possibility to go beyond the standard quantum electrodynamics of atoms in
vacuum, and explore the interplay between exotic states of light and matter.
Several theoretical studies [4, 5, 6, 7, 8, 9, 10, 11, 12] and a few
experimental works [13, 14] have already demonstrated interesting effects,
such as the generation of exotic coherent dipole-dipole interactions [4, 14,
10], or coupling to chiral edge states which are topologically protected [13,
8, 9]. These works trigger further questions, as there are many different
kinds of topological models that remain unexplored, and many different
topological concepts whose implications for quantum optics are not yet well
understood.

In one-dimensional (1D) topological photonic baths, the relation between the
underlying bath topology and the quantum optical effects that can be observed
is fairly well understood. For example, it has been shown that the bath
topology plays a role in the bound states that the emitters seed. For specific
values of the emitter transition frequency, the atoms can behave like a mirror
or, in other words, as if they produce a vacancy in the bath. Consequently,
the bound states can share many properties of the topological edge states that
the bath supports in the nontrivial phase [4, 7]. These vacancy-like bound
states were later recognized as a general feature that occurs in many photonic
lattices, whether topological or not [5].

In two-dimensional (2D) topological baths, the effects of a nontrivial
topology are much less clear. A straightforward way to observe topological
effects is to couple emitters to the boundaries of 2D topological baths. In
this case, the dynamics are sensitive to whether or not the bath supports edge
states, and this can be exploited for the development of interesting
applications [8, 9]. Another possibility is to couple emitters to the bulk of
a 2D topological bath, however, in this case it is hard to disentangle what
effects are simply due to the breaking of time-reversal symmetry, and what
effects are due to a nontrivial topology. This kind of systems have been
explored in Ref. [6] for a few emitters coupled to the Harper-Hofstadter model
in the Landau regime, and in Ref. [5] for a single emitter coupled to the
Haldane model. In the latter, it was shown that vacancy-like bound states
exist only when the bath is in a nontrivial topological phase, albeit only for
specific values of the parameters within these phases, where the model
recovers particle-hole symmetry. It is therefore an open question to what
extent the topology of the underlying bath can affect the dynamics of quantum
emitters coupled to it.

In this paper we aim to shed light on this issue. We consider a 2D lattice of
bosonic degrees of freedom that is a topological insulator, and a set of
emitters coupled to it. We investigate the dissipative and coherent dynamics
of emitters, paying especial emphasis on the way that topological properties
are reflected in the properties of polaritons (excitations trapped in atom-
photon bound states). We consider three different scenarios: i) few emitters;
ii) a 1D array of emitters coupled along one of the columns (or rows) of the
2D lattice; iii) 2D arrays, where each emitter is coupled to a single bath
mode in a periodic fashion. In case i), the physics is similar to that taking
place in other topologically trivial lattices, albeit with some differences,
such as the appearance of quasi-bound states around singular gaps which are
not vacancy-like, and are therefore, in a sense, fragile, as they require
finely tuned system parameters to exist. In case ii), the bath’s topology has
more dramatic effects. The effective polariton models have dispersion
relations whose complexity, or in which the presence of a gap, is linked to
the Chern numbers of the bath’s energy bands. This has clear observable
consequences, such as propagation with multiple wavefronts, or the chiral
dynamics of a single polariton in a 1D emitter array. In case iii), we show
how the polariton models can have different Chern numbers depending on the
array geometry and the light-matter coupling strength.

This paper is organized as follows. In Section II we describe the specific
models studied, and in Section III we recall the concept of vacancy-like bound
states, which is essential to understand the following sections. Then, in
Sections IV to VI we describe the phenomena that take place in different kinds
of emitter arrangements, in order of increasing complexity. We conclude in
Section VII by summarizing the main results and pointing out possible future
research directions.

##  II Models

We consider a set of $N$ identical two-level systems (quantum emitters), with
raising and lowering operators
$\\{\sigma^{+}_{n}\equiv\ket{e_{n}}\\!\bra{g_{n}},\sigma^{-}_{n}\equiv(\sigma^{+}_{n})^{\dagger}\\}_{n=1}^{N}$;
$\ket{g_{n}}$ and $\ket{e_{n}}$ denote the ground and excited state of the
$n$th emitter, respectively. Each one is coupled to a single local bosonic
mode of a 2D lattice of size $L_{x}\times L_{y}$ in units of the lattice
constant, which is modelled via a tight-binding Hamiltonian $H_{B}$. The
coupling, described by $H_{I}$, is of the Jaynes-Cummings type, with coupling
constant $g$. The total Hamiltonian of the system is $H=H_{E}+H_{B}+H_{I}$,
with

$\displaystyle H_{E}$
$\displaystyle=\Delta\sum_{n}\sigma^{+}_{n}\sigma^{-}_{n}\,,$ (1)
$\displaystyle H_{B}$
$\displaystyle=\sum_{\bm{k}}\Psi^{\dagger}_{\bm{k}}\mathcal{H}_{B}(\bm{k})\Psi_{\bm{k}}\,,$
(2) $\displaystyle H_{I}$
$\displaystyle=\frac{g}{\sqrt{L_{x}L_{y}}}\sum_{n,\,\bm{k}}e^{i\bm{k}\bm{r}_{n}}\sigma^{+}_{n}a_{\alpha_{n},\bm{k}}+\mathrm{H.c.}$
(3)

$\Psi_{\bm{k}}\equiv(a_{1,\bm{k}},\dots,a_{M,\bm{k}})^{T}$ is a column vector
of bosonic annihilation operators, $a_{\alpha,\bm{k}}$ annihilates a boson
with quasimomentum $\bm{k}$ in the $\alpha$ sublattice;
$\bm{r}_{n}\in\mathbb{Z}_{L}^{2}$ and $\alpha_{n}\in\mathbb{Z}_{M}$ denote the
unit cell and sublattice to which the $n$th emitter is coupled, respectively.

Figure 1: Schematic of the system under consideration. The bath can be split
in two parts: a bath with vacancies, described by the Hamiltonian
$H_{B}^{\circ}$, and the sites to which the emitters are coupled, described by
$H_{B}^{\bullet}$.

For the bath, we will consider two paradigmatic examples of two-dimensional
Chern insulators:

_1\. Harper-Hofstadter model_ [15]. This model corresponds to a square lattice
with nearest-neighbor hoppings of strength $J$, threaded by a magnetic flux.
For a magnetic flux per plaquette $\phi=2\pi p/q$ ($p,q\in\mathbb{N}$,
coprime), the momentum-resolved bath Hamiltonian $\mathcal{H}_{B}(\bm{k})$,
$\bm{k}=(k_{x},k_{y})$, can be represented by a $q\times q$ matrix, with
elements

$\frac{1}{J}\left[\mathcal{H}_{B}(\bm{k})\right]_{mn}=-2\cos(k_{x}+n\phi)\delta_{m,n}-\delta_{m,n+1}\\\
-\delta_{m,n-1}-e^{-ik_{y}}\delta_{n,m+q-1}-e^{ik_{y}}\delta_{m,n+q-1}\,.$ (4)

_2\. Haldane honeycomb model_ [16]. This model corresponds to a honeycomb
lattice with nearest-neighbor hoppings of strength $t_{1}$, and complex next-
to-nearest-neighbor hoppings of strength $t_{2}e^{\pm i\phi}$. The complex
phases are chosen such that there is no net magnetic flux threading the
lattice, yet the model is not time-reversal invariant. In addition, there is a
different on-site potential $\pm M$ for each sublattice. The momentum-resolved
bath Hamiltonian can be expressed in terms of the Pauli matrices
$\\{\sigma^{\nu}\\}_{\nu=x,y,z}$, and the identity matrix $\sigma^{0}$, as
$\mathcal{H}_{B}(\bm{k})=h_{0}(\bm{k})\sigma^{0}+\sum_{\nu}h_{\nu}(\bm{k})\sigma^{\nu}$,
with

$\begin{split}h_{0}(\bm{k})&=-2t_{2}\operatorname{c}(\phi)\left[\operatorname{c}(k_{x})+\operatorname{c}(k_{y})+\operatorname{c}(k_{y}-k_{x})\right]\,,\\\
h_{x}(\bm{k})&=-t_{1}\left[1+\operatorname{c}(k_{x})+\operatorname{c}(k_{y})\right]\,,\\\
h_{y}(\bm{k})&=-t_{1}\left[\operatorname{s}(k_{x})+\operatorname{s}(k_{y})\right]\,,\\\
h_{z}(\bm{k})&=M+2t_{2}\operatorname{s}(\phi)\left[\operatorname{s}(k_{x})-\operatorname{s}(k_{y})+\operatorname{s}(k_{y}-k_{x})\right]\,,\end{split}$
(5)

where we have used the abbreviations $\operatorname{s}(x)\equiv\sin(x)$ and
$\operatorname{c}(x)\equiv\cos(x)$.

##  III Vacancy-like bound states

Since many of the effects discussed in this paper are due to the presence of
vacany-like bound states (VBS) [5], in this section we briefly review their
properties and translate them into properties of the emitter self-energy (also
known as the level-shift operator [17]). Furthermore, we show their relevance
in understanding the whole system’s spectrum in the large-$g$ limit.

To define them, it is useful to split the bath Hamiltonian in two parts,
$H_{B}=H_{B}^{\circ}+H_{B}^{\bullet}$, such that $H_{B}^{\bullet}$ contains
only the terms involving the sites to which the quantum emitters are coupled,
while $H_{B}^{\circ}$ contains the rest of the terms, see Fig. 1 for a
schematic representation.

######  Theorem 1 (Existence of VBS, Ref. [5])

if $\ket{\Psi_{\mathrm{VM}}}$ is a single-particle eigenstate of
$H^{\circ}_{B}$ with energy $E_{\mathrm{VM}}$, the state
$\ket{\Psi_{\mathrm{VBS}}}=\sum_{n}c_{e,n}\sigma^{+}_{n}\ket{\mathrm{vac}}+\alpha\ket{\Psi_{\mathrm{VM}}}$
is an eigenstate of the whole Hamiltonian $H=H_{E}+H_{B}+H_{I}$ with the same
energy $E_{\rm BS}=E_{\rm VM}$, provided $\Delta=E_{\mathrm{VM}}$.

We can identify vacancy modes with those energies $E_{\rm VM}$ for which the
self-energy of the emitter subsystem $\Sigma$, is singular, $\det\Sigma(E_{\rm
VM})=0$. Note that for a system in which $N$ emitters are coupled locally,
each to a single site of a $D$-dimensional photonic lattice, the self-energy
is a $N\times N$ matrix, whose elements are given by (see Appendix A)

$\Sigma_{mn}(z)=g^{2}\left[\Phi(z,\bm{r}_{m}-\bm{r}_{n})\right]_{\alpha_{m},\alpha_{n}}\,,$
(6)

with $\Phi$ an $M\times M$ matrix given by

$\Phi(z,\bm{r})=\int_{\mathrm{BZ}}\frac{d^{D}\bm{k}}{(2\pi)^{D}}\,e^{i\bm{k}\cdot\bm{r}}\left[z-\mathcal{H}_{B}(\bm{k})\right]^{-1}\,,$
(7)

where the integral is over the first Brillouin zone. In particular, for a
single emitter, the self-energy is a complex function, and the energies of the
vacancy modes are its zeros on the real line, $\Sigma(E_{\rm VM})=0$. Theorem
1 follows from the fact that the bound-state energies satisfy the “pole
equation” [18]

$\det[E_{\rm BS}-\Delta-\Sigma(E_{\rm BS})]=0\,.$ (8)

Furthermore, one can show that the emitter weights of the bound states satisfy
(Appendix A)

$\bm{c}_{e}^{\dagger}\left[1-\Sigma^{\prime}(E_{\mathrm{BS}})\right]\bm{c}_{e}=1\,,$
(9)

where $\Sigma^{\prime}(z)\equiv\partial_{z}\Sigma(z)$, and
$\bm{c}_{e}\equiv(c_{e,1},\dots,c_{e,N})^{T}$.

Theorem 1 holds regardless the value of the light-matter coupling constant
$g$. As a consequence:

  1. (i)

The shape, in particular, the decay length, or the inverse participation
ration (IPR), of the photonic component of the bound state, does not change as
a function of $g$—it is fixed by $\ket{\Psi_{\mathrm{VM}}}$.

  2. (ii)

The populations of the emitter excited states $\lvert c_{e,n}\rvert^{2}$,
decay monotonically like $\sim g^{-2}$ as $g\to\infty$.

It should be noted that even if the condition $\Delta=E_{\mathrm{VM}}$ is not
met, the presence of vacancy modes is relevant in the strongly-interacting
regime $g\gg\Delta,\lVert H_{B}\rVert$. As shown in Fig. 2(a), the bound-state
energies may show two qualitatively different behaviors as the light-matter
coupling increases: either they diverge or they converge to the energy of a
vacancy mode. The former bound states, which we shall call _strongly-
hybridized polaritons_ , are, in this regime, approximate eigenstates of the
interaction Hamiltonian $H_{I}$. For our particular Hamiltonian, the
eigenstates of $H_{I}$ are the symmetric and antisymmetric superpositions
$\left(\sigma^{+}_{n}\pm
a_{\bm{r}_{n},\alpha_{n}}^{\dagger}\right)\ket{\mathrm{vac}}/\sqrt{2}$. The
latter kind is, in the same regime, approximately equal to the vacancy modes.
This explains the observed behavior of the emitter occupation $\lvert
c_{e}\rvert^{2}$, and the localization of the photonic component of the bound
states, shown in Fig. 2(b). Namely, for vacancy-like bound states, or bound
states that tend to a vacancy-mode, $\lvert c_{e}\rvert^{2}$ and the IPR show
the behavior described in points (i) and (ii), whereas for strongly-hybridiced
polaritons, $\lvert c_{e}\rvert^{2}\to 1/2$ and ${\rm IPR}\to 1$ as
$g\to\infty$.

Figure 2: (a) Bound-state energies for a single emitter coupled to a single
site of the Harper-Hofstadter model, with $\phi=2\pi/3$, as a function of the
light-matter coupling strength. Grey areas mark the range of the bath’s energy
bands, while the ticks on the right vertical axis mark the energies of the
vacancy modes; $\Delta=E_{\rm VM,2}$ (blue, solid line). (b) Inverse
participation ratio (IPR), and emitter’s excited state occupation of the same
bound states. The IPR has been computed as
$\mathrm{IPR}:=\sum_{j}\lvert\psi_{j}\rvert^{4}$, where the sum runs over all
bath sites, and $\psi_{j}$ denotes the probability amplitude to find the
photon in the $j$th bath site, normalized such that
$\sum_{j}\lvert\psi_{j}\rvert^{2}=1$; $1/(L_{x}L_{y})\leq\mathrm{IPR}\leq 1$.
The lowest value is attained for a homogeneous state (constant probability),
while the maximum is attained for a fully localized state (at a single site).
A lattice of $(L_{x},L_{y})=(3\cdot 2^{4},2^{4})$ cells along each dimension,
and periodic boundary conditions, has been used in this exact-diagonalization
calculation.

##  IV Few emitters

Although the connection to the underlying bath topology in systems with few
emitters is less straightforward than in the case of extended emitter
arrangements, it is nevertheless interesting to look at the effects induced by
the breaking of time-reversal symmetry and point out similarities and
differences with other 2D baths. As we show in this section, for the Harper-
Hofstadter model, the major difference as compared with the square lattice
model is the non-monotonicity of the bound states that appear in the gaps
opened in the spectrum by the presence of a magnetic flux. For the Haldane
model, similarly to the honeycomb lattice, a single emitter can seed a
quasibound state that gives rise to true long-range interactions, which also
break time-reversal symmetry at short distances.

To show all this, it is instrumental to compute the emitter self-energy.
Knowing $\Sigma_{mn}(\Delta)$ for arbitrary distances
$\bm{r}_{mn}\equiv\bm{r}_{m}-\bm{r}_{n}$ between the emitters and values of
$\Delta$ outside the band ranges is equivalent to knowing the shape of the
single-emitter bound states (see Appendix A), and thus also provides
information about the effective coherent interactions mediated by the bath.
Unfortunately, obtaining analytical expressions for the self-energy is rather
difficult for both the Harper-Hofstadter and Haldane models. A workaround is
to integrate numerically its expression in terms of the momentum-resolved bath
Hamiltonian, a task that can be facilitated by using as a starting point the
exact expressions for the partial integrals shown in Appendix C.

###  IV.1 Harper-Hofstadter model

The Harper-Hofstadter model has been partially covered in Ref. [6], where they
consider quantum emitters coupled to a square lattice threaded by a relatively
small magnetic flux, one such that the magnetic length
$l_{B}=l_{0}/\sqrt{\phi}$ is larger than the lattice constant $l_{0}$, but
smaller than the lattice size, $1<l_{B}/l_{0}<\min(L_{x},L_{y})$. In this
limit the bath spectrum splits into many sub-bands with a vanishingly small
width (Landau levels); a single emitter resonant with one of the Landau levels
does not decay, but its excited-state population oscillates as the photon hops
back and forth between the emitter and the Landau orbital located around it.

Here, instead, we consider the limit of large magnetic flux $l_{B}/l_{0}<1$,
in which the bath spectrum consists only of a few sub-bands with a finite
width. In this case, the emission dynamics is very similar to that found for
the square lattice [1, 2], provided the coupling strength is sufficiently
small compared to the width of the sub-bands. The density of states for each
sub-band has a similar shape, with a van-Hove singularity giving rise to the
same long-term algebraic decay and emission profile peaked along the main
diagonals of the lattice (see Fig. 3). Regarding the bound states, the main
difference with respect to the square lattice is that they do not decay
monotonically away from the emitter. This has important consequences for
extended emitter arrays that we will review in the next section.

Figure 3: Self-energy (a), $\Sigma(\Delta+i0^{+})=\delta\omega-i\Gamma/2$, and
emission dynamics (b), for a single emitter coupled to the Harper-Hofstadter
model with $\phi=2\pi/3$. In (b) we show the emitter excited state occupation
as a function of time for two different values of the emitter transition
frequency $\Delta_{j}$ ($j=1,2$), tuned to different van-Hove singularities of
the bath’s density of states. Next to it we show respective snapshots of the
bath occupation at time $t_{0}=50J^{-1}$. The coupling strength is in both
cases $g=0.1J$.

###  IV.2 Haldane model

Haldane’s model has been less explored in the quantum-emitter literature.
Here, we focus on the regime where both $\lvert M\rvert,\lvert
t_{2}\rvert\ll\lvert t_{1}\rvert$, in which the Haldane model can be regarded
as a small perturbation to the honeycomb lattice [19]. We shall call this
regime the “honeycomb limit”. In the following, we derive approximate
expressions for the self-energy, valid in the honeycomb limit for values of
$\Delta$ outside the band ranges, and analyze their implications for the
effective bath-induced interactions between the emitters.

Particularizing Eqs. (6) and (7) to the case of a 2D, two-band bath, we can
express it in integral form as

$\Sigma_{mn}(z)=g^{2}\int_{\rm
BZ}\frac{d^{2}\bm{k}}{(2\pi)^{2}}\frac{e^{i\bm{k}\cdot\bm{r}_{mn}}f_{\alpha_{m}\alpha_{n}}(z,\bm{k})}{\left[z-h_{0}(\bm{k})\right]^{2}-h^{2}(\bm{k})}\,.$
(10)

Here, $h^{2}(\bm{k})\equiv
h_{x}^{2}(\bm{k})+h_{y}^{2}(\bm{k})+h_{z}^{2}(\bm{k})$, and the numerator can
take one of the following forms:

$\displaystyle f_{AA}(z,\bm{k})$
$\displaystyle=z-h_{0}(\bm{k})+h_{z}(\bm{k})\,,$ (11) $\displaystyle
f_{BB}(z,\bm{k})$ $\displaystyle=z-h_{0}(\bm{k})-h_{z}(\bm{k})\,,$ (12)
$\displaystyle f_{AB}(z,\bm{k})$
$\displaystyle=h_{x}(\bm{k})-ih_{y}(\bm{k})\,,$ (13) $\displaystyle
f_{BA}(z,\bm{k})$ $\displaystyle=h_{x}(\bm{k})+ih_{y}(\bm{k})\,,$ (14)

depending on the sublattice to which the emitters are coupled. The largest
contribution to the integral appearing in Eq. (10) comes from the bath modes
whose frequencies are closest to that of the emitters. It is possible to
obtain an analytical approximation expanding the integrand around the
quasimomenta of such modes, if the distance between the emitters is large
$\lvert\bm{r}_{mn}\rvert\gg 1$ (see Appendix B). When $\Delta$ lies above or
below the bath’s bands and close to the band edges we obtain

$\Sigma_{mn}\propto\operatorname{K}_{0}\left(\frac{r^{\prime}_{mn}}{\xi_{0}}\right)\,.$
(15)

Here, $r^{\prime}_{mn}$ is the Euclidean norm of $\bm{r}^{\prime}_{mn}$, which
is simply $\bm{r}_{mn}$ expressed in a different coordinate basis, and
$\operatorname{K}_{0}$ denotes the zeroth-order Bessel function of the second
kind. If instead, $\Delta$ lies in the gap between the two bands,

$\Sigma_{mn}(\Delta)\simeq I_{+}(\Delta)+I_{-}(\Delta)\,,$ (16)

where each contribution comes from a different K-point,
$\pm\bm{K}\equiv\pm(1,-1)2\pi/3$. For two emitters coupled to the same
sublattice,

$I^{\pm}_{mn}(\Delta)\propto e^{\pm
i\bm{K}\cdot\bm{r}_{mn}}\operatorname{K}_{0}\left(\frac{r^{\prime}_{mn}}{\xi_{\pm}}\right)\,,$
(17)

while for two emitters coupled to different sublattices, say
$(\alpha_{m},\alpha_{n})=(A,B)$,

$I^{\pm}_{mn}(\Delta)\propto\frac{e^{\pm
i(\bm{K}\cdot\bm{r}_{mn}-\theta_{mn}-\pi/2)}}{\xi_{\pm}}\operatorname{K}_{1}\left(\frac{r^{\prime}_{mn}}{\xi_{\pm}}\right)\,.$
(18)

Here, $\theta_{mn}$ denotes the angle between $\bm{r}^{\prime}_{mn}$ and the
horizontal axis, and $\operatorname{K}_{1}$ denotes the first-order modified
Bessel function of the second kind.

In light of these expressions we conclude that for emitters tuned above and
below the bands the effective interactions mediated by the bath are
approximately real. In the middle band gap, $I^{+}_{mn}(\Delta)$ and
$I^{-}_{mn}(\Delta)$ may have different magnitude, resulting in complex
effective interactions. This is a consequence of breaking the symmetry between
the two different K-points, which, by contrast, is always present in the
honeycomb lattice. Individual phases do not determine by themselves whether
the polariton dynamics is time-reversal invariant or not. Instead, a gauge-
invariant quantity quantifying this is the total phase accumulated when
hopping from one emitter to another in a closed path, or, in other words, the
flux threading the path. In Fig. 4, we present this calculation for three
emitters forming a triangle, for different values of the detuning: above the
bath’s bands, and in the middle band gap. As it turns out, the flux tends to
zero as the distance between the emitters increases. Only when the symmetry
between the two K-points is restored ($M=0$), the interaction ranges are equal
$\xi_{+}=\xi_{-}$, and the flux does not decay, but it oscillates between
$\pm\pi/2$, which corresponds to perfectly chiral dynamics [6].

Figure 4: Strength of the effective interactions (right column) and total
magnetic flux (left column) threading a triangular emitter plaquette, as shown
in the schematics above, as a function of the distance $\ell$ between the
emitters. The dots have been obtained by numerical integration of Eq. (6),
while the crosses correspond to the approximate analytic formulas (detailed in
Appendix B). The remaining parameters are $t_{2}=0.1t_{1}$ and $\phi=0.1$, the
same in all the cases shown.

For a critical bath, i.e., when the gap closes at either $\pm\bm{K}$, the
respective interaction range diverges, ${\xi_{\pm}\to\infty}$. From Eq. (18),
using the asymptotic relation $\operatorname{K}_{1}(x)\sim 1/x$ for $x\to 0$,
we realize that the effective interaction between two emitters coupled to
different sublattices decays as $\sim 1/r^{\prime}_{mn}$ when they are
spectrally tuned to the vanishing gap, so it is truly long-range. The same
thing happens in the honeycomb lattice, which presents two Dirac cones at both
$\pm\bm{K}$, and has therefore a singular gap. When the emitter transition is
tuned to the frequency of the singular gap a quasibound state appears, which
scales as $\sim 1/r$, where $r$ is the distance to the emitter, giving rise to
long-range effective interactions [19]. This quasibound state is, in the case
of the honeycomb lattice, actually a vacancy-like quasibound state, meaning
that its energy and photonic component remain the same, regardless the value
of $g$. In the Haldane model, the origin of the long-range interactions can be
attributed also to a quasibound state, although in this case it is not of the
vacancy-like type. In the Haldane model the gap does not close simultaneously
at both K-points, and as a consequence, the self-energy for a single emitter
does not vanish at the frequency of the single Dirac point; while the density
of states indeed vanishes, the Lamb shift is different from zero, so there is
no longer a vacancy mode at that frequency. Instead, there will be a
quasibound state whenever the detuning $\Delta$ compensates precisely for the
Lamb shift. In such a case, an emitter initially excited will not decay
completely, as shown in Fig. 5, even though its frequency is resonant with
some guided modes in the bath. This quasibound state is fragile, in the sense
that it only exists for particular values of $(\Delta,g)$. From the
approximate expressions for the self-energy, we can see that at long distances
it also scales like $\sim 1/r$, and has support mostly on the opposite
sublattice to which the emitter is coupled.

Figure 5: Emitter excited state occupation as a function of time for three
different configurations, showcasing the sensitivity of the dynamics to the
emitter detuning and light-matter coupling strength. The bath parameters are
the same in all cases and correspond to a critical bath with a Dirac point at
$-\bm{K}$: $t_{2}=0.1t_{1}$, $\phi=0.5$, and $M=3\sqrt{3}t_{2}\sin(\phi)$. The
Dirac-point frequency is denoted by $\omega_{0}$, while $\delta\omega_{0}$
denotes the Lamb shift at $\omega_{0}$ for $g=0.5t_{1}$. The simulation has
been carried out in a lattice with $L_{x}=L_{y}=2^{8}$ unit cells along each
dimension and periodic boundary conditions. The right-hand side plot shows the
same data in logarithmic scale.

##  V 1D emitter arrangements

The main characteristic of 2D Chern insulators is the presence of chiral edge
states in lattices with open boundary conditions (OBC) [20]. Equivalently, we
can consider lattices with vacancies forming a 1D (line) defect. Since the
energies of the chiral edge states span the whole range of a bandgap, we are
guaranteed to find vacancy-like bound states for arrangements of emitters
forming a line defect, provided the emitter transition frequency lies in one
of the topological bandgaps. Furthermore, since the number of edge states is
linked to the Chern numbers of the bands above and below the gap, the number
of vacancy-like bound states that can be found is also determined by the Chern
numbers in general. In the following, we demonstrate this effect in the
Harper-Hofstadter model and the Haldane model.

###  V.1 Harper-Hofstadter model

Let us consider first a set of quantum emitters coupled to consecutive sites
of the Harper-Hofstadter model along the $x$-direction, i.e., their
coordinates in the basis of the lattice vectors are
$\bm{r}_{n}\equiv(x_{n},y_{n})=(n,0)$, see Fig. 6(a) for a schematic
representation. The presence of emitters breaks the bath’s translational
symmetry in the direction perpendicular to the line defect, but we can still
Fourier-transform in the parallel direction, obtaining a collection of 1D
single-emitter systems parametrized by the quasimomentum $k_{x}$,
$H=\sum_{k_{x}}H^{\mathrm{1D}}(k_{x})$, with

$H^{\mathrm{1D}}(k_{x})=\Delta\sigma^{+}_{k_{x}}\sigma^{-}_{k_{x}}+g(\sigma^{+}_{k_{x}}a_{k_{x},0}+\mathrm{H.c.})+H_{B}^{\mathrm{1D}}(k_{x})\,.$
(19)

In this case, $H_{B}^{\mathrm{1D}}$ corresponds to a chain with nearest-
neighbor hopping of amplitude $J$, and periodic on-site potentials
$\mu_{n}=-2J\cos(k_{x}+n\phi)$ (see Appendix C.1). Due to their proximity, the
bound states of neighboring emitters will hybridize forming a band. Its
dispersion relation $\omega_{\mathrm{BS}}(k_{x})$, can be obtained solving the
corresponding pole equation for each 1D system,

$\omega_{\mathrm{BS}}-\Delta-\Sigma(\omega_{\mathrm{BS}};k_{x})=0\,.$ (20)

Here, $\Sigma(z;k_{x})$ is the single-emitter self-energy corresponding to the
1D system described by Eq. (19), which can be computed analytically (see
Appendix C.1).

Figure 6: (a) Schematic representation of an emitter line defect in the
Harper-Hofstadter model. The dots represent the bath sites, while lines
connecting them represent hopping processes (marked with an arrow if they are
complex). Blue crosses mark the sites to which the quantum emitters are
coupled. The big, red and blue arrows represent chiral edge states that are
present in the model with vacancies. (b) Spectrum of $H_{B}^{\circ}$ for the
same system, with $\phi=2\pi/7$. The spectrum is color-coded according to the
localization of the eigenstates with respect to the line defect. On it, we
show the polariton bands (orange and green lines) obtained for two different
values of $\Delta$ (dashed horizontal lines). Vacancy-like bound states
correspond to the crossings between the edge-state energies and $\Delta$. The
light-matter coupling constant is $g=0.5J$. (c) Photonic component of the
bound state with energy closest to $\Delta$ for a finite emitter line defect.
The parameters are the same as in panel (b), with $\Delta=-1.3J$. In addition
to the probability of finding the photon in a given bath site (color plot), we
also show the probability current [21].

According to Theorem 1, whenever $\Delta$ matches the energy of an edge state
(for a given $k_{x}$), there will be a vacancy-like bound state. Thus, the
number of vacancy-like bound states for any $\Delta$ inside the $n$th gap is,
in general, equal to the number of edge modes crossing that gap, see Fig.
6(b). In other words, for $\max\omega_{n-1}<\Delta<\min\omega_{n}$, the number
of zeros of $\Sigma(\Delta;k_{x})$ as a function of $k_{x}$, $\zeta(\Delta)$,
is lower bounded by $\zeta(\Delta)\geq 2(c_{n}-c_{n-1})$. Here $\omega_{n}$
and $\omega_{n-1}$ denote the dispersion relations of the bath bands above and
below the $n$th gap, respectively, and $c_{n}$ and $c_{n-1}$ denote their
corresponding Chern numbers. The factor 2 is due to the fact that we have to
take into account the edge states on both sides of the line defect, which have
opposite chirality. The only cases in which this bound is not satisfied is
when $\Delta$ is tuned precisely to a crossing between edge modes.

For an infinitely long emitter line defect, or for a system with periodic
boundary conditions (PBC), the photonic component of these vacancy-like bound
states is actually a chiral edge state of the bath. And therefore, it has
support only in one of the two half-planes into which the line defect splits
the bath. However, for sufficiently large but finite line defects, the true
bound states are linear combinations of the (degenerate) bound states found
for PBC, combined in such a way that they meet the OBC. For the Harper-
Hofstadter model, degenerate vacancy-like bound states happen to have exactly
the same decay length, but opposite chirality and directionality, so the
vacancy-like bound state ends up decaying symmetrically away from the emitter
line defect, see Fig. 6(c).

As we can appreciate in Fig. 6(b), not only the number of zeros of
$\Sigma(\Delta,k_{x})$ depends on the topology, but the overall “bumpiness” of
the bound-state dispersion relation increases for increasing number of edge
states. We can write an effective tight-binding model for bound states
localized around the emitter positions as
$H_{\mathrm{eff}}=\sum_{n,\,d}J^{\mathrm{eff}}_{d}\beta^{\dagger}_{n+d}\beta_{n}$,
where $\beta^{\dagger}_{n}$ and $\beta_{n}$ create and annihilate a bound
state (polariton) localized around the $n$th emitter, respectively. The
hopping constants can be obtained from the dispersion relation as
$J^{\mathrm{eff}}_{d}=\sum_{k}\omega_{\mathrm{BS}}(k)e^{ikd}/N$. Since $H^{\rm
1D}_{B}$ is symmetric with respect to the change $k_{x}-m\phi\to-k_{x}-m\phi$,
where $m\in\\{1,2,\dots,q\\}$ labels the sublattice to which the emitters are
coupled, the same kind of inversion symmetry is present in all the spectrum,
which, in turn, implies that the effective hopping constants are real up to a
phase, $\operatorname{Im}(e^{im\phi d}J^{\rm eff}_{d})=0$, which can be
absorbed in the definition of the polariton creation and annihilation
operators. The “bumpiness” of the dispersion relation translates into exotic
hopping constants, which, although decaying with an overall exponential
factor, may be larger for longer distances than for shorter distances, see
Fig. 7(a). This can be observed, e.g., by letting an initially localized
polariton, $\ket{\psi(0)}=\beta^{\dagger}_{0}\ket{\mathrm{vac}}$, evolve for a
short period of time. The resulting wave function, $\ket{\psi(\delta
t)}\simeq(1-i\delta tH_{\mathrm{eff}})\ket{\psi(0)}$, essentially maps the
hopping amplitudes to occupations of the neighboring sites. The peculiar shape
of the effective hoppings can be traced back to the particular shape of the
single-emitter bound states [22], which in the continuum limit ($\phi\to 0$)
have the nodes typical of Landau orbitals [6].

The particular shape of $\omega_{\mathrm{BS}}$ also has visible consequences
in the ballistic transport of a localized excitation. For a localized initial
state, the wave function develops several wavefronts, which propagate linearly
in time with a critical group velocity $v_{c}$, see Fig. 7(b). This can be
understood from the behavior of the probability
$p_{x}(t):=\bra{\psi(t)}\beta^{\dagger}_{x}\beta_{x}\ket{\psi(t)}$ at large
space-time scales. In the limit $x,t\to\infty$, keeping constant the ratio
$x/t=\xi$, it is possible to approximate it using the stationary-phase
approximation method [23],

$\begin{split}p_{\xi}(t)&=\left\lvert\int_{\rm
BZ}\frac{dk}{2\pi}\,e^{i[k\xi-\omega_{\rm BS}(k)]t}\right\rvert^{2}\\\
&\sim\left\lvert\sum_{k_{0}}\frac{e^{if(k_{0})t+\operatorname{sign}\left(f^{\prime\prime}(k_{0})\right)i\pi/4}}{\sqrt{2\pi
t\lvert
f^{\prime\prime}(k_{0})\rvert}}+O(t^{-1})\right\rvert^{2}\,.\end{split}$ (21)

Here, the sum runs over all stationary points of the phase
$f(k)=k\xi-\omega_{\rm BS}(k)$, $f^{\prime}(k_{0})=0$, within the integration
range. Clearly, the asymptotic expansion diverges at $\xi=v_{c}$, since there
the second derivative of $f$ vanishes, $f^{\prime\prime}(k_{0})=\omega_{\rm
BS}^{\prime\prime}(k_{0})=0$.

Figure 7: (a) Effective hopping amplitudes corresponding to the polariton
bands shown in Fig. 6(b). (b) Ballistic diffusion of an initially localized
polariton for the two different polariton bands.

These features of 1D arrangements of emitters only appear in the Markovian
regime. For $\lVert H_{I}\rVert\gg\lVert H_{B}\rVert,\lVert H_{S}\rVert$, the
effective hopping amplitudes decay monotonically as a function of distance,
and initially-localized excitations will diffuse with just a single wavefront
(see Appendix D). This can be understood noting that the single-emitter bound
states appearing below (above) the lowest (highest) bath band decay
monotonically as we move away from the emitter, and the strongly-hybridized
polaritons have energies in those ranges.

###  V.2 Haldane model

Now, we turn our attention to the Haldane model. The setup we analyze is
essentially the same as for the Harper-Hofstadter model, shown in Fig. 6(a),
albeit in this case we consider emitters coupled to both sites of each unit
cell along a given direction, in order to make $H_{B}^{\circ}$ correspond to a
bath with open boundary conditions. Thus, $H^{\mathrm{1D}}(k_{x})$ corresponds
in this case to a ladder with two emitters, each coupled to one of the sites
in a given unit cell. As it happened for the Harper-Hofstadter model, the
self-energy matrix can also be computed analytically in this case (see
Appendix C.2). Since there are two emitters per unit cell of the line defect,
we obtain two bands for the polariton modes, which are the solutions of

$\det\left[\omega_{\mathrm{BS},\pm}-\Delta-\Sigma(\omega_{\mathrm{BS},\pm};k_{x})\right]=0\,,$
(22)

or equivalently,
$\omega_{\mathrm{BS},\pm}-\Delta-\lambda_{\pm}(\omega_{\mathrm{BS},\pm};k_{x})=0$,
where $\lambda_{\pm}(\omega;k_{x})$ denote the eigenvalues of the self-energy
matrix $\Sigma(\omega;k_{x})$,

$\lambda_{\pm}=\frac{1}{2}\left[\Sigma_{11}+\Sigma_{22}\pm\sqrt{(\Sigma_{11}-\Sigma_{22})^{2}-4\lvert\Sigma_{12}\rvert^{2}}\right]\,.$
(23)

Again, for any crossing of $\Delta$ with the edge-mode frequencies, we obtain
a perfect vacancy-like bound state, see Fig. 8(a). Note that degeneracies of
the edge states translate into degeneracies of the vacancy-like bound states.
Thus, for topologically non-trivial Haldane baths, there exist a value of
$\Delta$—tuned to the crossing of the edge states—such that the resulting
polariton bands are degenerate for some $k_{x}$. In general, we cannot expect
degeneracies of the polariton bands to occur for a range of values of $\Delta$
and $k_{x}$, since, according to Eq. (23), they require that both
$\Sigma_{11}=\Sigma_{22}$ and $\Sigma_{12}=0$ simultaneously. In practice,
however, we observe that the effective polariton bands remain almost gapless
regardless the value of $\Delta$, as long as it lies in a topologically non-
trivial gap, see Fig. 8(b). By contrast, for the same emitter arrangement, the
polariton bands are always gapped if the bath is in the trivial phase. A
difference with respect to the Harper-Hofstadter model is that now the bound
states for finite line defects do not have to be symmetric with respect to the
emitter line, as shown in Fig. 8(c).

Figure 8: (a) Energy spectrum of $H_{B}^{\circ}$, and polariton bands
$\omega_{\mathrm{BS},\pm}$ (black lines) for a set of emitters coupled to the
Haldane model forming a line defect. The parameters of the bath are
$t_{2}=0.1t_{1}$, $\phi=0.8$ and $M\simeq 0.11t_{1}$ (topological) or $M\simeq
0.63t_{1}$ (trivial). The spectrum of $H_{B}^{\circ}$ is color-coded according
to the localization of the eigenstates with respect to the line defect. Below
we show a zoom to the polariton bands, which are color-coded according to the
weight of the eigenstates in each sublattice,
$\sigma_{z}\equiv\sigma^{+}_{A}\sigma^{-}_{A}-\sigma^{+}_{B}\sigma^{-}_{B}$.
The dashed horizontal line marks the value of $\Delta\simeq 0.25t_{1}$; and
$g=0.5t_{1}$. (b)
$\mathrm{Gap}:=\min_{k_{x}}\lvert\omega_{\mathrm{BS},+}(k_{x})-\omega_{\mathrm{BS},-}(k_{x})\rvert$,
and average chiral group velocity, Eq. (24), for the same values of the bath
parameters and varying emitter detuning; $g=0.1t_{1}$. The inset in the upper
plot is a zoom to the values in the topological case. (c) Bound state with
energy closest to $\Delta$, for a finite emitter line defect in a topological
bath with the same parameters as in panel (a). Figure 9: Emitter dynamics for
an initial state consisting of a single excitation localized in sublattice
$A$, $\ket{\psi(0)}=\sigma^{+}_{L/2,A}\ket{\rm vac}$. The parameters of the
system are $\Delta\simeq 0.25t_{1}$, $g=0.1t_{1}$, $t_{2}=0.1t_{1}$,
$\phi=0.8$ and $M\simeq 0.11t_{1}$ (topological) or $M\simeq 0.63t_{1}$
(trivial).

Another difference with respect to the case studied in the Harper-Hofstadter
model is the sublattice degree of freedom present in the emitter line defect,
which allows for truly complex effective hoppings. One of the consequences of
this, is the chiral dynamics that occur when an excitation is initially
localized in one of the emitter sublattices, as demonstrated in Fig. 9. We can
understand this effect computing the average chiral group velocity
$\overline{v}_{C}$ [24], which quantifies the average group velocity of the
polaritons in the system, weighted by their relative occupation in each
emitter sublattice. Specifically,

$\overline{v}_{C}=\frac{1}{2}\sum_{\alpha=\pm}\int_{-\pi}^{\pi}\frac{dk}{2\pi}v_{\alpha}(k)\langle\sigma_{z}\rangle_{k,\alpha}\,,$
(24)

where

$\langle\sigma_{z}\rangle_{k,\alpha}\equiv\bra{\psi_{\mathrm{BS},\alpha}(k)}\sigma^{+}_{A,k}\sigma^{-}_{A,k}-\sigma^{+}_{B,k}\sigma^{-}_{B,k}\ket{\psi_{\mathrm{BS},\alpha}(k)}\,,$
(25)

and $v_{\alpha}(k)=\partial_{k}\omega_{\mathrm{BS},\alpha}(k)$. As can be seen
in Fig. 8(b), the value of $\overline{v}_{C}$ for line defects in topological
baths is much larger than in the case of trivial baths. This chiral dynamics
is lost if the light-matter coupling constant is too large (see Appendix D).

##  VI 2D emitter arrangements

Going beyond 1D line defects, we can arrange the emitters periodically in the
two dimensions of the bath forming what we may call an “emitter superlattice”.
Such systems have been considered in the past as a platform where one could
generate strongly-correlated states of light and matter [25], fractional
quantum Hall states [26], or more generally, use them as quantum simulators
[27, 28, 29, 30]. One might expect the properties of the resulting polariton
models, in particular their topology, to depend strongly on the properties of
the underlying bath. In this section we show how this is not the only aspect
to take into account, but the geometry of the emitter arrangement and the
coupling strength $g$ also play a crucial role.

First, we consider the case in which there is an emitter coupled to every site
of every unit cell in some region of the bath. If $U_{\bm{k}}$ is a unitary
that diagonalizes the bath Hamiltonian,
$H_{B}=\sum_{\bm{k}}\tilde{\Psi}_{\bm{k}}^{\dagger}\Omega(\bm{k})\tilde{\Psi}_{\bm{k}}$,
with

$\displaystyle\Omega(\bm{k})$ $\displaystyle\equiv
U_{\bm{k}}\mathcal{H}_{B}(\bm{k})U^{\dagger}_{\bm{k}}=\operatorname{diag}\left(\omega_{1}(\bm{k}),\dots,\omega_{M}(\bm{k})\right)\,,$
(26) $\displaystyle\tilde{\Psi}_{\bm{k}}$ $\displaystyle\equiv
U_{\bm{k}}\Psi_{\bm{k}}=[\tilde{a}_{\bm{k},1},\dots,\tilde{a}_{\bm{k},M}]^{T}\,,$
(27)

we can apply the same transformation to the emitter raising and lowering
operators, defining
$[\tilde{\sigma}^{-}_{\bm{k},1},\dots,\tilde{\sigma}^{-}_{\bm{k},M}]^{T}=U_{\bm{k}}[\sigma^{-}_{\bm{k},1},\dots,\sigma^{-}_{\bm{k},M}]^{T}$,
such that the whole system Hamiltonian can be expressed as
$H=\sum_{j=1}^{M}\sum_{\bm{k}}H_{j}(\bm{k})$, with

$H_{j}(\bm{k})=\Delta\tilde{\sigma}^{+}_{\bm{k},j}\tilde{\sigma}^{-}_{\bm{k},j}+\omega_{j}(\bm{k})\tilde{a}^{\dagger}_{\bm{k},j}\tilde{a}_{\bm{k},j}\\\
+g(\tilde{\sigma}^{+}_{\bm{k},j}\tilde{a}_{\bm{k},j}+\mathrm{H.c.})\,.$ (28)

Since $H_{j}(\bm{k})$ only involves a single mode and a single emitter
raising/lowering operator, it can be easily diagonalized, yielding the
following eigenmodes and eigenenergies:

$\displaystyle\lambda_{j,\pm}(\bm{k})=\frac{\Delta+\omega_{j}(\bm{k})}{2}\pm\sqrt{\left(\frac{\Delta-\omega_{j}(\bm{k})}{2}\right)^{2}+g^{2}}\,,$
(29)
$\displaystyle\beta^{\dagger}_{\bm{k},j,+}=\cos\frac{\theta_{\bm{k}}}{2}\tilde{\sigma}^{+}_{\bm{k},j}+\sin\frac{\theta_{\bm{k}}}{2}\tilde{a}^{\dagger}_{\bm{k},j}\,,$
(30)
$\displaystyle\beta^{\dagger}_{\bm{k},j,-}=\sin\frac{\theta_{\bm{k}}}{2}\tilde{\sigma}^{+}_{\bm{k},j}-\cos\frac{\theta_{\bm{k}}}{2}\tilde{a}^{\dagger}_{\bm{k},j}\,,$
(31)
$\displaystyle\cos\theta_{\bm{k}}=\frac{\Delta-\omega_{j}(\bm{k})}{\sqrt{[\Delta-\omega_{j}(\bm{k})]^{2}+4g^{2}}}\,.$
(32)

We can see that in the Markovian regime,
$g\ll\lvert\Delta-\omega_{j}(\bm{k})\rvert$, $\lambda_{j,+}\simeq\Delta$,
$\lambda_{j,-}\simeq\omega_{j}(\bm{k})$,
$\beta_{\bm{k},j,+}\simeq\tilde{\sigma}_{\bm{k},j}$ and
$\beta_{\bm{k},j,-}\simeq\tilde{a}_{\bm{k},j}$, assuming
$\Delta>\omega_{j}(\bm{k})$ (the same is true also for
$\Delta<\omega_{j}(\bm{k})$, interchanging the labels $+\leftrightarrow-$).
Since the transformation $U_{\bm{k}}$ is the same for both the bath modes and
the emitter modes, the Chern numbers of the polariton bands will be the same
as those of the original bath’s bands. On the other hand, in the strong-
coupling limit, $g\gg\lvert\Delta-\omega_{j}(\bm{k})\rvert$, each bath’s band
gives rise to two strongly-hybridized polariton bands with dispersions
$\lambda_{j,\pm}(\bm{k})\simeq[\Delta+\omega_{j}(\bm{k})]/2\pm g$, and
eigenmodes
$\beta_{\bm{k},j,\pm}\simeq(\tilde{\sigma}^{+}_{\bm{k},j}\pm\tilde{a}^{\dagger}_{\bm{k},j})/\sqrt{2}$.
It is easy to see that the Chern number of these strongly-hybridized-polariton
bands is the same as the Chern number of the free bath band that originates
them. In fact, from Eq. (29) it is clear that there are no band touchings in
the spectrum of the whole system, provided the free bath spectrum is gapped,
$\omega_{i}(\bm{k})\neq\omega_{j}(\bm{k})$ for all $i\neq j$. Thus, the Chern
numbers of all the bands are constant and do not depend on the coupling
strength $g$.

A very different situation happens if we only couple emitters to a single
sublattice. In the Markovian regime, it is easy to see that the Chern number
of the single polariton band must be 0, since the sum of Chern numbers of all
the bands must vanish [31], and in the limit $g\to 0$ the presence of emitters
does not change the Chern numbers of the original bath’s bands. On the other
hand, the strongly-hybridized-polariton bands should also have a zero Chern
number, since the corresponding eigenmodes are approximately given by
$\beta^{\dagger}_{\pm}\simeq(\sigma^{\dagger}_{\bm{k}}\pm
a^{\dagger}_{\bm{k},1})/\sqrt{2}$ (assuming, w.l.o.g. that the emitters are
coupled to the first sublattice). This behavior is exemplified in Fig. 10,
where we show two topological phase diagrams for a system of quantum emitters
coupled to the $A$ sublattice of the Haldane model. For small coupling
strengths the phase diagram is similar to that of the bare Haldane model, with
new phases emerging near the original phase transition points. The main two
lobes (red and blue regions) correspond to models where the center band is
trivial (Markovian polariton band) and the bath’s bands remain unaffected by
the presence of emitters. As we increase $g$, the area of these two lobes is
reduced, until they completely disappear for $g>3\sqrt{6}t_{2}$. In the limit
$g\to\infty$ all the bands become trivial (fixing the values of the rest of
parameters). This can be understood noting that, in this limit, the emitters
effectively decouple the two sublattices of the Haldane model. As this case
shows, the presence of emitters can also modify the topology of the guided
modes in the bath.

Figure 10: Topological phase diagram of the Haldane model with quantum
emitters coupled to the $A$ sublattice in the weak (a) and strong (b) coupling
regime. Different phases are marked using different colors. Each color
corresponds to a different set of Chern numbers as indicated in the legend
(from the lowest- to the highest-energy band). The emitter frequency is tuned
in all the cases to the middle of the inner bandgap $\Delta=3t_{2}\cos(\phi)$.
The dashed red and black curves, with equation $M=\pm
3\sqrt{3}t_{2}\sin(\phi)$, correspond to the phase boundaries of the bare
Haldane model. In the presence of emitters, the boundaries are given by $M=\pm
3\sqrt{3}t_{2}\sin(\phi)\pm g/\sqrt{2}$ (solid red and black curves).

We can argue intuitively that for strongly-hybridized polariton bands it is
necessary to have a dense emitter superlattice to have non-trivial Chern
numbers, since the photonic component of these bound states is very localized.
For sparse emitter supperlatices, topological polariton bands only occur in
the Markovian regime. Furthermore, for dense superlattices the bands of guided
modes in the bath can become trivial in the strong-coupling limit, while for
sufficiently sparse lattices, the topology of the bath’s bands remains
unaffected by the presence of emitters.

##  VII Conclusion

To sum up, in this work we have explored the physics of polaritons appearing
in systems of quantum emitters coupled to two-dimensional topological baths.
In particular, we have focused on two paradigmatic examples of topological
insulators: the Haldane model and the Harper-Hofstadter model. In the Haldane
model, we have shown how, in a certain parameter regime, the effective
couplings resemble the ones expected for the honeycomb lattice, and for
emitter detunings in the bandgap they break time-reversal symmetry. A single
emitter coupled to the Haldane model can seed a quasibound state, which gives
rise to truly long-range interactions scaling like $\sim 1/r$, where $r$ is
the distance between the emitters. In constrast with the quasibound states
found for a single emitter in other lattices (e.g., the honeycomb lattice), it
is not a vacancy-like bound state, and therefore is not robust to changes in
the light-matter coupling constant. For one-dimensional emitter arrangements
(emitter line defects) coupled to the Haldane model in the Markovian regime,
the gap of the effective polariton bands in the topological phase is orders of
magnitude smaller than the one found in the trivial phase, and in the former
case the polariton dynamics is chiral. In the Harper-Hofstadter model, the
Chern numbers of the surrounding bands dictate the non-monotonicity of the
effective couplings, which give rise to a characteristic ballistic diffusion
of polaritons with several wavefronts. For two-dimensional emitter
arrangements (emitter superlattices), we have shown how, in the case of one
emitter coupled to every bath site, the resulting polariton bands always have
the same Chern numbers of the original bath’s bands. For other geometries,
however, it is possible to have different Chern numbers, and even change the
topology of the bath’s guided modes.

One aspect that we have not addressed in this work is the effect of disorder
and imperfections in the emitter arrangements. For example, it could be
interesting to consider systems where the emitters are coupled at random
locations in the bath, and see whether the resulting amorphous effective
models are topologically non-trivial [32]. Another interesting future research
direction is to go beyond single-excitation physics. Within the Markovian
regime, the emitter dynamics can be described with an effective XX-type spin
Hamiltonian. Thus, for an emitter line-defect in the Harper-Hofstadter model,
the non-monotonic long-range couplings could result in multiple-component
Luttinger-liquid phases. In the Haldane model, the possibility of having
complex, truly long-range couplings could result in interesting many-body
phases that have not been described previously.

## Acknowledgements

M.B. and J.I.C. acknowledge funding from the ERC Advanced Grant QUENOCOBA
under the EU Horizon 2020 program (Grant Agreement No. 742102).

##  Appendix A Single-particle bound states

In this section we review how to compute the single-particle bound states for
a system of emitters coupled locally, i.e., each one to a single site of a
$D$-dimensional bath, in the thermodynamic limit ($L\to\infty$). Their shapes
and energies are intimately related to the properties of the level-shift
operator (aka self-energy) $\Sigma(z)$ [17]. For such systems, the matrix
elements of the level-shift operator in the basis of the single-excitation
emitter subspace $\\{\ket{n}\equiv\sigma_{n}^{+}\ket{\Omega}\\}_{n=1}^{N}$,
where $\ket{\Omega}\equiv\ket{g}^{\otimes N}\ket{\mathrm{vac}}$, are given by

$\begin{split}\left[\Sigma(z)\right]_{mn}&\equiv\bra{m}\Sigma(z)\ket{n}\\\
&=g^{2}\left[\Phi(z,\bm{r}_{m}-\bm{r}_{n})\right]_{\alpha_{m},\alpha_{n}}\,,\end{split}$
(33)

with (the integral is over the first Brillouin zone)

$\Phi(z,\bm{r})=\int_{\mathrm{BZ}}\frac{d^{D}\bm{k}}{(2\pi)^{D}}\,e^{i\bm{k}\cdot\bm{r}}\left[z-\mathcal{H}_{B}(\bm{k})\right]^{-1}\,.$
(34)

With a slight abuse of notation, we can denote by $\Sigma(z)$ not the level-
shift operator itself, but the above mentioned matrix.

Bound states are solutions to the eigenvalue equation
$H\ket{\psi_{\mathrm{BS}}}=E_{\mathrm{BS}}\ket{\psi_{\mathrm{BS}}}$ with
$E_{\mathrm{BS}}$ belonging to the discrete spectrum of $H$. Any state within
the single-excitation subspace can be written as

$\ket{\psi_{\mathrm{BS}}}=\left[\sum_{n}c_{e,n}\sigma_{n}^{+}+\sum_{\alpha,\,\bm{k}}c_{\alpha,\bm{k}}a_{\alpha,\bm{k}}^{\dagger}\right]\ket{\Omega}\,.$
(35)

Denoting $\bm{c}_{e}\equiv\left(c_{e,1},\dots,c_{e,N_{e}}\right)^{T}$,
$\bm{c}_{\bm{k}}\equiv\left(c_{1,\bm{k}},\dots,c_{M,\bm{k}}\right)^{T}$, and
$G_{\bm{k}}$ a $N_{e}\times M$ matrix with elements
$\left(G_{\bm{k}}\right)_{nm}=gL^{-D/2}e^{i\bm{k}\bm{r}_{n}}\delta_{m,\alpha_{n}}$,
the eigenvalue equation translates as

$\displaystyle\Delta\bm{c}_{e}+\sum_{\bm{k}}G_{\bm{k}}\bm{c}_{\bm{k}}=E_{\mathrm{BS}}\bm{c}_{e}\,,$
(36) $\displaystyle
G_{\bm{k}}^{\dagger}\bm{c}_{e}+\mathcal{H}_{B}(\bm{k})\bm{c}_{\bm{k}}=E_{\mathrm{BS}}\bm{c}_{\bm{k}}\,,\
\forall\bm{k}\in\mathrm{BZ}\,.$ (37)

From Eq. (37),

$\bm{c}_{\bm{k}}=\left[E_{\mathrm{BS}}-\mathcal{H}_{B}(\bm{k})\right]^{-1}G_{\bm{k}}^{\dagger}\bm{c}_{e}\,.$
(38)

Substituting back in Eq. (36), one can see that the amplitudes of the emitter
excited states are nontrivial solutions of

$\left[E_{\mathrm{BS}}-\Delta-\Sigma(E_{\mathrm{BS}})\right]\bm{c}_{e}=0\,,$
(39)

therefore, the energies of all bound states (in the thermodynamic limit) are
the solutions to the nonlinear equation

$\det\left[E_{\mathrm{BS}}-\Delta-\Sigma(E_{\mathrm{BS}})\right]=0\,.$ (40)

For a given $E_{\mathrm{BS}}$, the emitter amplitudes can be obtained solving
Eq. (39), and with them, the photonic part can be computed using Eq. (38).
Importantly, Eq. (39) does not fully determine $\bm{c}_{e}$. In order to do
so, one has to use the normalization condition
$\braket{\psi_{\mathrm{BS}}}{\psi_{\mathrm{BS}}}=1$, which implies

$\bm{c}_{e}^{\dagger}\left[1-\Sigma^{\prime}(E_{\mathrm{BS}})\right]\bm{c}_{e}=1\,,$
(41)

where $\Sigma^{\prime}(z)\equiv\partial_{z}\Sigma(z)$. The bath amplitudes in
position space can be obtained by an inverse Fourier transform of Eq. (38),
obtaining

$c_{\alpha,\bm{r}}=g\sum_{n}c_{e,n}\left[\Phi(E_{\mathrm{BS}},\bm{r}-\bm{r}_{n})\right]_{\alpha,\alpha_{n}}\,.$
(42)

Thus, knowing the self-energy, we can readily compute the photonic component
of any bound state of a given energy.

##  Appendix B Approximate expressions for the self-energy in the Haldane
model

The self-energy is essentially a Fourier transform of the kernel function

$\mathcal{K}(\bm{k})\equiv\frac{f_{\alpha_{m}\alpha_{n}}(z,\bm{k})}{[z-h_{0}(\bm{k})]^{2}-h^{2}(\bm{k})}\equiv\frac{f_{\alpha_{m}\alpha_{n}}(z,\bm{k})}{d(z,\bm{k})}\,.$
(43)

Note that $d(z,\bm{k})=\det[z-\mathcal{H}_{B}(\bm{k})]$. For
$\lvert\bm{r}_{mn}\rvert\gg 1$, the integrand is highly oscillatory, so it is
the sharpest features of $\mathcal{K}$ that give non-negligible contributions
to the integral. If $z=\Delta$ lies above or below the bands, in the regime
$\lvert M\rvert,\lvert t_{2}\rvert\ll\lvert t_{1}\rvert$, the magnitude of the
kernel is peaked around the $\Gamma$-point, $\bm{\Gamma}\equiv(0,0)$, which is
actually an extremum of $d$. Thus, we can obtain an approximation to the self-
energy expanding

$d(\Delta,\bm{\Gamma}+\bm{q})\simeq
d(\Delta,\bm{\Gamma})+\frac{1}{2}\bm{q}^{T}\,{\rm H}\,\bm{q}\,,$ (44)

where $\rm H$ denotes the Hessian of $d$ evaluated at the $\Gamma$-point,
which turns out to be of the form

$\rm H=2c\begin{pmatrix}1&-1/2\\\ -1/2&1\end{pmatrix}\,.$ (45)

With a linear transformation

$\bm{p}=\frac{1}{2}\begin{pmatrix}\sqrt{3}&-\sqrt{3}\\\
1&1\end{pmatrix}\bm{q}\equiv L\bm{q}\,,$ (46)

we can express

$\Sigma_{mn}(\Delta)\simeq\frac{g^{2}}{(2\pi)^{2}}\int_{\mathbb{R}^{2}}d^{2}\bm{p}\,\frac{ae^{i\bm{p}\cdot\bm{r}_{mn}^{\prime}}}{b+cp^{2}}\,,$
(47)

where $\bm{r}^{\prime}_{mn}=(L^{-1})^{T}\bm{r}_{mn}$,
$p\equiv\lVert\bm{p}\rVert$, and real constants

$\begin{split}a&=f_{\alpha_{m}\alpha_{n}}(\Delta,\bm{\Gamma})\lvert\det
L\rvert^{-1}\\\
&=\frac{2}{\sqrt{3}}\times\begin{cases}\Delta+6t_{2}\cos(\phi)\pm
M\,,&\text{if }\alpha_{m}=\alpha_{n}\\\ -3t_{1}\,,&\text{if
}\alpha_{m}\neq\alpha_{n}\end{cases}\end{split}$ (48) $\displaystyle b$
$\displaystyle=d(\Delta,\bm{\Gamma})=[\Delta+6t_{2}\cos(\phi)]^{2}-M^{2}-9t_{1}^{2}\,,$
(49) $\displaystyle c$
$\displaystyle=2t_{1}^{2}-4t_{2}\cos(\phi)[\Delta+6t_{2}\cos(\phi)]\,.$ (50)

In Eq. (48), the term $+M$ or $-M$ has to be chosen if $\alpha_{m}=A$ or
$\alpha_{m}=B$, respectively. Using polar coordinates, integrating the angular
variable using the Jacobi-Anger expansion, we obtain

$\Sigma_{mn}(\Delta)\simeq\frac{g^{2}a}{2\pi}\int_{0}^{\infty}dp\,\frac{p\operatorname{J}_{0}(pr^{\prime}_{mn})}{b+cp^{2}}\,.$
(51)

Last, integrating the radial variable we obtain

$\Sigma_{mn}(\Delta)\simeq\frac{g^{2}a}{2\pi
c}\operatorname{K}_{0}\left(\sqrt{\frac{b}{c}}r^{\prime}_{mn}\right)\,.$ (52)

In these expressions, $\operatorname{J}_{n}$ and $\operatorname{K}_{n}$ denote
the $n$th-order Bessel function of the first kind and modified Bessel function
of the second kind, respectively.

When $\Delta$ lies in between the two bands, the sharpest features of the
integrand are instead located at the K-points,
$\pm\bm{K}\equiv\pm(1,-1)2\pi/3$. Thus, following an analogous procedure for
each of these points, we can approximate the self-energy by the sum of two
different contributions $\Sigma_{mn}(\Delta)\simeq
I_{mn}^{+}(\Delta)+I_{mn}^{-}(\Delta)$. When $\alpha_{m}=\alpha_{n}$,
$I^{\pm}_{mn}(\Delta)$ has the same form as shown in the Eq. (52), with

$\begin{split}&a_{\pm}=\frac{2e^{\pm i\bm{K}\cdot\bm{r}_{mn}}}{\sqrt{3}}\\\
&\times\begin{cases}\Delta+M-3t_{2}[\operatorname{c}(\phi)\mp\sqrt{3}\operatorname{s}(\phi)]\,,&\text{if
}\alpha_{m}=A\\\
\Delta-M-3t_{2}[\operatorname{c}(\phi)\pm\sqrt{3}\operatorname{s}(\phi)]\,,&\text{if
}\alpha_{m}=B\\\ \end{cases}\end{split}$ (53) $\displaystyle b_{\pm}$
$\displaystyle=[\Delta-3t_{2}\operatorname{c}(\phi)]^{2}-[M\pm
3\sqrt{3}t_{2}\operatorname{s}(\phi)]^{2}\,,$ (54)
$\displaystyle\begin{split}c_{\pm}&=-t_{1}^{2}+6t_{2}^{2}[1-2\operatorname{c}(2\phi)]\\\
&\qquad+2t_{2}[\Delta\operatorname{c}(\phi)\pm\sqrt{3}M\operatorname{s}(\phi)]\,,\end{split}$
(55)

where we have used the shorthand notation
$\operatorname{c}(\phi)\equiv\cos(\phi)$ and
$\operatorname{s}(\phi)\equiv\sin(\phi)$. For $\alpha_{m}\neq\alpha_{n}$, the
numerator in the kernel vanishes at zeroth order,
$f_{\alpha_{m}\alpha_{n}}(\Delta,\pm\bm{K})=0$, so we have to consider the
next order to obtain a good approximation. Doing so, using polar coordinates
and integrating the angular variable, we obtain

$\begin{split}I^{\pm}_{mn}(\Delta)&=\frac{g^{2}a_{\pm}}{2\pi}\int_{0}^{\infty}dp\,\frac{p^{2}\operatorname{J}_{1}(pr^{\prime}_{mn})}{b_{\pm}+c_{\pm}p^{2}}\\\
&=\frac{g^{2}a_{\pm}}{2\pi
c_{\pm}}\sqrt{\frac{b_{\pm}}{c_{\pm}}}\operatorname{K}_{1}\left(\sqrt{\frac{b_{\pm}}{c_{\pm}}}r^{\prime}_{mn}\right)\,,\end{split}$
(56)

where $b_{\pm}$ and $c_{\pm}$ are the same as in Eqs. (54) and (55), and

$a_{\pm}=\frac{2t_{1}e^{\pm
i\bm{K}\cdot\bm{r}_{mn}}}{\sqrt{3}}\times\begin{cases}e^{\mp
i(\theta_{mn}-\pi/2)}\,,&\text{if }\alpha_{m}=A\\\ e^{\pm
i(\theta_{mn}+\pi/2)}\,,&\text{if }\alpha_{m}=B\end{cases}\,.$ (57)

Here, $\theta_{mn}$ is the angle between $\bm{r}^{\prime}_{mn}$ and the
horizontal axis.

##  Appendix C One-dimensional self-energies

When viewing 2D systems with an emitter line-defect as a collection of 1D
systems, the self-energy for each one of the 1D systems can be obtained using
Eq. (7), integrating only the quasimomentum perpendicular to the emitter line-
defect.

###  C.1 Harper-Hofstadter model

In this case, $H_{B}^{\mathrm{1D}}(k)$ corresponds to a chain with nearest-
neighbor hopping of amplitude $J$, and on-site chemical potentials
$\\{\mu_{n}=-2J\cos(k+n\phi):n\in\mathbb{Z}\\}$,

$H_{B}^{\rm
1D}=-J\sum_{n}\left(c^{\dagger}_{n+1}c_{n}+c^{\dagger}_{n}c_{n+1}\right)+\sum_{n}\mu_{n}c^{\dagger}_{n}c_{n}\,,$
(58)

where we assume that the bosonic creation and annihilation operators
$\\{c^{\dagger}_{n},c_{n}\\}_{n\in\mathbb{Z}}$ depend implicitly on the
parallel quasimomentum $k$. Since $\phi=2\pi p/q$ with co-prime integers $p$
and $q$, the chemical potentials are periodic with period $q$ (there are $q$
sites per unit cell). We may assume that the emitter is coupled at $n=1$,
since a redefinition of the quasimomentum $k_{x}\to k_{x}+m\phi$ (equivalent
to choosing a different gauge for the magnetic vector potential) corresponds
to a cyclic permutation of the chemical potentials $\mu_{n}\to\mu_{n+m}$. In
addition, the system is invariant under the combined action of time-reversal
and space inversion about the emitter position,
$\mu_{n}(k-\phi)=\mu_{q-n+2}(-k-\phi)$.

The matrix $\left[z-\mathcal{H}_{B}(k_{y})\right]/J$ is of the form

$A=\begin{pmatrix}a_{1}&1&&&&e^{-ik_{y}}\\\ 1&a_{2}&1&&&\\\
&1&\ddots&\ddots&&\\\ &&\ddots&&&\\\ &&&&&1\\\
e^{ik_{y}}&&&&1&a_{q}\end{pmatrix}\,,$ (59)

with $a_{n}=(z-\mu_{n})/J$. To compute the self-energy, we have to compute its
inverse, whose matrix elements are given by $(A^{-1})_{mn}=C_{nm}/\det A$.
Here $C$ is the cofactor matrix associated to $A$. Using Laplace’s formula to
compute determinants, we find (for $q\geq 3$ and $2<p<q$)

$\displaystyle C_{11}$ $\displaystyle=K(a_{2},\dots,a_{q})\,,$ (60)
$\displaystyle C_{12}$
$\displaystyle=-\left[K(a_{3},\dots,a_{q})+(-1)^{q}e^{ik_{y}}\right]\,,$ (61)
$\displaystyle\begin{split}C_{1p}&=(-1)^{p+1}\big{[}K(a_{p+1},\dots,a_{q})\\\
&\qquad\qquad\qquad+(-1)^{q}e^{ik_{y}}K(a_{2},\dots,a_{p-1})\big{]}\,,\end{split}$
(62) $\displaystyle C_{1q}$
$\displaystyle=(-1)^{q+1}\left[1+(-1)^{q}e^{ik_{y}}K(a_{2},\dots,a_{q-1})\right]\,,$
(63) $\displaystyle\begin{split}\det
A&=K(a_{1},\dots,a_{q})-K(a_{2},\dots,a_{q-1})\\\ &\qquad+(-1)^{q+1}2\cos
k_{y}\,,\end{split}$ (64)

where we have expressed the different quantities using the continuant $K$,
defined as

$K(x_{1},\dots,x_{r})\equiv\det\begin{pmatrix}x_{1}&1&&&&\\\ 1&x_{2}&1&&&\\\
&1&\ddots&\ddots&&\\\ &&\ddots&&&\\\ &&&&&1\\\ &&&&1&x_{r}\end{pmatrix}\,.$
(65)

Note that this continuant satisfies the recurrence relation:

$K_{n}=x_{n}K_{n-1}-K_{n-2}\,,\quad 1\leq n\leq r\,,$ (66)

where $K_{n}\equiv K(x_{1},\dots,x_{n})$, and we set $K_{0}=1$, $K_{-1}=0$.

Finally, the integral giving the emitter self-energy is straightforward to
compute using residue integration. For example:

$\displaystyle\left[\Phi(z,y)\right]_{11}$ $\displaystyle=\frac{1}{2\pi
J}\int_{-\pi}^{\pi}dk_{y}\,\frac{K(a_{2},\dots,a_{q})e^{ik_{y}y}}{K(a_{1},\dots,a_{q})-K(a_{2},\dots,a_{q-1})+(-1)^{q+1}2\cos
k_{y}}$ (67) $\displaystyle=\frac{(-1)^{q+1}K(a_{2},\dots,a_{q})}{2\pi
iJ}\oint dw\,\frac{w^{\lvert y\rvert}}{w^{2}+bw+1}$ (68)
$\displaystyle=\frac{(-1)^{q+1}K(a_{2},\dots,a_{q})}{J\sqrt{b^{2}-4}}\left[w_{+}^{\lvert
y\rvert}\Theta(1-\lvert w_{+}\rvert)-w_{-}^{\lvert y\rvert}\Theta(1-\lvert
w_{-}\rvert)\right]\,.$ (69)

To go from Eq. (67) to Eq. (68) we perform a change of variable
$w=\exp[ik_{y}\operatorname{sign}(y)]$, such that the integral becomes a
contour integral in the complex plane along the unit circumference
(anticlockwise). Last, in Eq. (69) we express the result in terms of the poles
$w_{\pm}=(-b\pm\sqrt{b^{2}-4})/2$, which are the roots of the second-order
polynomial $p(w)=w^{2}+bw+1$, with

$b=(-1)^{q+1}\left[K(a_{1},\dots,a_{q})-K(a_{2},\dots,a_{q-1})\right]\,;$ (70)

$\Theta$ denotes Heaviside’s step function. Since $p$ is a palindromic
polynomial (its coefficients form a palindrome), its roots are one the inverse
of the other, $w_{+}w_{-}=1$. As a consequence, only one of the poles
contributes to the integral in general, except when $\lvert w_{+}\rvert=\lvert
w_{-}\rvert=1$, which occurs whenever $z$ belongs to the spectrum of
$H^{\mathrm{1D}}_{B}$. In that case, the integral in Eq. (68) is not well
defined, and the self-energy is discontinuous. For real $z$, $b$ is real, and
if $z$ is not in the spectrum of $H^{\mathrm{1D}}_{B}$, then the poles
$w_{\pm}$ are also real. In this case, the pole inside the unit circumference
is $w_{\mathrm{in}}=w_{\operatorname{sign}(b)}$, and Eq. (69) can be
simplified as

$\left[\Phi(z,y)\right]_{11}=\frac{\operatorname{sign}(b)(-1)^{q+1}K(a_{2},\dots,a_{q})}{J\sqrt{b^{2}-4}}w_{\mathrm{in}}^{\lvert
y\rvert}\,.$ (71)

The self-energy appearing in the pole equation, Eq. (20), used to compute the
polariton band is simply
$\Sigma(\omega_{\mathrm{BS}})=g^{2}\left[\Phi(\omega_{\mathrm{BS}},0)\right]_{11}$.
It depends implicitly on the quasimomentum $k_{x}$ through the on-site
energies $\\{a_{n}\\}_{n=1}^{q}$. The zeroes of this function, which
correspond to the vacancy-like bound states, are the zeroes of the continuant
$K(a_{2},\dots,a_{q})$.

Once we know the energy of a bound state, we can compute analytically its
wavefunction using Eq. (42). For this, we need to compute
$\left[\Phi(z,y)\right]_{n1}$, for real $z$ not in the spectrum of
$H^{\mathrm{1D}}_{B}$. From Eqs. (60) to (63), we realize that the cofactors
have the form $C_{1n}=F_{n}+G_{n}e^{ik_{y}}$. Thus,

$\left[\Phi(z,y)\right]_{n1}=\frac{\operatorname{sign}(b)(-1)^{q+1}\left[F_{n}+G_{n}w_{\mathrm{in}}^{\operatorname{sign}(y)}\right]}{J\sqrt{b^{2}-4}}w_{\mathrm{in}}^{\lvert
y\rvert}\,.$ (72)

Eqs. (71) and (72) allow us to extract the following conclusions regarding the
shape of the bound states:

  1. 1.

On the sublattice to which the emitter is coupled, the amplitudes are
symmetric with respect to the emitter position. This is not the case for the
other sublattices.

  2. 2.

In general, the bound state is localized exponentially around the emitter
position, with the same decay length in all sublattices. It may decay
monotonically or alternate sign, depending on the sign of $w_{\mathrm{in}}$.

  3. 3.

For a vacancy-like bound state, the amplitude of the wavefunction on the
sublattice to which the emitter is coupled vanishes everywhere.

Except in the case of accidental degeneracies, the vacancy-like bound states
are completely localized on one side of the emitter. This can be understood,
since the vacancy modes of the system are the eigenstates of the two semi-
infinite chains in which the emitter divides the bath. As each chain is
characterized by a different order of the on-site potentials, the vacancy
modes will not be degenerate in general.

We remark that the conclusions drawn here are valid for arbitrary 1D systems
of the form shown in Eq. (58) and not just for the emitter-line defects in the
Harper-Hofstadter model. That is, Eqs. (71) and (72) are valid regardless the
specific values of the chemical potentials, as long as they are $q$-periodic,
$\mu_{n}=\mu_{n+q}$.

###  C.2 Haldane model

In this case, $H_{B}^{\mathrm{1D}}(k)$ corresponds to a ladder with complex
hopping amplitudes, such as the one shown in Fig. 11.

Figure 11: (a) Schematics of the 1D model relevant for the study of the bound
states in the Haldane model. The on-site potentials depend on the original
model’s parameters as $\mu_{A}=-2t_{2}\cos(k_{x}+\phi)+M$, and
$\mu_{B}(\phi,M)=\mu_{A}(-\phi,-M)$; the hoppings, which may be complex, are
given by $t_{AA}=-t_{2}e^{i\phi}-t_{2}e^{i(k_{x}-\phi)}$,
$t_{BB}=t_{AA}(-\phi)$, $t_{AB}=-t_{1}$, $t_{BA}=-t_{1}(1+e^{ik_{x}})$. The
arrows indicate the hopping direction, with the opposite direction having the
complex-conjugated hopping amplitude.

Instead of focusing on the Haldane model, in this section we describe a
general procedure to compute the self-energy, valid for any two-band, one-
dimensional bath with finite hopping range $R$. The momentum-resolved
Hamiltonian for such models can always be expressed in terms of the Pauli
matrices $\\{\sigma^{\nu}\\}_{\nu=x,y,z}$ and the identity matrix
$\sigma^{0}$, as
$\mathcal{H}_{B}(k)=\sum_{\alpha}h_{\alpha}(k)\sigma^{\alpha}$. The functions
$h_{\alpha}$ are finite Laurent series of the variable $w=e^{ik}$, i.e., they
are functions of the form

$h_{\alpha}(w)=a_{0}w^{-R}+a_{1}w^{-R+1}+\dots+a_{2R}w^{R}\,,$ (73)

so they are completely determined by the vector of coefficients
$[a_{0},\dots,a_{2R}]$, which are themselves simple linear combinations of the
hopping amplitudes of the model. We define the polynomials
$p_{\alpha}(w)=w^{R}h_{\alpha}(w)$, which have the same vector of coefficients
in the standard polynomial basis. Hermiticity of $\mathcal{H}_{B}$ then
implies that these polynomials are self-reciprocal, i.e., their coefficients
satisfy $a_{n}^{*}=a_{2R-n}$, for $0\leq n\leq R$.

For a 2-by-2 matrix, the inverse
$\left(z-\mathcal{H}_{B}\right)^{-1}=C^{T}/\det\left(z-\mathcal{H}_{B}\right)$
can be computed straightforwardly,

$\displaystyle\det\left(z-\mathcal{H}_{B}\right)=(z-h_{0})^{2}-h_{x}^{2}-h_{y}^{2}-h_{z}^{2}$
(74) $\displaystyle C^{T}=\begin{pmatrix}z-h_{0}+h_{z}&h_{x}-ih_{y}\\\
h_{x}+ih_{y}&z-h_{0}-h_{z}\end{pmatrix}\,.$ (75)

Thus, doing a change of variable $w=\exp[\operatorname{sign}(y)ik]$, we can
express the self-energy as an integral of a rational function over the unit
cricumference in the complex plane,

$\left[\Phi(z,y)\right]_{\alpha\beta}=\frac{1}{2\pi i}\oint dw\,w^{\lvert
y\rvert+R-1}\frac{P_{\alpha\beta}(w)}{Q(w)}\,.$ (76)

For $y\geq 0$,

$Q(w)=\left[zw^{R}-p_{0}(w)\right]^{2}-p_{x}(w)-p_{y}(w)-p_{z}(w)\,,$ (77)

while the numerator is, depending on the sublattice indices, given by

$\displaystyle P_{AA}(w)=zw^{R}-p_{0}(w)+p_{z}(w)\,,$ (78) $\displaystyle
P_{BB}(w)=zw^{R}-p_{0}(w)-p_{z}(w)\,,$ (79) $\displaystyle
P_{AB}(w)=p_{x}(w)-ip_{y}(w)\,,$ (80) $\displaystyle
P_{BA}(w)=p_{x}(w)+ip_{y}(w)\,.$ (81)

For $y<0$, the same expressions for both $P_{\alpha\beta}$ and $Q$ can be
used, replacing $p_{\alpha}$ by
$\overline{p}_{\alpha}(w)=w^{R}h_{\alpha}(w^{-1})$, which is the same
polynomial as $p_{\alpha}$ with the coefficients in reverse order. While $Q$
is not palindromic (its coefficients are complex in general), for real $z$ it
is self-reciprocal. This implies that its roots come in pairs of the form:
$w_{j}$, $1/w_{j}^{*}$. So, when $z$ is real and it lies outside the band
regions, there will be at most $R$ poles $\\{w_{\mathrm{in},j}\\}_{j=1}^{S}$,
$S\leq R$, inside the unit circumference that contribute to the integral.
Therefore, we can express

$\left[\Phi(z,y)\right]_{\alpha\beta}=\sum_{j}\left(w_{\mathrm{in},j}\right)^{\lvert
y\rvert+R-1}\operatorname{Res}\left(\frac{P_{\alpha\beta}}{Q},w_{\mathrm{in},j}\right)\,.$
(82)

As the hopping range increases, computing the poles and residues manually can
become challenging. However, this task can be carried out using well-
established numerical algorithms [33].

The self energy matrix appearing in the “pole equation”, Eq. (22), is given by
$\Sigma(\omega_{\mathrm{BS}})=g^{2}\Phi(\omega_{\mathrm{BS}},0)$.

##  Appendix D Strong light-matter coupling

When the light-matter interaction is the leading energy scale, $g\gg\Delta,J$,
it is convenient to work in a basis that diagonalizes $H_{I}$. So, we consider
the hybrid modes $\beta_{n,\pm}\equiv\left(\sigma^{-}_{n}\pm
a_{\bm{r}_{n},\alpha_{n}}\right)/\sqrt{2}$, with which we can write

$H_{I}=g\sum_{n}(\beta^{\dagger}_{n,+}\beta_{n,+}-\beta^{\dagger}_{n,-}\beta_{n,-})\,.$
(83)

$H_{S}$ and $H_{B}^{\bullet}$ now couple these hybrid modes between themselves
and with the rest of bath modes. For example, for the 1D system corresponding
to a line-defect in the Harper-Hofstadter model, we have

$\displaystyle H_{S}$
$\displaystyle=\frac{\Delta}{2}(\beta^{\dagger}_{+}+\beta^{\dagger}_{-})(\beta_{+}+\beta_{-})\,,$
(84)
$\displaystyle\begin{split}H_{B}^{\bullet}&=\frac{\mu_{1}}{2}(\beta^{\dagger}_{+}-\beta^{\dagger}_{-})(\beta_{+}-\beta_{-})\\\
&\quad-\frac{J}{\sqrt{2}}\sum_{j=0,2}\left[(\beta^{\dagger}_{+}-\beta^{\dagger}_{-})b_{j}+b^{\dagger}_{j}(\beta_{+}-\beta_{-})\right]\,.\end{split}$
(85)

Here we have assumed that the emitter is coupled to a site with chemical
potential $\mu_{1}$. The total Hamiltonian can be represented graphically as
shown in Fig. 12 below.

Figure 12: Schematic representation of the total Hamiltonian as a graph,
written in the basis $\beta_{+}$ (blue), $\beta_{-}$ (red) and
$\\{b_{n}\\}_{n\neq 1}$ (gray). The hopping strength is indicated next to each
link, while on-site energies are indicated above or below each site.

Using standard perturbative methods [34], we can obtain effective Hamiltonians
for these strongly-hybridized polariton modes. Within the single-particle
sector, the Hamiltonian can be written as

$H=\sum_{m,\,n}J_{mn}\sigma_{mn}\,,\quad\sigma_{mn}\equiv\ket{m}\\!\bra{n}\,,$
(86)

where $\\{\ket{n}\\}_{n\geq 0}$ is a basis of the single-particle subspace of
states. If $\lvert J_{00}-J_{nn}\rvert\gg\lvert J_{0n}\rvert$ for $n>0$, then,
a unitary transformation

$\tilde{H}\equiv e^{-S}He^{S}\simeq
H+[H,S]+\frac{1}{2!}\left[[H,S],S\right]+\cdots$ (87)

with

$S=\sum_{n}\frac{1}{J_{nn}-J_{00}}\left(J_{0n}\sigma_{0n}-J_{n0}\sigma_{n0}\right)\,,$
(88)

can remove all zeroth-order off-diagonal terms coupling the state
$\ket{\tilde{0}}$ with the rest of states in the new basis,
$\\{\ket{\tilde{n}}\equiv e^{-S}\ket{n}\\}_{n}$, leaving only terms
$\propto(J_{00}-J_{nn})^{-\alpha}$ with $\alpha\geq 1$. Thus, to first order
in perturbation series, the dispersion relation of the strongly-hybridized
polaritons appearing in the model shown in Fig. 12 can be computed as

$\begin{split}\omega_{\pm}&=\pm
g+\frac{\Delta+\mu_{1}}{2}\pm\frac{(\Delta-\mu_{1})^{2}}{8g}\\\
&\qquad+\sum_{\alpha=2,q}\frac{J^{2}}{\pm
2g+\Delta+\mu_{1}-2\mu_{\alpha}}\,.\end{split}$ (89)

While in the weak coupling limit the emitter dynamics can be modelled with an
effective single-band model, in the strongly-interacting regime the effective
emitter Hamiltonian has two bands, one for the “symmetric” and another for the
“antisymmetric” strongly-hybridized polaritons, that appear above and below
the bath’s energy bands, respectively. This symmetry with respect to the
exchange of the excitation between the emitters and the bath is an approximate
conserved quantity. From Eq. (89) we can see that, to leading order, the
dispersion of each kind of polariton has the usual cosine shape of a 1D tight-
binding model with nearest-neighbor hopping of strength $J/2$. As a
consequence, an initially localized excitation has the usual ballistic
diffusion along the emitter line with a single wavefront.

A similar analysis can be carried out for the Haldane model. In this model,
there are four polariton bands that appear above and below the bath’s bands.
Interestingly, while the minimal gap of the polaritons’ spectrum still changes
depending on the topology of the underlying bath, it remains more or less of
the same order. In addition, the average chiral group velocity is much smaller
than in the weak-coupling regime (relative to the effective hopping).

## References

  * González-Tudela and Cirac [2017a] A. González-Tudela and J. I. Cirac, Quantum emitters in two-dimensional structured reservoirs in the nonperturbative regime, Phys. Rev. Lett. 119, 143602 (2017a).
  * González-Tudela and Cirac [2017b] A. González-Tudela and J. I. Cirac, Markovian and non-markovian dynamics of quantum emitters coupled to two-dimensional structured reservoirs, Phys. Rev. A 96, 043811 (2017b).
  * Roos  _et al._ [2020] J. Roos, J. I. Cirac, and M. C. Bañuls, Markovianity of an emitter coupled to a structured spin-chain bath, Phys. Rev. A 101, 042114 (2020).
  * Bello  _et al._ [2019] M. Bello, G. Platero, J. I. Cirac, and A. González-Tudela, Unconventional quantum optics in topological waveguide qed, Science Advances 5, eaaw0297 (2019).
  * Leonforte  _et al._ [2021a] L. Leonforte, A. Carollo, and F. Ciccarello, Vacancy-like dressed states in topological waveguide qed, Phys. Rev. Lett. 126, 063601 (2021a).
  * De Bernardis  _et al._ [2021] D. De Bernardis, Z.-P. Cian, I. Carusotto, M. Hafezi, and P. Rabl, Light-matter interactions in synthetic magnetic fields: Landau-photon polaritons, Phys. Rev. Lett. 126, 103603 (2021).
  * Vega  _et al._ [2021] C. Vega, M. Bello, D. Porras, and A. González-Tudela, Qubit-photon bound states in topological waveguides with long-range hoppings, Phys. Rev. A 104, 053522 (2021).
  * Zhang  _et al._ [2022] H. Zhang, J. Li, H. Jiang, N. Li, J. Wang, J. Xu, C. Zhu, and Y. Yang, Coherent interaction of a quantum emitter and the edge states in two-dimensional optical topological insulators, Phys. Rev. A 105, 053703 (2022).
  * Vega  _et al._ [2022] C. Vega, D. Porras, and A. González-Tudela, Topological multi-mode waveguide qed (2022).
  * Roccati  _et al._ [2022] F. Roccati, S. Lorenzo, G. Calajò, G. M. Palma, A. Carollo, and F. Ciccarello, Exotic interactions mediated by a non-hermitian photonic bath, Optica 9, 565 (2022).
  * Gong  _et al._ [2022a] Z. Gong, M. Bello, D. Malz, and F. K. Kunst, Bound states and photon emission in non-hermitian nanophotonics (2022a).
  * Gong  _et al._ [2022b] Z. Gong, M. Bello, D. Malz, and F. K. Kunst, Anomalous behaviors of quantum emitters in non-hermitian baths (2022b).
  * Barik  _et al._ [2018] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, A topological quantum optics interface, Science 359, 666 (2018).
  * Kim  _et al._ [2021] E. Kim, X. Zhang, V. S. Ferreira, J. Banker, J. K. Iverson, A. Sipahigil, M. Bello, A. González-Tudela, M. Mirhosseini, and O. Painter, Quantum electrodynamics in a topological waveguide, Phys. Rev. X 11, 011015 (2021).
  * Hofstadter [1976] D. R. Hofstadter, Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * Haldane [1988] F. D. M. Haldane, Model for a quantum hall effect without landau levels: Condensed-matter realization of the "parity anomaly", Phys. Rev. Lett. 61, 2015 (1988).
  * Cohen-Tannoudji  _et al._ [1992] C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, and P. Thickstun,  _Atom-Photon Interactions: Basic Processes and Applications_, A Wiley-Interscience publication (Wiley, 1992).
  * Leonforte  _et al._ [2021b] L. Leonforte, D. Valenti, B. Spagnolo, A. Carollo, and F. Ciccarello, Dressed emitters as impurities, Nanophotonics 10, 4251 (2021b).
  * González-Tudela and Cirac [2018] A. González-Tudela and J. I. Cirac, Exotic quantum dynamics and purely long-range coherent interactions in dirac conelike baths, Phys. Rev. A 97, 043831 (2018).
  * Hatsugai [1993] Y. Hatsugai, Chern number and edge states in the integer quantum hall effect, Phys. Rev. Lett. 71, 3697 (1993).
  * Boykin  _et al._ [2010] T. B. Boykin, M. Luisier, and G. Klimeck, Current density and continuity in discretized models, European Journal of Physics 31, 1077 (2010).
  * Matsuki  _et al._ [2021] Y. Matsuki, K. Ikeda, and M. Koshino, Fractal defect states in the hofstadter butterfly, Phys. Rev. B 104, 035305 (2021).
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I_ (Springer-Verlag, New York, 1999).
  * Hügel and Paredes [2014] D. Hügel and B. Paredes, Chiral ladders and the edges of quantum hall insulators, Phys. Rev. A 89, 023619 (2014).
  * Greentree  _et al._ [2006] A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. Hollenberg, Quantum phase transitions of light, Nature Physics 2, 856 (2006).
  * Hayward  _et al._ [2012] A. L. C. Hayward, A. M. Martin, and A. D. Greentree, Fractional quantum hall physics in jaynes-cummings-hubbard lattices, Phys. Rev. Lett. 108, 223602 (2012).
  * Hartmann  _et al._ [2007] M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, Effective spin systems in coupled microcavities, Phys. Rev. Lett. 99, 160501 (2007).
  * Hung  _et al._ [2016] C.-L. Hung, A. González-Tudela, J. I. Cirac, and H. J. Kimble, Quantum spin dynamics with pairwise-tunable, long-range interactions, Proceedings of the National Academy of Sciences 113, E4946 (2016).
  * Manzoni  _et al._ [2017] M. T. Manzoni, L. Mathey, and D. E. Chang, Designing exotic many-body states of atomic spin and motion in photonic crystals, Nature Communications 8, 14696 (2017).
  * Bello  _et al._ [2022] M. Bello, G. Platero, and A. González-Tudela, Spin many-body phases in standard- and topological-waveguide qed simulators, PRX Quantum 3, 010336 (2022).
  * Avron  _et al._ [1983] J. E. Avron, R. Seiler, and B. Simon, Homotopy and quantization in condensed matter physics, Phys. Rev. Lett. 51, 51 (1983).
  * Mitchell  _et al._ [2018] N. P. Mitchell, L. M. Nash, D. Hexner, A. M. Turner, and W. T. M. Irvine, Amorphous topological insulators constructed from random point sets, Nature Physics 14, 380 (2018).
  * Mahoney and Sivazlian [1983] J. Mahoney and B. Sivazlian, Partial fractions expansion: a review of computational methodology and efficiency, Journal of Computational and Applied Mathematics 9, 247 (1983).
  * Bravyi  _et al._ [2011] S. Bravyi, D. P. DiVincenzo, and D. Loss, Schrieffer–wolff transformation for quantum many-body systems, Annals of Physics 326, 2793 (2011).