Many-body localization in waveguide QED

At the quantum many-body level, atom-light interfaces generally remain challenging to solve for or understand in a non-perturbative fashion. Here, we consider a waveguide quantum electrodynamics model, where two-level atoms interact with and via propagating photons in a one-dimensional waveguide, and specifically investigate the interplay of atomic position disorder, multiple scattering of light, quantum nonlinear interactions and dissipation. We develop qualitative arguments and present numerical evidence that such a system exhibits a many-body localized~(MBL) phase, provided that atoms are less than half excited. Interestingly, while MBL is usually formulated with respect to closed systems, this system is intrinsically open. However, as dissipation originates from transport of energy to the system boundaries and the subsequent radiative loss, the lack of transport in the MBL phase makes the waveguide QED system look essentially closed and makes applicable the notions of MBL. Conversely, we show that if the system is initially in a delocalized phase due to a large excitation density, rapid initial dissipation can leave the system unable to efficiently transport energy at later times, resulting in a dynamical transition to an MBL phase. These phenomena can be feasibly realized in state-of-the-art experimental setups.


I. INTRODUCTION
Quantum light-matter interfaces are being actively investigated, for their many possibilities to explore fundamental science and for applications. However, the full quantum many-body problem of interacting atoms and photons remains highly challenging to understand, not only due to the large Hilbert space, but also its out-of-equilibrium and open nature. Most of our theoretical knowledge is thus restricted to certain regimes. For example, much of our quantum theory, such as to model applications, treats the atomic medium as being a smooth quantum field [1], thus ignoring the potential complexities associated with granularity, disorder averaging, and multiple scattering of light. On the other hand, there has been significant work to deal with granularity and multiple scattering [2][3][4][5][6][7][8][9], but largely restricted to classical or mean-field limits, where potentially strong non-linear optical interactions and the buildup of quantum correlations are ignored.
It is hence important to more fully develop insights into what opportunities and many-body quantum phenomena might arise from the combination of multiple scattering and nonlinear interactions between atoms and light [10][11][12][13][14]. For example, it has been suggested that multiple scattering can be enhanced using ordered atomic arrays and become an important resource for applications, which can lead to spin squeezing interactions [15], improved photon-photon gates [16,17], or prolonged interrogation times for optical lattice clocks through the suppression of spontaneous emission [18,19]. Separately, in the regime of dilute, disordered media, perturbative diagrammatic techniques to account for nonlinear interactions and multiple scattering have only recently been developed [20]. Aside from being a possible resource, the effects of multiple scattering could also impose new, unexpected limits on the performance of disordered quantum light-matter interfaces [21,22].
Here, we propose the existence of and theoretically investigate a non-perturbative, quantum many-body phenomenon in a minimal "waveguide quantum electrodynamics (QED)" model, consisting of disordered two-level atoms or qubits interacting with photons in an infinite one-dimensional (1D) waveguide. In particular, we provide qualitative arguments and quantitative numerical simulations to argue that the system can support a many-body localized phase. The topic of many-body localization (MBL) has generated intense interest recently within condensed matter physics [23][24][25][26][27]. A simple phenomenological Hamiltonian is hypothesized to capture the novel MBL phase [28,29], which allows for predictions of exotic many-body properties such as the failure of a system to thermalize, the absence of transport, and the slow growth of entanglement entropy following a quench. Interestingly, while the concept of MBL is usually formulated with respect to closed systems, our waveguide QED system is intrinsically open, as light emitted by any finite set of excited atoms will be irreversibly lost once passing the boundaries of the atomic ensemble. This results in an additional rich interplay between dissipation and essentially closed, MBL-like dynamics. In particular, we show that an MBL phase occurs at low atomic excitation densities (approximately less than half excited in the large atom limit). The absence of transport then suppresses propagation of energy to the system boundaries and the resulting dissipation, which in turn makes the notion of MBL applicable in the first place. On the other hand, for systems that are initially highly excited, transport can be responsible for a rapid dissipation at early times. This causes the atomic excitation density to drop until transport is no longer efficient, resulting in a dynamical transition into an MBL phase. This physics should be feasible to explore in existing state-of-the-art waveguide QED systems, such as superconducting quantum bits coupled to microwave transmission lines [30][31][32] or structured waveguides [33,34].
The rest of the paper is organized as follows. First, in Sec. II, we introduce the waveguide QED Hamiltonian and present a qualitative argument that an MBL phase should exist in the thermodynamic limit, provided that the atoms are less than half excited. Given the large Hilbert space of the system, consisting of the exponentially large space associated with the atoms and the infinite space associated with the continuum of photon modes, we then restrict ourselves to the regime of near-resonant photon interactions. Then, in Sec. III, we show how the photonic modes can be integrated out, reducing the system to an open, interacting "spin model" describing photon-mediated interactions between the atomic degrees of freedom. Before going to MBL, and to gain intuition, we show how Anderson localization manifests itself in the spin model formalism, within the linear optics (single-excitation) limit. This example illustrates that although the system is formally open, localization implies that dynamics of a large system is well-described just by the Hermitian Hamiltonian component of the spin model. We briefly introduce the phenomenology of MBL phases in Sec. IV. As MBL is typically studied with respect to closed systems, and given the observation above regarding the connection between localization and Hermitianity, in Sec. V we present numerical evidence of MBL for the Hermitian Hamiltonian, namely the absence of transport and logarithmic growth of entanglement entropy for sufficiently low excitation densities. In Sec. VI we then consider the physical open system, and propose and numerically investigate a number of realistic observables for MBL. At low excitation densities, we again find an MBL phase. Surprisingly, however, for higher excitation densities, where MBL is not observed for the closed system, we find that transport to the system boundaries and the corresponding dissipation enable the system to initially and rapidly lower the excitation density, and thus dynamically transition into an MBL phase. Finally, in Sec. VII, we provide an outlook of possible future interesting directions to explore, and discuss the prospects of experimentally observing such physics in state-of-the-art waveguide QED platforms.

II. DESCRIPTION OF THE SYSTEM
Our system of interest is represented in Fig. 1(a). It consists of a one-dimensional, bi-directional waveguide, whose photons couple to the ground-excited (|g -|e ) transitions of N two-level atoms located at random positions z i . The atomic transition frequency is given by ω 0 . A minimal model for such a system is given by the following Hamiltonian (with = 1): Here, H a describes the excited-state energy of the atoms, and we adopt the notationσ i αβ = |α i β i | with {α, β} ∈ {g, e} for the operators of atom i. H ph describes the energy of the photons, characterized by a continuum of modes with two possible propagation directions (labeled ν = ±) and wavevector k. We assume that within the bandwidth of modes to which the atoms significantly couple, the dispersion relation for the guided modes can be linearized as ω k = c|k|. The interactions, as given by H int , describe the processes by which an excited atom can emit a photon, or a ground-state atom can absorb a photon, with a coupling constant g assumed to be identical for all atoms, and an interaction phase e iνkzi encoding the photon propagation. In the following, we consider disorder in the positions of the atoms: z i = (i + i )d, with d the average distance between two neighboring atoms and i a random variable. For numerics, we will focus on the regime of full disorder, where i is drawn between −1/2 and 1/2, in order to minimize the effective localization length of the system and the number of atoms needed in simulations, although we believe our conclusions are general for any amount of disorder (see below). The Hamiltonian of Eq. (1) can be an excellent approximation for a number of systems such as superconducting qubits coupled to a transmission line [30][31][32], or atoms coupled to a photonic crystal waveguide [35,36], where additional sources of dissipation (not included in this model) can potentially be much smaller than the coherent atom-waveguide interactions.
To gain some basic insight into Eq. (1), we begin by relating it to well-known results regarding Anderson localization of light in 1D, in the linear optical or single-excitation limit with disorder in the position of the atoms. (a) Schematic illustration of a spatially disordered atomic chain coupled to photons in a one-dimensional, bidirectional waveguide. The magnitude of atom-light interactions is set by the rate Γ1D by which a single, isolated excited atom irreversibly emits a photon into the waveguide. With multiple atoms and excitations, the system dynamics becomes complex and interacting due to the following possibilities. First, a photon can be reflected or transmitted with non-trivial, frequency-dependent amplitudes r(ω), t(ω) upon scattering off of an atom, which then leads to multiple scattering in the presence of multiple atoms. Second, multiple photons (e.g., two photons with frequencies ω1,2 illustrated) scattering off of the same atom can lead to the generation of entangled, frequency-mixed photons (ω3 + ω4 = ω1 + ω2), which also become multiply scattered. (b) Proposed phase diagram in the thermodynamic limit of large atom number N → ∞, versus amount of disorder and fraction of atoms excited, nexc/N . For a single excitation, nexc/N = 1/N , the system is non-interacting, and arbitrarily small disorder gives rise to Anderson localization. Beyond a single excitation, the system becomes interacting due to the nonlinearity of two-level atoms. Up to nexc/N = 1/2, we hypothesize that the system exhibits MBL. For nexc/N > 1/2, photon transport is allowed and leads to delocalization. Going beyond the thermodynamic limit, and considering a finite atomic system, transport in the delocalized phase leads to dissipation, in the form of emission of light beyond the atomic system boundaries. The subsequent loss of excitation density can lead to a dynamical transition into the MBL phase, as depicted by the dashed arrow showing time evolution. (c) Schematic illustration of a spatially disordered atomic chain coupled to photons in a one-dimensional, "half" waveguide, with one end closed by a mirror at z = 0. The atomic system can thus only dissipate by emitting photons from one side of the chain.
The single-photon, single-atom scattering dynamics of Eq. (1) can be exactly solved [37,38]. In particular, the response of a ground-state atom to an incoming photon of well-defined frequency ω is characterized by the reflection and transmission coefficients r(ω) = −Γ 1D /(Γ 1D − 2i∆) and t(ω) = 1 + r(ω) ( Fig. 1(a), with ∆ = ω − ω 0 being the detuning between the photon and the atomic transition, and Γ 1D = 4πg 2 /c being the spontaneous emission rate of a single, isolated excited atom into the waveguide. When many atoms are coupled together via the waveguide, the multiple scattering arising from the succession of transmission and reflection events ( Fig. 1(a) by different atoms can be solved by the transfer matrix formalism [39]. Performing a disorder average over the position fluctuations i of the N atoms, the 1D waveguide is always Anderson localized for arbitrarily small disorder (arbitrarily narrow distributions of i , see Fig. 1(b) [40][41][42]. The transmittance exhibits exponential attenuation with large fluctuations log(T tot ) = −N/N loc and var(log(T tot )) = 2N/N loc , with (...) denoting the average value over the different disorder realizations and var(...) the variance of the random variable. N loc is the localization length (expressed in terms of number of atoms), and is related to the single-atom linear transmittance T (∆) = |t(∆)| 2 by N loc (∆) = 1/| log T (∆)| = | log(4∆ 2 /(Γ 2 1D + 4∆ 2 ))| −1 in the dilute regime (d > λ = 2πc/ω) and for large disorder [39]. This quantity goes to zero at resonance, and increases with increasing detuning |∆| as the atomic response becomes weaker away from resonance. For systems larger than the localization length N > N loc , strong destructive interference suppresses the propagation of light through the medium.
While rigorously establishing an MBL transition beyond the single-excitation limit is more complicated, we can nonetheless make a qualitative argument, by calculating a nonlinear single-atom transmittance T (∆, ρ ee ) that depends on the excited state population ρ ee , and then substituting into the Anderson localization result N loc (∆, ρ ee ) ∼ 1/| log T (∆, ρ ee )| to estimate a length scale over which transport becomes prohibited. As discussed in the Appendix, this transmittance can be exactly calculated for an incident coherent state field of arbitrary amplitude, with corresponding Rabi frequency Ω. One finds the relations ρ ee = Notably, for large photon input (Ω → ∞), one finds that ρ ee → 1/2 and T → 1. This expresses the well-known result that an atom becomes saturated at very high intensities and is no longer able to respond to light. In turn, reflection and multiple scattering are suppressed, and the localization length N loc (∆, ρ ee → 1/2) → ∞ diverges. We thus hypothesize that an atomic excitation density of n exc /N = 1/2 sets the MBL-delocalization transition, as illustrated in Fig 1(b).

III. SPIN MODEL
Beyond a single excitation, directly solving Eq. (1) is difficult due to the infinite Hilbert space of the photons and frequency mixing induced by the atoms. Furthermore, the localization length N loc (∆, ρ ee ) increases both with increasing detuning and atomic population, and can rapidly exceed numerically tractable system sizes (see Appendix).
To partially mitigate these challenges, we will restrict ourselves to the regime in which the photons involved in the dynamics are near resonance with the atoms. As discussed in Ref. [43], since the field part of the Hamiltonian (1) is quadratic, it can be formally integrated out, resulting in photon-mediated interactions between atoms. Furthermore, near resonance, the delay time of these interactions due to the speed of light can be ignored to good approximation. This is because the atoms in typical situations have very large ratios of resonance frequencies to linewidths, ω 0 Γ 1D . Thus, the atoms are very dispersive, for example, as characterized by the large, frequency-dependent phase shifts in r(ω), t(ω), and this causes the propagation delay of near-resonant photons to be dominated by interaction with the atoms, rather than the speed of light itself. Equivalently, the hybrid atom-photon polaritons that diagonalize the system are in fact almost entirely atomic in nature. Ignoring retardation, one obtains a reduced master equation describing instantaneous photon-mediated interactions, between only the atomic "spin" degrees of freedom [43][44][45][46][47] with and Γ i,j = Γ 1D cos(k 1D |z i − z j |) and k 1D = ω 0 /c. Here, we have transformed to a rotating frame, so that the trivial phase evolution associated with H a can be ignored. We note that in one dimension, the photon-mediated interactions between atoms are infinite in range. Furthermore, integrating out the photons results in a dissipative (open) system, as atoms can physically emit a photon that goes beyond the atomic system boundaries, resulting in irreversible loss of excited population. If loss is indeed associated with the system boundaries, then conceptually its effects (for a fixed initial number of excitations) can be reduced by considering progressively larger system sizes, where the ratio of "boundary" (e.g., atoms within a distance ∼ N loc of an edge) to "bulk" regions decreases. While the maximum number of atoms in numerics is constrained, we can nonetheless decrease this ratio by introducing a closely related "half-waveguide" model, as illustrated in Fig. 1(c). Here, one end of the waveguide (say at z = 0) is terminated by a perfect mirror, which results in only one open boundary where photons can escape. Applying similar considerations as the derivation of Eq. (2), it is straightforward to show that the master equation iṡ with and While Eqs. (2) and (4) have now eliminated the infinite Hilbert space of photons, an interesting new feature emerges: dissipation. In particular, MBL is usually formulated with respect to closed systems. Furthermore, in experimental platforms to investigate MBL thus far, any degree of dissipation is viewed as an unwanted imperfection on top of an idealized closed system. Here, in contrast, the dissipative part of the master equation a priori has the same strength (∼ Γ 1D ) as the coherent interactions. A key idea that we aim to establish is that within an MBL regime, the bulk region of a many-body system behaves as if it is increasingly closed as it becomes further in distance from a boundary, as suppressed transport makes the emission of a photon past the system edges improbable. These ideas can be better clarified by returning to the simpler problem of Anderson localization, but now studied from the spin model perspective.
A. Disorder in the single-photon limit The key single-excitation properties are encoded in the N single-excitation eigenmodes and eigenvalues of the non-Hermitian Hamiltonian H 1D . The eigenvalues themselves are complex, with the real and imaginary parts accounting for the shift in resonance frequency ω ξ of the collective mode with respect to the bare atomic frequency, and half of the collective decay rate Γ ξ /2, respectively.
As a concrete example, we first take one single realization of disorder, and sort the eigenstates in increasing order of their resonance frequencies, labeled by 1 ≤ ξ ≤ N . Denoting the wave functions by |ψ ξ = j c ξ jσ j eg |g ⊗N , in Figs. 2(a,b) we plot |c ξ j |, the square root of the probability for atom j to be excited, for the atom numbers of N = 50 and N = 200, respectively. For all numerical calculations, we take an average atomic spacing of d = 2.7π/k 1D (to be in the dilute regime) and full disorder. We see that states in the middle of the spectrum (ξ ∼ N/2) are localized, while states at the edges of the spectrum (ξ ∼ 1 and ξ ∼ N ) are extended. Furthermore, in Fig. 2(c), we fix the atom number N = 50, and plot the average value of the resonance frequency ω ξ /Γ 1D versus ξ over a large number of disorder realizations (blue curve). We simultaneously plot (orange curve) the previously established Anderson localization length N loc (ω ξ ) corresponding to these frequencies. Comparing with Fig. 2(a), we see that the transition from localized to extended eigenstates occurs when |ω ξ | is detuned enough that the localization length N loc (ω ξ ) > N exceeds the number of atoms inside the chain. In particular, the existence of extended states of the spin model (such as found in Refs. [48][49][50]) does not contradict the fact that Anderson localization always exists in the 1D disordered system. Specifically, for any detuning ∆ there always exists some sufficiently large N > N loc where one will eventually observe Anderson localization.
In Fig. 2(d), we plot the average value over a large number of disorder realizations of the decay rate Γ ξ /Γ 1D , normalized by the single-atom decay rate, as a function of ξ/N for two different system sizes N = 50 and N = 200. We see that the extended modes in the edges of the spectrum have a large decay rate, while the localized modes in the middle have an exponentially small decay rate. This confirms that localized modes, with populations away from the boundaries, cannot efficiently decay by spontaneous emission. Furthermore, as N increases, we find that the percentage of delocalized modes [here quantified by the percentage of modes with a decay rate greater than Γ 1D /2, above the green line in Fig. 2(d)] decreases approximately as ∼ 1/ √ N , as shown in the inset. While this analysis was restricted to the single-excitation manifold, it naturally also follows that for our particular model, the existence of MBL cannot be inferred purely by looking for the localized nature of eigenstates of the spin model with large number of excitations. In particular, a study along those lines would have to determine whether an apparently delocalized state might become localized as the system size is increased. As that is difficult within numerical capabilities, we must demonstrate the existence of MBL using other measures. Separately, we note that while localized and extended single-excitation eigenstates can be identified by spectral filtering (i.e. looking at eigenstates at the center and edges of the spectrum, respectively), this procedure cannot be extended to more excitations. For example, assuming that interactions can be considered as a perturbation, two singleexcitation extended eigenstates with energies ±ω could combine to give an extended two-excitation eigenstate with approximately zero energy.
We can again easily confirm this in the single-excitation limit. Specifically, for a single disorder realization and N = 150, we calculate the eigenstate amplitudes |c ξ j | for the non-Hermitian and Hermitian Hamiltonians, which we plot in Figs. 3(a,b), respectively. Then, in Fig. 3(c), we plot the overlap D(ξ, ξ ) = | ψ (1D) ξ |ψ (hermi) ξ | between the non-Hermitian and Hermitian eigenstates. We find that this matrix is nearly the identity, particularly in the middle of the spectrum where the modes are well-localized [51].
This suggests that we can establish the existence of an MBL in two complementary ways. First, we can take the Hermitian Hamiltonian H hermi , and look for well-established mathematical signatures for closed systems, such as logarithmic growth of entanglement entropy following a quench. Then, assuming that the existence of MBL is not affected by losses for the reasons above, we can simultaneously look for realistic observables in the physical, open system.

IV. INTRODUCTION TO CONVENTIONAL MBL
The key properties of MBL can be understood from a "canonical" Hamiltonian hypothesized for all MBL systems [28,29]. For spin systems like ours, it reads Here, theτ i z operators are so-called local integrals of motion, pseudo-σ z operators that are quasi-local in space. Because H LIOM only involves products ofτ z operators, eachτ j z commutes with the Hamiltonian and thus their occupancies are conserved quantities. Furthermore, as these operators only have support on a few sites, there will be no transport of energy, and an MBL system prepared initially out of equilibrium will never thermalize.
The interacting terms (e.g., J i,j for two-body interactions) have exponentially decreasing amplitude with the distance between modes [28,52,53]. The interactions differentiate MBL from Anderson localization, and cause each local integral of motionτ i z to acquire different phases in evolution depending on the occupancy of otherτ j z . For closed systems, this results in a dephasing for any local subsystem, due to the gradual entanglement of these degrees of freedom with others further away. Likewise, starting from a product state in the physical basis, these interactions cause a subsystemρ A consisting of half the entire system to exhibit a logarithmic growth of entanglement entropy S(t) = −Tr[ρ A (t) log(ρ A (t))] in time in the MBL phase [54]. As a comparison, systems that do not exhibit MBL experience a ballistic growth of entanglement entropy (linear with time).

V. MANY-BODY LOCALIZATION IN THE HERMITIAN WAVEGUIDE MODEL
In this section, we study the Hermitian component H hermi of the effective Hamiltonian, looking at transport and entanglement entropy. Specifically, we consider initial states consisting of product states of n exc excitations, |ψ init = j∈Einitσ j eg |g ⊗N with E init being the set of indices of the initially excited atoms. Enforcing a fixed number of excitations enables us to study modestly larger system sizes, by throwing out the Hilbert subspace with higher excitation numbers.
We first consider transport, in a system of N = 30 atoms. In Fig. 4(a), we show the time-dependent excitedstate population σ i ee (t) per site, averaged over ∼ 20 − 100 disordered configurations, for n exc = {1, 3, 5} roughly equidistantly spaced excitations. In comparison, in Fig. 4(b), we show the evolution for the same initial states, but in an ordered chain of lattice constant d = 2.7π/k 1D . In Fig. 4(c) we plot the site-dependent population σ i ee (t max ) for the final time of the simulation, t max = 100/Γ 1D , both for the disordered (blue) and ordered cases (red). In the ordered system, the population rapidly becomes distributed evenly among all the spins, so that σ i ee (t max ) ≈ n exc /N (dashed black line). (The only exception is for n exc = 1, where the single excitation evolves according to a well-defined band dispersion relation). In contrast, in the presence of disorder, the system retains a clear memory. This can be quantified by defining an observable that captures the memory of the initial state: This quantity would be equal to 1 in the case of an initial population distribution conserved through time, and to 0 in the case of an equally shared population between all atoms. In Fig. 4(d), we plot the disorder-averaged M(t) and observe a clear memory of the initial state (or, equivalently, an absence of transport at long times) in all the cases with disorder.
In Fig. 5 we study the evolution of the half chain entanglement entropy S(t) in time, for the same disordered system and initial conditions. A clear region of logarithmic growth is observed when the number of excitations is large enough (up to the largest number n exc = 5 we can simulate in a system of N = 30 atoms), before saturating due to the limited system size. Figures 4 and 5 provide evidence that the Hermitian system, as defined by H hermi , exhibits a many-body localized phase up to an excitation density of at least ∼ 1/6.

A. Delocalization transition
As discussed in Sec. II, based on qualitative arguments, we expect a disordered waveguide QED system to exhibit an MBL phase up to an excitation density of ρ ee = n exc /N = 1/2 in the thermodynamic limit. However, we also expect that the localization length N loc (∆, ρ ee ) should gradually increase as the excitation density increases, eventually diverging as ρ ee → 1/2. Thus, for a finite system of moderate size as can be simulated, this would lead to a smooth crossover to delocalization as n exc is increased.
To observe this, we consider a smaller chain of N = 20, and repeat the same calculations as before for the Hermitian waveguide, but now for 1 ≤ n exc ≤ 7 to reach a higher excitation density. In Fig. 6(a) we plot the transport properties for n exc = 3, 5 and 7 excitations. It can be seen that the memory of the initial state becomes negligible already for n exc = 5, and essentially vanishes for n exc = 7, as the initial population is eventually equally shared among all atoms in the waveguide. In Fig. 7 we plot the half-chain entanglement entropy, and see a smooth transition toward ballistic growth as the number of excitations exceeds n exc 4 − 5.
While a delocalization transition that occurs at a lower density of excitations, n exc /N ∼ 1/4, cannot be ruled out from these numerics, the results are also consistent with our previous hypothesis of a smooth crossover. In particular, assuming that the dynamics have a characteristic bandwidth of ∼ Γ 1D , one finds a heuristic many-body localization length N loc (∆ = Γ 1D , ρ ee ) that exceeds the system size of N = 20 for n exc 4.

A. Quantum transport
We now investigate practical signatures of MBL in the physical, open system. Anticipating that regions within the excitation density-dependent localization length ∼ N loc (∆, ρ ee ) of the boundaries are strongly affected by dissipation, we reduce their relative contribution in numerics by going to the open half-waveguide, as described by Eq. (4), which only has one radiating boundary located to the right.
We begin by studying transport in the open system. Based on the intuition that the system will resemble a closed system (and thus exhibit stronger features of MBL) provided that excitations remain far from the boundary, we slightly alter the initial conditions compared to Sec. V, and take initial states consisting of n exc excitations all located on the left half of the chain. To eliminate any bias originating from a specific choice of location within the left half, we consider equal superpositions of all possible basis states, with each basis state multiplied by a random phase. We have checked, however, that other choices do not affect the final conclusions. Numerically, we directly solve for the master equation dynamics of Eq. (4), for n exc ≤ 2, and by a quantum jump approach for higher excitation number, with a minimum of 150 trajectories for each configuration, and with at least 50 random configurations in all cases to calculate disorder-averaged observables.
We then define the time-dependent imbalance M half (t) between the left and right halves, where P e (t) = i σ i ee (t) is the total excited population that now decays in time. In Fig. 8 (left panel) we plot the disorder averaged, time-dependent imbalance for various initial excitation numbers n exc = 1, 2, 4, 6 and for different system sizes (solid lines). Alongside the imbalance, we also plot P e (t).
Beginning with a low density of excitations (n exc = 1, 2) in the disordered system, we see that both transport and dissipation are strongly suppressed, with the asymptotic behavior apparently trending toward perfect suppression, M half → 1 and P e → n exc , as N → ∞. This confirms the increasingly closed nature of the system. Similar trends also are observed for a higher number of excitations, n exc = 4. Here, however, the localization length N loc (∆ = Γ 1D , ρ ee ) (again assuming a characteristic bandwidth ∼ Γ 1D for dynamics) becomes comparable to the system sizes simulated, and dissipation is no longer negligible. Interestingly, one sees that the imbalance does not monotonically decay, but rather reaches a minimum over a short time scale, while partially recovering at later times. This behavior is more prominent for smaller N , due to the higher initial excitation density. Physically, part of the initial population on the left experiences transport and quickly propagates toward the right boundary of the chain (decreasing the imbalance), where it can be dissipated. The lowering of density on the left half subsequently suppresses the transport, and the system dynamically enters into an MBL phase, while the imbalance partially recovers as the remaining excitations on the right half then slowly but preferentially dissipate away. Given this observation, we are then motivated to investigate even higher initial densities (n exc = 6 and N = 16), where now the localization length N loc (∆ = Γ 1D , ρ ee = 3/8) ≈ 65 N far exceeds the system size. Despite clearly starting in a transport and dissipation-allowed phase, the same dynamical behavior is clearly and even more prominently observed.
We now return to the hypothesized phase diagram of Fig. 1(b), which we initially considered to be for a system in the thermodynamic limit. Considering that any physical system will consist of finite atom number, we propose that its dynamical behavior, if it begins in the delocalized phase, will consist of transport-facilitated dissipation, until it reaches an MBL phase, as illustrated by the arrow showing evolution in time. Furthermore, this dissipation process should be quite rapid and dominate the initial dynamics, due to the lack of subradiant states (with decay rate Γ 1D ) for excitation densities n exc /N 1/2, and by the fact that the average dissipation rate of eigenstates of the non-Hermitian Hamiltonian for such high excitation densities must grow extensively with system size, ∼ N Γ 1D .
Finally, for comparison, we plot with dashed curves the same observables for all the system sizes and number of excitations previously considered, but for an ordered system. It should be noted that the total population in the ordered system also decays slowly for long times, albeit faster than the disordered case. This can be attributed to the large number of single and multi-excitation subradiant states in a 1D waveguide at low excitation density [47]. In the case of a full waveguide (dissipation at both boundaries) and sufficient number of excitations, it has previously been shown that this leads to a power-law decay for ordered systems at long times [18].
FIG. 8. Time-dependent, disorder-averaged population imbalance M half (t) (left panels) and total population Pe(t) (right panels) for the open half-waveguide, for nexc = 1, 2, 4, 6 initial excitations, and for various system sizes N (indicated by different colors). We plot these quantities both for the disordered (solid curves) and ordered (dashed) systems.

B. Quantum revivals
As entanglement entropy is difficult to measure in a closed system of large size, and has no immediate analogue in an open system, capturing some essence of the slow entanglement entropy growth is generally one of the most challenging aspects of experimentally studying MBL [55]. When the local integral of motion closely resembles the physical spin, it has been theoretically proposed that the slow growth of entanglement can be probed by generalized spin echo techniques [56]. In our system, however, we find that this technique is ineffective, as localized modes have a non-negligible extent over at least a few sites (see the single-excitation dynamics in Fig. 4, for example). Physically, this is because suppressing the emission of light by one excited atom arises from the excitation of at least the neighboring atoms and their destructive interference in radiation.
We thus adopt an alternative proposal by Vasseur et al. [57], which is readily implementable within waveguide QED. In particular, we consider the addition of an extra ancilla atom A that does not directly couple or radiate into the waveguide, but is allowed to coherently exchange excitations with the waveguide atom with index n c [see Fig. 9(a)]. The waveguide atoms thus serve as an exotic bath for the ancilla. The corresponding Hamiltonian for this system (studied for the half-waveguide) is with the non-Hermitian Hamiltonian H half given in Eq. (5) and H c = C σ A egσ nc ge + h.c. coupling the ancilla to atom n c in the chain. In the following we use C = Γ 1D /2 and we couple the ancilla to the atom number n c = 3 in the chain of N − 1 atoms.
If the system was governed purely by H c , an ancilla initially excited would indefinitely and reversibly exchange its excitation with atom n c . In the Anderson localized regime, where the waveguide interactions are turned on but all waveguide atoms are initially unexcited, the ancilla couples to a few localized modes (those whose support contain atom n c ). This leads to a more complex structure of collapse and revival dynamics in the ancilla population over time, as illustrated in Fig. 9(b) (calculated here for the Hermitian half-waveguide, with no dissipation). Although the details of the revivals depend on the disorder configuration, since the number of localized modes to which the ancilla couples is finite and fixed, the revivals would continue with the same frequency and amplitude indefinitely in the absence of dissipation, and for an increasingly long time with increasing system size in the presence of dissipation. In contrast, in the MBL phase where multiple waveguide atoms are initially excited: with |ψ wg containing n exc − 1 excitations, interactions between excitations effectively cause the ancilla to couple to an increasing number of degrees of freedom in time, which makes revivals less likely. This is illustrated in Fig. 9(c), for the case of n exc − 1 = 2 excitations initialized in the Hermitian half-waveguide as well, at positions |e 1 , e 4 . For closed systems, the revival rate should approximately be inversely proportional to the effective size of the Hilbert space to which the ancilla is coupled. We can thus numerically study the revival rate R(t) = N rev (t)/t (see Appendix B for how the number of revivals N rev (t) in the time window [0, t] is determined). To do this, we again directly solve the master equation for low excitation numbers n exc = 1, 2, and by the quantum jump formalism for higher numbers (with at least 500 trajectories for each configuration). We perform averages over at least 50 configurations for each choice of n exc , N , and the initial waveguide state is chosen to be in an equal superposition of all basis states with n exc − 1 excitations, multiplied by random phases.
To quantify the evolution of the average number of degrees of freedom N (t) to which the ancilla is coupled by interactions, we take the ansatz that where α is a proportionality constant. We further split N (t) = N 0 (t) + N int (t) with N 0 (t) the average number of degrees of freedom that interacts with the ancilla without interactions inside the waveguide and N int (t) the average number of degrees of freedom that comes from interactions between the different qubits. N int (t) should grow logarithmically in time for closed, large MBL systems and linearly in time for delocalized systems, while N 0 (t) saturates to a constant value that depends on the localization length N loc . Thus, using the unloaded waveguide as a reference to obtain R 0 (t), we can extract from the revival rates the average number of interacting degrees of freedom (up to a proportionality factor), by considering the quantity We plot O(t) for different initial excitation numbers n exc = 2, 3, 4, 6 and atom number N = 14, 17, 24 in Fig. 9, both for the Hermitian (d) and open systems (e) up to Γ 1D t = 300. Starting with the Hermitian half-waveguide system, one clearly observes a transition from logarithmic to linear growth of O(t) as the system goes from a MBL to a delocalized phase when the density of excitations increases. Interestingly, for the open system, we observe a slow logarithmic growth of O(t) in all the cases considered over this time range. We again attribute this to a dynamical transition, where transport results in rapid dissipation of large excitation densities, until the waveguide atoms reach an MBL phase. This transition occurs faster than the characteristic interaction rate of the ancilla with the bath (taken to be comparable to the single atom-waveguide coupling strength Γ 1D ), resulting in suppression of quantum revivals reminiscent of a closed MBL system. Separately, we have seen that for sufficiently long times Γ 1D t 500, a deviation from logarithmic scaling can be observed, due to the very slow dissipation of remaining excitations left in the system.

VII. CONCLUSION
In summary, we have proposed and presented numerical evidence that a system of disordered atoms coupled to a waveguide exhibits an MBL phase, provided that the density of atomic excitations is less than 1/2. Compared with many other MBL systems already studied, this system has a number of interesting features. In particular, the system contains two types of particles (atomic spins, and a continuum of photons), and furthermore, the continuum nature of the photon modes intrinsically causes the atomic dynamics to appear open, with a dissipation strength that a priori is equal to the coherent interaction strength, as seen in Eq. (2). Thus, beyond typical MBL systems where dissipation is simply an unwanted effect, our system presents interesting opportunities to study the apparent transition toward Hermitianity [51] when the system is localized, and a dynamical transition toward MBL induced by loss.
State-of-the-art waveguide QED systems involving superconducting qubits coupled to transmission lines should already be capable of studying the proposed phenomena. In particular, it has been demonstrated that systems consisting of at least N = 7 controllable, individually measurable, and identical qubits can be realized [32], with ratios of atom-waveguide interaction strengths to additional, unwanted dissipation of Γ 1D /Γ ∼ 10 2 -10 3 [32,58] that enable the predicted dynamics to be observed before additional effects set in. Furthermore, since our proposed scheme relies on disorder, the use of identical qubits is not necessary, and up to N = 72 qubits have been realized in such an instance in similar systems [59]. Coupling between superconducting qubits and microwave photons can also be realized in two dimensions [60], thus offering promising opportunities to investigate MBL, including in regimes beyond what can be studied directly with numerics. We note that superconducting qubit systems are already being used to investigate MBL [53,61,62], albeit in regimes where photons are not a central degree of freedom. Moreover, although still in their infancy, systems consisting of atoms coupled to photonic crystal waveguides can also potentially reach the desired combinations of large atom-waveguide interaction strengths [35,36,63,64] and large atom number to investigate MBL.
Beyond our initial theoretical investigations, our work also stimulates other theoretical questions to explore. First, while we have focused solely on position disorder (and full disorder in the numerics, to minimize the localization length), our qualitative arguments about the existence of MBL seem quite general. We thus envision future efforts to confirm a thermodynamic phase diagram similar to Fig. 1(b), for arbitrary amounts and types of disorder (e.g., in resonance frequencies). Furthermore, thus far, we have reduced the complexity of our system by integrating out the photons and focusing on the atomic "spin" degrees of freedom. While this is an excellent approximation near resonance, it would be interesting to more fully explore the system from the photonic standpoint. For example, we anticipate that the MBL phase is reflected in interesting quantum correlations of light, either generated through the excited atoms themselves, or explicitly via quantum transport by sending in optical pulses. Including the photons might also provide an avenue to develop diagrammatic techniques [20,65,66] to understand MBL and the delocalization transition, and in a way that is not as limited by system size as with pure numerics. Finally, while we have considered the most basic continuum of photon modes here, consisting of a linear dispersion and infinite bandwidth, current waveguide QED systems also offer excellent potential for dispersion engineering, such as through the introduction of band edges and gaps [67] or even its global shape [60,68], and other features such as realizing some degree of chirality in interactions [69]. These can dramatically alter the nature and the range of the photon-mediated interactions, and result in non-trivial boundaries between MBL and delocalized phases.
with L[ρ] = 2σ ge ρσ eg −σ ee ρ − ρσ ee , one obtains the optical Bloch equations whose steady state solutions are: and ρ eg = iΩ Γ 1D + 2i∆ 4∆ 2 + Γ 2 1D + 2Ω 2 . (A4) The input-output equation for a waveguide [43] allows one to express the transmitted field aŝ whereÊ in represents the (coherent state) input field. This leads to the expression of the transmittance T = Ê †Ê / Ê † inÊ in : From this single-atom transmittance, we can estimate the saturation and frequency dependent localization length using N loc (∆, ρ ee ) = 1/| log(T )|, with N loc being the localization length in units of number of atoms. In the equation above, the excited-state population ρ ee parametrically depends on the Rabi frequency via Eq. (A3). In a waveguide QED system, one can estimate that the characteristic bandwidth of MBL dynamics is given by Γ 1D . One can then evaluate the excited-state population ρ ee and transmittance T at a detuning ∆ = Γ 1D (illustrated in Fig. 10(a), as a function of Ω/Γ 1D ), and subsequently calculate the excited-state population-dependent MBL localization length ∼ N loc (∆ = Γ 1D , ρ ee ), which we plot in Fig. 10(b). Note that the localization length exceeds the maximum system sizes we can numerically study N ∼ 20-30 for excitation densities of n exc /N ∼ 0.3.

Revival counting
In the quantum revival simulations, one has to count the revival rates of the ancilla in order to extract the number of interacting degrees of freedom in the waveguide. For each realization of the disorder, we obtain the population of the ancilla σ A ee (t) as a function of time from which we extract all the times t n corresponding to an extremum of σ A ee (t) . Our algorithm is designed to prevent local maxima with arbitrarily small visibility from qualifying as true revivals. To this end, for each local maximum, we compute Q = σ A ee (tn) − σ A ee min,n−1 σ A ee min,n−1 , which compares the value at the position of the maximum with the minimum value of the population of the ancilla, during the time since the last revival maximum. We count the local maximum at t n as a true revival if Q exceeds a prescribed value Q min , and if the population of the ancilla σ A ee (t n ) exceeds 0.25. For each disorder realization, we then define N rev (t) as the cumulative number of revivals of the population in the time window [0, t]. Then we average N rev (t) over the disorder in order to obtain R(t) = N rev (t)/t (note that the details of the revival dynamics vary for each configuration, so one cannot obtain R(t) from the disorder-averaged population of the ancilla σ A ee (t) ). All the curves presented in this work have been obtained taking Q min = 0.4.