Emergence of stationary many-body entanglement in driven-dissipative Rydberg lattice gases

Non-equilibrium quantum dynamics represents an emerging paradigm for condensed matter physics, quantum information science, and statistical mechanics. Strongly interacting Rydberg atoms offer an attractive platform to study driven-dissipative dynamics of quantum spin models with long-range order. Here, we explore the conditions under which stationary many-body entanglement persists with near-unit fidelity and high scalability. In our approach, coherent many-body dynamics is driven by Rydberg-mediated laser transitions, while atoms at the lattice boundary reduce the entropy of the many-body state. Surprisingly, the many-body entanglement is established by continuously evolving a locally dissipative Rydberg system towards the steady-state, as with optical pumping. We characterize the dynamics of multipartite entanglement in a 1D lattice by way of quantum uncertainty relations, and demonstrate the long-range behavior of the stationary entanglement with finite-size scaling, reaching"hectapartite"entanglement under experimental conditions. Our work opens a route towards dissipative preparation of many-body entanglement with unprecedented scaling behavior.

Non-equilibrium quantum dynamics represents an emerging paradigm for condensed matter physics and statistical mechanics. Strongly interacting Rydberg atoms offer an attractive platform to study the drivendissipative dynamics of quantum spin models with long-range order. We explore conditions under which stationary many-body entangled phases persist with near-unit fidelity in the presence of optically-driven dipolar interactions and site-localized dissipations, thereby providing a scalable and robust scheme for the production of many-body entangled steady-states. We characterize the dynamics of multipartite entanglement in a 1D lattice, and demonstrate a long-range behavior for the stationary entanglement with finite-size scaling. Quantum control of open many-body systems has become a major theme in the quest to explore new physics at the interface between condensed matter physics, quantum information science, and statistical mechanics [1][2][3][4][5]. The ability to control and tailor the many-body interactions and their dissipative jump processes has been identified as a powerful resource for the preparation of stationary entangled states [6] and the investigations of noise-driven quantum phase transitions [2]. Indeed, quantum reservoir engineering provides the frameworks for dissipative quantum computation [3,4] and communication [7] with built-in fault-tolerance, and for dissipative quantum simulation [2], with tremendous experimental progresses [8]. Furthermore, continuous evolution of open quantum systems offers intriguing prospectives to the underlying relationship between entanglement and nanoscale quantum thermodynamics [5].
Laser-driven Rydberg atoms offer unique possibilities for creating and manipulating open quantum systemsρ (r) mb of dipolar interacting spin models [9]. By exciting atoms to high-lying Rydberg states, strong and long-range interactions between the Rydberg atoms can be exploited to induce spinspin interactions, whereas atoms comprising the many-body state can couple to their radiative environments by spontaneous emission [7]. The competition between the coherent and incoherent dynamics can drive the system to bipartite entangled states for two atoms [11] and novel states of matter for a mesoscopic number of atoms, exhibiting topological order, glassiness, and crystallization dynamics [12][13][14][15][16][17][18]. Dissipative light propagation in Rydberg medium can show highly nonclassical optical behaviours [19]. Remarkably, the basic primitives behind such a principle have been demonstrated in the laboratory by several groups [20].
Another frontier is the characterization of entanglement in the many-body stateρ (r) mb (t) under dynamical transitions [1,21]. For interacting spin systems, spin wavesŜ m are the quasiparticle excitations describing the "collective" state of the atomic mode. Entanglement in such a system is then defined by the quantum correlations among the atomic modes [1]. Hence, verification protocols for number-state entanglement can be extended to properly extract the many-body entanglement stored in the collective mode of these quantum spin systems [22,23]. For instance, uncertainty relations, as defined in Refs. [23], can be applied to access the genuine multipartite entanglement in the subspaces ofρ (r) mb containing a fixed number of total excitations n = Ŝ † mŜm [6]. In this Letter, we explore such non-equilibrium many-body entangled states persisting with high fidelity in the stationary limit for dipolar interacting atoms in a lattice. As illustrated in Fig. 1, the main idea is to globally pump the atoms with a driving field Ω to its Rydberg state |r and to induce site-selective decoherence on "reservoir" atoms. While the dissipative jump operators are strictly local, the resulting dynamics amounts to an optical pumping for the atomic sample to one of the entangled eigenstates of the many-body HamiltonianĤ (n) xy in the single-excitation subspace n = 1. We apply our method to generate bipartite entanglement for system-size N = 4, and then extend our work for N = 6 to investigate the drivendissipative dynamics of the many-body entanglement via uncertainty relations [23]. We find that quadripartite W -state can persist indefinitely with fidelity F ≥ 0.99 for N = 6, and that entanglement depth k shows favorable scaling relative to its system size, reaching "hectapartite" entanglement (k = 100) for N = 126 atoms.
We consider the many-body states of N atoms configured in a lattice [See Fig. 1(a)], irradiated by a uniform driving field that couples the atomic ground state |g to the highly excited Rydberg state |r with Rabi frequency Ω and detuning δ. A pair of atoms i, j in the Rydberg state at lattice sites x i , x j couple each other via the potential ∆ with power-law scaling, for which we take p = 6 for the van der Waals regime of blockade shifts [7]. In a frame rotating with the laser frequency, the system Hamiltonian is given bŷ mb }) for the atomic coupling to their local radiative reservoirs. Here, we can set arbitrarily the decay rates Γ i by dressing the atoms i with a rapidly decaying excited state |e , as shown in the inset of Fig. 1 (See Ref. [25]).
By working in the off-resonant limit |δ − ∆ Fig. 1(b)], we decompose the system Hamiltonian into an effective spin Hamiltonian by adiabatic elimination [25], where w d = Γ 2 /4 + 2Ω 2 is the power-broadened linewidth andδ (i) is the effective detuning respect to the blockade shifted level for site i. The (anti)ferromagnetic XY Hamil-tonianĤ xy with C p > 0 (C p < 0) delocalizes the Rydberg excitations between sites i, j at a rate J ij constrained by the excitation manifold n, whereas the "squeezing" termĤ 2 provides the two-photon coupling between the manifolds l, l + 2.
The expressions forδ (i) , J ij and S ij are given in Ref. [25]. The dissipative many-body entanglement forρ (r) mb is prepared as follows. (1) A set of lasers with detunings δ l is used to resonantly drive transitions between the most shifted states of the manifolds l, l with l ≤ l max except for a particular subspace n < l max , where n is the subspace of the target eigenstate [25]. (2) The combination of squeezing termĤ 2 (l → l ) and dissipation Lρ    {p} is the basis state of the subspace with n Rydberg excitations on the sites p i ∈ {1, · · · , N } (e.g., |r (2) ij = |g 1 · · · r i · · · r j · · · g N ) and the bipartition A (B) is a cut of Hilbert space for entangled atoms ("reservoir" atoms). Hence, the overall atomic dynamics produces a stationary entangled eigenstate |ψ eig of the many-body HamiltonianĤ (n) xy for t → ∞ by way of the dipole-mediated optical pumping of "system" atoms A with "reservoir" atoms B [28]. In the quantum-jump picture, as the environment projects B to |g · · · g B , atoms A collapse to |ψ eig in an interaction-free manner due to the correlations between A and B viaĤ (n) xy . As illustrated in Fig. 1(b), we consider the case of enhanced decay rates Γ 1,N = Γ Γ r for "edge" atoms (number of defect sites N d = 2) and pumping to single-excitation subspace (n = 1). In order to achieve this condition, we selectively dress atoms {1, N } at each end of the lattice boundary with fields Ω d to induce Γ [25], which would otherwise have a natural decay rate Γ r as for atoms {2, · · · , N − 1}. For the 1D staggered triangular lattice in Fig. 1(a), ξ is a parameter controlling the relative interaction strength between nearest (a 0 ) and next-nearest (a 1 ) neighbor sites via a 1 = ξa 0 .
As we tune the driving laser to δ = ∆ (i,i+1) p /2, we set the target subspace n = 1 by driving excitations between manifolds l = 0 and l = 2. At the early stage of Liouvillian dynamics (0 ≤ t 1 ≤ 1/Γ), atoms in |G = |g · · · g are rapidly pumped to the n = 1 subspace. The Rydberg excitation in the subspace then evolve underĤ sipatively pumped to a W -like entangled state |ψ eig by way of coherent population trapping for an effective Λ-level consisting of a "radiative" state |r (1) 2 [25]. In the following, we first perform a numerical analysis of the relaxation behavior of the Rydberg gas to a stationary bipartite entanglement for system-size N = 4 and N d = 2. Fig. 2(a) displays the contour map of entanglement fidelity ss ]|ψ 2 for the stationary stateρ as a function of interaction parameter ξ and distance a 0 (in units of blockade radius d B = 6 C 6 /w d ). The profile of fidelity along a 0 depicts the requirement of Rydberg-blockade regime a 0 < d B to provide sufficient nonlinearity in l [ Fig. 1(b)] to selectively drive transitions between |G and |r (2) i,i+1 , and to adiabatically eliminate subspace l = 0, 2 [25]. Atoms in the strongly interacting region 0.2 ≤ a 0 /d B ≤ 0.5 are efficiently pumped to the single-excitation subspace. The interaction parameter ξ is tuned to numerically maximize the steady-state entanglement fidelity up to F 2 = 0.9982 at ξ 1 = 6 √ 3 and ξ 2 = 0.36 for a 0 /d B = 0.26. To validate our entanglement pumping scheme, we also show the dissipative dynamics of bipartite entanglement by way of concurrence C [21] for ξ 1 in the inset to Fig. 2(b). The atomic sample A is driven to a maximally entangled steady-state |ψ 2 (F 2 = 0.9965) within tΓ = 200 (i.e., t ∼ 0.02/Γ r ). Now, let us treat the case of many-body entanglement with N = 6 atoms in the 1D lattice, as an example of multipartite system N 2. With same parameter set Ω and Γ, we simulate the dissipative dynamics of entanglement fi- xy by way of quantumtrajectory method [See Fig. 3(a)]. Here, we have optimized the steady-state fidelity max(F 4 ) = 0.9912 for the parameters {ξ, a 0 /d B } = {1.1996, 0.285}, thereby setting the eigenstate ofĤ (1) xy to the symmetric quadripartite W -state |ψ 4 . The dissipative transitions of many-body entanglement is detected by the quantum uncertainty relations [6,23], which serves as the collective entanglement witness {∆(t), y c (t)} for the low-lying excitations n ofρ detects the amount of higher-order spin-waves (e.g., represents the minimum uncertainty for (k − 1)-partite entangled states for a given y c . Violation of the uncertainty bound ∆(ρ , then, signals the presence of genuine kpartite entanglement stored inρ . Experimentally, {∆, y c } can be measured by the collective transverse spin fluctuations and by the excitation statistics [25]. By applying the uncertainty witness {∆, y c }, we observe that atoms initially in ground state are dissipatively driven to the quadripartite entangled W state by sequentially crossing the boundaries ∆ Fig. 3(b). The dissipative transitions of many-body entanglement are indicated by black, purple, green, and red lines of Fig. 3(b) for the average trajectoryρ (r) mb (t) of the Monte-Carlo wave-function. In particular, for pumping time t ∼ 100/Γ, the many-body system can exhibit a full quadripartite entanglement even with a moderate system-size N = 6 (N d = 2). In the stationary limit, the entanglement parameters reaches the minimum values {∆, y c }| ss → {1.5 × 10 −2 , 2 × 10 −4 } as the system is pumped to the desired eigenstate |ψ 4 . Next, we move on to the question of finite-size scaling behavior of the many-body entanglement. Here, we define the entanglement depth k in accord with the concept of kproducibility for qubits [1,21], thereby identifying the minimal depth for genuine k m -partite entanglement to produce the purported state Tr B [ρ 2), to tripartite entanglement (green, t3Γ 65), and to stationary quadripartite entanglement (red, t4Γ 100).
nalize the many-body HamiltonianĤ (1) xy for ξ 1 truncated up to next-nearest-neighbor interactions, and characterize the resulting many-body entanglement depth k stored in the stationary eigenstate |ψ eig (y c → 0) up to N s → 126. Fig.  4 captures our result of {∆, y c } for the eigenstates |ψ eig = |W k A ⊗ |g · · · g B , where |W k A is the k-partite symmetric W state. Due to the nonlinear sensitivity of our witness for some region k, we characterize the scaling of the minimal entanglement depth k m ≤ k [25]. The shaded area represents the physical region, whereby k m -partite entanglement could be defined for a given N s , and the dashed lines are the uncertainty bounds for 0 ∼ 100-partite entanglement (with 20partite increments). Remarkably, as we increase the system size N s , we observe a favorable scaling behavior of entanglement depth k m up to genuine "hectapartite" entanglement (k m = 100) for N s = 126.
We believe that our entanglement pumping scheme is feasible with realistic parameters. For example, by exciting 85 Rb atoms to Rydberg state |r = |62D 3/2 , quadripartite entan- ss is probed with quantum uncertainty witness ∆ for yc → 0 by way of direct diagonalization ofĤ (1) xy as a function of system-size Ns. We obtain stationary eigenstatesρ gled states could be prepared within t 4 = 3 µs. The limit for any driven-dissipative approach with Rydberg lattice gases will be the photoionization rate t ph ∼ 40 µs [25]. Since the pumping time to reach |ψ eig scales with N due to the quantum walk in n = 1 subspace, our method can be applied to generate stationary 30-partite entanglement with t p < t ph .
For further improvements, one could explore emergent atom-field systems embedded in photonic crystals (Ph.C). Here, dispersive optical interactions near band edges can induce dipole-dipole oscillationsĤ xy and effective "Rydberg" blockadesĤ 2 with tailored scaling ∆ (ij) p ∼ 1/x p ij between low-lying excited atoms in the band gap by exchanging virtual photons [29]. Decay rates Γ i = Γ can be controlled by the density of states of the Ph.C [30].
In addition, the coherent delocalization dynamicsĤ (n) xy in the restricted decoherence-free subspace n can be extended to examine non-equilibrium dynamics of random walks and to simulate diverse sparse Hamiltonian matricesĤ (n) xy for quantum algorithms [31]. The dissipative many-body pumping scheme proposed here is the starting point of a long-term research program of exploring dynamical behaviors of entanglement in open many-body systems and its underlying relation to quantum criticality out of equilibrium [1].
In conclusion, we have examined the condition under which driven-dissipative dynamics displays a rich family of manybody entangled states. The purported stationary many-body entanglement shows a favorable long-range scaling behavior k m relative to its system size N s up to k m = 100 for N s = 126. Our work, thereby, opens the door towards a new class of open quantum system simulator with well-controlled coherent and dissipative dynamics suited to the study of quantum many-body phenomena out of equilibrium, monitored by information-based quantities [1][2][3].
Funding is provided by the KIST Institutional Programs, and, in part, by the Canada Foundation for Innovation, the Ontario Ministry of Research & Innovation and Industry Canada.
As discussed in Ref. [1], for defect sites i ∈ B, atoms initially in the Rydberg state |r with decay rate Γ r radiatively couple to a highly decohering state |e with decay rate Γ e Γ r so that atoms in bipartition B can behave as an effective "reservoir" channel for the "system" atoms in partition A. In this section, we discuss how we could manipulate the spontaneous emission rate Γ i of the Rydberg state |r for the "defect" atoms.
As illustrated in Fig. S1(a), we consider a ladder-type energy level diagram, where |r is dressed with |e by auxilary field Ω d . In the rotating-wave frame of the dressing laser Ω d , the Hamiltonian is given bŷ The resulting optical Bloch equations are, then, where γ e,r = Γ e,r /2 and ∆ d is the detuning for the dressing field Ω d relative to the transition |e ↔ |r . In writing Eqs. S2-S3, we have neglected the Langevin noise forcesF µν and assumed c-number counterparts forσ µν . Hence, we find that the atomic coherence σ gr e −γrt and∆ = iγ e + ∆ d . Fig. S1(b) shows the dynamics of Rydberg population σ  [1] with Γ e = 10 4 Γ r . The black solid (dashed) line is the atomic dynamics for Ω d = 10Γ r (Ω d ∈ {10 2 Γ r , · · · , 9×10 2 Γ r } with 10 2 Γ r increments). The red line is the result of atomic decay Γ 10 3 Γ r with Ω d = 10 3 Γ r . As we increase Ω d → Γ e , we find that the effective decay rate for the defect atoms scales with Γ ∼ |Ω d | 2 /Γ e up to Ω d ∼ 0.1Γ e .

EFFECTIVE SPIN HAMILTONIAN
In the off-resonant limit |δ − ∆ (ij) p | w d , we obtain the effective HamiltonianĤ eff (Eq. 2 in Ref. [1]) by truncating the perturbative expansion to the second order and by timeaveraging highly oscillating terms [2], with the interaction Hamiltonian given bŷ where In particular, we usê withσ (i) µν = |µ i ν| for µ, ν ∈ {g, r} and blockade shift ∆ For the case of δ = ∆ (i.i+1) p /2, we obtain the following effective spin Hamiltonian (Eq. 2 in Ref. [1]) The coefficients for the n = 1 subspace to Eqs. S10a-S10c areδ with the shape function f ij = 1 − 2∆ and Kronecker delta δ kd . Now, let us discuss the possibility of optically pumping the systemρ (r) mb to the n-excitation subspace with n ≥ 2. As discussed in Ref. [1], this is achieved by a set of lasers resonantly driving the transitions l → l + 2 and n − 2 → n + 1, and by Lρ For the case of 1D spin chain with lattice constant a 0 , the most shifted state containing l excitations is obtained by degenerate Rydberg configurations with l-nearest neighbor excitations (e.g., |l = |r 1 , · · · , r l , g l+1 , · · · , g N ). The Rydberg spectrum is then given by where β (k) p = |k| −p is the power-law scaling parameter, and δ s (l − 1) is a step function. The transition energy for l → l + 2 is then so that the anharmonicity for the higher order l ≥ 2 is FIG. S3: The power law scaling behavior of spin-spin coupling strength Jij with 1D staggered triangular lattices for ξ1 = 6 √ 3 (blue) and ξ2 0.36 (red). For ξ1, the spatial range of Jij depicts a monotonic power-law decay, whereas jig-jag oscillatory pattern exists for ξ2.
As shown in Fig. S2, for a given target subspace n, we terminate the two-photon excitations to l max = n + 1. All subspaces with l ∈ {0, · · · , n − 3} are connected by way of twophoton transitions Ω 2 with detunings δ Fig. S2(a)], and subspaces with l ∈ {n − 2, n − 1, n + 1} are connected by a pair of two-photon (Ω 2 ) and three-photon (Ω 3 ) couplings with detunings δ Fig. S2(b)], thereby setting a decoherence-free subspace for n. The efficacy of our approach is then given by the Rydberg blockade condition δV l+2,l / > w d , whereby we find d B /a 0 > l max . For our analysis of Figs. 2-4 in Ref. [1], d B /a 0 4 with l max = n + 1 = 2, so that the strong blockade condition d B /a 0 > l max is met for n = 1.
For the 1D staggered triangular lattice in Fig. S1(a) ( Fig.  1(a) of Ref. [1]), the position vectors are x i = {(k − 1)a 1 , 0} for odd sites (i = 2k + 1) and x i = {a 0 cos θ + (k − 1)a 1 , a 0 sin θ} for even sites (i = 2k), where θ = arccos(ξ/2) and ξ = a 1 /a 0 . Fig. S3 shows the finite-range behavior of nonlocal coupling strength J ij between |r (1) i and |r (1) j states for atoms interacting in the van der Waals regime (p = 6). For the sufficiently large ξ > 1, we find that the exchange rate J ij is significantly diminished for sites |i − j| > 2 due to the ∼ 1/r 6 range of the van der Waals interaction. For the following discussion on the spectrum of effective spin HamiltonianĤ (1) xy , we truncate our analysis up to next-nearest-neighbor interactions.
The effective HamiltonianĤ (1) xy =Ĥ ls +Ĥ xy is fully parametrized by the parameter f ij , hence ξ = a 1 /a 0 [See Fig. S1(a)]. The eigenstate in the n = 1 subspace can be written as |ψ =  (α i = 0) can occur by way of the destructive quantum interference between excitation pathways connecting "bright" state |r (1) i (decay rate Γ 10 4 Γ r ) to "metastable" states |r (1) j =i (decay rate Γ r ) via coherent population trapping. The emergence of "dark state" for such a toy model provides an insight on the choice of interaction parameter ξ 1 → 6 √ 3, whereby J 12 = −J 13 , resulting perfect destructive interference for the edge atoms 1, N . This is indeed confirmed by solving the full spectrum of the sparse Hamiltonian matrixĤ (n) xy with J |i−j|>2 → 0 and by numerically simulating the stationary state of the master equation.

N-PARTITE UNCERTAINTY WITNESS
In this section, we describe our method of constructing the N -partite uncertainty witness from Ref. [4]. Our entanglement witness {∆, y c } consists of identifying the boundaries ∆ as well as their mixed siblings with less entanglement depth. As shown in Refs. [3,4], the lower bound of ∆ Fig. 3 of Ref.
[1], we depict the boundaries ∆ for all possible realizations of fully separable states, bipartite entangled and tripartite entangled states, respectively, by following the procedures of Refs. [5,6].
Generally, we can determine the projection operatorsΠ i = |W i W i | with i ∈ {1, · · · , 2 m } for arbitrary number of systems N m = 2 m with the recursive relationship, i |r i and |G (m) = |g · · · g for 2 m atoms. As discussed in [4], we then construct the uncertainty witness ∆ = i δ 2Π i to identify the bounds {∆ For convenience, we set the maximal N m ≥ N s to be larger than the systemsize N s , so that we could distinguish the entanglement depth k for most k ≤ N s .
For Fig. 4 of Ref.
[1], we assumed the stationary limit, so that y c → 0. In order to verify the minimum bounds ∆ (k−1) b , we only need to optimize the overlap of pure (k − 1)partite entangled states of the form |ψ with one of the projectors |W i . This is achieved when the test state is a balanced (k − 1)-partite Wstate (i.e., |α i | = 1/ √ k − 1). Fig. S4(a) depicts the uncertainty bounds ∆ (k−1) b with k ∈ {1, · · · , N m } calculated for y c = 0 and N m = 2 7 = 128. The shaded regions represent the parameter spaces for which ambiguity exists for ∆ . This is caused by the nonlinear structure for k ∈ {1, · · · , 128} and Nm = 128. Shaded regions indicate the parameter spaces for which ambiguity exists for ∆ , due to the nonlinear sensitivity of ∆(|ψ . For such regions, we conservatively quote the minimum k-value for the genuine k-partite entanglement with ∆(ρ (r) The minimum entanglement depth km certified by {∆, yc} (red dots), and the entangelment depth k in the purported eigenstate |ψ eig = |W k A ⊗ |g · · · g B (blue dots), with fully balanced k-partite W state |W k A.
For such regions, we conservatively quote the minimum value of k m to certify the presence of genuine entanglement depth k m + 1 stored in the purported stateρ Fig. S4(b)]. Hence, the entanglement depth k m (red dot) is a conservative estimate, which can be detected in an experiment, as opposed to the model-dependent analysis of k (blue dot) for the pure state form |ψ eig = |W k A ⊗ |g · · · g B (e.g., by counting the number of non-zero probability amplitudes in |W k A ).
Experimentally, the witness {∆(ρ x , sin θ dσ (i) y } and the excitation statistics {p 0 , p 1 , p ≥2 }, where θ d is the detection angle in the transverse plane x − y. As discussed in Ref. [6], ij |d ij | is the average off-diagonal coherence d ij = |g i r| ⊗ |r j g| for the reduced density matrix ρ 1 in the single-excitation subspace. Since min δ 2 S t θ d = 2 ij |d ij |, we find the following upper bound of the measured variancẽ (S17) The quantum statistics y c = 2N can be detected by the total excitation statistics {p 0 , p 1 , p ≥2 } with MCP ionization signals. Hence, our entanglement witness can be readily implemented even for low-resolution Rydberg experiments without the capability to locally detect the state of single atoms in the lattice.

EXPERIMENTAL PARAMETERS WITH ALKALI ATOMS
Here, we consider 85 Rb atoms interacting with optical field near the transition between |g = |5S 1/2 and |r = |nD 3/2 with two-photon Rabi frequency Ω drive = Ω 1 Ω 2 /∆ and detuning δ. As shown by Fig. S5, this is achieved by a twophoton transition with Rabi frequencies Ω 1 , Ω 2 via the intermediate state |e = |5P 1/2 with one-photon detuning ∆ . For the one-dimensional lattice with ξ > 1 in Ref. [1], as d B /a 0 > 1, the Rydberg excitation spectrum displays a highly nonlinear ladder l due to the dipole-dipole interaction ∆ (ij) p = C p | r i − r j | −p , with the most shifted level given by configuration states consisting of nearest-neighbor excitations |r In order to achieve the parameter sets of Figs. 2-4 in Ref.
The blockade shift ∆ (ij) p is determined by Rydberg coefficient C p , for which we take C 6 = 740 GHz·µm 6 for the van der Waals interaction between two Rydberg atoms in |r = |62D 3/2 , m j = 3/2 [8]. For Ω drive = 30 GHz and δ = 0.65 GHz, the blockade shift for nearest-neighbors is ∆ (i,i+1) 6 = 1.3 GHz, while the power-broadened linewidth for the transition |g ↔ |r is w d √ 2Ω drive . The resulting blockade radius is r B 6 C 6 /w d = 2.9µm. Hence, Rydberg atoms interacting in the strong blockade regime with the lattice constants a 0 ∼ 1µm (Figs. 2-4 of Ref. [1]) can be spatially resolved by way of quantum gas microscope [10,11], so that Ω d can be locally addressed to the defect sites. Therefore, the pumping time for F > 0.95 for Figs. 2-3 of Ref. [1] is then τ p ∼ 10 2 /Γ = 3µs, which is not limited by the ionization time-scale τ π ∼ 40µs. In addition, since the quantum jumps in the n = 1 subspace occur on a time-scale of t j ∼ O(N 2 ) due to the random walk for |r (1) i until it reaches the "defect" sites with Γ i = Γ Γ r , we expect that the pumping time to reach stationarity also scales as t p ∼ O(N 2 ). By extrapolation, we estimate that our method could be extended to generate 30-partite entangled steady states with the parameters {Ω 1,2 , δ, ∆ , |g , |e , |r }. Further improvement is possible by using Rydberg states |nS 1/2 , which have an order-of-magnitude smaller photo-ionization cross-sections σ π and longer radiative lifetimes 1/Γ r [9]. Alternative strategies, including the use of photonic crystals with excited atoms in low-lying electronic states n p , will be discussed elsewhere.