# Topological multi-mode waveguide QED

C. Vega  Institute of Fundamental Physics IFF-CSIC, Calle Serrano 113b, 28006
Madrid, Spain.  D. Porras  Institute of Fundamental Physics IFF-CSIC, Calle
Serrano 113b, 28006 Madrid, Spain.  A. González-Tudela  Institute of
Fundamental Physics IFF-CSIC, Calle Serrano 113b, 28006 Madrid, Spain.

###### Abstract

Two-dimensional topological insulators with quantized large Chern numbers
feature a number of protected, chiral edge modes linked to the value of the
invariant. When implemented in photonic setups, these systems then naturally
display a topologically-protected multi-mode waveguide at their edges. Here,
we show how to take advantage of these setups by interfacing them with quantum
emitters. Using a Harper-Hofstadter lattice as a particular model, we find
situations in which the emitters feature quasi-quantized decay rates due to
the increasing number of edge modes in the different emergent band-gaps, and
where their spontaneous emission spatially separates in different modes,
enabling the generation of non-local time-bin entangled states already with a
single $\pi$-pulse. Besides, we show how the emitters can selectively interact
with the different channels using non-local light-matter couplings as the ones
that can be obtained with giant atoms. Such capabilities pave the way for
generating quantum gates among topologically-protected photons as well as
generating more complex entangled states of light in topological channels.

Topological photonics [1, 2, 3, 4] is a burgeoning field aiming at exporting
topological concepts into photonics to bring novel, and more robust, ways of
controlling the properties of light. At the classical level, one of the
original motivations of the field was to exploit the chiral, edge modes
appearing in two-dimensional topological insulators to obtain unidirectional
and robust-to-disorder photon flows [5, 6]. In their simplest instance, that
is, when only a single edge mode appears, it was soon realized in several
platforms [7, 8, 9]. Remarkably, these realizations opened up applications
beyond their initial motivation, such as the realization of “topological
lasers”[10, 11, 12, 13] and chiral light-matter couplings [14, 15, 16, 17,
18], or the generation of gaussian quantum correlations with parametric
drivings [19, 20, 21, 22, 23, 24]. Besides, they have been proposed as robust
_quantum buses_ to mediate quantum state transfer [25, 26, 27, 28], without
the exponential time dependence on their distance of their zero-dimensional
counterparts [29, 30, 31, 32, 33, 34, 35].

The more complex scenario when the two-dimensional system features several
edge modes, as proposed and observed in Refs. [36, 37], opens up new
opportunities, e.g., to increase quantum communication capacity using several
topologically-protected channels. Remarkably, this situation has been much
more scarcely explored in the literature, likely because it is still unclear
how to profit from these extra modes. Coupling quantum emitters to such
structures provide a natural route to harness these additional degrees of
freedom. On the one hand, local light-matter coupling enables the coupling to
all photonic modes at the emitter frequencies, thus being able to interact
with several modes simultaneously. This has already been shown to lead to
unconventional emitter-emitter interactions when coupled to the bulk modes of
one-[38, 39, 40], two-[39, 41], and three-[42, 43] dimensional topological
photonic insulators. Besides, the strongly interacting character of the
emitters can induce (non-gaussian) quantum correlations beyond the ones that
can be obtained with parametric drivings [19, 20, 21, 22, 23, 24], opening a
path to observe exotic quantum many-body states [44, 45, 46, 47, 48, 49]. All
these reasons are motivating the development of such topological light-matter
interfaces in various platforms, ranging from superconducting qubits coupled
to microwave resonators [50, 18] to solid-state emitters coupled to
topological photonic crystals [14, 15, 16, 17].

Figure 1: Scheme of the topological light-matter interface: Each blue globe
represents a lattice site, with associated annihilation (creation) operators
of local excitations $a_{x,y}^{(\dagger)}$. The yellow-shaded region denotes a
single lattice plaquette. We also represent the spontaneous emission of an
emitter (in red) coupled to the edge of the lattice. The excitation is
radiated through all edge modes ($2$ in the figure), that propagate at
different group velocities (denoted by $v_{g,0}$ and $v_{g,1}$), which will
cause a spatial separation of the emitted pulse.

In this Letter, we develop a theory for the topological multi-mode waveguide
QED scenario that appears when quantum emitters couple to the edges of
topological photonic insulators with large Chern numbers, and show how to
exploit it for the generation of non-trivial states of light. Motivated by
recent experiments [18], we particularize for a photonic lattice with an
effective magnetic flux–the so-called Harper-Hofstadter (HH) lattice [51]–,
where bands with large Chern number appear for small magnetic fluxes [52, 53].
We derive an effective theory for the edge modes, and study the most relevant
quantum optical observables, unveiling, for example, the quasi-quantization of
the emitter’s decay rate at the different topological band-gaps. Additionally,
we propose a protocol to generate time-bin-like entanglement between the
different topologically protected channels by combining a $\pi$-pulse driving
of the emitters, and the spontaneous spatial separation of the modes that
occur in these systems. Finally, we also explain how to make the emitters
selectively interact with the different channels through the use of the non-
local couplings [54, 55, 56, 57, 58, 59].

_System_. The model that we consider along this manuscript is depicted in Fig.
1(a): a quantum emitter interacts locally with one of the edges of a HH
lattice [51, 52, 53], that is, a square photonic lattice with nearest-
neighbour hoppings of rate $J$, subject to an effective magnetic flux, $\phi$.
The bath Hamiltonian then reads (we set $\hbar=1$ along this manuscript):

$H_{B}=-J\left(\sum_{x,y}a_{x+1,y}^{\dagger}a_{x,y}+e^{-2\pi i\phi
x}a_{x,y+1}^{\dagger}a_{x,y}\right)+\text{H.c.}\;,$ (1)

where $a^{(\dagger)}_{x,y}$ represent the annihilation (creation) operator at
the $(x,y)$ position, and where we take the cavity energy, $\omega_{a}=0$, as
the energy reference of the problem. The emitter is assumed to have a single
optical transition between its ground ($g$) and excited state ($e)$ with
frequency $\omega_{e}$, that couples to one of the photonic lattice sites at
the edges through the standard light-matter Hamiltonian,
$H_{I}=(ga_{x_{e},y_{e}}\sigma_{eg}+\mathrm{H.c.})$, with $g$ being its
coupling strength, $(x_{e},y_{e})$ the position where it couples, and using
the notation $\sigma_{\alpha\beta}=\ket{\alpha}\bra{\beta}$ for the emitter’s
operators. The emitter’s Hamiltonian then reads,
$H_{S}=\omega_{e}\sigma_{ee}$, such that the full topological light-matter
interface is described by a Hamiltonian comprising the three terms,
$H=H_{S}+H_{B}+H_{I}$.

Figure 2: (a) HH spectrum with periodic boundary conditions in the Y direction
- open for X - for $\phi=1/9$ and a lattice size of $65\times 65$ sites. Each
dot is coloured according to the localization parameter $\eta$ defined in the
main text: $\eta=0$ corresponds to a delocalized state, while $\eta=(-)1$
depicts complete localization at the (left) right boundary. (b) Emitter
spontaneous emission rate as a function of its transition frequency
$\omega_{e}$, computed from Eq. (3) approximating $\delta(\omega_{e}-E)$ by a
Gaussian function with mean $\omega_{e}-E_{B}$ and width $\theta=0.07J$ (see
Ref. [60] for details). The emitter is coupled to the edge of a HH lattice of
size $150\times 150$ with flux $\phi=1/25$. Yellow vertical fringes are
centered at Landau levels and have a width equal to $\theta$.

_Photonic bath as topological multi-mode waveguide: effective description._
Despite its apparent simplicity, the HH model displays a very rich behaviour
depending on the value of $\phi$ [51, 52, 53]. For concreteness, along this
work we take $\phi=1/q$, with $q\in\mathbb{N}$, that will be enough to
illustrate the behaviour we are interested in. With this parametrization, the
spectrum of the system with periodic boundary conditions features $q$ bands,
labeled as Landau-levels [52], separated by $l=0,\dots,q-1$ band-gaps.
Besides, it can be shown that the first $\left[(q-1)/2\right]$-bands (being
$\left[\cdot\right]$ the floor function) have an associated quantized Chern
number $C=-1$ (see Ref. [52] and Sup material [60] for more details). Thus,
with open boundary conditions, it is expected that the $l$-th band-gap
features $2(l+1)$ gapless edge states, associated to the $l+1$ Landau-levels
below that energy 111As shown in Sup. Material [60], this situation is
qualitatively equivalent to the case where the multiple edge states emerge
from a single band with $|C|>1$ [36, 37]..

To illustrate this emergence of multiple edge modes in the different band-
gaps, we consider a cylinder geometry for the bath with periodic (open)
boundary conditions in the Y (X) direction. This allows one to write the
spatial wavefunction of the bath eigen-states as
$\psi_{l}(x,y)=e^{ik_{y}y}\Psi_{l}(x)$, and calculate their eigen-energies
$\omega(k_{y})$ numerically. In Fig. 2(a), we plot an example of the bath
spectrum for a bath with $\phi=1/9$ and $65\times 65$ sites, showing the
emergence of the gapless modes between the bulk flat bands. Besides, we define
a localization parameter
$\eta=\sum_{x=0}^{L-1}\left(-1+2\frac{x}{L-1}\right)|\Psi(x)|^{2}$, that
ranges $\eta\in[-1,1]$, achieving the extremes $(-)1$ when the modes are
maximally localized at the left (right) edge, encoded in purple (yellow)
color, respectively, in Fig. 2(a). Like this, one can see how, for energies
below $\omega_{a}$, the edge modes at the left (right) have always negative
(positive) group velocity along Y, thus, having a chiral character. Thus, when
a coupling an emitter to one of the edges, it will only interact with the
modes of certain chirality.

An important observation from Fig. 2(a) is that the edge modes dispersion
deviate significantly from the linear behaviour which is typically assumed in
the literature to deal with such topological channels [25, 26, 27]. Since this
can have important consequences in the quantum optical behaviour, we derive a
more accurate effective description of these modes for $\phi\ll 1$. In
particular, we show that the edge-mode dispersion for the left-localized
states emerging from the $l$-th Landau Level is well approximated by [60]:

$\displaystyle\omega_{l}(k_{y})$
$\displaystyle\approx\omega_{\mathrm{LL}}(\phi)+a_{l}(\phi)(k_{y}-k_{l}(\phi))^{2}\,,$
(2)

for $(k_{y}-k_{l}(\phi))\in(-\pi,0)$, where
$\omega_{\mathrm{LL}}(\phi)/J\approx-4+4\pi\phi\left(l+\frac{1}{2}\right)-\left(\pi\phi\right)^{2}(l^{2}+l+1/2)$
is the energy of the $l^{\text{th}}$-Landau level [52], $a_{l}(\phi)$ is the
effective curvature of the edge modes which we extract from numerical fittings
[60], and which converges to $a_{l}(\phi\rightarrow 0)\sim 0.6$, and
$k_{l}(\phi)$ is the momentum resonant to the minimum edge mode energy.
Besides, the spatial distribution of these modes along the $X$ direction is
simply given by
$\Psi_{l}(x)=\sqrt{2/\lambda_{l}(\phi)}e^{-x/\lambda_{l}(\phi)}$, where
$\lambda_{l}$ grows as $\phi\rightarrow 0$ as expected, since in $\phi=0$ we
should recover the delocalized bath eigenstates of the standard square lattice
model. Furthermore, for a fixed energy, resonant edge modes with larger $l$
tend to be less localized (see Ref. [60] for more details).

Figure 3: Spontaneous emission features of a quantum emitter coupled to the
middle site of the left boundary of a HH lattice for $\phi=1/9$, with quantum
emitter frequencies of $\omega_{e}/J=-3.31$ (first row), $\omega_{e}/J=-2.16$
(second row) and $\omega_{e}=-1.25$ (last row), resonant to $1$, $2$ and $3$
edge modes respectively. (a, c, e) Snapshots at $TJ=200$ of the bath
population in real-space for a lattice size of $50\times 50$ and a coupling
constant of $g/J=0.1$, in the presence of a defect, depicted as a white-
colored region in the lowest part of the lattice. (b, d, f) Dynamics of the
emitted photon. For each $\omega_{e}$, we plot the evolution of the population
of the left boundary sites $|A(0,Y;t)|^{2}$, in arbitrary units. At the top,
we include a snapshot of the pulse shape at the final time instant, showing
$1$, $2$ and $3$ peaks respectively, which correspond to the number of
resonant edge modes in each case.

_Spontaneous emission features in topological multi-mode waveguide QED_ After
having characterized the multi-mode waveguide emerging at the HH lattice
edges, let us now see how coupling emitters to it leads to several unique
phenomena. First, let us note that since the emitter probes the system at
fixed frequency, $\omega_{e}$, one can control the number of modes that will
be relevant for its dynamics just by adjusting its relative detuning with the
bath energies. A magnitude that evidences that control is the Markovian decay
rate, $\Gamma(\omega_{e})$, defined by [62]:

$\displaystyle\Gamma(\omega_{e})=\sum_{E_{B}}\left|\bra{e}H_{I}\ket{E_{B}}\right|^{2}\cdot\delta(\omega_{e}-E_{B})\;,$
(3)

being $\ket{E_{B}}$ are the bath eigenstates, i.e.,
$H_{B}\ket{E_{B}}=E_{B}\ket{E_{B}}$, for the considered configuration. In Fig.
2(b), we plot in blue solid line the expected $\Gamma(\omega_{e})$ for an
emitter coupled to the edge of a HH lattice for $\phi=1/25$ as a function of
$\omega_{e}$. There, we see how the expected decay rate abruptly increases
from one band-gap to the other as the emitter’s energies is varied. The jumps
occur when the emitter’s energy $\omega_{e}$ starts crossing the Landau level
energies, indicated in shaded yellow region in the figure, due to the
emergence of another edge mode that couples to the emitter. Note also that the
decay rate remains almost constant along the whole band-gap region, except for
a deviation that occurs due to the non-linear energy dispersion of the modes,
that becomes more prominent for the larger band-gaps. As shown in Ref. [60],
this quasi-quantized behaviour is well captured by our effective model, which
gives a semi-analytical approximation for the decay rates into the different
topological channels:

$\displaystyle\Gamma_{l}(\omega_{e})\approx\frac{g^{2}|\Psi_{l}(0)|^{2}}{|v_{g,l}(\omega_{e})|}=\frac{g^{2}}{\lambda_{l}(\phi)|a_{l}(k_{e}^{(l)}-k_{l}(\phi))|}$
(4)

being $k_{e}^{(l)}$ is the momentum resonant with the emitter’s energy, i.e.,
$\omega_{l}(k_{e}^{(l)})=\omega_{e}$. In the second equality, we explicitly
calculate the group velocity
$v_{g,l}(\omega_{e})=\partial_{k_{y}}\omega_{l}(k_{y})$ using Eq. (2). Then,
the total decay can be obtained by summing the contributions of the active
channels, i.e., $\Gamma=\sum_{l}\Gamma_{l}$, which explains why when
$\omega_{e}$ crosses a Landau level, the decay rate of an emitter abruptly
increases.

Beyond this quasi-quantization of the emitter’s decay rate, a more interesting
aspect is what occurs with the spontaneously emitted photons when the emitters
are driven. Let us first assume a perfect $\pi$-pulse driving which prepare
the system in the state $\ket{\Psi_{0}}=\ket{e}\otimes\ket{\mathrm{vac}}$,
with $\ket{\mathrm{vac}}$ being the bath state with no photons. When the laser
is switched off, the whole system evolves according to the total Hamiltonian
$\ket{\Psi(t)}=e^{-iHt}\ket{\Psi_{0}}$, leading eventually to a single-photon
wavepacket state

$\ket{\Psi(t\rightarrow\infty)}=\ket{g}\otimes\sum_{x,y}A(x,y)a_{x,y}^{\dagger}\ket{\mathrm{vac}}\;,$
(5)

as it occurs in other quantum optical setups [63]. However, in this case such
wavepackets have several unique features which we illustrate in Fig. 3 for an
emitter coupled to the edge of a HH lattice with $\phi=1/9$. First,
irrespective of the band-gap that the emitters are resonant to, the photons
are emitted in a chiral and robust fashion due to their topological origin.
This is illustrated in Figs. 3(a, c, e), where we plot the full photonic bath
population, $|A(x,y;t)|^{2}$, at a time $TJ=200$. There, we observe that,
irrespective of the band-gap chosen, the single-photon wavepacket can overcome
the defect introduced in one of the edges without altering significantly its
propagation. Besides, we also observe how the situations with more than one
edge-state, Figs. 3(c, e), display a localized, but more complex,
wavefunction. To appreciate better the inner structure of these wavepackets,
we plot in Figs. 3(b, d, f) their temporal dynamics focusing only at the edge
population, $|A(0,Y;t)|^{2}$. Like this, we observe a unique effect of these
multi-mode waveguides, that is, that after certain time, the emission into the
different topological channels becomes spatially separated due to the
different group velocities of the modes at the emitter frequency. Intuitively,
this separation starts to occur for times such that
$|v_{g,l}-v_{g,l^{\prime}}|T\gtrapprox\Gamma^{-1}_{l}+\Gamma^{-1}_{l^{\prime}}$,
that is, that the separation between the wavepackets is larger than their
intrinsic linewidth ($\Gamma^{-1}_{l^{\prime}}$), as given by Eq. (4). When
this separation occurs, one can say that our emitter has generated single-
particle entanglement [64] between orthogonal time-bins $T_{l}$. Defining
$\ket{1}_{l}$ as the presence of a photon in the $T_{l}$ time-bin and $0$ in
the rest, the photonic state created in the asymptotic limit can be written as
$|\Psi(t\to\infty)\rangle_{\mathrm{ph}}\approx\sum_{l}c_{l}\ket{1}_{l}$. Note
that in more conventional quantum optical setups [65, 66, 67, 68, 69, 70, 71,
72, 73, 74, 75, 76, 77, 78, 79, 80, 81] where such time-bin entanglement is
generated, it is required to combine superpositions in multi-level emitters
and multiple-drivings, while here already with a single $\pi$-pulse, a
$W$-type [82] entangled structure appears due to the multi-mode nature of the
waveguide. We foresee that combining these topological multi-mode waveguide
QED setups, with more complicated emitters, e.g., with an anharmonic level
structure like transmons [76, 78, 77, 81], one can generate more complex
entangled states of light, with the added value of their topological
protection [83, 84, 85, 21, 23, 86, 24, 87].

Finally, let us note that one can make the quantum emitters interact
selectively with one of the resonant channels [60]. The key idea is to couple
the emitter with more than a single cavity mode, as it can be achieved with
giant atoms [54, 55, 56, 57, 58]. Let us illustrate it in the simplest case
where we want to cancel only one resonant momenta. This requires that the
emitter couples to two adjacent cavities with the same strength, and relative
phase $e^{i\varphi}$ i.e. $g_{(0,y)}=g$ and $g_{(0,y+1)}=ge^{i\varphi}$. In
that case, the k-dependent light-matter coupling reads

$|G_{\varphi}(k_{y})|^{2}\propto\left(1+\cos(k_{y}+\varphi)\right)\;,$ (6)

and thus, vanishes at $k_{y}=k_{e}$ if we choose $\varphi=\pi-k_{e}$. This was
the key idea introduced in simple waveguide QED setups [88, 89, 90] to obtain
chiral emission. Here, the emission is already chiral, but we can still use it
to cancel the emission into the resonant momenta of the undesired resonant
channels. A proof-of-principle demonstration of this effect is shown in Fig.
4, where we compare the photonic profiles of the spontaneous emission for an
emitter with local coupling vs. tailored non-local ones to suppress emission
into one of the two resonant channels. In the Sup. Material [60], we also
prove that the number of non-local couplings required to cancel $N_{k}$
resonant momenta scales only linearly with $N_{k}$. This mode selectivity is
an interesting tool in this scenario because when a photon propagates in a
chiral channel and interacts with an emitter, it acquires a $\pi$-phase [91].
Thus, if multiple photons are sent in the different channels, such mode
selectivity can lead to different phases between the topologically-protected
photons.

Figure 4: Snapshot of photonic population after the decay of an emitter
coupled to the left edge of a HH lattice of $600\times 600$ sites for
$\phi=1/9$, $g/J=0.2$ and $\omega_{e}/J=-1.86$, resonant to $2$ modes. We
compare the emission for of a locally-coupled emitter (orange markers) vs. a
giant atom with tailored non-local couplings, chosen to suppress $k_{e}^{(0)}$
(purple solid line) and $k_{e}^{(1)}$ (blue dashed line).

_Conclusions & outlook._ Summing up, we have characterized the topological
multi-mode waveguide QED setup that appears when quantum emitters couple to
the edges of a Harper-Hofstadter lattice. We find several unique features such
as the quasi-quantization of the decay rates and the _spontaneous_ generation
of entanglement in the different topological channels, as well as a way to
make the emitters interact selectively with some of the channels. We foresee
that the combination of these setups with multi-level emitters and/or complex
time-dependent [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
80, 81] or parametric [19, 20, 21, 22, 23, 24] drivings can be used to
generate more complex states-of-light in these topologically-protected
channels either in a transient [83, 84, 85, 21, 23, 86, 24, 87] or a steady-
state fashion [92, 93, 94, 95, 96, 97, 98, 99, 100, 101], as well as to induce
gates between topologically-protected photons. Although we find our results
for the Harper-Hofstadter model, we expect our findings can be of interest to
other system where such multiple edge states appear [36, 37, 102, 103, 104].

###### Acknowledgements.

The authors acknowledge support from the Proyecto Sinérgico CAM 2020
Y2020/TCS-6545 (NanoQuCo-CM), the CSIC Research Platform on Quantum
Technologies PTI-001 and from Spanish project
PGC2018-094792-B-100(MCIU/AEI/FEDER, EU). During the writing process of this
article, another work studying the coupling of emitters to the edges of two-
dimensional photonic insulators appeared [105], although focused only in the
single-edge mode scenario.

## References

  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Rev. Mod. Phys. 91, 15006 (2019).
  * Rider  _et al._ [2019] M. S. Rider, S. J. Palmer, S. R. Pocock, X. Xiao, P. Arroyo Huidobro, and V. Giannini, Journal of Applied Physics 125, 120901 (2019).
  * Smirnova  _et al._ [2020] D. Smirnova, D. Leykam, Y. Chong, and Y. Kivshar, Appl. Phys. Rev. 7, 021306 (2020).
  * Price  _et al._ [2022] H. Price, Y. Chong, A. Khanikaev, H. Schomerus, L. J. Maczewsky, M. Kremer, M. Heinrich, A. Szameit, O. Zilberberg, Y. Yang, B. Zhang, A. Alù, R. Thomale, I. Carusotto, P. St-Jean, A. Amo, A. Dutt, L. Yuan, S. Fan, X. Yin, C. Peng, T. Ozawa, and A. Blanco-Redondo, arXiv:2201.06315 10.1088/2515-7647/ac4ee4 (2022).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, Phys. Rev. Lett. 100, 13904 (2008).
  * Raghu  _et al._ [2008] S. Raghu, X. L. Qi, C. Honerkamp, and S. C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Nature 461, 772 (2009).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, Y. Plotnik, J. M. Zeuner, D. Song, Z. Chen, A. Szameit, and M. Segev, Phys. Rev. Lett. 111, 103901 (2013).
  * Hafezi  _et al._ [2013] M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Nature Photonics 7, 1001 (2013).
  * Bahari  _et al._ [2017] B. Bahari, A. Ndao, F. Vallini, A. El Amili, Y. Fainman, and B. Kanté, Science 358, 636 (2017).
  * Bandres  _et al._ [2018] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Science 359, aar4005 (2018).
  * Klembt  _et al._ [2018] S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. Liew, M. Segev, C. Schneider, and S. Höfling, Nature 562, 552 (2018).
  * Amelio and Carusotto [2020] I. Amelio and I. Carusotto, Phys. Rev. X 10, 041060 (2020).
  * Barik  _et al._ [2018] S. Barik, A. Karasahin, C. Flower, T. Cai, H. Miyake, W. DeGottardi, M. Hafezi, and E. Waks, Science 359, 666 (2018).
  * Barik  _et al._ [2020] S. Barik, A. Karasahin, S. Mittal, E. Waks, and M. Hafezi, Phys. Rev. B 101, 205303 (2020).
  * Hallett  _et al._ [2022] D. Hallett, A. P. Foster, D. Whittaker, M. S. Skolnick, and L. R. Wilson, ACS Photonics 9, 706 (2022).
  * Mehrabad  _et al._ [2020] M. J. Mehrabad, A. P. Foster, R. Dost, E. Clarke, P. K. Patil, A. M. Fox, M. S. Skolnick, and L. R. Wilson, Optica 7, 1690 (2020).
  * Owens  _et al._ [2021] J. C. Owens, M. G. Panetta, B. Saxberg, G. Roberts, S. Chakram, R. Ma, A. Vrajitoarea, J. Simon, and D. Schuster, arXiv: 2109.06033 (2021).
  * Mittal  _et al._ [2016a] S. Mittal, S. Ganeshan, J. Fan, A. Vaezi, and M. Hafezi, Nat. Photon. 10, 180 (2016a).
  * Mittal  _et al._ [2021] S. Mittal, V. V. Orre, E. A. Goldschmidt, and M. Hafezi, Nat. Photon. 15, 542 (2021).
  * Blanco-Redondo  _et al._ [2018] A. Blanco-Redondo, B. Bell, D. Oren, B. J. Eggleton, and M. Segev, Science 362, 568 (2018).
  * Wang  _et al._ [2019] H. Wang, J. Qin, X. Ding, M. C. Chen, S. Chen, X. You, Y. M. He, X. Jiang, L. You, Z. Wang, C. Schneider, J. J. Renema, S. Höfling, C. Y. Lu, and J. W. Pan, Phys. Rev. Lett. 123, 250503 (2019).
  * Blanco-Redondo [2020] A. Blanco-Redondo, Proceedings of the IEEE 108, 837 (2020).
  * Doyle  _et al._ [2022] C. Doyle, W.-W. Zhang, M. Wang, B. A. Bell, S. D. Bartlett, and A. Blanco-Redondo, arXiv: 2202.01873 (2022).
  * Yao  _et al._ [2013] N. Y. Yao, C. R. Laumann, A. V. Gorshkov, H. Weimer, L. Jiang, J. I. Cirac, P. Zoller, and M. D. Lukin, Nat. Commun. 4, 1 (2013).
  * Lemonde  _et al._ [2019] M. A. Lemonde, V. Peano, P. Rabl, and D. G. Angelakis, New Journal of Physics 21, 113030 (2019).
  * Cao  _et al._ [2021] J. Cao, W. X. Cui, X. Yi, and H. F. Wang, Annalen der Physik 533, 2100120 (2021).
  * Dlaska  _et al._ [2017] C. Dlaska, B. Vermersch, and P. Zoller, Quantum Science and Technology 2, 015001 (2017).
  * Almeida  _et al._ [2016] G. M. Almeida, F. Ciccarello, T. J. Apollaro, and A. M. Souza, Phys. Rev. A 93, 032310 (2016).
  * Mei  _et al._ [2018] F. Mei, G. Chen, L. Tian, S. L. Zhu, and S. Jia, Phys. Rev. A 98, 012331 (2018).
  * Lang and Büchler [2017] N. Lang and H. P. Büchler, npj Quantum Information 3, 1 (2017).
  * Longhi  _et al._ [2019] S. Longhi, G. L. Giorgi, and R. Zambrini, Advanced Quantum Technologies 2, 1800090 (2019).
  * Tran  _et al._ [2020] M. C. Tran, A. Deshpande, A. Y. Guo, A. Lucas, and A. V. Gorshkov, arXiv:2010.02930 (2020).
  * Wang  _et al._ [2021a] C. Wang, L. Li, J. Gong, and Y.-x. Liu, arXiv: 2111.05561 (2021a).
  * Han  _et al._ [2021] J. X. Han, J. L. Wu, Y. Wang, Y. Xia, Y. Y. Jiang, and J. Song, Phys. Rev. A 103, 032402 (2021).
  * Skirlo  _et al._ [2014] S. A. Skirlo, L. Lu, and M. Soljačić, Phys. Rev. Lett. 113, 113904 (2014).
  * Skirlo  _et al._ [2015] S. A. Skirlo, L. Lu, Y. Igarashi, Q. Yan, J. Joannopoulos, and M. Soljačić, Phys. Rev. Lett. 115, 253901 (2015).
  * Bello  _et al._ [2019] M. Bello, G. Platero, J. I. Cirac, and A. González-Tudela, Science Advances 5, eaaw0297 (2019).
  * Leonforte  _et al._ [2021] L. Leonforte, A. Carollo, and F. Ciccarello, Phys. Rev. Lett. 126, 063601 (2021).
  * Vega  _et al._ [2021] C. Vega, M. Bello, D. Porras, and A. González-Tudela, Phys. Rev. A 104, 053522 (2021).
  * De Bernardis  _et al._ [2021] D. De Bernardis, Z.-P. Cian, I. Carusotto, M. Hafezi, and P. Rabl, Phys. Rev. Lett. 126, 103603 (2021).
  * García-Elcano  _et al._ [2020] I. n. García-Elcano, A. González-Tudela, and J. Bravo-Abad, Phys. Rev. Lett. 125, 163602 (2020).
  * García-Elcano  _et al._ [2021] I. n. García-Elcano, J. Bravo-Abad, and A. González-Tudela, Phys. Rev. A 103, 033511 (2021).
  * Bello  _et al._ [2022] M. Bello, G. Platero, and A. González-Tudela, Phys. Rev. X Quantum 3, 010336 (2022).
  * Grusdt  _et al._ [2013] F. Grusdt, M. Höning, and M. Fleischhauer, Physical Review Letters 110, 260405 (2013).
  * Grusdt  _et al._ [2014] F. Grusdt, F. Letscher, M. Hafezi, and M. Fleischhauer, Phys. Rev. Lett. 113, 155301 (2014).
  * Maghrebi  _et al._ [2015] M. F. Maghrebi, N. Y. Yao, M. Hafezi, T. Pohl, O. Firstenberg, and A. V. Gorshkov, Phys. Rev. A 91, 033838 (2015).
  * Schine  _et al._ [2016] N. Schine, A. Ryou, A. Gromov, A. Sommer, and J. Simon, Nature 534, 671 (2016).
  * Clark  _et al._ [2020] L. W. Clark, N. Schine, C. Baum, N. Jia, and J. Simon, Nature 582, 41 (2020).
  * 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, Phys. Rev. X 11, 011015 (2021).
  * Hofstadter [1976] D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * Harper  _et al._ [2014] F. Harper, S. H. Simon, and R. Roy, Phys. Rev. B 90, 075104 (2014).
  * Goldman [2009] N. Goldman, Journal of Physics B: Atomic, Molecular and Optical Physics 42, 055302 (2009).
  * Kockum  _et al._ [2018] A. F. Kockum, G. Johansson, and F. Nori, Phys. Rev. Lett. 120, 140404 (2018).
  * González-Tudela  _et al._ [2019] A. González-Tudela, C. Sánchez Muñoz, and J. Cirac, Phys. Rev. Lett. 122, 10.1103/PhysRevLett.122.203603 (2019).
  * Frisk Kockum [2021] A. Frisk Kockum (Springer, Singapore, 2021) pp. 125–146.
  * Kannan  _et al._ [2020a] B. Kannan, M. J. Ruckriegel, D. L. Campbell, A. Frisk Kockum, J. Braumüller, D. K. Kim, M. Kjaergaard, P. Krantz, A. Melville, B. M. Niedzielski, A. Vepsäläinen, R. Winik, J. L. Yoder, F. Nori, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Nature 583, 775 (2020a).
  * Wang  _et al._ [2021b] X. Wang, T. Liu, A. F. Kockum, H. R. Li, and F. Nori, Phys. Rev. Lett. 126, 043602 (2021b).
  * Giorgi  _et al._ [2019] G. L. Giorgi, S. Longhi, A. Cabot, and R. Zambrini, Annalen der Physik 531, 1900307 (2019).
  * [60] In this Supplementary Material, we provide more details on: i) Numerical and effective mode description of the topological multi-mode channels; ii) Calculation of the local density of states; iii) Characterization of multi-mode entanglement.
  * Note [1] As shown in Sup. Material [60], this situation is qualitatively equivalent to the case where the multiple edge states emerge from a single band with $|C|>1$ [36, 37].
  * Cohen-Tannoudji  _et al._ [1998] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Atom—Photon Interactions_ (1998).
  * Law and Kimble [1997] C. K. Law and H. J. Kimble, Journal of Modern Optics 44-11, 2067 (1997).
  * Van Enk [2005] S. J. Van Enk, Phys. Rev. A 72, 064306 (2005).
  * Gheri  _et al._ [1998] K. M. Gheri, C. Saavedra, P. Törmä, J. I. Cirac, and P. Zoller, Phys. Rev. A 58, R2627 (1998).
  * Saavedra  _et al._ [2000] C. Saavedra, K. M. Gheri, P. Törmä, J. I. Cirac, and P. Zoller, Phys. Rev. A 61, 062311 (2000).
  * Schön  _et al._ [2005] C. Schön, E. Solano, F. Verstraete, J. I. Cirac, and M. M. Wolf, Phys. Rev. Lett. 95, 110503 (2005).
  * Lindner and Rudolph [2009] N. H. Lindner and T. Rudolph, Phys. Rev. Lett. 103, 113602 (2009).
  * Economou  _et al._ [2010] S. E. Economou, N. Lindner, and T. Rudolph, Phys. Rev. Lett. 105, 093601 (2010).
  * Schwartz  _et al._ [2016] I. Schwartz, D. Cogan, E. R. Schmidgall, Y. Don, L. Gantz, O. Kenneth, N. H. Lindner, and D. Gershoni, Science 354, 434 (2016).
  * Pichler  _et al._ [2017] H. Pichler, S. Choi, P. Zoller, and M. D. Lukin, Proceedings of the National Academy of Sciences , 201711003 (2017).
  * Borregaard  _et al._ [2020] J. Borregaard, H. Pichler, T. Schröder, M. D. Lukin, P. Lodahl, and A. S. Sørensen, Phys. Rev. X 10, 021071 (2020).
  * Wein  _et al._ [2022] S. C. Wein, J. C. Loredo, M. Maffei, P. Hilaire, A. Harouri, N. Somaschi, A. Lemaître, I. Sagnes, L. Lanco, O. Krebs, A. Auffèves, C. Simon, P. Senellart, and C. Antón-Solanas, Nat. Photon. 16, 374 (2022).
  * Tiurev  _et al._ [2022] K. Tiurev, M. H. Appel, P. L. Mirambell, M. B. Lauritzen, A. Tiranov, P. Lodahl, and A. S. Sørensen, Phys. Rev. A 105, L030601 (2022).
  * Gimeno-Segovia  _et al._ [2019] M. Gimeno-Segovia, T. Rudolph, and S. E. Economou, Phys. Rev. Lett. 123, 070501 (2019).
  * Kurpiers  _et al._ [2019] P. Kurpiers, M. Pechal, B. Royer, P. Magnard, T. Walter, J. Heinsoo, Y. Salathé, A. Akin, S. Storz, J. C. Besse, S. Gasparinetti, A. Blais, and A. Wallraff, Phys. Rev. Appl. 12, 044067 (2019).
  * Besse  _et al._ [2020] J. C. Besse, K. Reuer, M. C. Collodo, A. Wulff, L. Wernli, A. Copetudo, D. Malz, P. Magnard, A. Akin, M. Gabureac, G. J. Norris, J. I. Cirac, A. Wallraff, and C. Eichler, Nat. Commun. 11, 1 (2020).
  * Kannan  _et al._ [2020b] B. Kannan, D. L. Campbell, F. Vasconcelos, R. Winik, D. K. Kim, M. Kjaergaard, P. Krantz, A. Melville, B. M. Niedzielski, J. L. Yoder, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Science Advances 6, 10.1126/SCIADV.ABB8780 (2020b).
  * Wei  _et al._ [2021] Z. Y. Wei, D. Malz, A. González-Tudela, and J. I. Cirac, Phys. Rev. Res. 3, 023021 (2021).
  * Wei  _et al._ [2022] Z. Y. Wei, D. Malz, and J. I. Cirac, Phys. Rev. Lett. 128, 010607 (2022).
  * Ferreira  _et al._ [2022] V. S. Ferreira, G. Kim, A. Butler, H. Pichler, and O. Painter, arXiv: 2206.10076 10.48550/arxiv.2206.10076 (2022).
  * Dur  _et al._ [2000] W. Dur, G. Vidal, and J. I. Cirac, Phys. Rev. A 62, 062314 (2000).
  * Mittal  _et al._ [2016b] S. Mittal, V. V. Orre, and M. Hafezi, Optics Express 24, 15631 (2016b).
  * Rechtsman  _et al._ [2016] M. C. Rechtsman, Y. Lumer, Y. Plotnik, A. Perez-Leija, A. Szameit, and M. Segev, Optica 3, 925 (2016).
  * Mittal  _et al._ [2018] S. Mittal, E. A. Goldschmidt, and M. Hafezi, Nature 561, 502 (2018).
  * Dai  _et al._ [2022] T. Dai, Y. Ao, J. Bao, J. Mao, Y. Chi, Z. Fu, Y. You, X. Chen, C. Zhai, B. Tang, Y. Yang, Z. Li, L. Yuan, F. Gao, X. Lin, M. G. Thompson, J. L. O’Brien, Y. Li, X. Hu, Q. Gong, and J. Wang, Nat. Photon. 16, 248 (2022).
  * Tschernig  _et al._ [2021] K. Tschernig, . Jimenez-Galán, D. N. Christodoulides, M. Ivanov, K. Busch, M. A. Bandres, and A. Perez-Leija, Nat. Commun. 12, 1 (2021).
  * Ramos  _et al._ [2016] T. Ramos, B. Vermersch, P. Hauke, H. Pichler, and P. Zoller, Phys. Rev. A 93, 62104 (2016).
  * Guimond  _et al._ [2020] P. O. Guimond, B. Vermersch, M. L. Juan, A. Sharafiev, G. Kirchmair, and P. Zoller, npj Quantum Information 6, 1 (2020).
  * Soro and Kockum [2022] A. Soro and A. F. Kockum, Phys. Rev. A 105, 023712 (2022).
  * Lodahl  _et al._ [2017] P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, Nature 541, 473 (2017).
  * Porras and Fernández-Lorenzo [2019] D. Porras and S. Fernández-Lorenzo, Phys. Rev. Lett. 122, 143901 (2019).
  * Wanjura  _et al._ [2020] C. C. Wanjura, M. Brunelli, and A. Nunnenkamp, Nature Comm. 11, 1 (2020).
  * Ramos  _et al._ [2021] T. Ramos, J. J. Garcia-Ripoll, and D. Porras, Phys. Rev. A 103, 033513 (2021).
  * Gómez-León  _et al._ [2021] A. Gómez-León, T. Ramos, A. González-Tudela, and D. Porras, arXiv:2109.10930 (2021).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Phys. Rev. X 9, 041015 (2019).
  * Zirnstein and Rosenow [2021] H. G. Zirnstein and B. Rosenow, Phys. Rev. B 103, 195157 (2021).
  * McDonald  _et al._ [2018] A. McDonald, T. Pereg-Barnea, and A. A. Clerk, Phys. Rev. X 8, 041031 (2018).
  * McDonald and Clerk [2020] A. McDonald and A. A. Clerk, Nature Comm. 11, 1 (2020).
  * McDonald  _et al._ [2021] A. McDonald, R. Hanai, and A. A. Clerk, arXiv: 2103.01941 (2021).
  * Trescher and Bergholtz [2012] M. Trescher and E. J. Bergholtz, Phys. Rev. B 86, 241111 (2012).
  * Udagawa and Bergholtz [2014] M. Udagawa and E. J. Bergholtz, Journal of Statistical Mechanics: Theory and Experiment 2014, P10012 (2014).
  * Bergholtz  _et al._ [2015] E. J. Bergholtz, Z. Liu, M. Trescher, R. Moessner, and M. Udagawa, Phys. Rev. Lett. 114, 016806 (2015).
  * Zhang  _et al._ [2022] H. Zhang, J. Li, H. Jiang, N. Li, J. Wang, J. Xu, C. Zhu, and Y. Yang, Phys. Rev. A 105, 053703 (2022).
  * Fermi [1932] E. Fermi, Rev. Mod. Phys. 4, 87 (1932).
  * Shi  _et al._ [2015] T. Shi, D. E. Chang, and J. I. Cirac, Phys. Rev. A 92, 53834 (2015).

Supplemental Material: Topological multi-mode waveguide QED

In this Supplementary material, we provide the details of the calculations
supporting the main manuscript. In Section SM1, we describe the main
characteristic of the spectrum of the topological photonic lattice we
consider, explaining how we calculate the spectrum for open and periodic
boundary conditions. In Section SM2, we focus on the edge modes, and show how
their features can be captured by a simple phenomenological model. In Section
SM3, we give more details of the distinctive features of the spontaneous
emission of emitters coupled to the edge of these systems, such as their
expected decay rates and the spontaneous spatial separation of the photonic
emission patterns.

##  SM1 Characterization of the Harper-Hofstadter lattice model

The Harper-Hofstadter (HH) model that we consider along this manuscript is a
two-dimensional bosonic lattice where time-reversal symmetry is broken by an
artificial gauge field introduced through a Peierls phase [51, 52, 53].
Denoting by $a^{(\dagger)}_{\textbf{r}}$, the annihilation (creation)
operators of the bosonic mode at site ${\textbf{r}}=(x,y)$ of the lattice, the
HH Hamiltonian can be written (setting $\hbar=1$ along the manuscript):

$H_{B}=\omega_{a}\sum_{\textbf{r}}a_{\textbf{r}}^{\dagger}a_{\textbf{r}}-J\left(\sum_{x,y}\left(a_{x+1,y}^{\dagger}a_{x,y}+e^{-2\pi
i\phi x}a_{x,y+1}^{\dagger}a_{x,y}\right)+\text{H.c.}\right)\;,$ (SM1)

where we have assumed that time-reversal symmetry is broken by an artificial
uniform magnetic field in the direction perpendicular to the direction of the
lattice, which strength is encoded in the value of complex phase $\phi$
acquired in the nearest-neighbour hopping $J$ along the Y direction. Note that
we have also assumed the cavities energy, $\omega_{a}$, to be the same along
the lattice, and therefore we will consider it as the energy reference of the
problem setting $\omega_{a}\equiv 0$. Then, there are three magnitudes that
determine the shape of the spectrum and eigenmodes of the system, that are,
$J,\phi$, and the system size $L_{x}\times L_{y}$.

Let us start analyzing the spectrum of the system imposing periodic boundary
conditions along the X and Y direction, so that momentum
${\textbf{k}}=(k_{x},k_{y})$ is a good quantum number running over the
following values
$k_{\alpha}=-\pi,-\pi+2\pi/L_{\alpha},\dots,\pi-2\pi/L_{\alpha}$. As we said
in the main text, we will restrict to rational values of $\phi=p/q$, being
$(p,q)$ co-primes. This simplifies the diagonalization of the Hamiltonian,
since one can write an effective unit cell that describe the bath lattice of
$1\times q$ sites. With that unit cell, and using a plane-wave expansion to
account for the periodicity of the lattice, we find that the bath Hamiltonian
is given by $q$-bands due to the degeneracy introduced by the super-cell:

$\displaystyle
H_{B}=\sum_{\alpha=1}^{q}\sum_{\textbf{k}}\omega_{\alpha}({\textbf{k}})c^{\dagger}_{\alpha,{\textbf{k}}}c_{\alpha,{\textbf{k}}}\,,$
(SM2)

where $\alpha$ is the index running over the different bands, and
$c^{(\dagger)}_{\alpha,{\textbf{k}}}$ are the operators describing the
eigenstates of the bath for a given band $\alpha$ and momentum k. In Fig.
SM1(a), we plot an example of that bulk spectrum for $\phi=1/12$ and system
size $L_{x}=L_{y}=48$. There, we see the emergence of $q$ almost flat bands
well-separated in energies. These are the so-called Landau levels that appear
in such model [51, 52, 53], and which describe a cyclic motion of lattice
excitations whose radius in lattice constant units is given by the magnetic
length, $l_{B}=1/\sqrt{2\pi\phi}$. Thus, in order for these Landau level
picture to survive, it is required that such orbit fits in the lattice, i.e.,
$l_{B}\ll L_{\alpha}$. When this condition is not satisfied, the spectrum
tends to the square-lattice tight-binding spectrum
$\omega(k)=\omega_{a}-2J(\cos(k_{x})+\cos(k_{y}))$. This transition can be
observed in the bath density of states (DoS), defined as

$\text{DoS}(E)=\sum_{E_{B}}\delta(E-E_{B})\;,$ (SM3)

where the sum is performed over all the bath eigenenergies $E_{B}$. Since the
Dirac-$\delta$ has only mathematical sense in the continuum limit, for the
finite systems we consider we approximate it by a Gaussian distribution of
width $\theta$, that is:

$\text{DoS}(E)\approx\sum_{E_{B}}f_{\theta}(E-E_{B})$ (SM4)

with $f_{\theta}(E-E_{B})$:

$f_{\theta}(E-E_{B})=\frac{1}{\sqrt{2\pi\theta^{2}}}\exp\left(\frac{(E-E_{B})^{2}}{2\theta^{2}}\right)\;.$
(SM5)

Using that trick, we plot the DoS of the Harper-Hofstadter Hamiltonian in Fig.
SM1(c), as a function of $\phi$ for a fixed system size, showing the
transition from $q$-separated band to a unique band with a Van-Hove
singularity at the central energy $\omega_{a}$, characteristic of the nearest-
neighbour two-dimensional tight-binding model.

An important characteristic of these energy bands that appear in the bulk
spectrum is that they can have a non-zero quantized topological invariant
associated to them, that is, the Chern number. In fact, it can be shown that
the Chern number of the $l$-th band is given by $C_{l}=t_{l}-t_{l-1}$, where
the $t_{l}$ are integer numbers obtained by solving Azbel-Hofstader
Diophantine equation [52, 53]:

$l+1=qs_{l}+pt_{l}\;,\text{ with }|t_{l}|<\frac{q}{2}\,,$ (SM6)

being $t_{-1}=0$. Notice that the relation $C_{l}=t_{l}-t_{l-1}$ implies that
$t_{l}=\sum_{l^{\prime}<l}C_{l^{\prime}}$, i.e. the topological index $t_{l}$
reveals the sum of the Chern numbers of the lowest $l$ bands. For example, in
the case considered in the main text of $\phi=1/q$, it results that the first
$\left[(q-1)/2\right]$ bands (letting $\left[\cdot\right]$ be the floor
function) will have a Chern number $C_{l}=-1$. Thus, this choice provides us a
way to explore situations with different topological invariants $t_{l}$ just
by probing different energies. Note that other fluxes can lead to different
Chern number combinations. For example, as we will see below $\phi=4/9$ leads
to a lowest energy band with $C_{0}=2$, and $\phi=5/14$ leads to $C_{0}=3$.

Due to the bulk-boundary correspondence, these quantized topological invariant
will have important consequences when we consider open boundary conditions as
they will give rise to a number of localized, gapless edge states. To
illustrate that, we consider now the system to be placed in a cylinder
geometry, such that the system is periodic in the Y direction, but open in the
X direction, defining two edges where localized states can appear. To confirm
that localized edge states appear, we plot in Fig. SM1(b) the spectrum of the
system as function of $k_{y}$, which is still a good quantum number, for the
same parameters than the one in Fig. SM1(a). There, we observe how on top of
the flat bands appearing in the bulk modes, several edge-mode dispersion
appears associated to each of the lowest Landau-levels. This makes that the
larger the energy of the band-gap, the largest the number of edge states.

Figure SM1:  (a, b) Spectrum of a HH lattice of with (a) toroidal and (b)
cylinder topology, for a magnetic flux $\phi=1/12$ and a system size of
$48\times 48$ sites. When periodic boundary conditions are imposed over both
spatial directions, panel (a), the lower part of the spectrum consists of a
series of flat bands. If a boundary is opened, edge-localized, gapless
excitations emerge at band-gap energies, panel (b). The color of each dot is
assigned according to the localization parameter $\eta$ defined in Eq. (SM8).
(c, d) Density of states (DoS) of a HH lattice as defined in Eq. (SM3) for the
same lattice size, for different values of $\phi$: $1/12$ (dash-dotted blue
line), $1/30$ (dashed orange line) and $0$ (solid green line). The topology of
the lattice is (c) toroidal and (d) cylindrical, respectively. For $\phi\neq
0$, the DoS exhibits divergences at Landau levels. As $\phi\rightarrow 0$, a
van Hove divergence is formed at the middle of the spectrum, which is a
typical feature of a standard tight-binding square lattice with nearest-
neighbor hoppings. The most relevant difference between both instances is
that, for cylinder topology, the DoS deviates from zero in the band-gaps due
to the presence of topological edge modes. For these fluxes and system size,
the ratio between the magnetic lengths $l_{B}$ and the system size $L=48$ are
$l_{B}/L\approx 0.03$, $l_{B}/L\approx 0.05$ and $l_{B}/L=\infty$
respectively. Both DoS have been obtained using an auxiliary width (as
expressed in Eq. (SM5)) of $\theta/J=0.1$.

To further characterize the properties of the edge modes appearing in such
band-gaps, let us note that the eigenmodes in these configuration can be
written: $H_{B}\ket{k_{y},\beta}=E_{\beta}(k_{y})\ket{k_{y},\beta}$.
Projecting their wavefunction into their spatial coordinates $\langle
x,y|\Psi_{\beta,k_{y}}\rangle=e^{ik_{y}y}\psi_{\beta,k_{y}}(x)$, one can find
that $\psi_{\beta,k_{y}}(x)$ satisfies the Harper equation:

$\displaystyle\psi_{\beta,k_{y}}(x+1)+\psi_{\beta,k_{y}}(x-1)+2\cos\left(2\pi\phi
x-k_{y}\right)=E_{\beta}(k_{y})\psi_{\beta,k_{y}}(x)\;,$ (SM7)

and it features a localized shape. To make it more evident, we define a
localization parameter for each eigenstate as follows:

$\eta=\sum_{x=0}^{L-1}\left(-1+2\frac{x}{L_{x}-1}\right)|\psi_{\beta,k_{y}}(x)|^{2}\;,$
(SM8)

which features a maximum $\pm 1$ value when localized in left/right edge, and
$0$, when it is delocalized. We codify the value of that parameter in Fig.
SM1(a) in a color scale where purple/yellow indicates a maximum localization
in the left/right edges, whereas blue indicates delocalization. There, we
observe another important property of the edge states, that is, that the modes
along one edge are perfectly chiral, since they feature a positive/negative
group velocity depending on the edge where they are localized. This will have
important consequences when an emitter couples to one of the edges, as we will
see in Section SM3.

##  SM2 Effective edge mode description as a multi-mode waveguide

As shown in Fig. SM1(b), the first band-gaps of the HH can host a controllable
increasing number of edge modes. Since these are effectively one-dimensional
modes and chiral, they can be seen as an effective multi-mode one-way
waveguide [36, 37]. Generally, such topological have been described within
linear approximations [25, 26, 27]. However, from Fig. SM1 it is clear that
this is not the case in this scenario. In what follows, we will derive a more
accurate effective theory that is able to analytically capture the behaviour
of these multi-mode waveguides. For concreteness, we will derive such
expressions for the left-localized edge states, although a similar description
can be found of the right-localized ones.

###  SM2.1 Situation with $\phi=1/q$

Figure SM2: (a) Fitted values of the effective model parameter $a_{l}(\phi)$
in Eq. (SM9), characterizing the curvature of the edge mode dispersion for the
three lowest-energy edge states. In the limit of small fluxes,
$\phi\rightarrow 0$, all curvatures tend to the same value $\sim 0.6$. (b)
Fitted values of $k_{l}(\phi)$ as a function of the Landau level energy
$\omega_{l}$, for different magnetic flux $\phi=1/q$ ranging from $q=12$ to
$q=60$. We observe a linear trend $k_{l}=\beta_{l}\omega_{l}(\text{mod}2\pi)$.
Note that the discontinuity between the lines occur because of the definition
of $k_{l}(\pi)$ over the $[-\pi,\pi]$ range. Both figures (a) and (b) are
obtained considering a HH lattice of $40\times 40$ sites. (c) Localization
length $\lambda_{l}$ of the $3$ lowest edge modes for different values of
$\phi$, encoded in the markers’ color. All localization lengths are computed
at the same energy, namely the middle of the third spectral band-gap, as
depicted in the figure inset.

Let us start with the situation $\phi=1/q$ that we consider along the main
text and in Fig. SM1. After extensive numerical analysis, we find that a good
empirical ansatz for the $l$-th eigen-mode dispersion for small magnetic
fluxes is given by:

$\displaystyle\omega_{\mathrm{eff},l}(k_{y})$
$\displaystyle=\omega_{l}(\phi)+a_{l}(\phi)(k_{y}-k_{l}(\phi))^{2}\,,$ (SM9)

where the $\omega_{l}(\phi)$ can be found approximately in the perturbative
limit [52] as:

$\displaystyle\frac{\omega_{l}(\phi)}{J}=-4+2\pi\phi\left(l+\frac{1}{2}\right)-\left(\pi\phi\right)^{2}\left(l^{2}+l+\frac{1}{2}\right)+\mathcal{O}\left(\phi^{3}\right)\,,$
(SM10)

whereas $a_{l}(\phi)$ and $k_{l}(\phi)$ are fitting parameters that depend on
both the effective flux $\phi$ and $l$-th edge mode considered, although not
on system size as long as $l_{B}\ll L$. Note such quadratic energy dispersions
are typical of other (topologically-trivial) waveguides, where the finite size
effects introduce energy cut-off for the modes that lead to that behaviour. To
further characterize this effective model, we start plotting in Fig. SM2(a)
the evolution of the curvature of the edge modes, $a_{l}(\phi)$, as a function
of $\phi$ for the three lowest-energy edge modes in different colors. There,
we see how for big fluxes, $\phi$, the curvature of the modes differ
significantly, whereas for small fluxes, they converge to a value
$a_{l}(\phi)\approx 0.6$.

Regarding the value of the momentum cut-off $k_{l}(\phi)$, we find that there
is a linear dependence with $\omega_{l}(\phi)$, i.e.,
$k_{l}(\phi)=\beta_{l}\omega_{l}(\phi)(\mathrm{mod}2\pi)$. To illustrate that,
in Fig. SM2(b) we plot one against each other for several fluxes ranging from
$\phi\in[\frac{1}{60},\frac{1}{12}]$ in the different markers, together with a
dotted line which indicates the result of the fitting.

Figure SM3: (a, b) Effective model vs. exact diagonalization for a magnetic
flux $\phi=1/19$ and two different system sizes. In both subfigures, blue dots
represent the bath spectrum obtained by exact diagonalization. Purple/yellow
dots represent the prediction of the effective model of the energies of edge
modes localized at the left/right edge. Notice that varying $L$, the
dispersion relation of right-localized states is shifted along $k_{y}$, which
can also be captured by the effective model although we did not explicitly
show that. (c) Fitting error, as defined in Eq. (SM11), for the $3$ lowest
edge modes and varying flux. Empty and filled markers correspond to linear and
quadratic fittings respectively. We observe that, for all cases, the quadratic
fit error is few order of magnitudes lower, showing that this approach is
significantly more realistic.

In Fig. SM3(a), we show the comparison between the exact diagonalization
results and our effective description for a particular value of $\phi$,
showing indeed an excellent agreement for the lowest edge states dispersion.
In Fig. SM3(b), we make a more quantitative assessment on the quality of the
model by defining an error parameter:

$\varepsilon_{l}=\frac{1}{N_{k}}\sum_{k_{y}\in\Omega}|\omega_{\mathrm{exact},l}(k_{y})-\omega_{\mathrm{eff},l}(k_{y})|^{2}\;,$
(SM11)

where $\Omega$ is some region in the Brillouin zone along Y where we are
interested in performing the approximation to the exactly numerically
calculated edge-mode dispersion, $\omega_{\mathrm{exact},l}(k_{y})$, and
$N_{k}=|\Omega|$ is the number $k_{y}$ modes within that region. In Fig.
SM3(b), we plot $\varepsilon_{l}$ as a function of $\phi$ and compare the
accuracy between a linear fit, i.e., $\omega_{\mathrm{eff},l}(k_{y})\propto
k_{y}$ (empty markers) and the quadratic fit of Eq. (SM9) (filled markers),
showing how indeed the later provides a much more accurate approximation of
the modes for all fluxes.

Apart from the energies, another magnitude of interest of the modes is the
localization parameter. In particular, we know from Eq. (SM7) that when
$E_{\beta}(k_{y})$ lies within a band-gap, the spatial wavefunction along X
will be exponentially localized, i.e.,
$\psi_{\beta}(x)=\sqrt{2/\lambda_{l}(\phi)}e^{-x/\lambda_{l}(\phi)}$, where
$\lambda_{l}(\phi)$ will depend on both the energy level $l$ and the flux
$\phi$. In Fig. SM2(c) we also plot its dependence, showing that higher-energy
edge modes are less localized, and also that increasing the value of $\phi$
yields to higher de-localization.

###  SM2.2 Other situations $\phi\neq 1/q$, $q\in\mathbb{N}$

Along this work, we have restricted to magnetic fluxes of the form $\phi=1/q$,
with $q\in\mathbb{N}$, due to the topological features thoroughly discussed
along the manuscript, namely the sequence of $\left[q/2\right]$ lowest bands
with Chern number $C=-1$. This structure does not prevail if $\phi$ does not
fit this form. To illustrate this, we consider the cases of $\phi=4/9$ and
$5/14$ in Fig. SM4(a, c) and (b, d), which features a Chern-number of
$C_{0}=2,3$, respectively. We start by plotting the spectrum for such values
of $\phi$ in panels (a) and (b), where we see that the lowest band-gap feature
4 and 6 edge-state dispersions, respectively, as expected from the value of
$C_{0}$. In general, through numerical inspection we found that in these
situations the energy dispersions of the modes tend to be more similar than in
the different band-gaps of the $\phi=1/q$ situation. This will result in
qualitatively different spontaneous emission patterns, as observed in panels
(c) and (d), respectively, where we plot snapshots of the emission in two
different situations, illustrating the richness of this setup to obtain
qualitatively different photonic wavepackets.

Figure SM4: (a, b) Lowest energy spectrum of a Harper-Hofstadter lattice of
$80\times 80$ sites with a magnetic flux of (a) $\phi=4/9$ and (b)
$\phi=5/14$. In each case, the lowest band has a Chern number of $2$ and $3$
respectively, which correspond to the number of edge modes per boundary
resonant to the lowest spectral band-gap. (c, d) Snapshots of the photonic
dynamics at time $TJ=300$ of spontaneous emission along the $Y$ axis of an
emitter coupled to the left boundary of the lattice with coupling $g/J=0.2$.
The frequency of the emitter in each case is indicated in the (a, b) plots as
a solid orange line, and is $\omega_{e}^{(a)}/J=-2.45$ and
$\omega_{e}^{(b)}/J=-2.57$. We can observe that the pulse shape differs from
the single mode scenario. However, it is not possible in this case to resolve
the different peaks of the pulse due to the similarity of the group velocities
of the different topological channels.

##  SM3 Spontaneous emission of emitters coupled to the edge of the photonic
lattice

In this section, we will consider what happens when a two-level emitter, with
Hamiltonian $H_{S}=\omega_{e}\sigma_{ee}$, couples to one of the edges of such
HH lattice that, for concreteness, we assume to be the left one of Fig. 1(a)
In general, we will consider the most standard local light-matter couplings
given by:

$\displaystyle H_{I}=g\sigma_{eg}a_{{\textbf{r}}_{e}}+\mathrm{H.c.}\,,$ (SM12)

where $g$ is the coupling strength of the bath, ${\textbf{r}}_{e}$ the
position of the cavity mode the emitter couples to, and
$\sigma_{\alpha\beta}=\ket{\alpha}\bra{\beta}$ the dipole operator of the
optical emitter transition that couples the photonic bath. In Section SM3.5,
however, we will consider the non-local couplings that can be engineered with
giant atoms [54, 55, 56, 57, 58], as a way of selectively coupling some of the
topological edge modes.

###  SM3.1 Expected Markovian decay rates or Local density of states

A single emitter can be prepared in its excited state with a classical
driving, e.g., using a $\pi$-pulse. If one assumes that the bath has initially
no excitations, this state,
$\ket{\Psi_{0}}=\ket{e}\otimes\ket{\mathrm{vac}}_{B}$, can only evolve into an
state of the form:

$\displaystyle\ket{\Psi(t)}=\left(C_{e}(t)\sigma_{eg}+\sum_{\textbf{r}}A_{\textbf{r}}(t)a^{\dagger}_{\textbf{r}}\right)\ket{g}\otimes\ket{\mathrm{vac}}_{B}\,,$
(SM13)

because the full light-matter Hamiltonian: $H=H_{S}+H_{B}+H_{I}$, conserves
the number of excitations
$N_{\mathrm{exc}}=\sigma_{ee}+\sum_{\textbf{r}}a_{\textbf{r}}^{\dagger}a_{\textbf{r}}$
since $[H,N_{\mathrm{exc}}]=0$. Using time-dependent perturbation theory or,
equivalently, a Markovian approximation for the system-bath coupling, the
emitter is expected to show an exponential decay of its excitation, i.e.,
$|C_{e}(t)|^{2}\approx e^{-\Gamma t}$, with $\Gamma$ being the expected
Markovian decay rate given by Fermi’s Golden rule [106]:

$\displaystyle\Gamma(\omega_{e})=2\pi
g^{2}\sum_{\beta}|\psi_{\beta}(0)|^{2}\delta\left(\omega_{e}-E_{\beta}(k_{y})\right)=-2g^{2}\text{Im}\left[\frac{1}{\omega_{e}+i0^{+}-H_{B}}\right]_{{\textbf{r}}_{e},{\textbf{r}}_{e}}=2g^{2}\;\text{LDoS}({\textbf{r}}_{e},\omega_{e})\,.\leavevmode\nobreak\
$ (SM14)

In the last equality we introduced the local density of states (LDoS) at the
emitter position. This quantity is defined similarly to the regular DoS, but
weighting the contribution of each bath eigenstate by its support on the
position of the emitter:

$\text{LDoS}({\textbf{r}},E)=\sum_{E_{B}}|\langle{\textbf{r}}|E_{B}\rangle|^{2}\cdot\delta(E-E_{B})$
(SM15)

Figure SM5: (a) Markovian decay rate of a quantum emitter coupled to the edge
of a HH lattice of $100\times 100$ sites, for different values of the magnetic
flux $\phi$, and computed by exact diagonalization (computed using an
auxiliary width of $\theta/J=0.15$ as defined in Eq. (SM5)). The decay rate
curves in the lowest part of the spectrum is zoomed in the figure, showing
that lower fluxes support more constant plateaus in the spectral band-gaps.
(b) For the same lattice size, the Markovian decay rate for $\phi=1/25$ is
depicted as a solid line. Square markers depict the prediction of our
effective model for the decay rate in this specific configuration. (c)
Spontaneous decay of the emitter population in the same configuration as in
(b) for different values of $\omega_{e}$ and $g/J=0.05$. Each curve is linked
to a value of $\omega_{e}$ and its corresponding color depicted by the square
markers in (b).

From this definition, the relation of the LDoS with the Markovian decay rate
becomes clear by comparing with Eq. (3) of the main text, since

$|\langle e|H_{I}|E_{B}\rangle|^{2}=|\langle
e|g\sigma_{eg}a_{{\textbf{r}}_{e}}|E_{B}\rangle|^{2}=g^{2}|\langle
0|a_{{\textbf{r}}_{e}}|E_{B}\rangle|^{2}=g^{2}|\langle{\textbf{r}}_{e}|E_{B}\rangle|^{2}\;,$
(SM16)

where $|0\rangle=|g\rangle|\text{vac}\rangle$. We numerically compute the LDoS
in a similar fashion as we did for the DoS, defining a ’smoothed’ Dirac delta
function $f_{\theta}(E-E_{B})$, that we take to be a Gaussian distribution as
expressed in Eq. (SM5), and computing the LDoS as

$\text{LDoS}({\textbf{r}}_{e},E)\approx\sum_{E_{B}}|\langle{\textbf{r}}_{e}|E_{B}\rangle|^{2}\cdot
f_{\theta}(E-E_{B})\;,$ (SM17)

From Eq. (SM14) we read that the shape of the LDoS determines the shape of the
expected Markovian decay rate $\Gamma(\omega_{e})$. In Fig. SM5(a) we plot the
Markovian decay rate at the emitter position calculated with exact
diagonalization, for several values of $\phi$. We observe that smaller
$\phi$’s are associated with narrow band-gaps, situation in which the LDoS
’quasi-quantized’ behaviour leads to plateaus in $\Gamma(\omega_{e})$. In Fig.
SM5(b), we represent with a solid line the expected decay rate $\Gamma_{e}$ as
a function of $\omega_{e}$ for $\phi=1/12$, computed using exact
diagonalization, and compare it with the markers that represent the expected
decay rate computed as follows:

$\Gamma_{l}(\omega_{e})\approx\frac{g^{2}|\psi_{l}(0)|^{2}}{|v_{g,l}(\omega_{e})|}=\frac{2g^{2}}{\lambda_{l}|v_{g,l}(\omega_{e})|}\,,$
(SM18)

where $v_{l,g}(\omega_{e})$ is the group velocity of the $l$-th mode at the
emitter energy, that is,
$v_{g,l}=\partial_{k_{y}}\omega_{\mathrm{eff},l}(k_{y})|_{k_{y}=k_{e}}$, with
$\omega_{\mathrm{eff},l}(k_{e})=\omega_{e}$. Thus, the total decay rate in a
given band-gap will be given by
$\Gamma(\omega_{e})=\sum_{l}\Gamma_{l}(\omega_{e})$, where the sum runs over
the number of edge modes that are present in that band-gap. We observe that
the decay rate displays a ladder-like structure, with jumps located at Landau
levels. This quasi-quantization behaviour can be probed by monitoring the
decay of the emitter population during spontaneous emission. In Fig. SM5(c) we
show the quantum emitter dynamics of an initially excited emitter for the
range of $\omega_{e}$ depicted in Fig. SM5(b). There, we observe how its
timescale remains approximately constant until it crosses the Landau level
energy and is able to interact with the new mode of the higher band-gap. We
note, however, that the steps of the ladder are not strictly constant due to
the non-linear dependence of the mode dispersion. On the contrary, they
display the typical $1/\sqrt{\omega-\omega_{\mathrm{edge}}}$ dependence
associated to one-dimensional quadratic band-edge dispersions.

Figure SM6: Markovian decay rate from exact diagonalization at the boundary of
a HH lattice of size $120\times 120$ and $\phi=1/15$. The figures are obtained
using an auxiliary function width of (a) $\theta_{1}/J=0.03$, (b)
$\theta_{2}/J=0.07$ and (c) $\theta_{3}/J=0.16$. The figure insets represent
as crossed markers the distribution of bath eigenvalues
$\left\\{E_{B}\right\\}$. Above these markers, we include a representation of
the auxiliary probe Gaussian function $f_{\theta}(E-E_{B})$ that is used as
described in Eq. (SM17) to probe the energy spectrum and compute the LDoS. In
each figure, the vertical yellow fringes are centered at Landau levels, and
their width is equal to the corresponding value of $\theta$.

For the lowest-energy edge modes, we observe the agreement between the our
analytical model and the numerical calculations, except near the Landau level
divergences. At these energies the numerical computation of the LDoS entails
limited resolution, due to the finite width $\theta$ associated to the
auxiliary function $f_{\theta}(E-E_{B})$: if we probe the LDoS of an energy
$|E-\omega_{l}|<\theta$ from below, the LDoS will count some states above
$\omega_{l}$, softening the transition. This issue is unavoidable: if we try
to have an arbitrarily small value of $\theta$, the probe function will
eventually observe the discreteness of the spectrum. In such case, the
computation of the LDoS would suffer from numerical instabilities. On the
contrary, an excessively value for $\theta$ will make the approximation states
in Eq. (SM17) increasingly worse. In our case, the information about the LDoS
jumps would be lost, softening the shape of the curve. All these numerical
issues are graphically represented in Fig. SM6.

###  SM3.2 Photonic spontaneous emission patterns

As we see in the Fig. 3 of the main text, the non-uniform group velocity of
the modes along a given band-gap favours an spatial separation of the photons
propagating into the different channels. This generates naturally single-
photon time-bin entangled states [64], that when combined with sequential
generation methods [65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78,
79, 80, 81] can be used to generate complex states of light in these
topologically-protected channels. Let us now analyze here the relevant
magnitudes that determine such spontaneous separation of the photons focusing
on the second band-gap where there are two edge-modes. If we neglect the
broadening introduced by the non-linear mode dispersion, which we will see
below it is a good approximation in our system, the spatial separation between
the different modes is determined by:

  * •

The different group velocities, $v_{g,l}$, which makes that after a time $T$,
the wavepacket fronts are separated by: $|v_{g,l}-v_{g,l^{\prime}}|T$. With
our choice of units, that length is already normalized to lattice constant
units.

  * •

However, the spatial modes have an intrinsic broadening (in lattice constant
units) of order $1/\Gamma_{l}$ in lattice constant units [107].

Thus, one can define a parameter $R_{l,l^{\prime}}$:

$\displaystyle
R_{l,l^{\prime}}=\frac{|v_{g,l}-v_{g,l^{\prime}}|}{\Gamma_{l}^{-1}+\Gamma_{l}^{\prime-1}}\,.$
(SM19)

which quantifies how favorable is a given configuration to observe the
separation of the modes. This quantity has dimensions of energy (i.e. inverse
time) and captures the competition between the spatial broadening
$\propto\Gamma_{l}^{-1}+\Gamma_{l^{\prime}}^{-1}$ and the pulse separation
induced by the difference between group velocities
$|v_{g,l}-v_{g,l^{\prime}}|$. When $R_{ll^{\prime}}T\sim 1$, it is expected
that the modes $l$ and $l^{\prime}$ are fully resolvable at time $T$. In Fig.
SM7 we plot the different group velocities of the modes (panel (a)),
broadenings (panel (b)), and $R_{0,1}$ (panel (c)) as a function of the
quantum emitter energy $\omega_{e}$ for a given bath configuration with
$\phi=1/40$. From Fig. SM7(c) we observe that the pulses associated with modes
$0$ and $1$ will be fully resolvable at $TJ\sim 10^{3}$. As expected, the best
conditions to achieve resolution will occur when the quantum emitter energy
$\omega_{e}$ is slightly above the Landau level energy $\omega_{1}$: in this
case, $|v_{g,0}-v_{g,1}|$ will be maximum, which favours resolution.

Figure SM7: (a) Group velocities of the lowest $4$ edge modes of a HH lattice
of $50\times 50$ sites and $\phi=1/40$. (b) Inverse decay rate
$1/\Gamma^{(l)}$ of the quantum emitter onto each topological channel $l$, in
units of $J/g^{2}$, and for the case of $g/J=0.1$. (c) Distinguishability
parameter $R_{01}$, as defined in Eq. (SM19), for the $2$ lowest edge modes in
the $\phi=1/40$ configuration.

Apart from the intrinsic broadening of the emitted wavepackets, let us also
note that there is an additional source of broadening coming from the
curvature of the edge-mode dispersion at the emitter’s frequency, that is,
$\gamma_{\mathrm{dis}}$ in
$\omega_{l}(k)\approx\omega_{e}+v_{g,l}(k-k_{e})+\gamma_{\mathrm{dis}}(k-k_{e})^{2}/2$.
In particular, it is well known that a wavepacket with an initial broadening
$\sigma_{0}$ propagating in such a non-linear dispersive channel will have an
increasing size growing with:

$\sigma(t)=\sigma_{0}\sqrt{1+\frac{\gamma_{\mathrm{dis}}^{2}t^{2}}{\sigma_{0}^{4}}}\;.$
(SM20)

In the case of spontaneous emission of a quantum emitter into the
$l^{\text{th}}$-topological channel of a HH lattice, we can take $\sigma_{0}$
as the inverse of the decay rate onto such mode, $\Gamma_{l}^{-1}$. The pulse
width will then evolve as:

$\sigma(t)\sim\frac{1}{\Gamma_{l}}\sqrt{1+\gamma_{\mathrm{dis}}^{2}\Gamma_{l}^{4}t^{2}}\;.$
(SM21)

From this evolution equation, it follows that the role of dispersion will be
negligible at a certain time $T$ as long as

$\gamma_{\mathrm{dis}}^{2}\Gamma_{l}^{4}T^{2}\ll 1$ (SM22)

Taking $g/J\sim 0.1$, we find that $\Gamma_{l}/J\sim 10^{-3}$. On the other
hand, we can estimate $\gamma_{\mathrm{dis}}$ from our effective theory: in
particular, we will have that $\gamma_{\mathrm{dis}}\sim
2a_{l}(\phi)(k_{y}-k_{l}(\phi))$, which will be at most of the order of the
unity. Thus, at this value of $g$, dispersion effects will be thus relevant at
times $TJ\sim 10^{6}$, which is several orders of magnitude above the
timescale where the pulses are separable due to different group velocities,
characterized by $R_{ll^{\prime}}T\sim 1$.

###  SM3.3 Robustness to disorder

In Figs. 3(a, c, e) of the main text, we have shown the robustness of single,
two, and three-edge mode propagation to an edge defect. This is a consequence
of the topological nature of the edge modes. Here, we discuss in greater
detail the protection of photon propagation against disorder by introducing
random perturbations in the energy of local lattice modes, as follows:

$H\rightarrow
H+\sum_{{\textbf{r}}}\delta\omega_{{\textbf{r}}}a_{{\textbf{r}}}^{\dagger}a_{{\textbf{r}}}\;,$
(SM23)

where $\delta\omega_{{\textbf{r}}}$ is a random variable uniformly distributed
along the interval $(-\sigma,\sigma)$, where $\sigma$ is the strength of the
applied disorder. Topological gapless modes spectrally located at a band-gap
of width $E_{W}$ are typically expected to be robust to disorder strenghts
$\sigma$ up to the order of $E_{W}$.

In Fig. SM8, we analyse the effect of disorder in two key features: the LDoS
defined in Eq. (SM15) and photon emission and propagation. In Fig. SM8(a) we
plot the Markovian decay rate, computed from exact diagonalization, for an
emitter coupled at the boundary of a HH lattice for $\phi=1/12$ for different
values of $\sigma$ averaged among different realizations of disorder. We
observe that the laddered structure is preserved for values of $\sigma$
comparable to the width of the lowest spectral band-gaps. This is a
consequence of the protection of the multi-mode spectrum in the topological
band-gaps. Furthermore, in Fig. SM8(b-d), we plot photon population emitted
for an emitter resonant to the second lowest band-gap of the bath spectrum and
for different values of $\sigma$. For weak disorder i.e. small $\sigma$
compared to $E_{W}$, we observe the same light-cone configuration as in the
$\sigma=0$ case. We only start to witness light-cone distortion for values of
$\sigma$ comparable to the size of the band-gap that is resonant to the
emitter energy.

Figure SM8: (a) Ensemble average of the Markovian decay rate of a quantum
emitter coupled at the boundary of a HH lattice of $30\times 30$ sites for
$\phi=1/12$ for different values of disorder strength: $\sigma/J=0$ (blue),
$\sigma/J=0.35$ (orange) and $\sigma/J=0.7$ (green). Each solid line represent
the mean LDoS for a collection of $100$ realizations of disorder, which
variance is represented as a shaded region around the mean. Every curve is
computed by exact diagonalization, using an auxiliary width of $\theta/J=0.15$
as defined in Eq. (SM5). (b)-(c)-(d) Photon transport from the spontaneous
emission of such quantum emitter with energy $\omega_{e}/J=-2.48$ along the
left boundary of a HH lattice of $70\times 70$ sites and $\phi=1/12$, for
several values of disorder strength. We observe that photon transport is
robust for small values of $\sigma$. The effects of disorder only start to
arise for $\sigma\sim 0.5J$, which is of the order of the band-gap width
$E_{W}\sim 0.77J$.

###  SM3.4 Photon-loss effects

Along all the manuscript, we have neglected the possibility that the emitter
or bath modes couples to other additional bath, generating additional decay
rates. In that case, the emitter-bath dynamics must be described by a density
matrix, $\rho(t)$, formalism. Assuming that the coupling to these baths are
Markovian, the dynamics of such density matrix can be described by the
following time-local master equation:

$\displaystyle\frac{\partial\rho}{\partial
t}=i[\rho,H]+\sum_{\textbf{r}}\frac{\kappa_{\textbf{r}}}{2}(2a_{\textbf{r}}\rho
a_{\textbf{r}}^{\dagger}-a^{\dagger}_{\textbf{r}}a_{\textbf{r}}\rho-\rho
a^{\dagger}_{\textbf{r}}a_{\textbf{r}})+\frac{\Gamma^{*}}{2}(2\sigma_{ge}\rho\sigma_{eg}-\sigma_{ee}\rho-\rho\sigma_{ee})\,,$
(SM24)

where $\kappa_{\textbf{r}}$ and $\Gamma^{*}$ are the Markovian decay rates
induced by these additional baths in the bath and emitter modes, respectively.
In the spontaneous emission configuration that we have considered along this
manuscript, all the effect of these baths can be captured by replacing the
full light-matter Hamiltonian $H$ by a non-Hermitian version which include the
effect of the losses:

$\displaystyle
H^{*}=H_{S}+H_{B}+H_{I}-i\sum_{\textbf{r}}\frac{\kappa_{\textbf{r}}}{2}a^{\dagger}_{\textbf{r}}a_{\textbf{r}}-i\frac{\Gamma^{*}}{2}\sigma_{ee}\,.$
(SM25)

Figure SM9: (a) Markovian decay rate of a quantum emitter coupled at the edge
of a HH lattice of $40\times 40$ and $\phi=1/12$, for several values of the
local mode loss rate $\kappa$. (b, c) Excitation dynamics along the left
boundary of a HH lattice with flux $\phi=1/12$, induced by the spontaneous
emission of a quantum emitter with energy $\omega_{e}/J=-2.48$, resonant to
two edge modes, and a local loss rate of (a) $\kappa/J=0$ and (b)
$\kappa/J=0.05$ on the lattice sites. We observe two light-cones, associated
to the two group velocities of both topological channels.

In Fig. SM9(a), we plot its effect in the predicted Markovian decay rates for
increasing values of $\kappa_{\textbf{r}}=\kappa$ and fixing $\Gamma^{*}$.
There, we observe a “softening” of the steps of the ladder, although the
quasi-quantized behaviour remains unaltered to a great extent. A similar
behaviour will occur with $\Gamma^{*}\neq 0$. The most significant effect of
$\kappa$ is to generate a finite propagation length of the photon modes, as
observed in Fig. SM9(b-c) where we plot an example of the spontaneously
emitted photons without and with dissipation, respectively..

###  SM3.5 Mode-selectivity using non-local couplings

Local light-matter Hamiltonians, as the ones considered in Eq. (SM12), lead to
completely delocalized couplings in momentum space. For example, assuming the
cylinder geometry that we considered in the previous section, one can rewrite
the light-matter Hamiltonian $H_{I}$ as follows:

$\displaystyle
H_{I}=g\sigma_{eg}a_{{\textbf{r}}_{e}}+\mathrm{H.c.}=\sum_{k_{y}}\left(\frac{g}{\sqrt{L_{y}}}\sigma_{eg}a_{(x_{e},k_{y})}e^{-ik_{y}y_{e}}+\mathrm{H.c.}\right)\,,$
(SM26)

where:

$\displaystyle
a_{(x_{e},k_{y})}=\frac{1}{\sqrt{L_{y}}}\sum_{y}a_{(x_{e},y)}e^{-ik_{y}y}\,.$
(SM27)

This makes that when coupling an emitter with an optical transition of a given
energy, $\omega_{e}$, it couples to all the edge modes of a given energy.
Thus, a quantum emitter locally coupled to the edge of a HH lattice will
couple to all topological channels. In this section, we address the problem of
selecting the channels onto the quantum emitter can decay by using non-local
couplings such that the emitted photon couples selectively to some of the
topological channels. The key point is that if one let the emitter to couple
several edge sites with different strength, $g_{j}$, the light-matter coupling
acquires a k-dependence depending on these couplings:

$\displaystyle
H_{I}=\sum_{y}g_{y}\sigma_{eg}a_{(x,y)}+\mathrm{H.c.}=\sum_{k_{y}}\left(G(k_{y})\sigma_{eg}a_{(x_{e},k_{y})}+\mathrm{H.c.}\right)\,,$
(SM28)

with:

$\displaystyle G(k_{y})=\frac{1}{\sqrt{L_{y}}}\sum_{y}g_{y}e^{-ik_{y}y}\,.$
(SM29)

Such k-dependent coupling can be used to cancel the coupling to certain
momentum and thus prevent the emission in these modes. For example, this has
been used in one [88] and two-dimensional [55] single-mode waveguides to
generate chiral, one-directional emission in baths with isotropic bath
dispersions. In our topological multi-mode waveguides, for fixed $\omega_{e}$,
the resonant modes feature a different resonant $k_{e}^{(l)}$. Thus, choosing
the $g_{y}$ such that $G(k_{e}^{(l)})=0$ to some of the modes, one can
selectively couple to the other ones.

Let us now illustrate a potential method to achieve that selectivity inspired
by the results in Refs. [88, 55] that only requires coupling to $N_{c}+1$
cavities if one wants to cancel the coupling to $N_{c}$ resonant momenta
$k_{e}^{(l)}$. Let us start by the case where we want cancel only the coupling
to a single resonant momenta $k_{e}^{(l)}$. For that, let us assume that the
emitter couples to neighbouring sites with the same amplitude, but a relative
complex phase $\varphi_{l}$, that is: $g_{y}=g$, $g_{y+1}=ge^{i\varphi_{l}}$.
Using these values, it can be shown from Eq. (SM29) that:

$\displaystyle|G_{\varphi_{l}}(k_{y})|^{2}=\frac{2g^{2}}{L_{y}}(1+\cos(k_{y}+\varphi_{l}))\,.$
(SM30)

Thus, if we want to make it zero for certain $k_{y}=k_{e}^{(l)}$, one just
have to choose a relative phase: $\varphi_{l,e}=\pi-k_{e}^{(l)}$. In Fig. SM10
we show a proof-of-principle realization of that idea for a situation when the
emitter is resonant to two edge modes. In Fig. SM10(a), we plot the energy
spectrum for a lattice with $\phi=1/12$ and $65\times 65$ sites. In solid
horizontal red line we indicate the energy of the emitter’s that is chosen in
the second band-gap to be resonant to two channels with resonant momenta
$k_{e}^{(0)}$ and $k_{e}^{(1)}$, respectively, indicated in vertical dotted
black lines. This means that if the emitter couples locally, it will couple to
the two $k$-channels, as shown in Fig. SM10(b), and emit in a two-mode fashion
in real space, as shown in black dots in Fig. SM10(e). On the contrary, if we
choose non-local couplings to cancel the coupling to either the momentum
$k_{e}^{(0)}$ or $k_{e}^{(1)}$, as depicted in orange and blue dotted lines in
Fig. SM10(a), one can see that the emitter selectively emits only in one of
the channels, as it is illustrated in Figs. SM10(c-e) by plotting snapshots of
the population in momentum and real space. In particular, in Fig. SM10(e)
where we plot a snapshot at time $TJ=150$ of the spatial profile of the
emitted pulse for the different coupling choices, we observe very clearly that
the designed non-local couplings suppress the emission onto the selected mode
compared to local light-matter coupling situation.

Figure SM10: (a) Spectrum of a HH lattice of $65\times 65$ sites for
$\phi=1/12$. Each dot is coloured according to the localization index $\eta$
as in Fig. 2 of the main text. The solid pink line depicts the quantum emitter
frequency $\omega_{e}/J=-1.86$, that is resonant to two left-localized edge
modes at momenta $k_{e}^{(0)}$ and $k_{e}^{(1)}$ respectively. We prove mode-
selectivity using non-local couplings that in momentum space are of the form
$|G_{\varphi_{0}}(k_{y})|^{2}$ and $|G_{\varphi_{1}}(k_{y})|^{2}$ represented
as dashdotted lines, which units are indicated in the right vertical axis, and
are analytically described in Eq. (SM30). Note that the orange (blue) line,
representing $|G_{\varphi_{0}}(k_{y})|^{2}$ ($|G_{\varphi_{1}}(k_{y})|^{2}$)
vanishes at $k_{e}^{(0)}$ ($k_{e}^{(1)}$). (b, c, d) Normalized photon
population in momentum space for the photonic state resulting from spontaneous
emission at time $TJ=300$ of a quantum emitter with energy
$\omega_{e}/J=-1.86$ coupled to the left edge of a HH lattice of size
$250\times 250$ and $\phi=1/12$. The emitter coupling is (b) local (c, d) non-
local, designed to cancel emission at $k_{e}^{(0)}$ and $k_{e}^{(1)}$
respectively. (e) Snapshot at time $TJ=150$ of the photonic state in the left
lattice boundary, coming from spontaneous emission in the three distinct
coupling configurations, with same values for $\omega_{e}$ and $\phi$. The
lattice size is $600\times 600$, and the light-matter coupling constant is
$g/J=0.2$.

In the case where there are more than two edge modes, one might be interested
in cancelling the coupling to two or more momenta. A way of doing that would
be to obtain an effective $G(k_{y})$ out of the product required to cancel
each momentum, separately. In the general case, where we are interested in
suppressing $N_{k}$ modes, such proposal implies taking $G(k_{y})$ as

$G(k_{y})=\prod_{\delta=0}^{N_{k}-1}G_{\varphi_{\delta}}(k_{y})\;,$ (SM31)

where $\varphi_{\delta}=\pi-k_{e}^{(\delta)}$. With this choice, we get a non-
local coupling in k-space of the form:

$G(k_{y})=\prod_{\delta=0}^{N_{k}-1}\left(1+e^{i(\varphi_{\delta}-k_{y})}\right)=\sum_{M=0}^{N_{k}-1}e^{-iMk_{y}}\sum_{|\Omega|=M}\exp\left(i\sum_{\delta\in\Omega}\varphi_{\delta}\right)\;,$
(SM32)

where the sum over $|\Omega|=M$ denotes a sum over all the
$N_{k}!/M!(N_{k}-M)!$ subsets of $\left\\{1,2,...,N_{k}\right\\}$ with $M$
elements. From this expansion in powers of $e^{-ik_{y}}$, we can directly read
the couplings in real-space:

$g_{(0,y+M)}=\sum_{|\Omega|=M}\exp\left(i\sum_{\delta\in\Omega}\varphi_{\delta}\right)\;,$
(SM33)

with $M=0,...,N_{k}-1$. This shows that the number of non-local couplings
required to cancel $N_{k}$ modes scales linearly with $N_{k}$, although this
method will be eventually limited by the capacity to resolve the modes. As an
example, for $N_{k}=2$, we have

$\displaystyle G_{\varphi_{0},\varphi_{1}}(k_{y})\propto
1+(e^{i\varphi_{0}}+e^{i\varphi_{1}})e^{ik_{y}}+e^{i(\varphi_{0}+\varphi_{1})}e^{i2k_{y}}\,,$
(SM34)

which can be achieved with real-space couplings of the form $g_{(0,y)}=g$,
$g_{(0,y+1)}=ge^{i\varphi_{0}}+ge^{i\varphi_{1}}$ and
$g_{(0,y+2)}=ge^{i(\varphi_{0}+\varphi_{1})}$.