Understanding 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 background
A 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.
a, 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 bounds
Our 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 laws
Our 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 space
We 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\).
An 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).
For 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.
Data 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 < 0.95; Supplementary Section B.1). The top axis shows the lattice constant exclusively for 1D arrays. The dashed lines represent analytical scaling.
Experimental implications
The scaling laws impact a broad range of areas, including quantum optics, out-of-equilibrium phase transitions, quantum simulation and fault-tolerant quantum computation.
(1) Transient superradiance beyond the Dicke limit. Our findings on the scaling of \(R_\star\) crucially address fundamental problems in quantum optics, such as transient superradiance in extended systems16,19,36. For instance, they set a rigorous upper bound on the scaling of the superradiant burst. Although this upper bound may be violated if light is collected only over a small solid angle, new scaling laws can be derived taking into account the detector aperture. Furthermore, in Supplementary Section C.1, we reveal a connection between the early time correlations that determine the appearance of a superradiant burst19 (achieved under dynamical evolution) and \(R_\star\) (which may be inaccessible by dynamics). Although our approach does not capture dynamical evolution, it allows us to predict that the timescale of a superradiant burst from an initially fully excited array of atoms in free space scales (up to a logarithmic factor) as \({T}_{\mathrm{R}} \sim 1/{\varGamma }_{\max }\) (Supplementary Section B.3). This timescale imposes constraints on the array size required for retardation effects to be negligible, thereby ensuring the Markovian evolution assumed in equation (1).
(2) Superradiant lasing in free space. Superradiant lasing, in which incoherently pumped atoms spontaneously radiate coherent light, is known to occur in cavities1,2. It is, however, unknown whether this phenomenon occurs in other environments. By means of our scaling laws, in Supplementary Section C.2, we derive an upper bound on the emitted intensity, I ≤ ℏω0\(R_\star\), and on the optimal pump rate, \(W_\star\) ≤ 2\(R_\star\)/N. This upper bound is tight for the Dicke limit (that is, for atoms in a cavity), yielding \(W_\star\) = NΓ0/2. These results rule out superradiant lasing for 1D arrays in free space (as there is never a superlinear scaling for the intensity). They suggest that a superradiant lasing phase transition might occur for 2D and 3D arrays, with a scaling of the lasing region determined by equation (7).
(3) Driven-dissipative Dicke phase transition. Coherently pumped atoms in a cavity exhibit a second-order phase transition (at a critical pumping rate NΓ0/2) from a magnetized phase characterized by a collective polarization to a paramagnetic phase3,4,5. Beyond the well-studied case of atoms in a cavity, this phenomenon—also called collective resonance fluorescence—has remained largely unexplored until recently5,52,53,54,55. Our scaling laws can be harnessed to readily show that \(R_\star\) sets the scaling of both the threshold drive intensity (\({\eta }_{c}{\sim}\sqrt{{\varGamma }_{0}{R}_{\star }/N}\)) and the emitted light intensity below the threshold (I ∼ ℏω0\(R_\star\)) (Supplementary Section C.3). These analytical results are in agreement with recent numerical studies on collective resonance fluorescence of free-space arrays54,55.
(4) Quantum simulators and processors with Rydberg atoms. For microwave transitions, such as between Rydberg states (where λ0 ≈ 10 mm), typical experimental implementations lie within the collective decay regime. Enhanced decay in dense Rydberg ensembles has been attributed to Dicke-like superradiance56. In Rydberg manifolds, collective decay may increase leakage error rates from the computational subspace in quantum processors, particularly in massively parallelized architectures with macroscopic Rydberg occupation. In quantum simulators, collective decay may become a relevant error source in large systems operating over long simulation times, potentially contributing to Rydberg avalanches whose microscopic origin remains unclear57,58.
As an example, we estimate the potential impact of collective decay using the setup in ref. 59, which uses a 2D array of tweezer-trapped 87Rb atoms for quantum computing. A similar platform has also been used to simulate complex many-body quantum dynamics60. The sources of leakage from the Rydberg state 53S1/2 include both collective (\({\varGamma }_{{\rm{col}}}\equiv {\varGamma }_{\max }/4\approx {R}_{\star }/N\)) and independent (that is, single-atom spontaneous emission, Γsp) decay to other levels, as well as black-body-induced transitions at a rate Γbb. The total leakage rate Γtot = Γcol + Γsp + Γbb is the difference between the (possibly black-body enhanced) transition rates from and to the 53S1/2 Rydberg level. Figure 4 shows the ratio between the total leakage rate Γtot and the coherent Rydberg–Rydberg interaction rate Vnn = C6/d6 between neighbouring atoms. The solid white lines in Fig. 4 mark the boundaries of the Dicke and dimensional scaling regimes (Methods and Supplementary Section C.4). In the transition region, the ratio Γtot/Vnn is obtained by interpolating between the two asymptotic scalings. To estimate Γcol, we adopt the prefactor derived from the product-state ansatz among the possible values consistent with the bounds in equation (4). This represents a realistic estimate, since such states are both physically well-motivated and experimentally accessible, and they reproduce \(R_\star\) exactly in the Dicke limit. Similar considerations apply to other highly excited states, which can evolve—over a timescale \({\varGamma }_{\,\text{col}\,}^{-1}\)—into configurations with leakage rates scaling as \(R_\star\) through transient superradiance (see point (1) above). For simplicity, although the Rydberg state can decay via numerous channels, we consider only the transition with the largest \({\varGamma }_{\max }\) when evaluating Γcol, as it dominates all others for large N (refs. 61,62).
Fig. 4: Implications of the scaling laws for quantum computing and simulation with Rydberg atom arrays.
Ratio between the largest total leakage Γtot (defined in the main text) and the Rydberg–Rydberg interaction Vnn, for the 53S1/2 state in 2D Rydberg arrays of 87Rb atoms. The solid white lines delimit the regions in which the collective decay rate follows Dicke (lower line) and dimensional (upper line) scalings. The transition region interpolates between the two regimes. Above the dashed line, collective decay exceeds the total independent leakage Γsp + Γbb originating from spontaneous emission and black-body processes at T = 300 K (Methods). The dash–dotted line indicates the state-of-the-art limit size, imposed by the objective’s field of view (~1 mm).
Collective decay remains sub-dominant in current platforms of a few hundred atoms, but next-generation arrays hosting thousands of atoms63 should enable the direct exploration of its detrimental effects. We estimate that a 40 × 200 array of Rydberg atom pairs (with 2-μm intrapair and d = 12 μm interpair spacing) performing parallel gates as in ref. 59 would reach a gate error of ε ≈ 0.3%, with collective decay contributing one-third (0.1%) of the total error budget (Supplementary Section C.4). The results shown in Fig. 4 for 2D arrays also hold in three dimensions, as the ratio Γtot/Vnn only slightly increases with dimensionality. This suggests that decoherence due to collective decay becomes more relevant in 3D arrays, where larger arrays can be built. Our analysis indicates that further investigation of collectively enhanced leakage errors in large-scale Rydberg atom arrays is critical to assess the scalability of these platforms for quantum computation and simulation purposes. In particular, the use of alternative atomic species or other Rydberg levels may allow for dimensional scalings (instead of the more detrimental Dicke scaling), thereby reducing collective effects in such platforms.
(5) Quantum error correction and typical states. The scaling law in equation (7) indicates that correlated decay may hamper quantum error correction7,64,65, as the error rate per qubit scales (in the worst case) as ~\(R_\star\)/N, which grows with N in two dimensions and above. Correcting for these errors requires acting on a timescale that shortens with N. Nevertheless, in Supplementary Section A.6, we prove that the decay rate for typical stabilizer states is close to NΓ0/2, implying that they do not experience correlated decay, due to random phases between qubits. This does not mean that the scaling laws for \(R_\star\) are irrelevant in practice, since even simple states like the product state \({\left\vert +\right\rangle }^{\otimes N}\) may be superradiant. Understanding the effect of collective decay on specific quantum error correcting codes warrants further investigation.
Outlook
Our results highlight a class of universal scaling laws that govern the fundamental limits of decay in many-body quantum systems. These findings open several avenues for both experimental exploration and theoretical development. One important direction is quantum metrology, where correlated decay may impose fundamental constraints. Recent experiments on lattice clocks30 and spin squeezing31,32,33 have investigated the role of Hamiltonian power-law dipole–dipole interactions. The dissipative counterpart of the interaction is typically neglected (as dephasing noise is currently the main source of error), although it sets a fundamental limit on the time available to generate and utilize metrologically useful states. Although the decay rate of typical (Haar-random) states is close to NΓ0/2, quantum metrology protocols often rely on a carefully selected set of states (for example, Dicke states or squeezed states), which provide metrological advantage66. Investigating the many-body decay of such states will prove essential to establish possible limitations on such schemes.
From a theoretical perspective, our approach illustrates the power of quantum approximation techniques to predict relevant properties of open many-body systems. Our formalism can be extended to yield upper bounds on the rates of change of higher-order observables such as k-point correlation functions (in this work, \({\hat{n}}_{{\rm{exc}}}\) corresponds to k = 1). Replacing \({\hat{n}}_{{\rm{exc}}}\) in equation (2) with a k-local observable, the ‘auxiliary’ Hamiltonian \({\hat{H}}_{\varGamma }\) will, in general, be (k + 1)-local, and the optimal product state provides an approximation with a multiplicative error of at most 3k+1, independent of system size15. Furthermore, our approach could be adapted to study the steady-state behaviour of driven-dissipative systems and to predict the scaling of correlation functions or other physical observables49. More advanced SDP relaxations such as the quantum Lasserre hierarchy48 could yield tighter bounds for generic systems. These methods are not restricted to spin models, and could be extended to study ensembles of interacting fermions or bosons, as well as to disordered systems. Moreover, invoking time-reversal arguments, our scalings apply to the maximal absorption rate, which may have implications for quantum batteries and light harvesting protocols. Given their generality, we anticipate that these ideas will become a powerful tool to investigate universal properties of large-scale many-body open quantum systems.