Topological multi-mode waveguide QED

Topological insulators feature a number of topologically protected boundary modes linked to the value of their bulk invariant. While in one-dimensional systems the boundary modes are zero dimensional and localized, in two-dimensional topological insulators the boundary modes are chiral, one-dimensional propagating modes along the edges of the system. Thus, topological photonic insulators with large Chern numbers naturally display a topologically protected multimode waveguide at their edges. Here, we show how to take advantage of these topologically protected propagating modes by interfacing them with quantum emitters. In particular, using a Harper-Hofstadter lattice, we find situations in which the emitters feature quasiquantized decay rates due to the increasing number of edge modes, and where their spontaneous emission spatially separates in different modes. We also show how using a single $\pi$-pulse the combination of such spatial separation and the interacting character of the emitters leads to the formation of a single-photon time-bin entangled state with no classical analog, which we characterize computing its entanglement entropy. Finally, we also show how the emitters can selectively interact with the different channels using nonlocal light-matter couplings such 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.

The more complex scenario when the two-dimensional system features several edge modes [37,38] 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- [39][40][41], two- [40,42], and three- [43,44] dimensional topological insulators. Besides, the strongly interacting character of the emitters can induce (non- x,y . The yellow-shaded region denotes a single lattice plaquette with flux ϕ. We also represent the spontaneous emission of an emitter (in red) coupled to the the lattice edge, which radiates through all edge modes (2 in the figure), that propagate at different group velocities (denoted by vg,0 and vg,1), causing a spatial separation of the emitted pulse. gaussian) quantum correlations beyond the ones that can be obtained with parametric drivings [19][20][21][22][23][24][25], opening a path to observe exotic quantum many-body states [45][46][47][48][49][50]. All these reasons are motivating the development of such topological light-matter interfaces in various platforms, ranging from superconducting qubits coupled to microwave resonators [18,51] to solid-state emitters coupled to topological photonic crystals [14][15][16][17].
In this work, 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. The physics emerging from this scenario is very different from the case of emitters coupled to one-dimensional topological insulators [39][40][41], where the localized nature of the boundary modes leads generally to coherent emitter dynamics/interactions rather than irreversible or collective decay dynamics. Using that theory, in Section II we un-arXiv:2207.02090v2 [quant-ph] 18 Aug 2023 veil, for the first time, the entanglement structure of the spontaneously emitted photons in such topological multimode waveguides, and show one can obtain almost maximally entangled W -type states between the edge channels. Besides, combining several π-pulse [52], we show one can obtain strongly-correlated multi-photon states. Finally, in Section IV we also devise a method that enable the emitters interact selectivily with the different topological channels using giant atoms [53][54][55][56][57][58].

II. MODEL
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 two-dimensional topological insulator. Motivated by recent experiments [18], we particularize for a square photonic lattice with nearestneighbour hoppings of rate J, subject to an effective magnetic flux, ϕ -the so-called Harper-Hofstadter (HH) lattice [59]-, where bands with large Chern number appear for small magnetic fluxes [60,61]. The bath Hamiltonian then reads (setting ℏ = 1): ,y a x,y + e −2πiϕx a † x,y+1 a x,y + H.c. , (1) where a ( †) x,y represent the annihilation (creation) operator at the (x, y) position, and where we take the cavity energy as the energy reference. The emitter is assumed to have a single optical transition between its ground (g) and excited state (e) with frequency ω e , that couples to one of the photonic lattice sites at the edges through the standard light-matter Hamiltonian, H I = (ga xe,ye σ eg + H.c.), with g being its coupling strength, (x e , y e ) the position where it couples, and σ αβ = |α⟩ ⟨β| the emitter's operators. The emitter's Hamiltonian then reads, H S = ω e σ ee , such that the full topological light-matter Hamiltonian reads H = H S + H B + H I .
The HH model displays a very rich behaviour depending on the value of ϕ [59][60][61]. here, we take ϕ = 1/q, with q ∈ N, that is 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 [60], separated by l = 0, . . . , q − 1 band-gaps. Besides, it can be shown that the first [(q − 1)/2]-bands (being [·] the floor function) have an associated quantized Chern number C = −1 [60,62]. 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 [63]. To illustrate this, 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- Figure 2. (a) HH spectrum with periodic boundary conditions in the Y direction -open for X -for ϕ = 1/9 and a lattice size of 65 × 65 sites. Each dot is coloured according to the localization parameter η defined in the main text: η = 0 corresponds to a delocalized state, while η = (−)1 depicts complete localization at the (left) right boundary. (b) Emitter spontaneous emission rate as a function of its transition frequency ωe approximating δ(ωe −E) by a Gaussian function with mean ωe − EB and width θ = 0.07J [62]. The emitter is coupled to the edge of a HH lattice of size 150 × 150 with flux ϕ = 1/25. Yellow vertical fringes are centered at Landau levels and have a width equal to θ. (c) Entanglement entropy, E(Ψ ph ), of the emitted single-photon state in the asymptotic limit [62], as a function of the emitter energy, for the same lattice parameters as panel (b).
states as ψ l (x, y) = e ikyy Ψ l (x), and calculate their eigenenergies ω(k y ) numerically. In Fig. 2(a), we plot an example of the bath spectrum for a bath with ϕ = 1/9 and 65×65 sites, showing the emergence of the gapless modes between the bulk flat bands. Besides, we define a localization parameter η = L−1 x=0 −1 + 2 x L−1 |Ψ(x)| 2 , that ranges η ∈ [−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 ω a , the edge modes at the left (right) have always negative (positive) group velocity along Y, thus, having a chiral character. Thus, when 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 typically assumed in the literature for such topological channels [26][27][28]. Since this can have important consequences in the quantum optical behaviour, we derive a more accurate effective description of these modes for ϕ ≪ 1, showing the edge-mode dispersion for the left-localized states emerging from the l-th Landau Level approximates by [62]: for (k y − k l (ϕ)) ∈ (−π, 0). Here, ω LL (ϕ)/J ≈ −4 + 4πϕ l + 1 2 − (πϕ) 2 (l 2 + l + 1/2) is the energy of the l th -Landau level [60], a l (ϕ) is the effective curvature of the edge modes which we extract from numerical fittings [62], and which converges to a l (ϕ → 0) ∼ 0.6, and k l (ϕ) is the momentum resonant to the minimum edge mode energy. Besides, the spatial distribution of these modes along the X direction is Ψ l (x) = 2/λ l (ϕ)e −x/λ l (ϕ) , where λ l grows as ϕ → 0 as expected, since in ϕ = 0 we should recover the delocalized bath eigenstates of the standard square lattice model.

III. SPONTANEOUS EMISSION FEATURES IN MULTI-MODE WAVEGUIDE QED
Let us now see how coupling emitters to the edge of the HH lattice leads to several unique phenomena. First, let us note that since the emitter probes the system at fixed frequency, ω 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 defined by [64]: being |E B ⟩ the bath eigenstates, i.e., H B |E B ⟩ = E B |E B ⟩, for the considered configuration. In Fig. 2(b), we plot in blue solid line the expected Γ(ω e ) for an emitter coupled to the edge of a HH lattice for ϕ = 1/25 as a function of ω 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 ω 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. As shown in Ref. [62], 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 Γ l (ω e ), that reproduces the non-linear dependence with the frequency, i.e.
As expected, the total decay is then obtained by summing the contributions of the active channels, i.e., Γ = l Γ l .
A more remarkable feature of these topological multimode waveguide scenario is what occurs with the spontaneously emitted photons when the emitters are driven. Let us first assume a perfect π-pulse driving which pre-pare the system in the state |Ψ 0 ⟩ = |e⟩ ⊗ |vac⟩, with |vac⟩ being the bath state with no photons. When the laser is switched off, the whole system evolves according to the total Hamiltonian |Ψ(t)⟩ = e −iHt |Ψ 0 ⟩, leading eventually to a single-photon wavepacket state as it occurs in other quantum optical setups [65]. However, in this case such wavepackets have unique features which we illustrate in Fig. 3 for an emitter coupled to the edge of a HH lattice with ϕ = 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 T J = 200. There, we observe that 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, 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 ′ |T ⪆ Γ −1 l + Γ −1 l ′ , that is, that the separation between the wavepackets is larger than their intrinsic linewidth (Γ −1 l ′ ).
When this separation occurs, one can say that our emitter has generated single-photon entangled states [66] between orthogonal time-bins T l . Defining |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 Using that, one can calculate the entanglement entropy [67], E(Ψ ph ), of the asymptotic state in the different band-gaps [62] whose result is shown in Fig. 2(c) and compared with the one of a perfect W -entangled state [68] in black dashed line. There, we observe how indeed the entanglement entropy of the emitted state indeed approximates that of maximally-entangled state. Note that in more conventional quantum optical setups [52,[69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84] where such time-bin entanglement is generated, it is required to combine superpositions in multilevel emitters and multiple-drivings, while here already with a single π-pulse, a W -type [68] entangled structure appears due to the multi-mode nature of the waveguide.
Applying several π-pulses one can obtain more complex multi-photon states. For example, if a second π-pulse is applied at a time T D before the excitation from the first π-pulse decays completely, the emission of the second photon is correlated with the first. As shown in Ref. [52], this ends up generating two-photon Bell-like states in the photon number basis where a † E(L) represents the photon operator emitted from the first (second) π-pulse. Compared to Ref. [52], in the topological multi-mode setups the single-photon wavepackets already have an internal superposition structure, a † E(L) ∝ l c E(L),l a † l |vac⟩, with a † l being the effective operator associated to the photon emitted in the l-th topological channels. Therefore, its multi-photon structure will be much richer, and depends on the interplay between the pulse delay T D , the global Γ and individual Γ l decay times, and the asymptotic time where it is measured. Let us emphasize that the non-gaussian character of these states is a consequence of the strongly interacting character of the emitters, and could never be obtained in classical setups.

IV. MODE SELECTIVITY VIA NON-LOCAL COUPLINGS
Finally, let us show how to make the emitters interact selectively with one of the resonant channels [62]. The key idea is to couple the emitter with more than a single lattice site, as it can be done with giant atoms [53][54][55][56][57]. 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φ i.e. g (0,y) = g and g (0,y+1) = ge iφ . In that case, the k-dependent lightmatter coupling reads and thus, vanishes at k y = k e if we choose φ = π − k e . This was the key idea introduced in single mode waveguide QED setups [85][86][87] 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. In Fig. 4 we show a proofof-principle realization of that idea for a situation when the emitter is resonant to two edge modes. In Fig. 4(a), we plot the energy spectrum for a lattice with ϕ = 1/12 and 65 × 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 , 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. 4(b), and emit in a two-mode fashion respectively. We prove mode-selectivity using non-local couplings that in momentum space are of the form |Gφ 0 (ky)| 2 and |Gφ 1 (ky)| 2 represented as dashdotted lines, which units are indicated in the right vertical axis, and are analytically described in Eq. (8). Note that the orange (blue) line, representing |Gφ 0 (ky)| 2 (|Gφ 1 (ky)| 2 ) vanishes at k  Fig. 4(e) where we plot a snapshot at time T J = 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 lightmatter coupling situation. In Sup. Material [62], 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 πphase [88]. Thus, if multiple photons are sent in the different channels, such mode selectivity can lead to different phases between the topologically-protected photons. This can be a resource for generating photon gates among topologically-protected photons by adapting existing protocols [89].
CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM), the CSIC Research Platform on Quantum Technologies PTI-001 and from Spanish project PID2021-127968NB-I00 (MCIU/AEI/FEDER, EU). AGT also acknowledges support from a 2022 Leonardo Grant for Researchers and Cultural Creators, BBVA Foundation. During the writing process of this article, another work studying the coupling of emitters to the edges of two-dimensional photonic insulators appeared [106], although focused only in the single-edge mode scenario.

VI. APPENDIX
In this Appendix, we provide the details of the calculations supporting the main manuscript. In Section A, 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 B, we focus on the edge modes, and show how their features can be captured by a simple phenomenological model. In Section C, 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.

Appendix A: 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 [59][60][61]. Denoting by a ( †) r , the annihilation (creation) operators of the bosonic mode at site r = (x, y) of the lattice, the HH Hamiltonian can be written (setting ℏ = 1 along the manuscript): (A1) 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 ϕ acquired in the nearest-neighbour hopping J along the Y direction. Note that we have also dismissed local cavity energy terms ω a r a † r a r , assuming all cavity energies to be the same along the lattice ω r = ω a , and therefore considering it as the energy reference of the problem setting ω a ≡ 0. Then, there are three magnitudes that determine the shape of the spectrum and eigenmodes of the system, that are, J, ϕ, and the system size L x ×L y .
Let us start analyzing the spectrum of the system imposing periodic boundary conditions along the X and Y direction, so that momentum k = (k x , k y ) is a good quantum number running over the following values k α = −π, −π + 2π/L α , . . . , π − 2π/L α . As we said in the main text, we will restrict to rational values of ϕ = 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 × 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: where α is the index running over the different bands, and c α,k are the operators describing the eigenstates of the bath for a given band α and momentum k. In Fig. A.1(a), we plot an example of that bulk spectrum for ϕ = 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 [59][60][61], and which describe a cyclic motion of lattice excitations whose radius in lattice constant units is given by the magnetic length, l B = 1/ √ 2πϕ. Thus, in order for these Landau level picture to survive, it is required that such orbit fits in the lattice, i.e., l B ≪ L α . When this condition is not satisfied, the spectrum tends to the square-lattice tight-binding spectrum ω(k) = ω a − 2J(cos(k x ) + cos(k y )). This transition can be observed in the bath density of states (DoS), defined as where the sum is performed over all the bath eigenenergies E B . Since the Dirac-δ has only mathematical sense in the continuum limit, for the finite systems we consider we approximate it by a Gaussian distribution of width θ, that is: Using that trick, we plot the DoS of the Harper-Hofstadter Hamiltonian in Fig. A.1(c), as a function of ϕ 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 ω 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 nonzero 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 [60,61]: 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 ϕ = 1/q, it results that the first [(q − 1)/2] bands (letting [·] 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 ϕ = 4/9 leads to a lowest energy band with C 0 = 2, and ϕ = 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. A.1(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. A.1(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.
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 |k y , β⟩ = E β (k y ) |k y , β⟩. Projecting their wavefunction into their spatial coordinates ⟨x, y|Ψ β,ky ⟩ = e ikyy ψ β,ky (x), one can find that ψ β,ky (x) satisfies the Harper equation: and it features a localized shape. To make it more evident, we define a localization parameter for each eigenstate as follows: which features a maximum ±1 value when localized in left/right edge, and 0, when it is delocalized. We codify the value of that parameter in Fig. A.1(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 C.
Appendix B: Effective edge mode description as a multi-mode waveguide As shown in Fig. A.1(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 [37,38]. Generally, such topological have been described within linear approximations [26][27][28]. However, from Fig. A.1 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.

Situation with ϕ = 1/q
Let us start with the situation ϕ = 1/q that we consider along the main text and in Fig. A.1. 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: where the ω l (ϕ) can be found approximately in the perturbative limit [60] as: whereas a l (ϕ) and k l (ϕ, L) are fitting parameters that depend on both the effective flux ϕ and l-th edge mode considered, although not on system size as long as l B ≪ 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. B.2(a) the evolution of the curvature of the edge modes, a l (ϕ), as a function of ϕ for the three lowest-energy edge modes in different colors. There, we see how for big fluxes, ϕ, the curvature of the modes differ significantly, whereas for small fluxes, they converge to a value a l (ϕ) ≈ 0.6. Regarding the value of the momentum cut-off k l (ϕ, L), we find that there is a linear dependence with ω l (ϕ), i.e., k l (ϕ) = β l ω l (ϕ)(mod2π). To illustrate that, in Fig. B.2(b) we plot one against each other for several fluxes ranging from ϕ ∈ [ 1 60 , 1 12 ] in the different markers, together with a dotted line which indicates the result of the fitting. To discuss the dependence of k l on the system size L, we may rewrite Eq. A7 as where the potential V ϕ (x) is given by: We observe that varying L modifies the boundary condition of the Harper equation in the right edge. However, if ϕ = 1/q and L is modified in increased or decreased in nq sites (n ∈ N), the boundary condition remains, leading to k l (ϕ, L) − k l (ϕ, L ′ ) = 2πϕ(L − L ′ ). In Fig. B.3(a), we show the comparison between the exact diagonalization results and our effective description for a particular value of ϕ, showing indeed an excellent agreement for the lowest edge states dispersion. In Fig. B.3(b), we make a more quantitative assessment on the quality of the model by defining an error parameter: where Ω 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, ω exact,l (k y ), and N k = |Ω| is the number k y modes within that region. In Fig. B.3(b), we plot ε l as a function of ϕ and compare the accuracy between a linear fit, i.e., ω eff,l (k y ) ∝ k y (empty markers) and the quadratic fit of Eq. (B1) (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. (A7) that when E β (k y ) lies within a band-gap, the spatial wavefunction along X will be exponentially localized, i.e., ψ β (x) = 2/λ l (ϕ)e −x/λ l (ϕ) , where λ l (ϕ) will depend on both the energy level l and the flux ϕ. In Fig. B.2(c) we also plot its dependence, showing that higher-energy edge modes are less localized, and also that increasing the value of ϕ yields to higher de-localization.
2. Other situations ϕ ̸ = 1/q, q ∈ N Along this work, we have restricted to magnetic fluxes of the form ϕ = 1/q, with q ∈ N, due to the topological features thoroughly discussed along the manuscript, namely the sequence of [q/2] lowest bands with Chern number C = −1. This structure does not prevail if ϕ does not fit this form. To illustrate this, we consider the cases of ϕ = 4/9 and 5/14 in Fig. B.4(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 ϕ in panels (a) and (b), where we see that the lowest bandgap 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 ϕ = 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 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 ky, which can also be captured by the effective model although we did not explicitly show that. (c) Fitting error, as defined in Eq. (B5), 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.
of this setup to obtain qualitatively different photonic wavepackets.

Appendix C: 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 = ω e σ 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) of the main text. In general, we will consider the most standard local e /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.
light-matter couplings given by: where g is the coupling strength of the bath, r e the position of the cavity mode the emitter couples to, and σ αβ = |α⟩ ⟨β| the dipole operator of the optical emitter transition that couples the photonic bath. In Section C 6, however, we will consider the non-local couplings that can be engineered with giant atoms [53][54][55][56][57], as a way of selectively coupling some of the topological edge modes.

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 π-pulse. If one assumes that the bath has initially no excitations, this state, |Ψ 0 ⟩ = |e⟩ ⊗ |vac⟩ B , can only evolve into an state of the form: because the full light-matter Hamiltonian: H = H S + H B + H I , conserves the number of excitations N exc = σ ee + r a † r a r since [H, N 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 ≈ e −Γt , with Γ being the expected Markovian decay rate given by Fermi's Golden rule [107]: 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: where |0⟩ = |g⟩|vac⟩. We numerically compute the LDoS in a similar fashion as we did for the DoS, defining a 'smoothed' Dirac delta function f θ (E −E B ), that we take to be a Gaussian distribution as expressed in Eq. (A5), and computing the LDoS as From Eq. (C3) we read that the shape of the LDoS determines the shape of the expected Markovian decay rate Γ(ω e ). In Fig. C.5(a) we plot the Markovian decay rate at the emitter position calculated with exact diagonalization, for several values of ϕ. We observe that smaller ϕ's are associated with narrow band-gaps, situation in which the LDoS 'quasi-quantized' behaviour leads to plateaus in Γ(ω e ). In Fig. C.5(b), we represent with a solid line the expected decay rate Γ e as a function of ω e for ϕ = 1/12, computed using exact diagonalization, and compare it with the markers that represent the expected decay rate computed as follows: where v l,g (ω e ) is the group velocity of the l-th mode at the emitter energy, that is, v g,l = ∂ ky ω eff,l (k y )| ky=ke , with ω eff,l (k e ) = ω e . Thus, the total decay rate in a given band-gap will be given by Γ(ω e ) = l Γ l (ω 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. C.5(c) we show the quantum emitter dynamics of an initially excited emitter for the range of ω e depicted in Fig. C.5(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/ √ ω − ω edge dependence associated to one-dimensional quadratic band-edge dispersions.
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 θ associated to the auxiliary function f θ (E − E B ): if we probe the LDoS of an energy |E − ω l | < θ from below, the LDoS will count some states above ω l , softening the transition. This issue is unavoidable: if we try to have an arbitrarily small value of θ, 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 θ will make the approximation states in Eq. (C6) 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. C.6.

Photonic spontaneous emission patterns
As we see in the Fig.3 of the main text, the nonuniform 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 [66], that when combined with sequential generation methods [52,[69][70][71][72][73][74][75][76][77][78][79][80][81][82][83][84] 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 ′ |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/Γ l in lattice constant units [108].
Thus, one can define a parameter R l,l ′ : 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 ∝ Γ −1 l +Γ −1 l ′ and the pulse separation induced by the difference between group velocities |v g,l − v g,l ′ |. When R ll ′ T ∼ 1, it is expected that the modes l and l ′ are fully resolvable at time T . In Fig. C.7 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 ω e for a given bath configuration with ϕ = 1/40. From Fig. C.7(c) we observe that the pulses associated with modes 0 and 1 will be fully resolvable at T J ∼ 10 3 . As expected, the best conditions to achieve resolution will occur when the quantum emitter energy ω e is slightly above the Landau level energy ω 1 : in this case, |v g,0 − v g,1 | will be maximum, which favours resolution.
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, γ dis in ω l (k) ≈ ω e + v g,l (k − k e ) + γ dis (k − k e ) 2 /2. In particular, it is well known that a wavepacket with an initial broadening σ 0 propagating in such a non-linear dispersive channel will have an increasing size growing with: In the case of spontaneous emission of a quantum emitter into the l th -topological channel of a HH lattice, we can take σ 0 as the inverse of the decay rate onto such mode, Γ −1 l . The pulse width will then evolve as: From this evolution equation, it follows that the role of dispersion will be negligible at a certain time T as long as Taking g/J ∼ 0.1, we find that Γ l /J ∼ 10 −3 . On the other hand, we can estimate γ dis from our effective theory: in particular, we will have that γ dis ∼ 2a l (ϕ)(k y − k l (ϕ)), 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 T J ∼ 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 ′ T ∼ 1.

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: where δω r is a random variable uniformly distributed along the interval (−σ, σ), where σ 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 strengths σ up to the order of E W . In Fig. C.8, we first analyse the effect of disorder qualitatively in two key features: the LDoS defined in Eq. (C4) and photon emission and propagation. In Fig. C.8(a) we plot the Markovian decay rate, computed from exact diagonalization, for an emitter coupled at the boundary of a HH lattice for ϕ = 1/12 for different values of σ averaged among different realizations of disorder. We observe that the laddered structure is preserved for values of σ 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. C.8(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 σ. For weak disorder i.e. small σ compared to E W , we observe the same light-cone configuration as in the σ = 0 case. We only start to witness light-cone distortion for values of σ comparable to the size of the band-gap that is resonant to the emitter energy. Now, let us make a more quantitative description of the impact of disorder in the propagation of the edge modes in such multi-mode scenario. For that, we use the method proposed in Ref. [92], where they quantify the robustness of photon transport by partitioning the Hilbert space into bulk and edge modes, and observing the edge mode content of the photonic state when defects or local disorder are present.
To distinguish between edge and bulk modes of the spectrum, we compute the localization properties of each eigenstate of the system using the Inverse Participation Ratio (IPR), defined as follows: If a wavefunction is spread in a discrete space of N sites, its IPR will be of the order of N . Then, we would expect an IPR of the order of L 2 for bulk modes (which spread along the whole lattice), and of L for edge states (which are localized along one dimension). In Fig. C.9(a) we observe these expected scaling relations. The implication of this is that larger lattices leads to larger distinction of the IPR of the bulk and edge modes, thus allowing a clear bi-partition by considering a given treshold ε in the IPR. After defining the partition of the Hilbert space, we quantify the edge mode content as the associated norm of the emitted wavefunction at a certain time, |ψ ph ⟩, projected into the edge-mode subspace: where the sum is performed over all lattice eigenstates |Ψ⟩ whose IPR is below a given treshold. In Fig. C.9(b) we consider a single square defect of increasing size |D|, and study how much the excitation spreads out of the edge modes after impinging with it. We compare the situation with two different system sizes and find that the edge mode content in both cases is practically constant. Furthermore, we see that this robustness measure is higher for larger lattices, where the protection of the edge modes is better. In the inset, we represent a snapshot of the spontaneous emission dynamics where we see that transport is highly protected even for large defects.
On top of that, we also consider a different disorder situation, that is, a random energy disorder over the lattice sites with a normal distribution of width σ (see Fig. C.9c). In this case, we find more instructive to project directly to the edge mode subspace of the clean system, since bulk modes undergo Anderson localization, and hence the IPR is no longer valid to differentiate bulk and edge states. In Eq. C14, this implies doing the sum over the edge eigenstates |Ψ⟩ of the pristine Hamiltonian. In this case, we see that the edge mode content decreases for higher disorder strengths, revealing lack of protection for values of σ of the order of ω e −ω bulk , which is ≈ 0.57J in the situation represented in the figure. Again, we also observe that larger sizes favour protection.

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, ρ(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: where κ r and Γ * 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: In Fig. C.10(a), we plot its effect in the predicted Markovian decay rates for increasing values of κ r = κ and fixing Γ * . 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 Γ * ̸ = 0. The most significant effect of κ is to generate a finite propagation length of the photon modes, as observed in Fig. C.10(b-c) where we plot an example of the spontaneously emitted photons without and with dissipation, respectively.

Single-photon entanglement generation
In the manuscript, we show that a spontaneously emitted photon from a locally-coupled quantum emitter will be distributed over the multiple boundary modes, and point out that this will lead to a single-particle entangled states between the different channels [66]. In this section, we give the analytical expression of such photonic state within the Markovian approximation, and quantify its entanglement.
Let us first determine the asymptotic photonic state. After a time t ≫ Γ −1 , the emitter population will be negligible and the state will be purely photonic i.e. e −iHt |e⟩ ⊗ |vac⟩ ≈ |g⟩ ⊗ |Ψ ph ⟩. Assuming a multimode waveguide of boundary modes and letting A l,k be the annihilation operator for a photon in the l th -mode with quasi-momentum k, we can rewrite the light-matter Hamiltonian as: with g l,k = Γ l /(2π), being Γ l the expected Markovian decay rate in the Γ l mode. Since the Hamiltonian preserves the number of excitations, if the emitter is initially excited, one can write the quantum state of the emitter+ bath system at any time as follows: |Ψ(t)⟩ = C e (t) |e⟩ ⊗ |vac⟩ + k C k (t) |g⟩ ⊗ a † k |vac⟩. Applying, the time-dependent Schrödinger equation to this wavefunction one can obtain the following set of equa- Formally integrating the second equation, inserting in the first equation and applying Markov approximation, we find that the emitter follows exponential decay C e (t) ≈ e −Γt/2 , where Γ = l Γ l . Then: where we assumed t ≫ Γ −1 . We conclude that the asymptotic photonic state |Ψ ph ⟩ can be written as For a single mode, we recover the well-known expression for the asymptotic photonic state of the emitted light from an emitter coupled to a single-mode waveguide [108]. From this expression, we observe that the contribution of the l th -mode to the asymptotic photonic state population is weighted by √ Γ l . Defining a multimode basis {|1 l ⟩}, where |1 l ⟩ contains one photon in the l th -mode, and defining C l = Γ l /Γ we can write the photonic state whose entanglement we want to characterize as: Here, we are consciously neglecting the contribution of the different phases between the modes, since in principle they can be corrected by local operations, and thus should not contribute to the entanglement measures. One way of characterizing the entanglement of this class of states is by calculating the entanglement entropy [67] between the different topological channels. The entanglement entropy of a given state Ψ for a bipartition of the Hilbert space H A ⊗ H B reads [67]: where S(ρ) = −Tr(ρ log ρ) and Tr A/B denotes the partial trace along the A/B subspace. In our case, we define as P l to the bi-partition separating the l th -mode from the rest. Assuming a bi-partition P l , we will take our Hilbert space as H l ⊕ (⊕ l ′ ̸ =l H l ′ ), and denote by |10⟩ (|01⟩) the state with one excitation in the l th -mode and zero in the rest (and viceversa). From Eq. (C21), and using that notation, we then know we can always write our state as: Denoting by p = |C l | 2 the probability of measuring the photon in the l th -mode, the entanglement entropy of this state can be readily computed as In Fig. C.11, we present the entanglement entropy for all possible bi-partitions in each spectral band-gap, in comparison with the entanglement entropy of a W -state with a number of qubits equal to the number of resonant edge modes in each band-gap (dashed lines). The latter can be calculated by imposing p = 1/N modes in Eq. (C24). In the figure, we observe that the entanglement of our generated state is very similar to the one obtained for a maximally-entangled W state in all the band-gaps. Let us note that this calculation is assuming that the emitter only decays through the edge modes. This will be a good approximation as long as the population of the bulk modes remains small. The latter can be perturbatively approximated by g 2 /(ω e − ω bulk ) 2 , which is what we plot in dashed orange curves in Fig. C.11. Thus, the calculation of the entanglement entropy will be valid outside of the yellow regions where the population into the Landau levels become negligible. To make a more quantitative estimation of the impact of the bulk-modes into the entanglement generated, one can use more sophisticated entanglement measures such as the negativity [109] that can be effectively computed for mixed states. However, we believe that for the sake of illustration the entanglement entropy between the different channels represent a more intuitive witness.

Mode-selectivity using non-local couplings
Local light-matter Hamiltonians, as the ones considered in Eq. (C1), 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: H I = gσ eg a re + H.c. = ky g L y σ eg a (xe,ky) e −ikyye + H.c. , (C25) where: a (xe,ky) = 1 L y y a (xe,y) e −ikyy .
This makes that when coupling an emitter with an optical transition of a given energy, ω 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: H I = y g y σ eg a (x,y) + H.c. = ky G(k y )σ eg a (xe,ky) + H.c. , with: G(k y ) = 1 L y y g y e −ikyy .
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 [85] and two-dimensional [54] singlemode waveguides to generate chiral, one-directional emission in baths with isotropic bath dispersions. In our topological multi-mode waveguides, for fixed ω e , Let us now illustrate a potential method to achieve that selectivity inspired by the results in Refs. [54,85] that only requires coupling to N c +1 cavities if one wants to cancel the coupling to N c resonant momenta k (l) e . Let us start by the case where we want cancel only the coupling to a single resonant momenta k (l) e . For that, let us assume that the emitter couples to neighbouring sites with the same amplitude, but a relative complex phase φ l , that is: g y = g, g y+1 = ge iφ l . Using these values, it can be shown from Eq. (C28) that: |G φ l (k y )| 2 = 2g 2 L y (1 + cos(k y + φ l )) .
Thus, if we want to make it zero for certain k y = k (l) e , one just have to choose a relative phase: φ l,e = π − k (l) e . 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 where φ δ = π − k 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 G φ0,φ1 (k y ) ∝ 1 + (e iφ0 + e iφ1 )e iky + e i(φ0+φ1) e i2ky , which can be achieved with real-space couplings of the form g (0,y) = g, g (0,y+1) = ge iφ0 + ge iφ1 and g (0,y+2) = ge i(φ0+φ1) .