Universal scaling laws for correlated decay of many-body quantum systems

MainUnderstanding the quantum dynamics of far-from-equilibrium open many-body systems is a major frontier in physics. From a fundamental perspective, the interplay between energy pumping and dissipation allows for the emergence of phases that transcend the paradigms established by equilibrium statistical physics. Examples in quantum optics include the superradiant laser1,2 and the driven Dicke phase transition3,4,5. From an applied standpoint, the full potential of quantum technologies—including quantum computing, quantum simulation and metrology—is realized only with large systems that remain coherent despite their coupling to a bath.In systems formed by many particles, the always-present vacuum fluctuations mediate long-range dissipative interactions that cannot be switched off, inducing correlated decay that may increase with system size. Such decay processes are collectively enhanced if the particles are tightly packed. Correlated decay may, thus, become the ultimate source of decoherence for many quantum technologies. For instance, it may alter the signal-to-noise ratio in metrology experiments such as atomic clocks or spin squeezing. Similarly, in large-scale quantum computers, it can lead to much shorter coherence times than the predicted timescales using independent noise models and may hinder quantum error correction6,7. On the other hand, correlated decay is a critical requirement for other applications, such as the development of new light sources1,2,8, the dissipative preparation of correlated many-body states9,10 or the protection of logical quantum information via dissipation11,12.Due to the exponential complexity associated with large quantum systems, exactly computing the largest decay rate is a formidable challenge. This problem remains unsolved except in trivial cases, such as permutationally symmetric models (for example, atoms coupled to a cavity) and non-interacting systems. In generic situations, finding the largest decay rate is as difficult as determining the ground state of a general 2-local Hamiltonian, which is known to be a quantum Merlin–Arthur-complete problem—believed to be hard even for a quantum computer13. This complexity is compounded by the diversity of experimental platforms, with many candidate qubits (neutral atoms, molecules, ions, superconducting qubits, quantum dots, vacancy centres among others) and various mediators of their interactions (electromagnetic field and other bosonic collective excitations such as phonons, magnons and so on).In this work, we find upper and lower bounds to the maximal decay rate by leveraging tools from Hamiltonian complexity theory14,15 and applying them in the context of out-of-equilibrium quantum dynamics. For a large class of systems, these bounds are asymptotically tight, thereby yielding scaling laws with system size that only depend on the spectral properties of the decoherence matrix Γ, whose dimension is linear in system size. The bounds are obtained by means of product-state ansatzes and, thus, imply that entanglement does not play any role in the scaling. Our results are formally rigorous and do not rely on mean-field approximations: although the bounds are derived using product states, they remain valid for any state (including entangled states), thereby going beyond previous approaches based on mean-field theory16 or approximate numerical methods17,18,19.We apply these tools to the specific case of ordered atomic arrays20,21,22 and lattices in free space23, which have become an all-around platform for different quantum technologies, ranging from quantum computing24 and quantum simulation25,26,27 to atomic clocks28,29,30 and spin squeezing31,32,33. In the physically relevant regime of lattice constants similar to the resonance wavelength, the maximal decay rate scales as \({\sim}{N}^{\frac{3}{2}-\frac{1}{2\text{D}}}\), where D is the array dimensionality. This scaling law is universal as it does not depend on specific details of the array such as lattice geometry or atomic polarization and has implications in a broad set of problems ranging from quantum dynamics to metrology and quantum computation.Theory backgroundA broad class of Markovian many-body open quantum systems of N qubits (Fig. 1a) is described by the Lindblad master equation$$\begin{array}{ll}\mathop{\hat{\rho}}\limits^{\cdot }=-\frac{{\rm{i}}}{\hslash }[\hat{H}+{\hat{H}}_{\mathrm{local}},\hat{\rho }]+{\mathcal {D}}_{\mathrm{local}}(\hat{\rho })\\\quad+\mathop{\sum}\limits_{i,\,j=1}^{N}{{\varGamma }}_{ij}\left({\hat{\sigma }}_{i}^{-}\hat{\rho }{\hat{\sigma }}_{j}^{+}-\frac{1}{2}\{{\hat{\sigma }}_{j}^{+}{\hat{\sigma }}_{i}^{-},\hat{\rho }\}\right),\end{array}$$ (1) where \({\hat{\sigma }}_{i}^{\pm }=({\hat{\sigma }}_{i}^{x}\pm {\rm{i}}{\hat{\sigma }}_{i}^{y})/2\) are the raising and lowering operators for qubit i. The first term describes coherent evolution, governed by an arbitrary qubit Hamiltonian \(\hat{H}\) that commutes with the total excitation operator \({\hat{n}}_{{\rm{exc}}}={\sum }_{i}{\hat{\sigma }}_{i}^{+}{\hat{\sigma }}_{i}^{-}\), and an arbitrary local Hamiltonian \({\hat{H}}_{{\rm{local}}}\) that can be written as a linear combination of local Pauli operators (that is, each Pauli operator acts non-trivially on a constant number of qubits). The remaining terms describe dissipation, either local—in the form of a linear combination of Lindbladians \({{\mathcal{D}}}_{{\rm{local}}}\) in which all Lindblad operators are arbitrary local Pauli operators—or collective, generated by Lindblad operators with non-local support that act on an extensive number of qubits. Precise definitions of \({\hat{H}}_{{\rm{local}}}\) and \({{\mathcal{D}}}_{{\rm{local}}}\), which model non-collective processes, are provided in the Methods.Fig. 1: Generic qubit ensemble described as an out-of-equilibrium, open, many-body quantum system.Full size imagea, In Markovian baths, integrating out the environment degrees of freedom yields a spin model with coherent and dissipative interactions. The dissipative couplings between N qubits are given by the decoherence matrix \({{\varGamma }}={({\varGamma }_{ij})}_{i,j = 1}^{N}\). b, For a closed system, the ground state is the state of minimal energy. For an open system, finding the state with maximal decay rate (\(R_\star\)) is analogous to finding the ground-state energy of a Hamiltonian.Collective dissipative interactions are represented by the Hermitian decoherence matrix \({{\varGamma }}={({\varGamma }_{ij})}_{i,j = 1}^{N}\). For the master equation to describe a physically valid evolution (that is, a completely positive and trace-preserving map), Γ must be positive semidefinite (that is, Γ ≽ 0). This ensures non-negative eigenvalues {Γμ}, which are physically interpreted as collective transition rates. Since Γ ≽ 0, the spectral norm \(\left\Vert {{\varGamma }}\right\Vert\) is equal to \({\varGamma }_{\max }\), the largest collective transition rate (that is, the largest eigenvalue of the matrix). We define \({\varGamma }_{0}=\mathop{\sum }\nolimits_{i = 1}^{N}{\varGamma }_{ii}/N\) to be the average individual decay rate.We note that the above master equation encompasses a wide variety of physical scenarios, including arbitrary Hamiltonian interactions, coherent driving, diverse decoherence channels (such as incoherent driving, coupling to a finite-temperature reservoir, and dephasing) and disorder in Γ (Supplementary Section A.1 provides a detailed discussion). However, for the sake of simplicity, we omit the local terms \({\hat{H}}_{{\rm{local}}}\) and \({{\mathcal{D}}}_{{\rm{local}}}\) in the following discussion, a choice that we justify later.The instantaneous decay rate of the many-body system, R, is computed as the expectation value of an ‘auxiliary’ (and Hermitian) Hamiltonian \({\hat{H}}_{\varGamma }\) (ref. 34), that is,$$R=-\frac{d}{dt}\langle {\hat{n}}_{{\rm{exc}}}\rangle \equiv \langle {\hat{H}}_{\varGamma }\rangle /\hslash$$ (2) where$${\hat{H}}_{\varGamma }=\hslash \mathop{\sum }\limits_{i,j=1}^{N}{\varGamma }_{ji}{\hat{\sigma }}_{i}^{+}{\hat{\sigma }}_{j}^{-}=\hslash \mathop{\sum }\limits_{\mu =1}^{N}{\varGamma }_{\mu }{\hat{c}}_{\mu }^{\dagger }{\hat{c}}_{\mu },$$ (3) and \({\varGamma }_{ji}={\varGamma }_{ij}^{* }\). The last equality is achieved by means of collective jump operators \({\hat{c}}_{\mu }=\mathop{\sum }\nolimits_{i = 1}^{N}{\alpha }_{i}^{(\mu )}{\hat{\sigma }}_{i}^{-}\), with \({\bf{\alpha }}^{(\mu )}\) being the normalized eigenvectors of Γ (in this notation, the largest eigenvalue is \({\varGamma }_{\max }\equiv {\varGamma }_{1}\)). Generically, the auxiliary Hamiltonian \({\hat{H}}_{\varGamma }\) describes an XY model defined on a weighted interaction graph with a local transverse field. In the specific case where the interactions are mediated by the electromagnetic field, Γ is proportional to the electromagnetic Green’s function35 and the decay rate is exactly equal to the photon emission rate integrated over all emission angles.Lower and upper boundsOur goal is to set the theoretical limits on the maximal decay rate \(R_\star\), the maximum value of R over all possible many-body quantum states. As demonstrated in Supplementary Section A.2, \(R_\star\) also provides upper bounds on the rates of change for general observables. For example, the average coherence \({N}^{-1}{\sum }_{j}\langle {\hat{\sigma }}_{j}^{-}\rangle\) changes at a rate at most \(\sqrt{2{R}_{\star }{\varGamma }_{0}}\), which captures the collectively enhanced decoherence due to correlated decay.Computing \(R_\star\) amounts to calculating the spectral radius of \({\hat{H}}_{\varGamma }\) (since \({\hat{H}}_{\varGamma }\succcurlyeq 0\)) or, equivalently, the ground-state energy of \(-{\hat{H}}_{\varGamma }\) (Fig. 1b). Finding the exact energy is expected to be hard13, except in two limiting cases. For non-interacting qubits (with Γij = Γ0δi,j), \(R_\star\) = NΓ0. In the Dicke limit (that is, with all-to-all interactions such that Γij = Γ0 ∀ i, j), \(R_\star\) = N(N + 2)Γ0/4 (refs. 16,36). These two cases serve as trivial lower and upper bounds, respectively, for \(R_\star\) in arbitrary environments. In general, for systems with non-negative dissipative couplings (Γij ≥ 0 ∀ i, j), we find \(R_\star\) ∼ ∑i≠jΓij (Supplementary Section A.3). In what follows, however, we do not make such a restrictive assumption.We establish bounds on \(R_\star\) using product states. Although other methods exist, our approach based on product states provides an important physical insight: entanglement is not necessary for a system to dissipate at a rate near the theoretical maximum scaling. This finding complements previous observations on the role of entanglement in spontaneous transient superradiance37. We stress that our use of product states here is fundamentally different from many previous analytical studies based on mean-field theory, since the state decaying at rate \(R_\star\) is, in general, highly entangled.We obtain the upper bound by harnessing well-established theoretical guarantees for product-state approximations15, and the lower bound from the variational product state ansatz \(| \psi \rangle {\propto \bigotimes }_{j=1}^{N}({| g\rangle }_{j}+{{\rm{e}}}^{-{\rm{i}}{\phi }_{j}}{| e\rangle }_{j})\), where {ϕj} are locally dependent relative phases between \(\left\vert e\right\rangle\) and \(\left\vert g\right\rangle\) determined by the dominant eigenvector of Γ. This yields the general bounds (Methods) as$$\max \left\{N{\varGamma }_{0},\frac{N{\varGamma }_{\max }}{4({\Delta }^{2}+1)}\right\}\le {R}_{\star }\le \frac{N}{2}\left(3{\varGamma }_{\max }-{\varGamma }_{0}\right),$$ (4) where \(0\le \Delta \le \sqrt{N-1}\) is the relative fluctuation of the entries of \({\bf{\alpha }}^{(1)}\), the dominant eigenvector of Γ. More specifically, we define \(\Delta =\sqrt{\,\text{Var}\,(\left\vert {\alpha }^{(1)}\right\vert )}/\overline{\left\vert {\alpha }^{(1)}\right\vert }\), where \(\overline{\left\vert {\alpha }^{(1)}\right\vert }\) and \(\,\text{Var}\,(\left\vert {\alpha }^{(1)}\right\vert )\) are the mean and variance of the absolute entries of \({\bf{\alpha }}^{(1)}\), respectively. The lower bound is tighter if decay is delocalized (that is, if the brightest collective jump operator has approximately uniform spatial support over all qubits), characterized by Δ = O(1). In particular, for a translationally invariant system, Δ = 0. For non-interacting qubits with identical decay rate Γii = Γ0 (such that \({\varGamma }_{\max }={\varGamma }_{0}\) and \(R_\star\) = NΓ0), both bounds are saturated. The lower bound is saturated for large N in the Dicke limit, where \({\varGamma }_{\max }=N{\varGamma }_{0}\), Δ = 0 and \(R_\star\) = N2Γ0/4 + O(NΓ0). Equation (4) also implies that for ‘sufficiently weak’ interactions (such that \({\varGamma }_{\max }\) is asymptotically independent of N), the maximal decay rate scales only linearly with system size. This generalizes some of the authors’ recent results on the impossibility of Dicke superradiance with nearest-neighbour interactions34, to systems with arbitrary interaction range and geometry.Universal scaling lawsOur bounds are tight for systems with delocalized decay (differing only by a constant factor), thereby yielding scaling laws for the maximal decay rate. Taking Δ = O(1) in equation (4), we find$${R}_{\star }{\sim} N{{\varGamma }}_{\max },$$ (5) which is one of the main results of this paper. Despite its apparent simplicity, the scaling law \({R}_{\star }\sim N{{\varGamma }}_{\max }\) in the delocalized regime is non-trivial and certainly not true for arbitrary systems. More broadly, in Supplementary Section A.5, we prove that there are no general scaling laws on \(R_\star\) that depend solely on the system size and the spectrum of Γ, which we show by constructing explicit, though non-physical, mathematical counterexamples. Since \({\varGamma }_{\max }\) can be computed numerically in O(N3) time, the scaling law provides an efficient scheme to approximate \(R_\star\) for large system sizes with quasi-translational invariance (that is, such that Δ = O(1)). Equation (5) remains valid for the general master equation (1), even in the presence of local terms \({\hat{H}}_{{\rm{local}}}\) and \({{\mathcal{D}}}_{{\rm{local}}}\). As we show in the Methods, these terms only shift \(R_\star\) by a subleading amount O(NΓ0). Physically, this implies that the largest decay rate is always dominated by collective dissipation. Moreover, as discussed in Supplementary Section A.1, the scaling law is also robust to disorder in Γ.Our scaling law reveals important insights about the N2 scaling in Dicke superradiance: one factor of N arises from the permutation symmetry, and the other one from the delocalized nature of the dominant decay channel together with a non-vanishing excitation density at large N. It may seem surprising that a product state yields the same asymptotic decay rate as the entangled Dicke state, but this can be thought of as an instance of the quantum de Finetti theorem38. Since \({\hat{H}}_{\varGamma }\) is 2-local, it suffices to only consider the two-body reduced density matrix of the permutationally symmetric Dicke state, which is close to a product state (with trace distance vanishing as 1/N). Our results show that the accuracy of the mean-field (product state) ansatz holds more generally, even when the permutation symmetry is broken.Maximal decay rate of atomic arrays in free spaceWe now focus on ordered lattices of two-level atoms in free space, whose interactions are described by the propagator of the electromagnetic field evaluated at the resonance frequency ω0, which is a long-ranged function with oscillating sign (Supplementary Section B.1). This makes the problem of finding the ground state of \(-{\hat{H}}_{\varGamma }\) non-trivial and is, thus, a perfect candidate to showcase the strength of our theoretical tools. Nevertheless, our formalism is not restricted to electric-dipole-mediated interactions in free space but can also describe magnetic-dipole or electric-quadrupole interactions in arbitrarily complex dielectric structures.In the large-N limit, and for a large range of lattice constants d, the functional dependence on system size of the largest transition rate is only determined by the dimensionality of the array39. One can relate the scaling with N to the presence of divergences of Γ(k) in reciprocal space as ∣k∣ approaches k0 ≡ ω0/c (c, the speed of light) by assuming that the single-excitation eigenstates of large arrays can be asymptotically well-approximated by spin waves40,41,42,43 and then performing a discrete sampling of the Brillouin zone39. Divergences do not occur for one-dimensional (1D) arrays. They appear for two-dimensional (2D) and three-dimensional (3D) lattices, as the number of atoms per volume increases, enhancing the constructive interference of photon emission for certain wavevectors. For a D-dimensional array, the largest transition rate scales as$${{\varGamma }}_{\max }^{(D)}/{{\varGamma }}_{0}{\sim}{({k}_{0}d)}^{-\frac{D+1}{2}}{N}^{\frac{D-1}{2D}}.$$ (6) This result can be easily generalized to compute the scaling of \({\varGamma }_{\max }\) for an arbitrary D-dimensional lattice in δ-dimensional free space, with D ≤ δ (Supplementary Section B.2). For this scaling to determine \(R_\star\) via equation (5), the jump operator associated with the largest transition rate must be delocalized. In three-dimensional free space (δ = 3), this is expected to hold. Although we explicitly verify this delocalization numerically for 1D and 2D arrays (Methods), an exact computation of the delocalization parameter Δ for large 3D arrays is numerically intractable. Nevertheless, as shown below, we independently validate the scaling law for all dimensions (including 3D) via semidefinite-program (SDP) relaxations.The asymptotic scaling of the maximal decay rate depends on the array dimensionality (D ∈ {1, 2, 3}) as$$\frac{{R}_{\star }^{(D)}}{{{\varGamma }}_{0}}{\sim}{N}^{\frac{3}{2}-\frac{1}{2D}}.$$ (7) This dimensional scaling, a main result of the paper, is universal. It is independent of microscopic details (such as lattice constant, geometry and polarization), which only appear as prefactors. The scaling in equation (7) differs substantially from that expected in the Dicke limit as N → ∞. The departure is the largest for 1D arrays, whose largest decay rate effectively scales as that of a collection of non-interacting atoms. In the Methods, we show that the upper bound on \(R_\star\) remains robust against common experimental imperfections, such as position disorder and finite temperature, and provide numerical evidence indicating that the scaling of the lower bound is also preserved.We benchmark our analytical scaling laws via a SDP relaxation44, which yields a rigorous upper bound on \(R_\star\) and, via approximation guarantees, a corresponding lower bound. SDP relaxations have been used to lower-bound different types of ground-state problems45,46,47,48 and, more recently, ground-state observables49. The SDP relaxation is formulated as$$\begin{array}{l}{R}_{\mathrm{SDP}}=\mathop{\max }\limits_{{X}\succcurlyeq 0}\quad \,{\mathrm{Tr}}\,({\it{\Gamma}}{X})\\ \,{\mathrm{subject}}\;\,{\mathrm{to}}\quad {{X}}_{ii}\,=\,1\quad \forall i=1,\ldots,N\end{array}$$ (8) which can be solved in polynomial time. This is a relaxation since we can write R = Tr(ΓM), where M is a positive matrix with elements \({{M}}_{ij}=\langle {\hat{\sigma }}_{i}^{+}{\hat{\sigma }}_{j}^{-}\rangle\). Note that the specific form of Mij imposes physical constraints on M. The SDP relaxation above omits these constraints and, thus, upper-bounds \(R_\star\). Using standard techniques for establishing approximation guarantees of SDP relaxations50, we show that RSDP also provides a lower bound to \(R_\star\) up to a constant factor (Methods).This establishes$${R}_{{\rm{SDP}}}{\sim}{R}_{\star },$$ (9) as a universal relation valid for arbitrary systems, regardless of whether decay is delocalized. Although RSDP is generally not analytically tractable, this result indicates that SDP relaxations offer a powerful numerical framework for extracting empirical scaling laws of \(R_\star\) at large system sizes, particularly for more complicated systems in which analytical approaches are unavailable. Using an SDP solver51, we obtain a good approximation of \(R_\star\) (Fig. 2) for arrays with lattice constant d = 0.4λ0, where λ0 = 2π/k0 is the wavelength associated with the resonant transition.Fig. 2: Scaling with system size of the numerical approximation for \(R_\star\).Full size imageAn approximation to the maximal decay rate \(R_\star\) is given by the SDP solution RSDP for different lattice dimensionalities with lattice constant d = 0.4λ0. For 1D (□) and 2D (△) arrays, the atoms are perpendicularly polarized; for the 3D (○) lattice, the atoms are polarized along one axis of the array. The dashed black line is a guide to the eye representing the analytical scaling law \({R}_{\star } {\sim} {N}^{\frac{3}{2}-\frac{1}{2\text{D}}}{\varGamma }_{0}\). Inset: comparison between the numerical approximation for the maximal decay rate \(R_\star\) given by RSDP and the largest eigenvalue of \({\hat{H}}_{\varGamma }\) for arrays of up to N = 25 atoms obtained by exact diagonalization. The dashed coloured lines indicate the upper and lower bounds, obtained by numerically evaluating equation (4).Source dataFor large N, the numerical approximations to \(R_\star\) given by the SDP follow the analytical scaling law of equation (7), thereby validating the dimensional scaling. For visualization purposes, we shift the dataset corresponding to each lattice dimension by a multiplicative factor (a constant upward shift in logarithmic scale). These shifts do not affect the scaling and highlight the excellent agreement between the numerical results and the analytical scaling laws. Since the dominant decay is delocalized, Fig. 2 validates our scaling law for \(R_\star\). We also compare the SDP results (without shifts) with the exact diagonalization of \({\hat{H}}_{\varGamma }\) for up to N = 25 atoms (inset).For sufficiently large atom numbers, the scalings hold regardless of the lattice constant. For finite N, however, they depend on both lattice constant and lateral size L = N1/Dd (Fig. 3). As discussed in Supplementary Section B.1, by taking the limits of the expression for Γ in the appropriate order (k0d → 0 before N → ∞), we confirm that for L ≪ λ0, one recovers Dicke’s scaling, that is, \({\varGamma }_{\max }=N{\varGamma }_{0}\). For arrays with a large lattice constant, d ≫ λ0, following a similar procedure yields the limit of non-interacting atoms, that is, \({\varGamma }_{\max }\simeq {\varGamma }_{0}\). These results indicate that asymptotic scaling is a robust bulk property, with boundary effects becoming negligible for sufficiently large arrays. We determine the crossover between ‘non-interacting’ and ‘collective’ behaviours by identifying the parameters for which there is an asymptotic change in the scaling of \(R_\star\), from linear to superlinear. For 2D and 3D arrays, the number of atoms required for such a crossover is \({N}^{(\text{crit})}\simeq \eta {({k}_{0}d)}^{6}\), where η = 0.02 and 5, respectively. As expected, for large interparticle distances, the number of atoms required for \(R_\star\) to be superlinear grows rapidly.Fig. 3: Finite-size effects in the scaling of the largest transition rate as a function of lattice constant.Full size imageData are obtained from a best fit to Γmax = βNαΓ0. These data are presented as fit values (solid lines) with the corresponding 1σ confidence interval (coloured regions). For small lattice constants, the scaling reaches the Dicke limit (α = 1), and for very large lattice constants, the ensemble becomes non-interacting (α = 0). For a region of physically relevant lattice constants, the fits agree with the dimensional scaling of equation (7), that is, α = 1/2 − 1/2D. The atoms form a square lattice and are polarized parallel to one axis of the array. The fits are done with n = {7, 6, 9} points equally distributed over a region \({N}_{1{\rm{D}}}\in [2,{N}_{1{\rm{D}}}^{\max }]\), where \({N}_{1{\rm{D}}}^{\max }=\{30000,250,40\}\) for one, two and three dimensions, respectively. The grey area shows the region in which the fit is not accurate (R2 {\varGamma }_{0}\)) for \(N\ge 16{({k}_{0}d)}^{6}/81{\pi }^{2}\).Generalizations of these results to the array regime for d > λ0/2, the case with parallel polarization, as well as 1D and 3D arrays, follow analogously, and are discussed in detail in Supplementary Section B.1.Delocalized decay in finite atomic arraysFor the scaling law \({R}_{\star } {\sim} N{\varGamma }_{\max }\) to hold for atomic arrays, it is necessary for the brightest jump operator \({\hat{c}}_{1}\) to be delocalized such that Δ ∼ O(1). As demonstrated in ref. 40, the collective jump operators of ordered atomic arrays are well-approximated by Fourier modes with quasi-momentum k in the limit of N ≫ 1 and are, therefore, delocalized. We confirm this in Extended Data Fig. 1a,b (top), where we show that the Fourier transform of the brightest collective jump operator is localized (in momentum space) for both 1D and 2D arrays.Specifically, we numerically obtain the principal eigenvector \({\mathbf{\alpha}}^{(1)}={({\alpha }_{1}^{(1)},\ldots,{\alpha }_{N}^{(1)})}^{\mathrm{T}}\) of Γ and compute its momentum-space distribution \(| {\alpha }_{{\bf{k}}}^{(1)}{| }^{2}\), where \({\alpha }_{{\bf{k}}}^{(1)}=\mathop{\sum }\nolimits_{j}^{N}\exp ({\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{j}){\alpha }_{j}^{(1)}/\sqrt{N}\). Furthermore, we show that the width of \(| {\alpha }_{{\bf{k}}}^{(1)}{| }^{2}\), quantified through the variance Δk2, decreases with the size of the array (Extended Data Fig. 1a,b, bottom). These results are obtained for an array polarized along one of the axis of the array. For perpendicular polarization, a similar narrowing of the distribution in k space is observed, although with some oscillation modulating the decreasing uncertainty Δk.Finally, to confirm that the upper and lower bounds to \(R_\star\) are asymptotically tight, we compute Δ for a 2D array (Extended Data Fig. 1c). The values of Δ are seen to saturate to a finite constant value as N increases. This is confirmed in Extended Data Fig. 1d, where for a few values of d in the relevant regime (0.1 ≤ d/λ0 ≤ 0.5), Δ saturates to a constant value for increasing system sizes.For 3D arrays, we expect a similar behaviour. However, due to computational constraints on N, we are unable to provide numerical evidence to support this expectation.Role of position disorder on the scaling laws for atomic arraysHere we analyse the robustness of the scaling law for arrays, that is, equation (7), to displacement in atomic positions. Position disorder occurs as either imperfections in the position of the trap (for example, in tweezer arrays) or temperature fluctuations of the atomic centre-of-mass position (for example, in optical lattices). In both cases, \(R_\star\) is obtained from equation (3), where Γ is replaced by \({{\varGamma }}{\prime}\) computed by sampling equation (32) according to atoms’ position probability distribution. We consider a Gaussian distribution centred around the lattice sites with a standard deviation σ = ηd. If η is kept constant, σ corresponds to different values of the temperature, or of the tweezers’ displacement depending on d. Different experiments use different atoms and atomic transitions, resulting in a wide range of values for λ0 and, hence, d. In the following, we fix η = 0.05, which corresponds to typical values reported in recent experiments63,69,70.We study the effects of disorder on the upper and lower bounds to \(R_\star\) in equation (4) separately.The upper bound is entirely determined by the scaling of \({\varGamma }_{\max }\). In the limit of weak position disorder (η ≪ 1), the largest eigenvalues of \({{\varGamma }}{\prime}\) can be analytically computed using perturbation theory. After averaging over the atom position distribution, we obtain (appendix C of ref. 71)$${\bar{\varGamma }}_{\max }^{{\prime} }\simeq (1-{\eta }^{2}){\varGamma }_{\max }+{\eta }^{2}{\varGamma }_{0}.$$ (37) For 1D arrays, the robustness of the upper bound is sufficient to conclude the validity of the scaling law in equation (7) because a lower bound to the decay rate is always given by NΓ0, which is trivially achieved by \({\left\vert e\right\rangle }^{\otimes N}\). For 2D and 3D arrays, the lower bound depends on the scaling of \({\varGamma }_{\max }\) and of the principal eigenvector of Γ through Δ (equation (4)). We numerically show that Δ ∼ O(1) for the case of 2D arrays. This is done in Extended Data Fig. 2a, where we plot Δ as a function of d for different array sizes. For fixed d, Δ converges to a constant value as N increases (Extended Data Fig. 2b). The numerical results shown in Extended Data Fig. 2 together with the robustness of the scaling of \({\varGamma }_{\max }\) allows us to conclude that lower and upper bounds are asymptotically tight, leading to the scaling law in equation (7) for 2D arrays. For 3D arrays, we expect the scaling law to hold, but computing Δ for large arrays is challenging and we do not have sufficient numerical evidence to support this expectation.Collective decay in Rydberg arraysAtoms excited to the 53S1/2 level can leak to all levels nP1/2 and nP3/2, with n ≥ 5. Out of all possible decay paths, \(R_\star\) is always larger for the 53S1/2 → 52P3/2 transition due to its long-wavelength nature. We simplify the multilevel nature of the system by assuming that for large atom numbers, the dominant collective channel will dominate over the rest61,62, so that \({{\varGamma }}_{\mathrm{col}}={{\varGamma }}_{\max }^{53{{\rm{S}}}_{1/2}\to 52{{\rm{P}}}_{3/2}}/4\). All other transitions are assumed to contribute only to the spontaneous (Γsp) and black-body (Γbb) leakage.The collective contribution to leakage (Γcol) is obtained from the lower bound in equation (4) by assuming Δ = 0. To estimate this quantity for large systems, we use the expansions in equations (35) and (36) when the expansion error ϵ is below 0.05. The error ϵ is defined as the contribution of the higher-order term in the expansions. The remaining region around L ≈ λ0 is calculated via a weighted smooth interpolation using the thin-plate spline method72. We look for a surface Γcol(d, N) that minimizes the function$$\sum _{i}{w}_{i}{\left({\varGamma }_{{\rm{col}}}-{\varGamma }_{\max }({d}_{i},{N}_{i})/4\right)}^{2}+\mu \iint \parallel {\nabla }^{2}{\varGamma }_{{\rm{col}}}\parallel ,$$ (38) where we take as input a 200 × 200 grid of points di, Ni with \({\varGamma }_{\max }({d}_{i},{N}_{i})\) given by the expansions, weights \({w}_{i}=1-{\epsilon }_{i}^{4}\), where ϵi is given by the expansion error ϵ evaluated at each point i, and smoothness parameter μ = 0.1. The grid is taken over a large parameter space, [0.05, 500] μm × [1, 1010] for Fig. 4 to minimize edge effects.