Subradiant edge states in an atom chain with waveguide-mediated hopping

We analyze the topological and dynamical properties of a system formed by two chains of identical emitters coupled to a waveguide, whose guided modes induce all-to-all excitation hopping. We find that, in the single excitation limit, the bulk topological properties of the Hamiltonian that describes the coherent dynamics of the system are identical to the ones of a one-dimensional Su-Schrieffer-Heeger (SSH) model. However, due to the long-range character of the exchange interactions, we find weakening of the bulk-boundary correspondence. This is illustrated by the variation of the localization length and mass gap of the edge states encountered as we vary the lattice constant and offset between the chains. Most interestingly, we analytically identify parameter regimes where edge states arise which are fully localized to the boundaries of the chain, independently of the system size. These edge states are shown to be not only robust against positional disorder of the atoms in the chain, but also subradiant, i.e., dynamically stable even in the presence of inevitable dissipation processes, establishing the capacity of waveguide QED systems for the realization of symmetry protected topological phases.


Introduction
A topological phase is a form of matter which can be characterized by a non-local order parameter such as the well celebrated Berry phase [1,2,3,4]. The properties of these topological phases are not altered by local perturbations, being inherently robust not only against local noise, but also disorder [5]. These properties have recently propelled systems that are capable of hosting such topological phases to front-runners for applications in the quantum technological context, such as in topological quantum computation [6,7,8], quantum state transfer [9,10,11,12], and quantum error correction [13,14], among others. Together with solid-state systems such as resonators or metamaterials [15,16,17,18,19], cold atom setups such as Rydberg atoms [20,21], or dense atomic lattice gases [22,23] have been shown to display symmetry-protected properties such as flat bands and robust edge states.
Systems governed by quadratic Hamiltonians and Lindbladians can be classified according to their topological properties via the so-called 'tenfold way classification' [24,25,26]. Whether or not topological phases can be hosted in the system is here determined by the dimension of the problem and the symmetries of the Hamiltonian or Lindbladian. Moreover, in these systems the bulk-boundary correspondence determines a relation between the bulk topological invariant and the amount of so-called edge states. These edge states are spatially localized on the boundaries of the system, decaying exponentially into the bulk [27,28]. Several recent works have shown (both theoretically and experimentally) that deviations from the tenfold way classification and weakening of the bulk-boundary correspondence may arise when one considers long-ranged interacting systems, e.g., with hopping constants that go beyond nearest neighbors [29,30,31,32,33,34,35,36,37,38,39,40,41,42]. Many of these works consider interactions that decay as a power law of the distance between sites. This assumption is particularly well-suited for describing cold atom physical systems, e.g. ions or Rydberg atoms.
In this work, we go a step further and analyze the topological properties of a system with infinite-range (or all-to-all) interactions. In par- ticular, we study the topological properties and presence of edge states in a generalized version of the celebrated Su-Schrieffer-Heeger (SSH) model [43] on a one-dimensional chain, where hopping may occur across all sites [see Fig. 1 (a)] [39]. Beyond shedding light on the topological properties of systems with long-range interactions, our motivation to study such a system stems from the latest developments in quantum hybrid systems such as atoms coupled to cavities [44,45] or waveguides, such as nanofibers or photonic crystals [46,19,47]. In these systems, the translationally symmetric light modes mediate all-to-all exchange interactions among the atoms [48,49,50,51]. In particular, the extended SSH model we study here describes approximately the physics of an array of two-level atoms coupled to a nanophotonic waveguide [see Fig. 1 We find that when the system conserves chiral symmetry, the bulk Hamiltonian belongs to the same topological class (BDI) as the original SSH model with short-range hopping. However, in general here the bulk-boundary correspondence is weakened, leading to the appearance of massive edge states with non-zero energy. We uncover parameter regimes where approximate zero-energy edge states, which are robust against local disorder, arise. Interestingly, we also reveal the presence of robust edge states in the case where the chiral symmetry is broken, and hence the tenfold way classification does not predict their existence.
Finally, we study the out-of-equilibrium dynamics of the system, which reveals that some of the uncovered edge states are subradiant, i.e., decay collectively emitting a photon into the waveguide in a much longer timescale than an isolated single atom would do.

Physical system and model Hamiltonian
Let us consider a chain of M identical atoms arranged in a bipartite lattice structure around a waveguide, e.g., a cylindrical optical nanofiber of refractive index n 1 . The two atomic chains, with lattice constant a and formed by M/2 atoms each, are placed close to the waveguide parallel to its longitudinal axis (here the z-axis), offset with respect to one another by a distance b [see Fig.  1(b)]. The atoms are modelled as two-level systems with ground and excited states |g i and |e i for the i-th atom, respectively.
Due to the coupling to the environment, the atoms incoherently emit photons and also interact with one another via dipole-dipole interactions due to the exchange of virtual photons mediated both by the guided and unguided modes of the waveguide [48,49,50]. Under the usual Born-Markov and secular approximations, one can describe the dynamics of this system through the quantum optical master equatioṅ Here, the Hamiltonian describes the coherent dipole-dipole interactions between the atoms in the system, the strength of which is determined by V ij . σ † i and σ i are the raising and lowering operators for the atom i, respectively, such that σ † i |g i = |e i and σ i |e i = |g i . The second term in (1) represents the dissipation via the emission of photons into the environment. The dissipation in this system has in general a collective character, determined by the form and size (relative to the single atom decay rate) of the non-diagonal elements of the matrix Γ ij [52]. We will elaborate further on this collective character, which leads to super-and subradiant emission, in the last section of the paper.
It is convenient to use cylindrical coordinates to represent the position (r i , ϕ i , z i ) of the i-th atom, assuming that the waveguide is concentric with the z axis [see Fig. 1(b)]. Both V ij and Γ ij can be split into the corresponding contributions from the guided modes and unguided (or radiation) modes as V ij = V g ij +V r ij and Γ ij = Γ g ij +Γ r ij . The contribution to the interactions due to the fundamental guided modes of the waveguide are given by Here, z ij = z i − z j , f = ±1 represents the propagation direction of the guided mode in the ±z direction and l = ±1 its polarisation in the clockwise or anticlockwise direction, respectively. The Green's function G ω 0 f li couples atom i to the guided mode of frequency ω 0 , with G * ω 0 f li denoting its complex conjugate. In the case of a nanofiber, this function reads where 0 is the vacuum permittivity, is the transition dipole moment of atom i, e (ωf l) (r) = e r , −le ϕ , f e z T is the electric field profile function of the guided mode evaluated at the distance r from the atom chain to the z axis (the center of the nanofiber), and β and β are the longitudinal propagation constant within the waveguide and its derivative, respectively. The specific choice of orientation and position of the atomic dipole moments with respect to the nanofiber and with respect to each other gives rise to different coherent and incoherent interactions. Moreover, the atomic dipole moment orientation determines whether the emission into the nanofiber is directional or not, i.e., whether the emission of photons from a single atom occurs in a preferential direction (i.e., to the rightor the left-propagating modes of the nanofiber). For simplicity, we will consider here the situation where all dipole moments are perpendicular to the nanofiber, i.e. d i r = d ≡ |d i |, which leads to a symmetric emission to the right and to the left. Moreover, all atoms are placed on a line, such that ϕ ij = ϕ i − ϕ j = 0 for all i, j = 1, . . . M .
The coefficients representing the interactions read here where is the single atom decay rate into the waveguide. Note, that the interactions induced by these guided modes do not decay with the distance. In this paper, we will be only considering the interactions mediated via these guided modes. This is a particularly good approximation when the distances a and b are comparable or larger than the wavelength λ 0 = 2πc/ω 0 of the |g → |e transition as, contrary to the guided case, both V r ij and Γ r ij decay as λ 0 /z 3 , with z being the distance between two atoms [53,54]. We start by considering the Hamiltonian (2) and consider the effects of dissipation in a later section. Hamiltonian (2) conserves the number of excitations in the system. Constraining our analysis to the single-excitation sector, it can be mapped via a Jordan-Wigner transformation into a quadratic Hamiltonian for spinless fermions where J |i−j| = V g ij and c j and c † j are the annihilation and creation operators of a spinless fermion on site j, respectively. Due to the bipartite form of the lattice, we can describe this system as a one-dimensional array formed by M sites distributed over M/2 unit cells, each hosting two sites belonging to two distinct sublattices, A and B [see Fig. 1(a)]. Hopping rates differ when the two sites involved belong to the same or different sublattices (i.e. |i − j| even and odd, respectively). Moreover, in the latter case we distinguish the case when i is odd or even, i.e., when the hopping occurs between sites separated by p − 1 or p unit cells, respectively. In summary, we denote the hopping rates [see Fig. 1

(a)] as
where we have introduced a q = c 2i−1 and b q = c 2i as the annihilation operators of a fermion in the q-th unit cell in sublattices A and B, respectively. Note, that when all hopping rates are set to zero except J 1 and J 1 , this Hamiltonian describes the well-known SSH model [43]. In our case, considering the form of (4), we obtain Hence, in this system the values of all the hopping rates are real and restricted to the interval

Bulk properties
We start by analyzing the bulk properties of Hamiltonian (7) assuming periodic boundary conditions. Here, the bulk Hamiltonian becomes a real symmetric block circulant matrix (2 × 2 blocks), allowing us to diagonalize the bulk Hamiltonian in the basis of Bloch states for the external degrees of freedom, i.e. cell index, as with k ∈ 0, 2π a . Here, we have introduced the eigenvectors It is convenient to write the bulk Hamiltonian where we have separated the parts containing the intra-sublattice (even, diagonal) hopping coefficients with and the inter-sublattice (odd, off-diagonal) hoppings, The two energy bands are then determined by The relevant symmetries to determine the topology of this system according to the tenfold way classification [24] are time-reversal, particlehole, and chiral symmetry. These are determined in our case by the operators T = K, C = σ z K and S = T · C = σ z , respectively, where K is the conjugation operator and σ z is the usual spin-1/2 Pauli matrix. In general, for real hopping However, h e (k) only possesses time-reversal symmetry and, most importantly, its presence makes the full Hamiltonian break chiral symmetry. This system is therefore to be attributed to the AI symmetry class, for which there is no topological invariant in one dimension. On the other hand, when n 0 (k) = 0 (e.g., when the all even hoppings J 2p = 0 ∀p) the system is attributed to the BDI symmetry class, which has a Z-type topological invariant in one dimension. This is, for example, the case for the usual SSH model in one dimension. The chiral symmetry is also reflected on the eigenenergies of the system, which lead to a spectrum that is symmetric about E = 0 when n 0 (k) = 0. In real space, the consequence of chiral symmetry is that if a pair of degenerate eigenstates have zero energy, each eigenstate can be rewritten such that they have support only on sublattice A and B, respectively. Conversely, if an eigenstate has support only on one sublattice, its energy must be zero [39]. In the following, unless stated otherwise, we will work on the assumption that chiral symmetry is indeed preserved.
The intrinsic geometric phase of the eigenstates of h(k) is given by the variation of the argument  (23) and (22) for the phase boundaries.
of n(k) = n x (k) + in y (k) through the First Brillouin Zone (FBZ). The winding number of n(k) can be directly mapped to the Berry phase and thus be used to classify the intrinsic geometry, and hence topology, of the system. The winding number of the system can be found as where ∆φ(k) is the variation of the argument of n(k) through the FBZ. Due to the presence of all-to-all hopping rates, the function n(k) depends on the number of atoms M even when considering periodic boundary conditions. However, we find when evaluating the winding number numerically that it only takes the values ν = 0 and 1 (see Fig. 2) and that it does not vary with the system size. The values of b and a for a fixed β at which the system changes from a topologically trivial to nontrivial phase and vice versa can be, moreover, obtained analytically by analyzing the symmetries of n(k): as a consequence of time-reversal symmetry, n x (k) and n y (k) are symmetric and antisymmetric, respectively, with respect to k = π/a. Hence, the closed path of n(k) in complex space within the FBZ is symmetric about the real axis. In order to determine the winding number, we need to verify whether the closed path of n(k) encloses the origin (ν = 1) or not (ν = 0). Since at the time reversal invariant momenta k = 0 and π/a, n(k) crosses the real axis, i.e., n y (k) = 0, it is enough to determine whether there is a change of sign in n x at the two momenta. In particular, if sgn( nx(π/a) nx(0) ) = 1, then ν = 0, and if sgn( nx(π/a) nx(0) ) = −1, then ν = 1, such that We can determine the phase boundaries analytically using that Hence, one easily extracts that the winding number changes at the values where λ = 2π β is the wavelength of the light propagating in the waveguide. Note, that due to the presence of all-to-all interactions, the gap between the two energy bands also closes whenever a/λ = n/M for any n ∈ Z. However, only when condition (22) is satisfied is this gap closing accompanied by a change in the topological invariant of the system. Finally, fixing a, we also find a change from ν = 0 to ν = 1 and viceversa when Note, that when only J 1 and J 1 are different from zero, the Hamiltonian (7) is precisely the one describing the SSH model, for which ν = 1 when |J 1 | > J 1 while ν = 0 when |J 1 | < J 1 [55]. Interestingly, we find that the winding number for our fully connected model (Fig. 2) coincides with the one of the SSH model, i.e., the topological invariant of this system does not change as the connectivity of the problem is increased from nearest neighbors to all-to-all interactions. valence bands. These states are expected to be exponentially localised on the edges of the chain [56]. For a semi-infinite chain, these eigenstates have exactly zero energy, hence the name 'zeroenergy edge states'. The bulk-boundary correspondence gets weakened in the presence of longranged hopping due to the correlations that are built between the bulk and the edges of the system, giving rise to massive edge states, i.e., edge states with energy ε different from zero (so-called mass gap), which may decay from the boundary into the bulk algebraically instead of exponentially [33].
We illustrate here the weakening of the bulkboundary correspondence for our system with allto-all interactions, when chiral symmetry is con-served (all even couplings J 2p = 0). The Hamiltonian is given by In Figure 3(a) we show the mass gap ε for a system of M = 10 sites as a function of a/λ and b/a. In order to calculate the mass gap, we obtain the eigenenergies ε n with n = 1, . . . , M of the Hamiltonian H for each value of a/λ and b/λ. Due to the imposed chiral symmetry, the resulting spectrum is symmetric around E = 0 [see, e.g., Fig.  3(b)]. The mass gap is thus calculated as ε = ε M/2+1 − ε M/2 /2. Note, that the mass gap becomes simply half of the gap between the two bulk bands in the topologically trivial area, i.e. where ν = 0 (c.f. Fig. 2). Since in this system the interactions are allto-all, the specific values of the couplings (determined in turn by a/λ and b/λ) give rise to a kaleidoscope of regimes with different effective geometries and edge state properties. In the following, we will analyze in detail some of these parameter regimes.

Flat states (|J
In this case, the first particle can hop to its nearest neighbor and to the last site of the chain with the same hopping rate. Moreover, due to the specific form of the rates (10-11), imposing that |J M −1 | = |J 1 | automatically implies that Here, not only the bulk but also the finite open system are translationally symmetric. Since all cells are equivalent, the concept of boundary becomes ill-defined. As a consequence, as one can observe in Fig. 3(d), all eigenstates of the Hamiltonian are such that all sites are uniformly occupied, reflecting the translational symmetry. Note, that the condition |J M −1 | = |J 1 | is satisfied when a/λ = m/M with m ∈ N. As we pointed out in the previous section, here the bulk gap in momentum space closes [see Fig. 3(e)]. Hence, even though this gap closing does not come accompanied by a change in the topological invariant, specifically in this parameter regime [horizontal dashed lines in Fig. 3(a)] the winding number is not well defined and hence the flat states are not topological edge states.

Localized states (J M −1 = 0)
The translational symmetry is broken whenever |J M −1 | = |J 1 |, as one can again distinguish the boundaries from the bulk. The extreme case here occurs when J M −1 = 0, i.e., when the first and the last atom in the chain are not coupled to each other [see Fig. 3(f)]. This occurs whenever the parameters satisfy In the topologically non-trivial region and along this line in parameter space [solid white lines in Fig. 3(a)], one finds two approximately zero energy [mass gap ε ≈ 0] eigenstates. While usual topologically protected edge states are exponentially localized on one of the two boundaries, we find here that the population is concentrated on both boundaries of the chain [cf. Fig. 3(g)], such that the occupation probability of the sites in the p-th cell in sublattice A and B, respectively, can be written as Here, we have introduced ξ as the localization length of the edge state for each sublattice. As we can observe in Fig. 3(h), the localization length is smaller (i.e. edge states more localized in the edges of the lattice) the smaller the ratio J 1 /J 1 . One can obtain that the localization length is approximately given by As one can see in Fig. 3(h), this simple expression yields indeed an excellent approximation to the localization length, particularly for small values of the ratio |J 1 /J 1 |. Note also that the edge states have support of both sublattices A and B, instead of only on one as it is the case in the SSH model. One recovers this feature only when the energy of the localized state is exactly zero. The edge states become completely localized on the edges when J 1 = 0. It can be easily shown that this occurs when with m and m being natural numbers. Here, the edge states can be analytically obtained, and they read (27) with n = M/4 , and |0 being the state with no fermions. Note, moreover, that at this point these edge states have exact zero energy, and that the spectrum becomes doubly degenerate [see Fig. 3(b)]. This structure is reminiscent of what occurs in the presence of a socalled strong zero mode [57,58,59,60,61,62,63]. Indeed, here there exists an operator Ψ = p i(−1) p a † p b p − b † p a p that anticommutes with a symmetry of the Hamiltonian D = Due to the long-range character of the hopping rates, though, in our case the strong zero mode Ψ is not concentrated on the edges of the lattice, as it is usually the case considering systems with short range interactions. Moreover, the Hamiltonian exactly commutes with Ψ for all system sizes (cf. [63]), instead of doing so only in the thermodynamic limit.

Robustness against disorder
Topologically protected edge states are robust against local disorder. In order to test whether the edge states described above are indeed robust, we introduce disorder as an uncertainty in the position of the atoms in the chain. We do this by letting the positions of the atoms z i in the chain have a Gaussian distribution around their nondisordered positions z where R is a number randomly generated from a normal distribution centred around zero and σ is the standard deviation of the normal distribution that is used to control the level of disorder applied. For each realization of disorder, we calculate the mass gap ε σ of the edge state |Ψ σ , and the fidelity between the measured edge state and the state with zero disorder. After many realizations, we extract the corresponding average fluctuations of the hoppings, J → J + δJ, and the average mass gap ε and fidelity F , shown in Fig. 4 (a) and (b), respectively, for the fully localized case, where both J 1 and J M −1 are equal to zero [where in the absence of disorder ε = 0 and the edge states are exactly given by (27)]. As one can observe, as we increase the average fluctuations up to more than a 10% (note that the maximum value δJ could take is γ/2) the mass gap increases linearly with the disorder, remaining much smaller than the gap between the bulk bands in all cases. Moreover, the fidelity of the edge state with respect to the no-disorder case remains above 0.99. In  5 Edge states with J 2p = 0 Up until now we have studied a system where the chiral symmetry is not broken, i.e., where the intra-sublattice hoppings J 2p have been put to zero. In a real atom-waveguide system, these cannot be generally put to zero unless a/λ = m/2 [see Eq. (9)]. Unfortunately, this regime is also precisely where the bulk gap closes (horizontal lines in Fig. 2), and hence no topological invariant can be defined. In this section, we will explore a situation where, even though the chiral symmetry is not strictly conserved, states with energy close to zero appear in the gap between the bulk bands which, moreover, are localized to the edges of the chain and are robust to local disorder. To find these states, we consider the number of atoms to be even but not a multiple of four, i.e. M = 4n + 2 for n ∈ N such that, at the point determined by condition (26), one can make exactly half of the even couplings zero, i.e.
Moreover, as it can be observed in Fig. 5(a), the spectrum becomes symmetric around E = 0.
Here, the eigenstates of the Hamiltonian can be obtained analytically. In particular, we find that there is a manifold of M/2 + 1 eigenstates with zero energy. Some of these eigenstates are indeed localized in the edge, such as These states are, moreover, robust in the presence of local disorder, in the sense that they stay localized to the edges, not spreading into the bulk [ Fig.  5(b)]. However, note that the mixing with the rest of eigenstates of the quasi-degenerate manifold lead to small variations in the populations of the edge sites involved.
To avoid this mixing, one may move away from this particular symmetric point. As we can observe in Fig. 5(a), both to the right and left of this point (larger and smaller b/a, respectively), the zero-energy degeneracy is indeed lifted. However, the eigenstates with energies now slightly different from zero are still localized on the edges and extremely robust to the addition of local disorder [see, e.g. Fig. 5(c)]. Moreover, again the structure of the spectrum and the eigenstates, including the existence of robust edge states, are independent of the number of atoms in the system, M . Note, however, that strictly speaking these localized states are not topologically protected states, as the chiral symmetry is here broken and hence there is no topological invariant associated with the corresponding Hamiltonian.

Non-equilibrium dynamics
Up until now, we have studied only the static properties of the atom-waveguide system in the absence of dissipation. While this has offered us the possibility to study the topological properties of the system and demonstrate the existence of edge states, one can see from the expressions of the coherent and dissipative coupling V g ij and Γ g ij , in Eqs. (4) and (5), respectively, that the two mechanisms are inevitably intertwined, i.e., the dissipation in general cannot be neglected. In this Section, we study the out-of-equilibrium dynamics of the system in the presence of dissipation, focusing in particular on the fate of the above found edge states.
Since the dissipation coefficients Γ g ij = γ cos (βz ij ) do not decay with the distance between the atoms, the dissipation in this system has a collective character. We establish this by calculating the dissipation channels and corresponding decay rates Γ n in the single-particle sector as the eigenstates and eigenvalues of the matrix Γ g ij , respectively. Independently of the parameter regime, there are only two superradiant modes with decay rates Γ 1,2 much larger than γ. The rest of the states forms a subradiant mani- fold with zero decay rate [64], as one can observe in Fig. 6(a) for a/λ = 1.25. Along this line in parameter space, one can see that one of these two superradiant eigenstates has no support on the lattice edges, while the other one does. However, the latter has an alternating phase profile, similarly to the state |A given in Eq. (32). Any initial state that one chooses to excite in this system that does not have a large overlap with these two superradiant states will be, thus, subradiant.
If, moreover, we are in a parameter regime where an edge state has been predicted to exist in the previous Section, we expect an initial state with large support on the edge to remain localized in the edge during the dynamics determined by (1). This is indeed what we observe in the upper panels of Figs. 6(b) and (c) by measuring the survival probability, P (t) = M/2 p=1 |A p (t)| 2 + |B p (t)| 2 for two values of the parameter b/a: the initial symmetric superposition localized state, (i.e., A 1 = B 1 = 1/ √ 2, A p = B p = 0 for all p = 2, . . . M/2) remains localized, and the excitation has an extremely long lifetime. In particular, for b/a = 0.4 the initial state is an eigenstate of both the coherent interactions (i.e. a zero-energy edge state), and of the dissipation (belonging to the subradiant manifold). Hence, as expected, the state remains completely unchanged throughout the dynamics. In the lower panels of Figs. 6(b) and (c) we show how this dynamics is also extremely robust against local disorder.

Conclusions and Outlook
We have studied an atom chain coupled to a waveguide and uncovered the existence of edge states, i.e. states localized in the boundary of the chain which are not only robust against the presence of local disorder, but also subradiant, i.e. they possess extremely long lifetimes.
The analysis presented here has been done in the absence of unguided or radiation modes, which has allowed us to obtain analytical results complementing the numerical ones. The contribution of the radiation modes to the dipole-dipole interactions V r ij is well justified when the atoms are far from each other (compared to the wavelength λ), and are expected to represent small perturbation that does not affect the robust edge states presented here. The unguided off-diagonal components of the decay matrix, Γ r ij , can be neglected on the same grounds. However, the diagonal elements of this matrix, representing the single atom decay into the radiation continuum, can be extremely destructive for the subradiant character of the edge states. This can be remedied under the right experimental implementations such as photonic crystal waveguides [65] or superconducting qubits embedded in 1D transmission lines [66], where the decay into the guided modes can be considered predominant.
Future work will entail a closer look into the emergence of strong zero modes in our system, and whether these can also be identified in the full model including dissipation. The fate of the topological properties and edge states depicted here when the emission is not symmetric (e.g. for different atomic dipole polarizations) will also be in the focus of our future research. Finally, a natural extension of the present system is the study of the topological properties and strong zero modes in a gas of atoms coupled to a 2D waveguide structure such as a 2D photonic crystal [67].