Research ArticlePHYSICS

Unconventional quantum optics in topological waveguide QED

See allHide authors and affiliations

Science Advances  26 Jul 2019:
Vol. 5, no. 7, eaaw0297
DOI: 10.1126/sciadv.aaw0297


The discovery of topological materials has motivated recent developments to export topological concepts into photonics to make light behave in exotic ways. Here, we predict several unconventional quantum optical phenomena that occur when quantum emitters interact with a topological waveguide quantum electrodynamics bath, namely, the photonic analog of the Su-Schrieffer-Heeger model. When the emitters’ frequency lies within the topological bandgap, a chiral bound state emerges, which is located on just one side (right or left) of the emitter. In the presence of several emitters, this bound state mediates topological, tunable interactions between them, which can give rise to exotic many-body phases such as double Néel ordered states. Furthermore, when the emitters’ optical transition is resonant with the bands, we find unconventional scattering properties and different super/subradiant states depending on the band topology. Last, we propose several implementations where these phenomena can be observed with state-of-the-art technology.


Even though the introduction of topology in condensed matter was originally motivated to explain the integer quantum Hall effect (1), its implications were more far-reaching than expected. On a fundamental level, topological order resulted in a large variety of new phenomena, as well as new paradigms for classifying matter phases (2). In practical terms, topological states can be harnessed to achieve more robust electronic devices or fault-tolerant quantum computation (3). This spectacular progress motivated the application of topological ideas to photonics, for example, to engineer unconventional light behaviors. The starting point of the field was the observation that topological bands also appear with electromagnetic waves (4). Soon after that, many experimental realizations followed using microwave photons (5), photonic crystals (6, 7), coupled waveguides (8) or resonators (911), exciton-polaritons (12), or metamaterials (13), to name a few [see (14) and references therein for an authoritative review]. Nowadays, topological photonics is a burgeoning field with many experimental and theoretical developments. Among them, one of the current frontiers of the field is the exploration of the interplay between topological photons and quantum emitters (QEs) (1517).

Here, we show that topological photonic systems cause a number of unprecedented phenomena in the field of quantum optics, namely, when they are coupled to QEs. We analyze the simplest model consisting of two-level QEs interacting with a one-dimensional (1D) topological photonic bath described by the Su-Schrieffer-Heeger (SSH) model (18) (see Fig. 1). When the QE frequency lies between the two bands (green region in Fig. 1B), we predict the emergence of chiral photon bound states (BSs), that is, BSs that localize to the left/right side of the QEs depending on the topology of the bath. In the many-body regime (i.e., with many emitters), those BSs mediate tunable, chiral, long-range interactions, leading to a rich phase diagram at zero temperature, e.g., with double Néel ordered phases. Furthermore, when the QEs are resonant with the bands (blue regions in Fig. 1B), we also find unusual dissipative dynamics. For example, for two equal QEs separated a given distance, we show that both the super/subradiance conditions (19) and the scattering properties depend on the parameter that governs the bath topology, even though the energy dispersion ω(k) is insensitive to it. This might open avenues to probe the topology of these systems in unconventional ways, e.g., through reflection/transmission experiments.

Fig. 1 System schematic.

(A) Schematic picture of the present setup: Ne two-level QEs interact with the photonic analog of the SSH model. This model is characterized by having alternating hopping amplitudes J(1 ± δ), where J defines their strength, while δ, the so-called dimerization parameter, controls the asymmetry between them. The interaction with photons (in transparent red) induces nontrivial dynamics between the emitters. (B) Bath’s energy bands for a system with a dimerization parameter ∣δ ∣ = 0.2. The main spectral regions of interest for this manuscript are the middle bandgap (green) and the two bands (blue).


The system that we study in this manuscript is shown in Fig. 1A: One or many QEs interact through a common bath, which behaves as the photonic analog of the SSH model (18). This bath model is described by two interspersed photonic lattices A/B of size N with alternating nearest-neighbor hoppings J(1 ± δ) between their photonic modes. Assuming periodic boundary conditions and defining V=(ak,bk), the bath Hamiltonian can be written in momentum space as HB=kVH˜B(k)V, with (setting ℏ = 1)H˜B(k)=(ωaf(k)f*(k)ωa)(1)where f(k) = − J[(1 + δ) + (1 − δ)eik] = ω(k)eiϕ(k) [with ω(k) > 0] is the coupling in momentum space between the A (B) modes, ak=jajeikj/N (bk=jbjeikj/N). Here, aj/bj (aj/bj) are the creation (annihilation) operators of the A/B photonic mode at the jth unit cell. We assume that the A/B modes have the same energy, ωa, that from now on will be the reference energy of the problem, i.e., ωa ≡ 0. This Hamiltonian can be easily diagonalized introducing the eigenoperators, uk/lk=[±ak+eiϕ(k)bk]/2, as HB=kωk(ukuklklk), leading to two bands with energy±ω(k)=±J2(1+δ2)+2(1δ2)cos(k)(2)Let us now summarize the main bath properties:

(1) The bath has sublattice (chiral) symmetry (18), such that all eigenmodes can be grouped in chiral symmetric pairs with opposite energies. Thus, the two bands are symmetric with respect to ωa, spanning [−2J, −2 ∣ δ ∣ J] (lower band) and [2 ∣ δ ∣ J,2J] (upper band). The middle gap is 4 ∣ δ ∣ J, such that it closes when δ = 0, recovering the normal 1D tight-binding model.

(2) This bath supports topologically nontrivial phases, belonging to the BDI class in the topological classification of phases (20). More concretely, both bands can be characterized by a topological invariant, the Zak phase Z (18), such that Z=0 corresponds to a trivial insulator, while Z=π implies a nontrivial insulator. For the parameterization we have chosen, this occurs for δ > 0 and δ < 0, respectively. Notice that for an infinite system (i.e., in the bulk), this definition depends on the choice of the unit cell, and the role of δ can be reversed by shifting the unit cell by one site. In the bulk, the band topology manifests in the fact that one cannot transform from one phase to the other without closing the gap (as long as the symmetry is preserved).

(3) With finite systems, however, the sign of δ determines whether the chain ends with weak/strong hoppings, which leads to the appearance (or not) of topologically robust edge states (21).

Now, let us lastly describe the rest of the elements of our setup. For the Ne QEs, we consider that they all have a single optical transition g-e, with a detuning Δ with respect to ωa, and they couple to the bath locally. Thus, their free and interaction Hamiltonian readHS=Δm=1Neσeem(3)HI=gm(σegmcxm+H.c.)(4)where cxm ∈ {axm, bxm} depends on the sublattice and the unit cell xm at which the mth QE couples to the bath. We use the notation σμνm=μmν, μ, ν ∈ {e, g}, for the mth QE operator. We highlight that we use a rotating-wave approximation, such that only excitation-conserving terms appear in HI.


In the next sections, we study the dynamics emerging from the global QE-bath Hamiltonian H = HS + HB + HI using several complementary approaches. When one is only interested in the QE dynamics and the bath can be effectively traced out, the following Born-Markov master equation (22) describes the evolution of the reduced density matrix ρ of the QEsρ̇=i[ρ,HS]+in,mJmnαβ[ρ,σegmσgen]+n,mΓmnαβ2[2σgenρσegmσegmσgenρρσegmσgen](5)

The functions Jmnαβ,Γmnαβ, which ultimately control the QE coherent and dissipative dynamics, respectively, are the real and imaginary parts of the collective self-energy Σmnαβ(Δ+i0+)=JmnαβiΓmnαβ/2. This collective self-energy depends on the sublattices α, β ∈ {A, B} to which the mth and nth QEs couple, respectively, as well as on their relative position xmn = xnxm. For our model, they can be calculated analytically in the thermodynamic limit (N → ∞) yieldingΣmnAA/BB(z)=g2z[y+xmnΘ+(y+)yxmnΘ(y+)]z44J2(1+δ2)z2+16J4δ2(6)ΣmnAB(z)=g2J[Fxmn(y+)Θ+(y+)Fxmn(y)Θ(y+)]z44J2(1+δ2)z2+16J4δ2(7)where Fn(z) = (1 + δ)zn + (1 − δ)zn + 1∣, Θ±(z) = Θ(±1 ∓ ∣ z ∣), Θ(z) is Heaviside’s step function, andy±=z22J2(1+δ2)±z44J2(1+δ2)z2+16J4δ22J2(1δ2)(8)

However, since we have a highly structured bath, this perturbative description will not be valid in certain regimes, e.g., close to band edges, and we will use resolvent operator techniques (23) or fully numerical approaches to solve the problem exactly for infinite/finite bath sizes, respectively. Since those methods were explained in detail in other works, here, we focus on the results and leave the details in the Supplementary Materials.


In this section, we assume that the QEs are in the bandgap regime; that is, their transition frequency lies outside of the two bands of the photonic bath. From here on, we only discuss results in the thermodynamic limit (when N → ∞) such that the edge states (21) play no role in the QE dynamics. We refer the interested reader to (24) and the Supplementary Materials to see some of the consequences the edge states have on the QE dynamics.

Single QE: Dynamics

Let us start considering the dynamics of a single excited QE, i.e., ∣ψ(0)〉 = ∣ e〉 ∣ vac〉, where ∣vac〉 denotes the vacuum state of the lattice of bosonic modes. Since H conserves the number of excitations, the global wave function at any time readsψ(t)=[Ce(t)σeg+j=1Nα=a,bCj,α(t)αj]gvac(9)

In both perturbative and exact treatments, the dynamics of Ce(t) can be shown [see (23) and the Supplementary Materials] to depend only on the single-QE self-energyΣe(z)=g2z sign(y+1)z44J2(1+δ2)z2+16J4δ2(10)obtained from Eq. 6 defining Σe(z)ΣnnAA(z). From here, we can already extract several conclusions: (i) Σe(z) is independent of the sign of δ, which means that the spontaneous emission dynamics is insensitive to the topology of the bands; (ii) perturbative approaches, like the Born-Markov approximation of Eq. 5, predict an exponential decay of excitations at a rate Γe(Δ) = −2ImΣe(Δ + i0+), which is strictly zero in the bandgap regime. Thus, one expects that the excitation remains localized in the QE at any time. However, in Fig. 2, we compute the exact dynamics Ce(t) for several δ’s and observe that this perturbative limit is only recovered in the limit of ∣δ ∣ → 1. On the contrary, when ∣δ ∣ ≪ 1 and δ ≠ 0, the dynamics displays fractional decay and oscillations. As it happens with other baths (25), the origin of this dynamics stems from the emergence of photon BSs, which localize around the QEs (2628). However, the BSs appearing in the present topological waveguide bath have some distinctive features with no analog in other systems and therefore deserve special attention.

Fig. 2 Single-QE dynamics.

Probability to find the emitter in the excited state, ∣Ce(t)∣2, for different values of ∣δ∣. The other parameters are Δ = 0 (middle of the bandgap) and g = 0.4J. As the bandgap closes, i.e., δ → 0, the decay becomes stronger. Dashed lines mark the value of ∣Ce(t → ∞ )∣2 = [1 + g2/(4J2 ∣ δ ∣)]−2.

Single QE: BSs

The energy and wave function of the BSs in the single-excitation subspace can be obtained by solving the secular equation H ∣ ΨBS〉 = EBS ∣ ΨBS〉, with EBS lying out of the bands, and ∣ΨBS〉 in the form of Eq. 9, but with time-independent coefficients. Without loss of generality, we assume that the QE couples to sublattice A at the j = 0 cell. After some algebra, one can find that the energy of the BS is given by the pole equation: EBS = Δ + Σe(EBS). Irrespective of Δ or g, there always exist three BS solutions of the pole equation (one for each bandgap region). This is because the self-energy diverges in all band edges, which guarantees finding a BS in each of the bandgaps (29, 30). The main difference with respect to other BSs (2630) appears in the wave function amplitudes, which readCj,a=gEBSCe2πππdkeikjEBS2ω2(k)(11)Cj,b=gCe2πππdkω(k)ei[kjϕ(k)]EBS2ω2(k)(12)where Ce is a constant obtained from the normalization condition that is directly related to the long-time population of the excited state in spontaneous emission. For example, in Fig. 2 where Δ = 0, it can be shown to be ∣Ce(t → ∞ )∣2 = ∣Ce4 = [1 + g2/(4J2 ∣ δ ∣)]−2.

From Eqs. 11 and 12, we can extract several properties of the spatial wave function distribution. On the one hand, above or below the bands (outer bandgaps), the largest contribution to the integrals is that of k = 0; thus, all the Cj have the same sign (see the left column of Fig. 3A, top and bottom row). In the lower (upper) bandgap, Cj of the different sublattices has the same (opposite) sign. On the other hand, in the inner bandgap, the main contribution to the integrals is that of k = π. This gives an extra factor (−1)j to the coefficients Cj (see Fig. 3A, middle row). Furthermore, the probability amplitudes of the sublattice where the QE couples to are symmetric with respect to the position of the QE, whereas they are asymmetric in the other sublattice; that is, the BSs are chiral. Changing δ from positive to negative results in a spatial inversion of the BS wave function. The asymmetry of the BS wave function is more extreme in the middle of the bandgap (Δ = 0). For example, if δ > 0, the BS wave function with EBS = 0 is given by Cj,a = 0 andCj,b={gCe(1)jJ(1+δ)(1δ1+δ)j,j00,j<0(13)whereas for δ < 0, the wave function decays for j < 0 while being strictly zero for j ≥ 0. At this point, the BS decay length diverges as λBS ∼ 1/(2 ∣ δ ∣) when the gap closes. Away from this point, the BS decay length shows the usual behavior for 1D baths λBS1/Δedge, with Δedge being the smallest detuning between the QE and the band-edge frequencies.

Fig. 3 BS properties.

(A) BS wave function for a QE placed at j = 0 that couples to the A sublattice; δ = 0.2 and g = 0.4J. Probability amplitudes Cj,a are shown in blue, while the Cj,b are shown in orange. The QE frequency is set to Δ = 2.2J (top row), Δ = 0 (middle row), and Δ = −2.2J (bottom row). The first column corresponds to the model without disorder, the second corresponds to the model with disorder in the couplings between cavities, and the third corresponds to the model with disorder in the cavities’ resonant frequencies. In both cases with disorder, the disorder strength is set to w = 0.5J. For each case, the value of the BS’s energy is shown at the bottom of the plots. (B) Inverse BS localization length for the two different models of disorder as a function of the disorder strength. Parameters: g = 0.4J and δ = 0.5. The dots correspond to the average value computed with a total of 104 instances of disorder, and the error bars mark the value of 1 SD above and below the average value (the blue curves are slightly displaced to the right for better visibility). Two cases are shown, which correspond to Δ ≃ 2.06J (triangles, outer bandgaps) and Δ = 0 (circles, inner bandgap). (C) Absolute value of the dipolar coupling for Δ = 0 and g = 0.4J; Markov, solid line; exact, dots. The insets show the shape of the BSs in the topological and the trivial phases. The situation for the BA configuration is the same, reversing the role of δ.

The physical intuition of the appearance of such chiral BS at EBS = 0 is that the QE with Δ = 0 acts as an effective edge in the middle of the chain or, equivalently, as a boundary between two semi-infinite chains with different topology. This picture provides us with an insight that is useful to understand other results of this manuscript: Despite considering the case of an infinite bath, the local QE-bath coupling inherits information about the underlying bath topology. One can show that this chiral BS has the same properties as the edge state, which appears in a semi-infinite SSH chain in the topologically nontrivial phase, for example, inheriting its robustness to disorder. To illustrate it, we study the effect of two types of disorder: one that appears in the cavities’ bare frequencies (diagonal), and another one that appears in the tunneling amplitudes between them (off-diagonal). The former corresponds to the addition of random diagonal terms to the bath’s Hamiltonian HBHB+j(εa,jajaj+εb,jbjbj) and breaks the chiral symmetry of the original model, while the latter corresponds to the addition of off-diagonal random terms HBHB+j(ε1,jbjaj+ε2,jaj+1bj+H.c.) and preserves it. We take the εν, j, ν = a, b,1,2, from a uniform distribution within the range [−w/2, w/2] for each jth unit cell. To prevent changing the sign of the coupling amplitudes between the cavities, w is restricted to w/2 < (1 − ∣ δ ∣) in the case of off-diagonal disorder.

In the middle (right) column of Fig. 3A, we plot the shape of the three BSs appearing in our problem for a situation with off-diagonal (diagonal) disorder with w = 0.5J. There, we observe that while the upper and the lower BS are modified for both types of disorder, the chiral BS has the same protection against off-diagonal disorder as a regular SSH edge state: Its energy is pinned at EBS = 0, and it keeps its shape with no amplitude in the sublattice to which the QE couples. On the contrary, for diagonal disorder, the middle BS is not protected anymore and may have weight in both sublattices.

Last, to make more explicit the different behavior with disorder of the middle BS compared to the other ones, we compute their localization length λBS as a function of the disorder strength w averaging for many realizations. In Fig. 3B, we plot both the average value (markers) of λBS1 and its standard deviation (SD) (bars) for the cases of the middle (blue circles) and upper (purple triangles) BSs. Generally, one expects that for weak disorder, states outside the band regions tend to delocalize, while for strong disorder, all eigenstates become localized [see, for example, (31)]. This is the behavior we observe for the upper BS for both types of disorder. However, the numerical results suggest that for off-diagonal disorder, the chiral BS never delocalizes (on average). Furthermore, the chiral BS localization length is less sensitive to the disorder strength w manifested in both the large initial plateau region and the smaller SDs compared to the upper BS results.

In summary, a QE coupled locally to an SSH bath (i) localizes a photon only on one side of the emitter depending on the sign of δ, (ii) with no amplitude in the sublattice where the QE couples to, and (iii) with the same properties as the topological edge states, e.g., robustness to disorder. As we discuss in more detail in the Supplementary Materials, the SSH bath is the simplest 1D bath that provides all these features simultaneously.

Two QEs

Let us now focus on the consequences of such exotic BS when two QEs are coupled to the bath. For concreteness, we focus on a parameter regime where the Born-Markov approximation is justified, although we have performed an exact analysis in the Supplementary Materials. From Eq. 5, it is easy to see that in the bandgap regime, the interaction with the bath leads to an effective unitary dynamics governed by the following HamiltonianHdd=J12αβ(σeg1σge2+H.c.)(14)

That is, the bath mediates dipole-dipole interactions between the QEs. One way to understand the origin of these interactions is that the emitters exchange virtual photons through the bath, which, in this case, are localized around the emitter. These virtual photons are nothing but the photon BS that we studied in the previous section. Thus, these interactions Jmnαβ inherit many properties of the BSs. For example, the interactions are exponentially localized in space, with a localization length that can be tuned and made large by setting Δ close to the band edges or fixing Δ = 0 and letting the middle bandgap close (δ → 0). Moreover, one can also qualitatively change the interactions by moving Δ to different bandgaps: For ∣Δ ∣ > 2J, all the Jmnαβ have the same sign, while for ∣Δ ∣ < 2 ∣ δ ∣ J, they alternate sign as xmn increases. In addition, changing Δ from positive to negative changes the sign of JmnAA/BB, but leaves JmnAB/BA unaltered. Furthermore, while JmnAA/BB are insensitive to the bath’s topology, the JmnAB/BA mimic the dimerization of the underlying bath, but allow for longer-range couplings. The most notable regime is again reached for Δ = 0. In that case, JmnAA/BB identically vanish, and thus, the QEs only interact if they are coupled to different sublattices. Furthermore, in such a situation, the interactions have a strong directional character; i.e., the QEs only interact if they are in some particular order. Assuming that the first QE at x1 couples to sublattice A, and the second one at x2 couples to B, we haveJ12AB={sign(δ)g2(1)x12J(1+δ)(1δ1+δ)x12if δx12>00if δx12<0Θ(δ)g2J(1+δ)if x12=0(15)

In Fig. 3C, we plot the absolute value of the coupling for this case computed exactly and compare it with the Markovian formula. Apart from small deviations at short distances, it is important to highlight that the directional character agrees perfectly in both cases.

Many QEs: Spin models with topological long-range interactions

One of the main interests of having a platform with BS-mediated interactions is to investigate spin models with long-range interactions (32, 33). The study of these models has become an attractive avenue in quantum simulation because long-range interactions are the source of nontrivial many-body phases (34) and dynamics (35), and are also very hard to treat classically.

Let us now investigate how the shape of the QE interactions inherited from the topological bath translates into different many-body phases at zero temperature as compared to those produced by long-range interactions appearing in other setups such as trapped ions (34, 35) or standard waveguide setups. For that, we consider having Ne emitters equally spaced and alternatively coupled to the A/B lattice sites. After eliminating the bath and adding a collective field with amplitude μ to control the number of spin excitations, the dynamics of the emitters (spins) is effectively given byHspin=m,n[JmnAB(σegm,Aσgen,B+H.c.)μ2(σzm,A+σzn,B)](16)denoting by σνn,α, ν = x, y, z, the corresponding Pauli matrix acting on the α ∈ {A, B} site in the nth unit cell. The Jmnαβ are the spin-spin interactions derived in the previous subsection, whose localization length, denoted by ξ, and functional form can be tuned through system parameters such as Δ.

For example, when the lower (upper) BS mediates the interaction, the Jmnαβ has negative (alternating) sign for all sites, similar to the ones appearing in standard waveguide setups. When the range of the interactions is short (nearest neighbor), the physics is well described by the ferromagnetic XY model with a transverse field (36), which goes from a fully polarized phase when ∣μ∣ dominates to a superfluid one in which spins start flipping as ∣μ∣ decreases. In the case where the interactions are long-ranged, the physics is similar to that explained in (34) for power-law interactions (∝1/r3). The longer range of the interactions tends to break the symmetry between the ferro/antiferromagnetic situations and leads to frustrated many-body phases. Since similar interactions also appear in other scenarios (standard waveguides or trapped ions), we now focus on the more different situation where the middle BS at Δ = 0 mediates the interactions, such that the coefficients JmnAB have the form of Eq. 15.

In that case, the Hamiltonian Hspin of Eq. 16 is very unusual: (i) Spins only interact if they are in different sublattices; i.e., the system is bipartite; (ii) the interaction is chiral in the sense that they interact only in case they are properly sorted: the one in lattice A to the left/right of that in lattice B, depending on the sign of δ. Note that δ also controls the interaction length ξ. In particular, for ∣δ ∣ = 1, the interaction only occurs between nearest neighbors, whereas for δ → 0, the interactions become of infinite range. These interactions translate into a rich phase diagram as a function of ξ and μ, which we plot in Fig. 4A for a small chain with Ne = 20 emitters (obtained with exact diagonalization). Let us guide the reader into the different parts:

Fig. 4 Spin models: Phase diagram and correlations.

(A) Ground state average polarization obtained by exact diagonalization for a chain with Ne = 20 emitters with frequency tuned to Δ = 0 as a function of the chemical potential μ and the decay length of the interactions ξ. The different phases discussed in the text, a valence-bond solid (VBS) and a double Néel ordered phase (DN), are shown schematically below, on the left and right, respectively. Interactions of different sign are marked with links of different color. For the VBS, we show two possible configurations corresponding to δ < 0 (top) and δ > 0 (bottom). In the topologically nontrivial phase (δ < 0), two spins are left uncoupled with the rest of the chain. (B) Correlations Cν(r)=σν9σν9+rσν9σν9+r, ν = x, y, z [Cx(r) = Cy(r)] for the same system as in (A) for different interaction lengths, fixing μ = 0 (left column). Correlations for different chemical potentials fixing ξ = 5; darker colors correspond to lower chemical potentials (right column). Note that we have defined a single index r that combines the unit cell position and the sublattice index. The yellow dashed line marks the value of 1/2 expected when the interactions are of infinite range.

(1) The region with maximum average magnetization (in white) corresponds to the places where μ dominates such that all spins are aligned upward.

(2) Now, if we decrease μ from this fully polarized phase in a region where the localization length is short, i.e., ξ ≈ 0.1, we observe a transition into a state with zero average magnetization. This behavior can be understood because in that short-range limit, JmnAB only couples nearest-neighbor AB sites, but not BA sites as shown in the scheme of the lower part of the diagram for δ > 0 (the opposite is true for δ < 0). Thus, the ground state is a product of nearest-neighbor singlets (for J > 0) or triplets (for J < 0). This state is usually referred to as valence-bond solid in the condensed matter literature (37). Note that the difference between δ ≷ 0 is the presence (or absence) of uncoupled spins at the edges.

(3) However, when the bath allows for longer-range interactions (ξ > 1), the transition from the fully polarized phase to the phase of zero magnetization does not occur abruptly but passes through all possible intermediate values of the magnetization. Besides, we also plot in Fig. 4B the spin-spin correlations along the x and z directions (note the symmetry in the xy plane) for the case of μ = 0 to evidence that a qualitatively different order appears as ξ increases. In particular, we show that the spins align along the x direction with a double periodicity, which we can pictorially represent by ∣ ↑ ↑ ↓ ↓ ↑ ↑ …〉x and label as double Néel ordered states. Such orders have been predicted as a consequence of frustration in classical and quantum spin chains with competing nearest-neighbor and next–nearest-neighbor interactions (3840), introduced to describe complex solid-state systems such as multiferroic materials (41). In our case, this order emerges in a system that has long-range interactions but no frustration as the system is always bipartite regardless the interaction length.

To gain analytical intuition of this regime, we take the limit ξ → ∞, where the Hamiltonian (16) reduces toHspin=UHspinUJ(SA+SB+H.c.)(17)where SA/B+=nσegn,A/B, and we have performed a unitary transformation U=noddσzn,Aσzn,B to cancel the alternating signs of JmnAB. Equality in Eq. 17 occurs for a system with periodic boundary conditions, while for finite systems with open boundary conditions, some corrections have to be taken into account due to the fact that not all spins in one sublattice couple to all spins in the other, but only to those to their right/left depending on the sign of δ. The ground state is symmetric under (independent) permutations in A and B. In the thermodynamic limit, we can apply mean field, which predicts symmetry breaking in the spin xy plane. For instance, if J < 0 and the symmetry is broken along the spin direction x, the spins will align so that (SAx)2=(SBx)2=SAxSBx=(Ne/2)2 and SAx2=SBx2=(Ne/2)2.

Since Ne is finite in our case, the symmetry is not broken, but it is still reflected in the correlations, so thatσνm,Aσνn,Aσνm,Aσνn,B1/2(18)with ν = x, y. In the original picture with respect to U, we obtain the double Néel order observed in Fig. 4B. As can be understood, the alternating nature of the interactions is crucial for obtaining this type of ordering. Last, let us mention that the topology of the bath translates into the topology of the spin chain in a straightforward manner: Regardless the range of the effective interactions, the ending spins of the chain will be uncoupled to the rest of the spins if the bath is topologically nontrivial.

This discussion shows the potential of the present setup to act as a quantum simulator of exotic many-body phases not possible to simulate with other known setups. The full characterization of such spin models with topological long-range interactions is interesting on its own and we will present it elsewhere.


Here, we study the situation when QEs are resonant with one of the bands. For concreteness, we only present two results where the unconventional nature of the bath plays a prominent role, namely, the emergence of unexpected super/subradiant states and their consequences when a single photon scatters into one or two QEs.

Dissipative dynamics: Super/subradiance

The band regime is generally characterized by inducing nonunitary dynamics in the QEs. However, when many QEs couple to the bath, there are situations in which the interference between their emission may enhance or diminish (even suppress) the decay of certain states. This phenomenon is known as super/subradiance (19), respectively, and it can be used, e.g., for efficient photon storage (42) or multiphoton generation (43). Let us illustrate this effect with two QEs: In that case, the decay rate of a symmetric/antisymmetric combination of excitations is Γe ± Γ12. When Γ12 = ± Γe, these states decay at a rate that is either twice the individual one or zero. In this latter case, they are called perfect subradiant or dark states.

In standard 1D baths, Γ12(Δ) = Γe(Δ) cos (k(Δ)∣xmn∣); thus, the dark states are such that the wavelength of the photons involved, k(Δ), allows for the formation of a standing wave between the QEs when both try to decay, i.e., when k(Δ) ∣ xmn ∣ = nπ, with n ∈ ℤ. Thus, the emergence of perfect super/subradiant states solely depends on the QE frequency Δ, bath energy dispersion ω(k), and their relative position xmn, which is the common intuition for this phenomenon.

This common wisdom is modified in the bath, where we find situations in which, for the same values of xmn, ω(k), and Δ, the induced dynamics is very different depending on the sign of δ. In particular, when two QEs couple to the A/B sublattice, respectively, the collective decay readsΓ12AB(Δ)=Γesign(Δ)cos(k(Δ)x12ϕ(Δ))(19)which depends both on the photon wavelength mediating the interactionk(Δ)=arccos[Δ22J2(1+δ2)2J2(1δ2)](20)an even function of δ, and on the phase ϕ(Δ) ≡ ϕ(k(Δ)), sensitive to the sign of δ. This ϕ dependence enters through the system-bath coupling when rewriting HI in Eq. 4 in terms of the eigenoperators uk, lk. The intuition behind it is that although the sign of δ does not play a role in the bath properties of an infinite system, when the QEs couple to it, the bath embedded between them is different for δ ≷ 0, making the two situations inequivalent.

Using Eq. 19, we find that to obtain a perfect super/subradiant state, the following conditions must be satisfied: ks)x12 − ϕ(Δs) = nπ, n ∈ ℕ. They come in pairs: If Δs is a superradiant (subradiant) state in the upper band, −Δs is a subradiant (superradiant) state in the lower band. In particular, it can be shown that when δ < 0, the super/subradiant equation has solutions for n = 0, …, x12, while if δ > 0, the equation has solutions for n = 0, …, x12 + 1. Besides, the detunings Δs at which the subradiant states appear also satisfy that J12AB(Δs)0, which guarantees that these subradiant states survive even in the non-Markovian regime [with a correction due to retardation, which is small as long as x12Γe(Δ)/(2 ∣ vg(Δ) ∣ ) ≪ 1]. Apart from inducing different decay dynamics, these different conditions for super/subradiance at fixed Δ also translate in different reflection/transmission coefficients when probing the system through photon scattering, as we show next.

Single-photon scattering

The scattering properties of a single photon impinging into one or several QEs in the ground state can be obtained by solving the secular equation with energies H ∣ Ψk〉 = ± ωk ∣ Ψk〉, with the ± sign depending on the band we are probing (44). Here, we focus on the study of the transmission amplitude t (see scheme of Fig. 5A) for two different situations: (i) a single QE coupled to both cavity A and cavity B in the same unit cell with coupling constants gα and g(1 − α), such that we can interpolate between the case where the QE couples only to sublattice A (α = 1) or B (α = 0), and (ii) a pair of emitters in the AB configuration separated x12 unit cells. After some algebra, we find the exact formulas for the transmission coefficients for the two situationst1QE=2iJ(1δ)sin(k)[J(1+δ)(±ωkΔ)g2α(1α)]2iJ2(1δ2)(±ωkΔ)sin(k)+g2ωk[2α(1α)(eiϕ1)±1](21)t2QE=[2J2(1δ2)(±ωkΔ)sin(k)]2g4ωk2ei2(kx12ϕ)[g2ωk±i2J2(1δ2)(±ωkΔ)sin(k)]2(22)

Fig. 5 Single-photon scattering.

(A) Pictorial representation of the scattering process: An incident photon impinges into a scatterer, part of which is reflected (transmitted) with probability amplitude r (t). Bottom row: Relevant level structure for the single-photon scattering for both scatterers considered: one and two QEs. ∣gg〉 ≡ ∣ g1g2 denotes the common ground state, while S,A=(e1g2±g1e2)/2 denotes the symmetric (antisymmetric) excited state combination of the two QEs. (B) Transmission probability for a single emitter coupled to both A and B cavities of the same unit cell (left) and two emitters in the AB configuration separated a total of x12 = 2 unit cells (right). The parameters in the single emitter case are g = 0.4J, δ = ±0.5, and Δ = 1.5J. The dashed line corresponds to the case where the emitter couples to a single sublattice (α = 0,1) (does not depend on the sign of δ). When the emitter couples to both sublattices (α = 0.3), the perfect reflection resonance experiences a shift that is different for δ > 0 (purple line) or δ < 0 (blue line). The parameters for the two-emitter case are g = 0.1J, δ = ±0.5, and Δ ≃ 1.65J, for which the two QEs are in a subradiant configuration if δ > 0.

In Fig. 5B, we plot the single-photon transmission probability ∣t2 as a function of the frequency of the incident photon for the single-QE (left) and two-QE (right) situations. Let us now explain the different features observed.

Single QE. We first plot in an orange dashed line (Fig. 5B) the results for α = 0,1, showing well-known features for this type of system (44), namely, a perfect transmission dip (∣t2 = 0) when the frequency of the incident photon matches exactly that of the QEs. This is because the Lamb shift induced by the bath in this situation is δωe = 0. The dip has a bandwidth defined by the individual decay rate Γe. Besides, it also shows ∣t2 = 0 at the band edges due to the divergent decay rate at these frequencies, also predicted for standard waveguide setups (44). The situation becomes more interesting for 0 < α < 1, since the QE energy is shifted by δωe = g2α(1 − α)/[J(1 + δ)], which is different for ±δ. This is why the dips in ∣t1QE2 appear at different frequencies for δ = ±0.3. Notice that t1QE is invariant under the transformation α → 1 − α (this is not true for the reflection coefficient, which acquires a δ-dependent phase shift for α = 0 but not for α = 1).

Two QEs. In the right panel of Fig. 5B, we plot ∣t2QE2 for two QEs coupled equally to a bath (same energy, distance, and coupling strength), and where the only difference is the sign of δ of the bath. The distance chosen is small such that retardation effects do not play a significant role. The differences between δ > 0 and δ < 0 in the ∣t2QE2 are even more pronounced than in the single-QE scenario since the responses are now also qualitatively different: While the case δ > 0 features a single transmission dip at the QE frequency, for δ < 0, the transmission dip is followed by a window of frequencies with perfect photon transmission, i.e., ∣t2QE2 = 1. A convenient picture to understand this behavior is depicted in Fig. 5A, where we show that a single photon only probes the symmetric/antisymmetric states in the single excitation subspace (S/A) with the following energies (linewidths) renormalized by the bath ωS,A=Δ±J12ABS, A = Γe ± Γ12). For the parameters chosen (see caption), it can be shown that for δ > 0, the QEs are in a perfect super/subradiant configuration in which one of the states decouples while the other one has a 2Γe decay rate. Thus, at this configuration, the two QEs behave like a single two-level system with an increased linewidth. On the other hand, when δ < 0, both the (anti)symmetric states are coupled to the bath, such that the system is analogous to a V-type system where perfect transmission occurs for an incident frequency ±ωk, EIT = (ωSΓA − ωAΓS)/(ΓA − ΓS) (45) (depicted in a black dashed line; Fig. 5B).

In both the single- and two-QE situations, the different response can be intuitively understood as the QEs couple locally to a different bath for δ ≷ 0. However, this different response of ∣t2 can be thought of as an indirect way of probing topology in these systems.


One of the attractive points of our predictions is that they can be potentially observed in several platforms by combining tools that, in most of the cases, have already been experimentally implemented independently. Some candidate platforms are as follows:

(1) Photonic crystals. The photonic analog of the SSH model has been implemented in several photonic platforms (6, 1012), including some recent photonic crystal realizations (7). The latter are particularly interesting due to the recent advances in their integration with solid-state and natural atomic emitters [see (46, 47) and references therein].

(2) Circuit QED. Superconducting metamaterials mimicking standard waveguide QED are now being routinely built and interfaced with one or many qubits in experiments (48, 49). The only missing piece is the periodic modulation of the couplings to obtain the SSH model, for which there are already proposals using circuit superlattices (50).

(3) Cold atoms. Quantum optical phenomena can be simulated in pure atomic scenarios by using state-dependent optical lattices. The idea is to have two different trapping potentials for two atomic metastable states, such that one state mostly localizes, playing the role of QEs, while the other state propagates as a matter wave. This proposal (51) was recently used (52) to explore the physics of standard waveguide baths. Replacing their potential by an optical superlattice made of two laser fields with different frequencies, one would be able to probe the physics of the topological SSH bath. These cold-atom superlattices have already been implemented in an independent experiment to measure the Zak phase of the SSH model (53).

Beyond these platforms, the bosonic analog of the SSH model has also been discussed in the context of metamaterials (54) or plasmonic and dielectric nanoparticles (55, 56), where the predicted phenomena could also be potentially observed.


In summary, we have presented several phenomena appearing in a topological waveguide QED system with no analog in other optical setups. When the QE frequencies are tuned to the middle bandgap, we predict the appearance of chiral photon BSs that inherit the topological robustness of the bath. Furthermore, we also show how these BSs mediate directional, long-range spin interactions, leading to exotic many-body phases, e.g., double Néel ordered states, which cannot be obtained, to our knowledge, with other bound-state mediated interactions. Besides, we study the scattering and super/subradiant behavior when one or two emitters are resonant with one of the bands, finding that transmission amplitudes can depend on the parameter that controls the topology even though the band energy dispersion is independent of it.

Except for the many-body physics, the rest of the phenomena discussed in this article, that is, the formation of chiral BSs and the peculiar scattering properties, could also be observed in classical setups, since these results are derived within the single-excitation regime. Given the simplicity of the model and the variety of platforms where it can be implemented, we foresee that our predictions can be tested in near-future experiments.

As an outlook, we believe that our work opens complementary research directions on topological photonics, which currently focuses more on the design of exotic light properties (1012, 57, 58). For example, the study of the emergent spin models with long-range topological interactions is interesting on its own and might lead to the discovery of novel many-body phases. Moreover, the scattering-dependent phenomena found in this manuscript can provide alternative paths for probing topology in photonic systems. On the fundamental level, the analytical understanding we develop for 1D systems provides a solid basis to understand quantum optical effects in higher-dimensional topological baths (59, 60).


Supplementary material for this article is available at

Section S1. Integration of the dynamics

Section S2. Subexponential decay

Section S3. Two QE dynamics in the non-Markovian regime

Section S4. Existence conditions of two QE BSs

Section S5. Finite-bath dynamics

Section S6. Middle BSs in 1D baths

Fig. S1. Schematics showing the contour of integration.

Fig. S2. Non-Markovian dynamics.

Fig. S3. Decaying part of the dynamics of a single emitter with parameters Δ = −2J, ∣δ∣ = 0.5, and g = 0.2J.

Fig. S4. Disappearance of the two-QEs BSs in the trivial and topological cases.

Fig. S5. Finite-size effects.

Table S1. Topological properties of several 1D baths, and their corresponding BS features A to C (see text for discussion) when an emitter couples to them.

References (61, 62)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by the Spanish Ministry of Economy and Competitiveness through grant nos. MAT2017-86717-P and BES-2015-071573. A.G.T. and J.I.C. acknowledge the ERC Advanced Grant QENOCOBA under the EU Horizon 2020 program (grant agreement 742102). M.B., G.P., and A.G.-T. acknowledge support from CSIC Research Platform on Quantum Technologies PTI-001. A.G.-T. acknowledges support from the project PGC2018-094792-B-100 (MCIU/AEI/FEDER, EU). Author contributions: M.B. and A.G.-T. conceived the original idea. M.B. did the analytical and numerical analysis under the supervision of A.G.-T. G.P. and J.I.C. provided very useful insights and guidance. All authors discussed and analyzed the results. M.B. and A.G.-T. wrote the manuscript with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article