Quantum Hall phases emerging from atom–photon interactions

We reveal the emergence of quantum Hall phases, topological edge states, spectral Landau levels, and Hofstadter butterfly spectra in the two-particle Hilbert space of an array of periodically spaced two-level atoms coupled to a waveguide (waveguide quantum electrodynamics). While the topological edge states of photons require fine-tuned spatial or temporal modulations of the parameters to generate synthetic magnetic fields and the quantum Hall effect, here we demonstrate that a synthetic magnetic field can be self-induced solely by atom–photon interactions. The fact that topological order can be self-induced in what is arguably the simplest possible quantum structure shows the richness of these waveguide quantum electrodynamics systems. We believe that our findings will advance several research disciplines including quantum optics, many-body physics, and nonlinear topological photonics, and that it will set an important reference point for the future experiments on qubit arrays and quantum simulators.


INTRODUCTION
Recent technological advances have underpinned the rapid development of cavity quantum electrodynamics (QED) and circuit QED, which allow to exploit quantum properties of light for applications in information processing [1][2][3] . A closely related subfield with growing theoretical and experimental interest is waveguide QED 4,5 , which studies one-dimensional arrays of natural or artificial atoms coupled to a waveguide. Existing platforms for waveguide QED systems include cold atoms 6 , defect centers 7 , superconducting qubits 2,3,[8][9][10][11][12] , and emerging structures based on exciton-polaritons 13 .
Waveguide QED is promising for many applications in quantum information processing. It can allow us to efficiently generate [14][15][16] , detect 17 , slow 18 , and store quantum light 19 . It is also useful as a platform for quantum simulators of complex many-mode physics 20,21 . A crucial advantage of waveguide QED systems is that they can exhibit long-range coupling between distant atoms mediated by light. This makes the dispersion of atomic excitations and their interactions markedly different from those in the typical case of nearest-neighbor coupling in circuit QED 1 or in conventional condensed matter systems. Thus, previously uncharted physical regimes can be naturally accessed in these wavequidebased quantum simulators.
The goal of this work is to explore the potential of waveguide QED to simulate quantum many-body topological phases of interacting matter. Topological phases underpin many concepts of modern matter physics. Topological edge states of electrons are usually created using magnetic fields 22 or spin-orbit interactions [23][24][25] . Studies of the quantum Hall effect and topological insulators have also inspired rapid progress in topological photonics where the central goal is to create robust edge states of light immune to disorder [26][27][28][29] . Since the effects of magnetic fields on light are weak, the realization of topological concepts in photonics requires artificial structures or metamaterials 30 .
Alternative approaches rely on time modulation of structure parameters [31][32][33][34] or engineered nonlinearities 35,36 . These techniques create effective gauge fields in real or synthetic dimensions, and mimic the effects of magnetic fields or spin-orbit couplings for photons.
Here, we uncover that the hallmarks of quantum Hall phases, Landau energy levels, topological edge states, and the Hofstadter butterfly spectrum [37][38][39] naturally arise in a finite one-dimensional array of closely spaced two-level atoms (qubits) coupled to photons in a waveguide, shown in Fig. 1a. The striking feature of our prediction is that the quantum Hall phase can emerge solely from interactions of indistinguishable particles without any external magnetic fields or special fine tuning or modulation of spatial or temporal parameters.
While the proposed idea requires highly coherent structures, it could potentially be tested in already available arrays of superconducting qubits coupled to waveguides 40 using the existing approaches to probe spatial profile of two-photon quantum states 41 . Our results open the way to engineer complex multiphoton states and realize their topological protection against disorder.

Interaction-induced topological states
Here, we show how to realize topological quantum Hall phases from atom-photon interactions. In waveguide QED setups (shown schematically in Fig. 1a), photons become strongly coupled to atoms and create light-matter quasi-particles called polaritons. These polaritons are not independent but are strongly interacting, because one atom cannot absorb two photons simultaneously 42 , leading to on-site repulsion. While the considered model is paradigmatic for quantum optics 4,5,43 , its two-particle Hilbert space was not analyzed until recently. As shown in Fig. 1c the polariton wave vector is comparable with that of light, a collective atomic state is easily excited optically, and it generally gets "darker" for larger wave vectors. In a two-particle "bright" state, the wave vectors of both excitations are small, which corresponds to Dicke superradiance 44,45 . Two-particle dark states where both wave vectors are large were predicted only last year. These states originate from the fermionization of strongly interacting polaritons 46 . It has also been suggested that interactions in the corner regions 47 of the diagram of Fig. 1d can localize one of the two polaritons in the center of the array 48 .
In this paper, we predict topological edge states driven by polariton-polariton interactions in the regions indicated by butterflies in Fig. 1d. In this region for a finite array, one polariton forms a standing wave with multiple nodes and a periodic potential for the other indistinguishable polariton, see Fig. 1b, e. As a result, the interaction is described by the self-induced Aubry-André-Harper 49 model that is mathematically equivalent to the quantum Hall problem on a lattice 50,51 .
The periodic modulation is an intrinsic feature that arises naturally due to the polariton-polariton interactions, in sharp contrast to previous studies 52, 53 where the modulation is imposed deliberately either by engineering the lattice 53,54 or applying external fields 29,33 . We show that the full Hofstadter butterfly-like spectrum could be obtained in a single shot from just a fixed atomic array, eliminating the need to continuously tune an external magnetic field in a conventional setup 33,39 .

Two-polariton states
The considered periodic one-dimensional array of N two-level atoms (qubits) coupled to light is described by an effective Dicketype Hamiltonian 46,47,55 H ¼ where c is the speed of light, σ y n is the operator creating excitation of the atom n with the resonance frequency ω 0 , ðσ y n Þ 2 ¼ 0 and Γ 0 is the radiative decay rate of a single atom. While for d = 0 the Hamiltonian (1) is equivalent to the conventional Dicke model 44 , even small interatomic spacings 0 < d ≪ 2πc/ω 0 make the model considerably richer. Single-particle eigenstates of Eq. (1) in the infinite array (N → ∞) are polaritons with the energy dispersion εðkÞ ¼ Γ 0 sin φ=ðcos k À cos φÞ 55,56 , which is shown schematically in Fig. 1e. The eigenstates are Bloch waves with real energies, characterized by continuous distribution of the wave vectors − π/d < k < π/d. The dispersion consists of two polaritonic branches, resulting from the avoided crossing of light with the atomic resonance where both branches are described by the same expression ε(k). The upper branch corresponds to |k| < ω 0 d/c and in this work we focus on the lower branch with large wave vectors, |k| > ω 0 d/c. In the finite array of N atoms, the wave vectors are quantized, k j d = πj/N, j = 1, 2, …N 57 , and the single-particle eigenstates are standing waves, as shown in Fig. 2a-d. Owing to the possibility of radiative decay into the waveguide, the eigenstates of the finite array (shown in Fig. 2a) have non-zero imaginary energies, Imε < 0 46 . Our main results below are obtained for strongly subradiant states where jImεj ( Γ 0 . They are not qualitatively affected by the radiative losses and remain valid even when the non-Hermitian part of the Hamiltonian Eq. (1) is neglected (see Supplementary Fig. 5). Crucially, the spectrum in Fig. 2a condenses near the resonance ε = ω 0 , where the group velocity of polaritons decreases.
Next, we proceed to the double-excited states Ψ ¼ P n;m ψ nm σ y n σ y m 0 j i. Their spectrum, obtained from the Schrödinger equation HΨ = 2εΨ, is shown in Fig. 2e-h in different energy scales. The spectrum consists of distinct groups of states with close energies. The distributions of the two-particle eigenvalue energies are plotted in the complex plane in Fig. 2e-g, and their shape resembles the sum of two single-particle eigenvalues (the single-particle spectra is plotted in Fig. 2a). This means that most states in Fig. 2e-g can be described by ε ≈ (ε j + ε i )/2, where ε j and ε i are the single-particle energies from Fig. 2a. However, the dense part of the group, which corresponds to ε i → ω 0 (see red arrows in Fig. 2f), is drastically transformed by the interaction. Three characteristic states from the group with j = 7 are presented in Fig. 2i-k. While the state in Fig. 2i is just a symmetrized product of two standing waves, weakly modified by interaction, the role of the interaction dramatically increases for Reε À ω 0 > À 0:66Γ 0 in Fig. 2h. The spectrum is split by interaction into relatively delocalized states with smaller radiative decay rate (yellow ellipse in Fig. 2h, k) and states with larger radiative losses where one of the two polaritons is localized at the edge of the structure (blue ellipse in Fig. 2h, j).
This interaction-induced transformation of the two-polariton spectrum is our central result. The delocalized states are almost (j − 1)-fold degenerate, where j is the group number and correspond to the Landau levels in the effective magnetic field. The states in Fig. 2j come in degenerate pairs corresponding to topological edge states localized at the opposite sides of the array.

Landau levels, topological edge states, and Hofstadter butterfly
We now present an analytical model explaining the topological origin behind the interaction-induced edge states in Fig. 2j. In the basis 58 , the following ansatz can be used for the two-polariton state where ψ ðjÞ x and χ x are the wavefunctions of the first and second polaritons. The former corresponds well to the standing wave ψ ðjÞ x ¼ cos k j ðx À 1 2 Þ. To determine the latter, we derive the Schrödinger equation that accounts for interaction between the  48 are also shown. e Single-particle polariton dispersion. Interaction of a lower-branch polariton with small k and that with large k is illustrated.
polaritons. It is given by (see "Methods" for more details) where ω j % ω 0 À 2φΓ 0 =k 2 j is the real part of the eigenfrequency of the single-polaritonic state ψ (j) and φ = ω 0 d/c. Equation (3) describes a motion of a particle on a lattice in an external potential of a standing wave with the period N/j. It has a striking similarity to the Harper equation for an electron moving in a square lattice subjected to the perpendicular magnetic field 59 : Here, α is the magnetic flux through the unit cell and k y is the wave vector in the perpendicular direction. For small magnetic fields α ≪ 1, the discreteness of the problem plays no role and the energy spectrum of Eq. (4) is a ladder of degenerate Landau levels for electrons moving along quantized cyclotron orbits. In the finite structure, the edge states of topological nature arise in the gaps between the Landau levels. Such states correspond to electrons moving along skipping orbits at the structure edge, and are the origin for the quantum Hall effect 22 .
In our system, the ratio j/N of the group number to the total number of atoms in Eq. The energy spectrum of the Harper Eq. (4) becomes very rich when the magnetic flux α increases. The Landau levels split and transform into a celebrated Hofstadter butterfly 38 , shown also in Supplementary Fig. 2. The butterfly has a self-similar structure with q allowed energy bands at the rational fluxes α = p/q 37 and a Cantor-set spectrum for irrational fluxes. Even though in our case the effective magnetic flux j/N is rational, we can still extract an analog of the Hofstadter butterfly from the two-polariton spectrum in Fig. 2e-g. We separate the groups of states in Fig. 2e-h formed by different standing waves (i.e., different effective magnetic fields) and align them horizontally, the details are presented in "Methods" and Supplementary Fig. 3. The resulting butterfly is shown in Fig. 3b and it qualitatively resembles the Hofstadter butterfly ( Supplementary Fig. 2).
In accordance with Figs 2 and 3a, for small magnetic fluxes j/N the butterfly in Fig. 3b features distinct Landau levels with edge states in the gaps between them. These edge states correspond to red points in Fig. 3b. At high magnetic fields the Landau levels split, but the spectrum still retains a surprisingly delicate structure.

Polariton-polariton entanglement
The internal structure of the two-polariton states is represented by their entanglement entropy 60 , S ¼ À P N ν¼1 jλ ν jln jλ ν j, obtained from the Schmidt expansion ψ nm ¼ P N  (2) also have an intrinsically low entropy, being just a product of a standing wave and a localized or an edge state. However, the states Eq. (2) can mix with each other resulting in larger entanglement entropy. This entangled mixing becomes especially prominent for the Landau level states, cf. points LL, ClS, and ES in Fig. 4b. When the real part of energy approaches ω 0 from the negative side, the mixing between different standing waves increases since the spectrum gets denser, and the states become chaotic-like, see also the top right corner of Fig. 1d. These chaotic-like states are characterized by an irregular wavefunction in the real space and a dense Fourier spectrum in the reciprocal space; they can not be reduced to a product of just two single-polariton states. The specific mechanism of chaotization requires a separate study. At very small negative energies ε − ω 0~− φΓ 0 the single-particle dispersion changes from ε ∝ −1/k 2 to ε ∝ − (k − π/d) 2 because both polaritons get closer to the Brillouin zone edge, and the fermionic correlations 46 emerge from chaos. The dense group of twopolariton states in the right panel of Fig. 4b, where Re ε > 0, is formed by the interaction with the quasi-superradiant mode with Re ε À ω 0 % 71Γ 0 in Fig. 2a. One of the states in this group has an entanglement entropy even higher than that of the chaotic-like states, as seen by the point labeled MES at Re ε À ω 0 % 35Γ 0 . Similar to the chaotic-like state, the most entangled state can not be reduced to a product of just several single-polariton states, but its wavefunction looks more regular.

DISCUSSION
We have discovered an interaction-induced internal topological order for the two-polariton states in a light-coupled onedimensional atomic array.
Our results indicate that the platform of waveguide quantum electrodynamics has tremendous uncovered potential for quantum simulators of many-body interacting systems. The underlying Dicke-type model demonstrates an incredible diversity of quantum states with different topologies, lifetimes, and entanglement. Its unexplored richness may mean it will become as celebrated as well-known many-body physics models such as the Heisenberg model, Bose-Hubbard model, or the Luttinger liquid. The waveguide-mediated long-ranged couplings, intrinsic for the waveguide QED setup are quite uncharacteristic for traditional quantum systems and there is much more to explore. For example, we have focused here only on the regime of extremely subwavelength distances between the atoms, where two-polariton bound states 58,61 play no role. Polariton-polariton interactions could be even more interesting in Bragg-spaced lattices, where the non-Markovian effects are drastically enhanced [62][63][64] . The ultra-strong coupling regime 43 is also unexamined for these quantum waveguides to the best of our knowledge.
On the more practical side, the discovered topological twopolariton states could be used to engineer quantum optical correlations and protect them against the disorder and decoherence. First, since one of the two polaritons is in a deeply subradiant state, the radiative losses for the considered states can be practically neglected (see more detailed analysis in the Supplementary Information). It may also be possible to suppress the radiative losses more strongly by exploring the quantum regime of recently proposed high-quality states 65 . Second, our calculation (see Supplementary Figs 6 and 7) indicates that the interaction-induced edge states of the polaritons are protected from certain kinds of disorder. While the inhomogeneous broadening of the qubit energies relatively easily destroys topological states, they are quite robust against the fluctuations of the qubit positions. Even though the fine structure of the Landau level states is smeared out by the fluctuations, the edge states remain well localized. It is thus very interesting to study whether the proposed concept can be used to generate complex multi-photon states with built-in nontrivial topology. In the considered twopolariton case, one polariton is localized at the edge and another one is in the standing wave state or vice versa. One can also envisage topological N11N states, where one polariton in the standing wave state induces topological localization of N > 1 polaritons. Such quantum states could find applications in quantum metrology like high-N00N states 66 .

Analytical model for polariton-polariton interactions
In this section, we start from the Hamiltonian Eq. (1) in the main text 47,55 and proceed to derive Eq. (3) that describes the interaction between the two polaritons. Substituting the ansatz Ψ j i ¼ P ψ mn σ y n σ y m 0 j i into the Schrödinger equation HΨ = 2εΨ we obtain the two-polariton Schrödinger equation in the form 47,48 H mn 0 ψ n 0 n þ ψ mn 0 H n 0 n À 2δ mn H nn 0 ψ n 0 n ¼ 2ðε À ω 0 Þψ mn : We note that Eq. (5) with the non-Hermitian Hamiltonian can be rigorously derived using conventional input-output approaches of quantum optics 67,68 . For example, another equivalent way to obtain the same results for double-excited states is to use the Linbdlad formalism for the density matrix of the system 68 . The double-excited eigenstates of the Hamiltonian Eq. (5) correspond to the resonances of the scattering matrix for two photons incident in the waveguide and interacting with the atomic array 47,67 . In order to obtain spatial resolution, the double-excited states could also be probed by addressing individual qubits 41 directly, rather than using the waveguide modes. In this case, it might be useful to explore the quantum tomography approach proposed in ref. 69 . This system is readily solved numerically after the wavefunction ψ is rewritten in the basis of N(N − 1)/2 localized states of the type Our next goal is to go beyond refs. 47,48,58 and obtain Eq. (3). To this end, we notice that 48,58 Here, the matrix ∂ 2 represents the one-dimensional discrete Laplacian (or the operator of discrete second-order derivative). This means that for a vector ψ n with a smooth dependence on n, one has ½∂ 2 ψ n ¼ ψ nþ1 þ ψ nÀ1 À 2ψ n % d 2 ψ n dn 2 : Thus, for a short-period array with φ ≪ 1 the operator K reduces to the second derivative operator. The inverted Hamiltonian K in Eq. (6) is a sparse matrix with only nearest-neighbor couplings. This fact inspires us to perform the transformation that means change of the basis to e iω0djxÀnj=c σ y n 0 j i; where x ¼ 1; 2; N: (8) This basis inherits the distribution of the electric field emitted by a given atom. Indeed, expðiω 0 djx À nj=cÞ is just the Green function for a photon in one dimension. Since the wave equations for the electric field are local, the transformed two-polariton Schrödinger equation will be local as well, i.e., it will involve only sparse matrices. Substituting Eq. (7) into Eq. (5) we find 48,58 where summation over the dummy indices x 0 and y 0 is assumed. Next, we look for the solution to the transformed equation (9) in the form corresponding to Eq. (2) in the main text. Here, one of the two excitations is a single-particle eigenstate of the matrix H with the eigenfrequency Ω j . The corresponding single-polariton eigenfrequency is ω j = ω 0 + Ω j . Using the definition K H À1 , we find The state is normalized as P x ½ψ ðjÞ x 2 ¼ 1. The normalization does not involve complex conjugation, because the original matrix H is not Hermitian but symmetric. As such, its eigenvectors ψ (j) satisfy the non-conjugated orthogonality condition jjj 0 h i P N x¼1 ψ ðjÞ x ψ ðj 0 Þ x ¼ δ jj 0 : Owing to the translational symmetry the vector ψ (j) is just a standing wave 56 : We note, that the ansatz (10) and (12), where the eigenstate ψ (j) does not take into account the interaction effects, works only for the transformed Schrödinger Eq.
We are going to consider strongly localized eigenstates that are orthogonal to the standing wave ψ (j) . Moreover, for relatively small j the function χ changes with x much faster than ψ (j) . Hence, we neglect the terms / ψ ðjÞ;2 x χ x and / ψ ðjÞ x ψ ðjÞ y χ y and find χ x Ω j þ K xy χ y À 2½ðψ ðjÞ;2 ÞK xy χ y ¼ 2ðε À ω 0 Þ Ω j K xy χ y : Taking Eq. (6) into account, we recover Eq. (3) from the main text.

Fourier analysis of the eigenstates
The calculation of the energy spectrum of the Hamiltonian Eq. (1), shown in Fig. 2, is relatively straightforward. The spectrum is found by standard linear algebra techniques, and is described in more detail in the Supplementary Information. However, it is more challenging to extract the butterfly spectrum in Fig. 3b from the spectrum in Fig. 2e. This task requires careful separation of the groups of two-polariton states. We start by performing the Schmidt decomposition of the twopolariton state for all the states that have Re ε < ω 0 with unconjugated orthogonality condition P x ψ ν x ψ μ x ¼ δ μν . The form of decomposition Eq. (15) where left and right singular vectors are the same enforces the symmetricity condition ψ xy = ψ yx . Our analysis of the Schmidt decomposition confirms that most states are well approximated using the two largest singular values λ 1 and λ 2 , that have close absolute values. Keeping only these two terms, we obtain linear combinations of the wavefunctions ψ 1 x and ψ 2 y as u ± n ¼ λ À1=4 1 ψ 1 n ± iλ À1=4 2 ψ 2 n . After that, the two-polariton state can be approximately presented as ψ xy / u þ x u À y þ u À x u þ y : Next, we select one of the two states u þ x , u À y that has lower inverse participation ratio, P ju x j 4 =½ P ju x j 2 2 , which means it is less localized in space. We designate this state as u (free) and the more localized one as u (loc) , perform the discrete Fourier transform u ðfreeÞ ðkÞ ¼ X N x¼1 e Àikx u ðfreeÞ x (16) and calculate the wave vector k max , corresponding to the maximum of the Fourier decomposition. The number of the group can be then determined from the quantization rule where square brackets indicate the rounding to the nearest integer. In order to improve the precision in Eq. (17) for large j, we also characterize the vectors χ by their mirror symmetry. Then we apply Eq. (17) separately for odd and even states with odd and even j, respectively. The results of Fourier transform for N = 200 atoms are shown in Supplementary Figs 3 and 4. Except for very large j close to N, the spectrum is clearly separated into well-defined steps of alternating parity. There is an exception of several outlier states with 26 k max N=π t 41 and smaller values of jRe ε À ω 0 < 10 À1 j that can be seen in lower left corner of Fig. S3. We exclude these states from our analysis and do not put them on the butterfly spectrum Fig. 3b. Another type of outlier states have k max N=π $ 50 and different distinct values of energies in the range 10 0 Γ 0 ÀRe ε þ ω 0 10 1 Γ 0 (lower right quarter of Fig. S3). They are out of the scope of the current manuscript but do not require any special exclusion procedure. Next, we assign each step to a different group of eigenvalues and align the groups with respect to each other. This is done by subtracting the energy with the largest (smallest negative) real part from the energies of the states of each group, ε j,max . In order to keep the points with the highest energy on the semilogarithmic plot after this subtraction, we also add a small value of 1.1 × 10 −4 Γ 0 to all the energies. The result is the butterfly spectrum, shown in Fig. 3b.