Subradiant states of quantum bits coupled to a one-dimensional waveguide

The properties of coupled emitters can differ dramatically from those of their individual constituents. Canonical examples include sub- and super-radiance, wherein the decay rate of a collective excitation is reduced or enhanced due to correlated interactions with the environment. Here, we systematically study the properties of collective excitations for regularly spaced arrays of quantum emitters coupled to a one-dimensional (1D) waveguide. We find that, for low excitation numbers, the modal properties are well-characterized by spin waves with a definite wavevector. Moreover, the decay rate of the most subradiant modes obeys a universal scaling with a cubic suppression in the number of emitters. Multi-excitation subradiant eigenstates can be built from fermionic combinations of single excitation eigenstates; such"fermionization"results in multiple excitations that spatially repel one another. We put forward a method to efficiently create and measure such subradiant states, which can be realized with superconducting qubits. These measurement protocols probe both real-space correlations (using on-site dispersive readout) and temporal correlations in the emitted field (using photon correlation techniques).


Superconducting qubits coupled to photons propagating in open transmission lines
offer a platform to realize and investigate the fascinating world of quantum light-matter interactions in one dimension -socalled "waveguide quantum electrodynamics (QED)" [4][5][6][7][8][9][10]. Such systems enable a number of exotic phenomena that are difficult to observe or have no obvious analogue in other settings, such as near-perfect reflection of light from a single resonant qubit [1,2,4,[11][12][13] or the dynamical Casimir effect [14]. One particularly interesting feature of these systems is that the interaction between multiple qubits, mediated by photon absorption and reemission, is of infinite range, as the mediating photon is unable to diffract energy into other directions in a onedimensional channel. This can give rise to strong collective effects in multi-qubit systems [8,15]. For example, it has been observed that two qubits separated by a substantial distance can exhibit super-or sub-radiance, wherein a single collective excitation can decay at a rate faster or slower than that of a single qubit alone [3].
The physics associated with collective effects in waveguide QED has attracted growing interest, and there have been a number of proposals that implicitly exploit suband super-radiant emission to realize atomic mirrors [16], photon Fock state synthesis [17], or quantum computation [18,19]. However, the fundamental properties of the collective qubit excitations themselves, such as their spatial character and decay spectrum, have been less systematically studied [20].
Here, we aim to provide a systematic description of subradiant states in ordered arrays by using a spin-model formalism, wherein emission and re-absorption of photons by qubits is exactly accounted for. Our study reveals a number of interesting characteristics. In particular, as the number of qubits N increases, we show that the Liouvillian "gap" closes, i.e., there exists a smooth distribution of decay rates associated with subradiant states whose value approaches zero. Furthermore, we find that the most subradiant multi-excitation states exhibit "fermionic" correlations in that the excitations obey an effective Pauli exclusion principle. We propose a realistic experimental protocol to measure these exotic spatial properties, and finally investigate the correlations in the corresponding emitted field. Interestingly, similar behavior can also be found in other settings, such as in arrays of atoms in three-dimensional space [21], which suggests a certain degree of "universality" to the phenomenon of subradiance. Taken together, these results show that the physics of subradiance is itself a rich many-body problem.
Spin model -We consider N regularly-spaced two-level transmon qubits [22] with ground and excited states |g  1. (a) Schematic of planar transmon qubits capacitively coupled to a coplanar waveguide with rate Γ1D. Photonmediated interactions couple the qubits together with an amplitude determined by the single-qubit emission rate Γ1D into the waveguide, and a phase determined by the phase velocity of the transmission line and the distance between qubits. Collective frequency shifts (b) and decay rates (c) for qubits coupled through a waveguide with k1Dd/π = 0.2. Blue circles correspond to the results for a finite system with N = 30 qubits. Black dotted lines correspond to k = ±k1D. The frequency shift for the infinite chain is denoted by the solid line (see Appendix A). The inset in (c) shows the scaling of the decay rate with qubit number Γ/Γ1D ∼ N −3 for the 4 most subradiant states.
tromagnetic environment in the Markovian regime, one finds that emission of photons into the waveguide leads to cooperative emission and exchange-type interactions between the qubits [16,18,21,23,24]. The dynamics of the qubit density matrix ρ can be described by a master equation of the formρ = −(i/h)[H eff ρ − ρH † eff ] + m,n Γ m,n σ m ge ρσ n eg , where the effective (non-Hermitian) Hamiltonian reads [16,23,24] with J m,n = Γ 1D sin(k 1D |z m − z n |)/2 and Γ m,n = Γ 1D cos(k 1D |z m − z n |) denoting the coherent and dissipative interaction rates, respectively. Here, Γ 1D is the single qubit emission rate into the transmission line, k 1D = ω eg /v is the resonant wavevector, σ m αβ = |α m β m | denotes an operator acting on the internal states {α, β} ∈ {g, e} of qubit m at position z m = md, with d the inter-qubit distance. The photonic degrees of freedom can be recovered after solving the qubit dynamics [16,24]. In particular, the positivefrequency component of the left-and right-going field emitted by the qubits readsÊ is measured directly beyond the first (last) qubit, at position z L = d (z R = N d). Here,Ê in L/R denotes the (quantum mechanical) input field. For a given number of excitations, the effective Hamiltonian H eff defines a complex symmetric matrix that can be diagonalized to find collective qubit modes with corresponding complex eigenvalues defining their resonance frequencies (relative to ω eg ) and decay rates.
In the simple case of k 1D d = 2nπ [(2n + 1)π], with n an integer, the coherent qubit-qubit interactions J m,n vanish and the effective Hamiltonian is purely dissipa- k is coupled superradiantly to the waveguide at a rate N Γ 1D , while all other modes are dark, with decay rate Γ = 0. This realizes the ideal Dicke model of superradiance [25]. Within the setting of a 1D waveguide, it has also been shown that this configuration has interesting quantum optical functionality. For example, the qubits act as a nearly perfect mirror for near-resonant photons [16,26,27] and can generate arbitrary photon Fock states on demand [28]. Away from this spacing, the system becomes multimode [20]; below, we provide the first comprehensive look at these modal properties.
Single-excitation modes -Numerical diagonalization of H eff in the single-excitation sector gives N distinct eigenstates |ψ , where J ξ and Γ ξ represent the frequency shift and decay rate associated with |ψ (1) ξ . Here, |e n = σ n eg |g ⊗N corresponds to having atom n excited, with all other atoms being in the ground state. We obtain a broad distribution of decay rates defining superradiant (Γ ξ > Γ 1D ) and subradiant (Γ ξ < Γ 1D ) states. Ordering the eigenstates by increasing decay rates, i.e., from ξ = 1 for the most subradiant to ξ = N for the most radiant, we find that strongly subradiant modes are characterized by a decay rate Γ ξ Γ 1D that is suppressed with qubit number as Γ ξ /Γ 1D ∝ ξ 2 /N 3 , similar to the case of a 1D chain of atoms in free space [21]. Thus, in the thermodynamic limit, the spectrum of decay rates becomes smooth and the "gap" of minimum decay rate closes. For an infinite chain, the eigenstates of H eff take the form of Bloch spin waves |ψ with k being a quantized wavevector within the first Brillouin zone (|k| ≤ π/d). In such states, the single excitation is delocalized and coherently shared among all the qubits. For finite chains, the eigenstates are instead described in momentum space by a wavepacket with a narrow distribution of wavevectors around a dominant wavevector k, which can thus serve as an unambiguous label of states. In Fig. 1 (b) and (c), we show the distribution of frequency shifts J k and decay rates Γ k with k for N =30 qubits and k 1D d/π = 0.2. We find large decay rates and frequency shifts for eigenstates with wavevectors k close to the resonant wavevec-tors ±k 1D . Conversely, we obtain decay rate minima and small frequency shifts around kd = 0 and |k|d = π. For k 1D d > 0.5π (k 1D d < 0.5π), wavevectors k d = 0 form the global (local) and |k| d = π the local (global) decay rate minimum, respectively. The k-dependence can be understood by considering the infinite lattice limit, where the qubits and waveguide generally hybridize to form two lossless polariton bands (see Appendix A). The dispersion relation is plotted in Fig. 1 (b) as a solid line, and matches well with the frequency shifts obtained for a finite system. For a finite system, the polaritons with wavevector around k ∼ k 1D (k = 0, π/d) are most (least) impedance-matched at their boundaries to the dispersion relation of propagating photons in the bare waveguide, thus giving rise to super-radiant (sub-radiant) emission.
Multi-excitation modes -A quadratic bosonic Hamiltonian would enable us to easily find the multi-excitation eigenstates of H eff from the single-excitation sector results. Here, however, the spin nature prevents multiple excitations of the same qubit. Specifically, two-excitation states |ϕ with N 2 a normalization factor, are not eigenstates of the effective Hamiltonian (1). Moreover, in the case of an index ξ corresponding to a subradiant single-excitation mode, the initial decay rate of |ϕ (2) ξ is significantly greater than twice the single-excitation decay rate Γ ξ . This discrepancy can be explained by noting that the spatial profile of |ϕ i.e. the probability p m,n = | e m , e n |ϕ (2) ξ | 2 for qubits m and n to be excited, contains a sharp cut along the diagonal m = n (p m,m ≡ 0). In reciprocal space, this corresponds to a broad distribution of wavevector components, including radiant contributions responsible for an increased decay rate. From this qualitative discussion, one expects the excitations forming a multi-excitation subradiant eigenstate to be smoothly repelled from one another.
We numerically find the existence of two-excitation subradiant eigenstates |ψ (2) ξ , with a decay rate scaling as Γ (2) ξ /Γ 1D ∼ N −3 -as in the single-excitation sector -for the most subradiant eigenstates. These eigenstates reveal interesting properties in real and momentum space. One example is illustrated in the top of Fig. 2 (a), where we consider the most subradiant two-excitation wavefunction |ψ (2) ξ=1 = m<n c mn |e m , e n for k 1D d/π = 0.2 and N = 20 qubits, and plot both the probability amplitude |c mn | 2 in real space (left) and |c k1,k2 | 2 in reciprocal space (right). Here, c k1,k2 refers to the two-dimensional discrete Fourier transform of c m,n . In real space, one sees that the maximum in |c mn | 2 occurs for m ≈ 15, n ≈ 6, revealing a tendency for the excitations to both repel each other, and avoid the system boundaries where they can be radiated. At the same time, in momentum space, a peak occurs around k 1,2 d/π = ±1, coinciding with the dominant wavevectors kd/π ≈ ±1 of the most subradiant single-excitation states [ Fig. 1 (c)]. A natural twoexcitation wavefunction ansatz that realizes both the real-and momentum-space properties consists of taking and in reciprocal space |c k 1 ,k 2 | 2 of the wavefunction profile of the most subradiant two-excitation eigenstate for k1Dd/π = 0.2 (top) and k1Dd/π = 0.5 (bottom), for N = 20 qubits. cmn is only defined for m < n, but for visual appeal here we symmetrize the plot by taking cmn = cnm. Dotted dashed circles are a guide to the eye to highlight the positions of the maximum momentum components. (b) Fidelity between the exact two-excitation eigenstate of the Hamiltonian with quasi-momentum values (k1,k2) and the fermionized ansatz for N =50 qubits and k1Dd/π = 0.2.
an anti-symmetric combination of single-excitation eigenstates, which enforces a Pauli-like exclusion ("fermionization"). In particular, starting from the wavefunctions of the two most subradiant single-excitation eigenstates, we find that we can construct an accurate approximation of the most subradiant two-excitation eigenstate, |ψ n |e m e n , with N a normalization factor. For k 1D d mod π = 0 and k 1D d/π away from 0.5, the ξ = 1, 2 single-excitation states have dominant wavevectors (k 1 , k 2 ) near the global decay rate minimum, e.g., at k = π/d for k 1D d/π = 0.2. For k 1D d/π = 0.5, the fermionic ansatz also works well to describe the most subradiant two-excitation eigenstate (bottom of Fig. 2 (a)). In this case, it is built from the most subradiant single-excitation eigenstates k 1 = π/d and k 2 = 0 (degenerate in decay rate), and results in the checkerboard pattern seen in the plot. To more generally examine the accuracy of the ansatz, we take the two-dimensional Fourier transform of each two-excitation eigenstate, and unambiguously assign a label of quasimomentum indices (k 1 , k 2 ) to each state |ψ (2) (k1,k2) based upon where the Fourier transform is peaked. We then compute the overlap fidelity (k1,k2) | 2 between the exact state and the fermionic ansatz composed of the single-excitation eigenstates (k 1 , k 2 ). As illustrated in Fig. 2 (b), the ansatz works well when the two single-excitation states composing the eigenstate are strongly subradiant. In this case we find that the infidelity 1 − F scales with the qubit number as 1/N 2 (see Appendix B 1). In the thermodynamic limit N → ∞, we find that the decay rate of such subradiant "fermionized" eigenstates approaches the sum of the decay rates of the single-excitation states they are composed of (see Appendix B 2). The conclusions made about the subradiant decay rate scaling and their fermionic nature -exem-plified here for two-excitations -are found to extend to higher excitation numbers provided that the density of excitations is dilute: m ex N .  Eigenstate preparation and measure of fermionic correlations. It would be interesting to observe the fermionic nature of subradiant states in an experiment. We note that deterministically generating a single subradiant eigenstate with near unity efficiency is likely unrealistic, given that these states form a continuum for large N . Furthermore, for a classical coherent state input that drives the qubit array, the array tends to inherit a photon number uncertainty. It can be shown [16] (see also Appendix C) that adding a single ancilla qubit to the array, which can be individually addressed, enables a collective Fock state with well-defined wavevector k to be generated, by alternately creating an excitation in the ancilla and coherently transferring it to the array.
Fock states |ϕ (mex) k ∼ (S † k ) mex |g ⊗N generated in this manner, for low numbers of excitations m ex and a k-vector corresponding to the decay rate minimum, are found to have a significant overlap with the most subradiant eigenstates. Hence, these states are natural starting points for subradiant eigenstate distillation. For instance, when k 1D d = 0.7π, the N = 10 two-excitation state |ϕ (2) k=0 is found to have an overlap F 0.58 with the most subradiant eigenstate (with only a weak dependence on the actual qubit number N ). An evolution of the aforementioned state in time t decreases the probability ℘ (2) for having two excitations in the system due to emission processes as depicted in Fig. 3 (a). However, conditioned on finding two excitations in the system, the state approaches the most subradiant two-excitation eigenstate, e.g., F 1D . This preparation is illustrated in Fig. 3 (b), where we show the emergence of fermionic spatial correlations in time. For completeness, we study the effect of independent dephasing and decay for each qubit at rates γ d and Γ , respectively, on this dynamical preparation of fermionic correlations in Appendix C 2. We find that a clear anti-bunching structure can be observed provided that Γ , γ d ≤ 10 −2 Γ 1D for the N = 10 qubit chain.
Fermionic spatial correlations can be probed by using in-parallel readout of two resonators which are each The red lines represent the delay times for which the amplitude of the oscillations in T (2) 1D for an initial state |ϕ (2) k=0 (solid curve) and |ψ  dispersively coupled to their own qubit [29,30]. While finite readout time adds experimental difficulty to taking precise snapshots of spatial correlations in time, practical readout times of 100 ns should be sufficient to capture dynamics on timescales of 5Γ −1 1D ≈ 800 ns while maintaining Γ /Γ 1D , γ d /Γ 1D ≈ 10 −2 . This assumes uncorrelated relaxation and dephasing rates (Γ and γ d , respectively) at the level of Γ , γ d ≈ 2π × 10 kHz, well within current experimental capabilities.
Correlations in the emitted field -A natural question is how the qubit fermionic correlations are mapped into the emitted photons, i.e., what kind of photon correlations can be observed in the radiation of the most subradiant two-excitation eigenstate.
We first analyze what happens to the most subradiant eigenstate in the two-excitation sector, |ψ (2) ξ=1 , once a photon is emitted and detected, for example, on the left side of the chain. We find that the new conditional state after detecting a photon, |ψ c ∼ E + L (t)|ψ is predominantly formed by a superposition of the two single-excitation states |ψ (2) ξ=1 it is composed of, i.e., |ψ c α 1 |ψ ξ=2 . More precisely, the projection of the conditional wavefunction onto any state besides the two most subradiant, ε = 1−|α 1 | 2 −|α 2 | 2 , scales as ε ∼ 1/N 2 for most lattice constants k 1D d = 0.5π.
After one photon is emitted at time t, the relative intensity of emission after a delay time τ , normalized by the intensity at time t, is given by the two-photon correlation function Prior experimental [31][32][33] and theoretical [34] work has demonstrated that such correlation functions can be measured in the microwave domain by amplifying the outgoing photon field and performing correlation measurements between two linear detectors. Figure 4 shows T (2) (t, τ ) for a chain of 10 qubits with lattice constant k 1D d = 0.7π. At t = 0, the qubits are prepared in the symmetrized state |ϕ (2) k=0 = N 2 (S † k=0 ) 2 |g ⊗N as discussed in the previous section. For short evolution times t, the overlap between the time-evolved symmetrized state and the most subradiant second-excitation eigenstate is not large, and the former still has significant superradiant components, leading to a sharp decay of T (2) (t, τ ) with τ . For longer times t, the radiant components of the initial state disappear, as illustrated on the upper plot of Fig. 4, and one observes beatings in T (2) (t, τ ) as a function of τ coming from the interference in emission of the two single-excitation subradiant components (see the right part of Fig. 4). The oscillation period is determined by the difference in frequencies of the two most subradiant single excitation eigenstates, J ξ=1 and J ξ=2 , respectively. In particular, the maxima in T (2) (t, τ ) occur at delay times τ max = nπ/|J ξ=1 − J ξ=2 |, with n an odd integer. While the observation of these oscillations is limited by the dephasing and spontaneous emission rates into other channels, they offer a glimpse into the decay structure of many body subradiant states.
In this Letter, we provided a comprehensive study of the subradiant properties of artificial atoms in waveguide QED. We have shown that this system represents an open quantum critical system with a closing of the Liouvillian gap in the thermodynamic limit, and a decay rate suppression scaling similar to that of atomic chains in free space. We have also shown that multi-excitation subradiant states exhibit "fermionic" spatial correlations, which can be probed in realistic experiments. More broadly, the investigation of subradiance overall has become of active interest in recent years [21,[35][36][37][38][39]. Waveguide QED systems appear to be an attractive platform to exploring this type of many-body dissipative dynamics. Exploring the degree of universality of these dynamics across different platforms remains as an interesting goal of future work.

APPENDICES Appendix A: Polaritons in the infinite lattice limit
Here, we describe the hybrid light-matter excitations of the system in the infinite lattice limit. This analysis also allows us to characterize the validity of the spin-model used in the main text, after the photonic degrees of freedom have been integrated out. For an infinite system, the total Hamiltonian describing both the qubits and photonic degrees of freedom is given by Here, S † k = 1/ √ N N j=1 e ikzj σ + j creates a collective spin excitation with k a quantized wavevector in the first Brillouin zone, and a † k is the creation operator of a propagating excitation with wave-vector k and frequency ω k = v|k| in the transmission line. The third term of Eq. (A1) describes the interaction between the qubits and the electromagnetic field, where the parameter g k quantifies the strength of the interaction. We take a light-matter coupling of the form [23] , where ω f > ω eg is a high-frequency cutoff and θ(.) is the Heaviside step function.
For each wavevector, and within the single-excitation sector, the Hamiltonian represents a 2x2 matrix that can be diagonalized to yield frequencies Ω ± k , as shown in Fig. 5. Physically, the two distinct solutions correspond to a qubit branch and a waveguide branch, with significant hybridization of the two around their intersection at k = ±k 1D . For a finite system, this implies that a collective excitation of qubits with wavevector close to ±k 1D efficiently radiates into the waveguide, as confirmed in Fig. 1 (c) of the main text. In the regions where |J ± k |/ω eg 1, with J ± k = Ω ± k − ω eg the frequency shift, we recover a good agreement with the expression obtained from the direct Bloch diagonalization of the effective spin-model Hamiltonian [Eq. (1) of the main text], which predicts J k ∼ Γ 1D [cot(k + k 1D )d/2 + cot(k 1D − k)d/2]/4 for k = ±k 1D and Γ k ∼ N Γ 1D δ k,±k1D /2, with Γ 1D = 2πg 2 ω eg (see Fig. 5).
Appendix B: Multi-excitation subradiant states 1. Fidelity scaling of the fermionic ansatz As illustrated in Fig. 2 (b) in the main text, two-excitation eigenstates |ψ (2) ξ around regions of decay rate minima can often be well-approximated by an ansatz |ψ (F ) ξ that constructs fermionic combinations of single-excitation states.
Here, we analyze the scaling of the infidelity of such a fermionic ansatz 1 − F with the qubit number N . The two-excitation eigenstates can be classified based on the associated wavevectors of their underlying single-excitation states. In particular, for the subradiant states considered here, these wavevectors take on values around the two decay rate minima at k = π/d and k = 0 shown in Fig. 1 (c) of the main text. Two types of fermionic combinations can be distinguished: (i) the ones for which both single-excitation states can be associated with wavevectors corresponding to the same minimum (i.e. with both wavevectors around the decay rate minimum of either k = 0 or k = π/d) and (ii) the ones that are combinations of both minima (i.e. one underlying single-excitation state associated with k = 0 and one with k = π/d). The former behavior can be seen in the most subradiant two-excitation eigenstate of k 1D d = 0.2π and the latter for k 1D d = 0.5π; both of these examples are pictured in Fig. 2 (a). In case (i), the wavevectors are localized within the same region in momentum space, and we find a scaling of 1 − F This difference in scaling is illustrated in Fig. 6, where in subfigure (a) the infidelity scaling for the most subradiant eigenstate is plotted for various k 1D d. In subfigure (b), the infidelity is plotted for the three most subradiant twoexcitation eigenstates of k 1D d = 0.5π. In these examples, only the most subradiant two-excitation eigenstate of k 1D d = 0.5π is of the type (ii) and thus shows a deviating infidelity scaling ∼ N −1 as opposed to the ∼ N −2 scaling of the remaining states.   ξ=1 (blue, scaling: s = 1) and the second-and third-most subradiant two-excitation eigenstates |ψ (2) ξ=2 , |ψ (2) ξ=3 (red, scaling: s = 2).

Scaling of decay rates
The most subradiant m ex -excitation eigenstates of the effective Hamiltonian are well described by the fermionic ansatz at low excitation densities (m ex N ) and can thus be written |ψ We find that r2 and r3 evolve as N −s with s 1.

Appendix C: Excitation transfer and subradiant eigenstate distillation
In the main text, we used that a collective Fock state of definite wavevector can be prepared in the qubit chain configuration. The preparation in such a state is achieved by adding an ancilla qubit, which can be individually addressed and enables the transfer of excitations to the chain qubits. Subsequent to the preparation step, this ancilla can be effectively decoupled from the dynamics of the rest of the chain by shifting its frequency far out of resonance. Here, we describe the preparation procedure and the distillation of subradiant eigenstates in more detail. Section C 1 introduces both the ancilla-chain configuration and the transfer protocol. The fidelities and success probabilities achieved by such a state preparation, in the presence of noise and imperfections, are discussed in Section C 2. A less universal preparation procedure, that however does not require an adjustment of phases, is outlined in Section C 3.   The transfer of excitations can be achieved in the so-called cavity configuration of waveguide QED [16] illustrated in Fig. 8 (a), where the N chain qubits (depicted in red) form "mirrors" and an ancilla "cavity" qubit (blue) is introduced at the midpoint. The chain qubits are equally spaced at a distance d and the ancilla qubit is separated by d c from the nearest chain qubits. Specifically, k 1D d = π and k 1D d c = π/2, where k 1D = ω eg /v is the wavevector of the qubit transition of frequency ω eg within the waveguide of group velocity v. The effective spin Hamiltonian for such a configuration, directly following from Eqn. (1) in the main text, is given by [16]

Excitation transfer setup and protocol
with S † mirr/rad = N/2 n=1 (−1) n (σ n eg ± σ −n eg )/ √ N the collective transfer and decay operators (respectively) of the chain of mirror qubits. Here, n and −n enumerate the mirror qubits to the right and left of the ancilla qubit a, respectively. The first term in Eq. (C1) allows for a coherent exchange of excitations between the ancilla a and the chain qubits, |e a ⊗ (S † mirr ) mex−1 |g ⊗N ⇔ |g a ⊗ (S † mirr ) mex |g ⊗N . The resulting chain state exhibits zero (for m ex = 1, ||S rad S † mirr |g ⊗N || 2 = 0) or low decay (∼ Γ 1D /N for 1 < m ex N ). The preparation protocol, ideally resulting in a state |ϕ (mex) mirr = N (S † mirr ) mex |g ⊗N , works as follows: Starting from all qubits in their ground state, a single-excitation Fock state in the mirrors can be prepared by applying a fast π-pulse to the ancilla qubit, |g a → |e a , and subsequently waiting for a time t for that m ex -th (here m ex = 1) excitation to be transfered to the mirrors. Higher number Fock states can be prepared by repeating the process. Eliminating the ancilla qubit subsequent to the transfer, e.g., by detuning its frequency, reduces the system to a qubit chain periodically spaced by distance d. Note that the resulting state, if ideally prepared, automatically consists of a collective state where each qubit is excited with the same amplitude, and with a well-defined phase.
Two further steps are needed to transform that Fock state to a state of definite wavevector on a qubit lattice of selected periodicity: First, a transformation to the state |ϕ (mex) k ∼ (S † k ) mex |g ⊗N of wavevector k, which can be realized by applying fast local phases to the qubits [29] (here S † k = n e ikdn σ n eg ). And second, a dynamical modification of k 1D d (in particular, via the resonant wavevector k 1D itself), which can be accomplished by flux tuning the qubit transition frequency [22].
As discussed in the main text, a state |ϕ with k a wavevector corresponding to the global decay minimum (e.g., k = 0 for k 1D d > 0.5π), is characterized by a significant overlap with the most subradiant m-excitation eigenstate. In particular, the state conditioned on m ex -excitations converges to the latter state in time, at the expense of a decreasing m ex -excitation probability.

Preparation fidelity and probability
In practice, both intrinsic decay through the waveguide in Eq. (C1), and excitation losses into channels other than the waveguide and dephasing, affect the probabilities and fidelities of the excitation transfer and eigenstate convergence process. We model the non-waveguide decay (of rate Γ ) and dephasing (of rate γ d ) as uncorrelated and Markovian, that is where {·, ·} denotes the anticommutator. We start by analyzing the transfer of m ex excitations, targeted towards ideally preparing the state |ϕ (mex) mirr , and the impact of imperfection mechanisms on its transfer probability and fidelity. For the simulations, the transfer times for each excitation t (mex) π ≈ π/(Γ 1D √ N m ex ) have been optimized to maximize the fidelity with the state |ψ (mex) ∼ |g a ⊗ (S † mirr ) mex |g ⊗N , where the ancilla qubit a is in the ground state and the chain qubits are in the target state. Moreover, we model the π-pulse to excite the ancilla qubit as an ideal gate that takes a negligible amount of time to perform. Incorporating loss and dephasing processes, the resulting chain state after the transfer, tracing out the ancilla qubit, is characterized by the reduced density matrixρ  Fig. 8 (b) for N = 10 qubits as a function of Γ and γ d . Note that even for Γ = γ d = 0, the two-excitation transfer probability (fidelity) is limited to ℘ (2) tra 0.45 (F (2) 0.99) by the loss processes through the waveguide. However, as the dominant loss mechanism stems from the emission of the ancilla qubit into the waveguide, decreasing the ancillawaveguide coupling compared to the waveguide coupling of the chain qubits allows for an increased fidelity and transfer probability at the expense of a longer transfer time. The fidelity only depends on the dephasing rate γ d and is independent of Γ .
The stateρ (mex) mirr resulting from the excitation transfer ideally excites the qubits with the well-defined phases illustrated in Fig. 8 (a) -the qubits in the left and right "mirrors" each have alternating phases corresponding to a spin wavevector of kd = π, while a "phase slip" occurs between the left and right mirrors as the qubits closest to the ancilla have the same phase. We want to subsequently convert this state to a state of well-defined k, here assumed to be k = 0, in the attempt to ideally prepare the m-excitation state |ϕ (mex) k=0 . We do so by the phase adjustment operation:ρ where S π = n exp(−i[1 − (−1) n ]π(σ n ee + σ −n ee )/2), with n and −n enumerating the left and right mirror qubits, respectively. As for the π-pulse operation, we assumed this gate to be instantaneous and of unit fidelity.
Subsequently, an evolution of the stateρ (mex) k=0 in time, ideally makes its m-excitation component converge to the most subradiant eigenstate (for k 1D d > 0.5π). This is illustrated in Fig. 9 (a) for m ex = 2, where both ℘ (2) (t) (the probability of finding two excitations in the system as a function of time t) and F (2) ξ=1 (t) = ψ (2) ξ=1 |ρ(t)|ψ (2) ξ=1 /℘ (2) (t) (the overlap fidelity with the most subradiant eigenstate as a function of time t) are shown for selected values of loss Γ and dephasing γ d . We find that the probability decay, for realistic parameters of γ d , Γ ∼ 10 −1 − 10 −3 Γ 1D , depends only on the sum Γ + γ d . This suggests that dephasing essentially destroys subradiance and thus results in fast losses. As in the transfer process, only dephasing γ d affects and limits the fidelity, which increases for times t < γ −1 d . Combining both the preparation and convergence state, the maximum preparation fidelity of the most subradiant two-excitation state in the presence of loss and dephasing is illustrated in Fig. 9 (b). Fidelities shown correspond to the maximum fidelity in convergence time t, conditioned on a probability of at least ℘ (2) (t) ≥ 0.2 to have two excitations in the system. A clear anti-bunching structure can be observed down to fidelities of around 75% (marked by the dashed line in Fig. 9 (b)), which limits Γ , γ d < ∼ 10 −2 Γ 1D for the N = 10 qubit chain.  k=0 , the state after an imperfect two-excitation population transfer and subsequent phase adjustment to a definite wavevector k = 0, for N = 10 qubits and k1D d = 0.7π. Left: Two-excitation probability ℘ (2) normalized by the initial t = 0 probability ℘ (2) tra as given in Fig. 8 (b). Right: overlap fidelity F (2) ξ=1 with the most subradiant eigenstate in time. The parameters (Γ ,γ d ) and (γ d ) in units of the waveguide coupling Γ1D, for the population and fidelity, respectively, are indicated in the figure. The fidelity is independent of Γ . Blue dashed lines follow for an initial t = 0 plane-wave state |ϕ (2) k=0 and Γ = γ d = 0. (b) Maximum fidelity maxt(F (2) ξ=1 (t)) for preparing the most subradiant two-excitation eigenstate in the presence of additional loss and dephasing at rates Γ and γ d , respectively. Here, N = 10 qubits and k1D d = 0.7π. The maximization over the convergence time t is conditioned on a probability of having two excitations in the system ℘ (2) (t) ≥ 0.2. The dashed line marks a fidelity of 75%.

Phase adjustment free k-state preparation
Above, we noted that states resulting from a transfer of excitations in the mirror configuration differ from plane wave excitations |ϕ (mex) k of fixed wavevector k by a π-phase discontinuity between left and right mirror qubits. Thus, a phase adjustment operation S π has been introduced for a transformation to the latter states. Here, we show that for k = π/d and specific values of k 1D d, where k 1D = ω eg /v is the wavevector of the qubit transition frequency ω eg within the waveguide of group velocity v and d the inter-qubit distance, such an adjustment step can be omitted. In the 'cavity configuration' depicted in Fig. 8 (a), the distance d c between the ancilla and chain qubits can be chosen as k 1D d c = (2η + 1) π/2 with η ∈ N 0 ; for the chain qubits k 1D d = π. The transfer Hamiltonian H C [Eq. (C1)], derived for η = 0, remains invariant under an arbitrary choice of η up to global phases in the collective operators. For η = 0, the chain becomes equally spaced once the ancilla is eliminated. For η = 0, and for arbitrary k 1D as obtained after adjusting the transition frequency, the phase factor (2η + 1) k 1D d between the two most central qubits can differ in a non-trivial way from k 1D d. This latter property can be exploited to obtain a configuration of the same spectral properties (e.g. eigenstate decay rates) as for a regularly spaced chain of lattice constant k 1D d, and where the eigenstates are identical to those of a regular chain apart from a π-phase discontinuity across the center. The Hamiltonian for that generalized setup, absent the ancilla qubit, follows as where the numbering of qubits is split into the subset of qubits to the left (N L ) and right (N R ) of the original ancilla qubit. The above Hamiltonian differs from that of a regularly spaced chain by the additional phase factor exp(i2ηk 1D d) between left and right qubits. However, for ϑ ≡ 2ηk 1D d = l π, where l ∈ N 0 , this Hamiltonian can be formally mapped onto the former one by re-defining the spin operators as σ n eg = e i ϑ σ n eg for n ∈ N R σ n eg for n ∈ N L . (C4) For ϑ ≡ 2ηk 1D d = (2l + 1)π (with l ∈ N 0 ) the left and right operators in Eqs. (C4) differ by a π-phase, and so do the eigenstates as compared to an equally spaced chain. Thus, the mirror excitation k = π/d with phase discontinuity in this setup takes the same role as that without such a discontinuity in a regularly spaced chain. The phase-adjustment step can therefore be omitted. Apart from the phase discontinuity in the eigenstates, all other properties are inherited from the regularly spaced chain. For η = 1 such a regime is reached for k 1D d = π/2, for η = 2 the distance is k 1D d = π/4 (and k 1D d = 3π/4). The phase adjustment free preparation comes with a restriction to specific k 1D d by fabrication. Moreover, the system must be probed in the experimentally more challenging regime k 1D d ≤ 0.5π where k = π/d corresponds to a global decay rate minimum, i.e. the transition frequency has to be changed by more than 50% subsequent to the transfer process.