Unconventional quantum optics in topological waveguide QED

Methods

In the next sections, we study the dynamics emerging from the global QE-bath Hamiltonian H = H_S + H_B + H_I 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]+i∑n,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 x_mn = x_n − x_m. For our model, they can be calculated analytically in the thermodynamic limit (N → ∞) yieldingΣmnAA/BB(z)=−g2z[y+∣xmn∣Θ+(y+)−y−∣xmn∣Θ−(y+)]z4−4J2(1+δ2)z2+16J4δ2(6)ΣmnAB(z)=g2J[Fxmn(y+)Θ+(y+)−Fxmn(y−)Θ−(y+)]z4−4J2(1+δ2)z2+16J4δ2(7)where F_n(z) = (1 + δ)z^∣n∣ + (1 − δ)z^{∣n + 1∣}, Θ_±(z) = Θ(±1 ∓ ∣ z ∣), Θ(z) is Heaviside’s step function, andy±=z2−2J2(1+δ2)±z4−4J2(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.

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†]∣g〉∣vac〉(9)

In both perturbative and exact treatments, the dynamics of C_e(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)z4−4J2(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 C_e(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 (26–28). 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.

• Download high-res image
• Open in new tab
• Download Powerpoint

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〉 = E_BS ∣ Ψ_BS〉, with E_BS 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: E_BS = Δ + Σ_e(E_BS). 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 (26–30) 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 C_e 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 ∣C_e(t → ∞ )∣² = ∣C_e∣⁴ = [1 + g²/(4J² ∣ δ ∣)]⁻².

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 C_j,α have the same sign (see the left column of Fig. 3A, top and bottom row). In the lower (upper) bandgap, C_j,α 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 C_j,α (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 E_BS = 0 is given by C_j,a = 0 andCj,b={gCe(−1)jJ(1+δ)(1−δ1+δ)j,j≥00,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 λBS∼1/∣Δedge∣, with Δ_edge being the smallest detuning between the QE and the band-edge frequencies.

• Download high-res image
• Open in new tab
• Download Powerpoint

The physical intuition of the appearance of such chiral BS at E_BS = 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 HB→HB+∑j(εa,jaj†aj+εb,jbj†bj) and breaks the chiral symmetry of the original model, while the latter corresponds to the addition of off-diagonal random terms HB→HB+∑j(ε1,jbj†aj+ε2,jaj+1†bj+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 E_BS = 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 λBS−1 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 x_mn 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 x₁ couples to sublattice A, and the second one at x₂ 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 N_e 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/r³). 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 H_spin 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 N_e = 20 emitters (obtained with exact diagonalization). Let us guide the reader into the different parts:

• Download high-res image
• Open in new tab
• Download Powerpoint

(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 (38–40), 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′=UHspinU†≃J(SA+SB−+H.c.)(17)where SA/B+=∑nσegn,A/B, and we have performed a unitary transformation U=∏n∈ℤoddσ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 〈SAx〉2=〈SBx〉2=(Ne/2)2.

Since N_e 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,B〉≃1/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.

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 ± Γ₁₂. When Γ₁₂ = ± Γ_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, Γ₁₂(Δ) = Γ_e(Δ) cos (k(Δ)∣x_mn∣); 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(Δ) ∣ x_mn ∣ = 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 x_mn, 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 x_mn, ω(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[Δ2−2J2(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 H_I in Eq. 4 in terms of the eigenoperators u_k, l_k. 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: k(Δ_s)x₁₂ − ϕ(Δ_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, …, x₁₂, while if δ > 0, the equation has solutions for n = 0, …, x₁₂ + 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 x₁₂Γ_e(Δ)/(2 ∣ v_g(Δ) ∣ ) ≪ 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 x₁₂ 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−α)(e−iϕ∓1)±1](21)t2QE=[2J2(1−δ2)(±ωk−Δ)sin(k)]2g4ωk2ei2(kx12−ϕ)−[g2ωk±i2J2(1−δ2)(±ωk−Δ)sin(k)]2(22)

• Download high-res image
• Open in new tab
• Download Powerpoint

In Fig. 5B, we plot the single-photon transmission probability ∣t∣² 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.