Ученые обнаружили квантовые кинетические ограничения в реляксации двухмерного газа Ридберга

Abstract We investigate the relaxation dynamics of a Rydberg gas in regimes where coherent processes and dissipation compete. In the strongly dissipative limit, the dynamics is known to be governed by an effective classical rate equation and to exhibit kinetically constrained, glassy relaxation towards a trivial stationary state. This behaviour originates from the Rydberg blockade, which prevents simultaneous excitations within a characteristic blockade radius. However, the fate of kinetic constraints in the weakly dissipative limit remains unexplored in large systems above one dimension. To access large system sizes and two-dimensional geometries, we employ the truncated Wigner approximation, a phase-space method that captures correlated many-body dynamics beyond classical rate equations. To probe the emergence of kinetic constraints on timescales where coherent and dissipative processes are comparable, we analyse the relaxation dynamics starting from two initial states: a fully polarised state and a Néel state, which belongs to a manifold of so-called quantum scars. In both cases, we observe a pronounced slowdown in the relaxation of the magnetisation towards the stationary state and identify transient signatures of quantum kinetically constrained dynamics in one and two dimensions. Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1. Introduction Understanding the relaxation dynamics of open quantum many-body systems remains a central topic of nonequilibrium many-body physics, where the interplay between coherent dynamics and dissipation gives rise to complex behaviour [1–6]. Programmable quantum simulators [7] based on ultracold atoms [8–10] and ions [11–13], particularly Rydberg arrays [14–17], provide a versatile setting to investigate such dynamics, as long-lived Rydberg states exhibit strong and tunable long-range interactions. Notably, the strong interactions between Rydberg excitations give rise to a blockade mechanism [18–21], where an excitation suppresses simultaneous excitations within a characteristic blockade radius. As a consequence, (de-)excitation events depend on the arrangement of neighbouring excitations, which imposes a kinetic constraint on the many-body dynamics. This limits the set of accessible configurations and therefore influences relaxation [22–24]. Kinetic constraints similar to those found in Rydberg gases have been studied in classical models of glass formers [25, 26]. Here, configuration-dependent transition rules give rise to slow relaxation accompanied by so-called dynamical heterogeneity, i.e. different spatial regions relaxing on substantially different timescales [27]. While these models provide valuable insight into the constrained dynamics of classical systems, the investigation of kinetic constraints in the quantum regime remains challenging. This is especially true for large systems and in higher dimensions, where the exponential growth of the Hilbert space and entanglement generation further complicate the simultaneous treatment of coherent and dissipative processes [1, 2, 28]. Approximate descriptions of the many-body dynamics are available in certain regimes. In the strong dephasing limit, quantum coherences are suppressed, and the evolution can be reduced to an effective classical stochastic process [29–33]. In this regime, relaxation towards the stationary state proceeds via the aforementioned glassy dynamics characterised by hierarchical relaxation and dynamical heterogeneity. In the weakly dissipative limit, when quantum coherences become potentially relevant, a description of the dynamics in terms of classical rate equations is generally not applicable. In this case, existing studies have been limited to comparatively small system sizes and in one dimension only [34–43], with matrix product state simulations enabling the study of larger one-dimensional systems [44–46]. However, the simulation of dynamics in two spatial dimensions remains difficult due to the substantial computational cost of tensor-network methods in open, long-range interacting systems. While recent advances in two-dimensional tensor-network and operator-based approaches offer promising routes forward, systematic studies of long-time relaxation remain difficult [47–50]. Therefore, understanding how blockade-induced kinetic constraints manifest in the regime where quantum coherent and dissipative dynamics compete remains an open problem in large, two-dimensional Rydberg gases. In this work, we investigate the relaxation dynamics in the regime where coherent driving is comparable to or exceeds dissipation. To access this parameter range for large system sizes and in two dimensions, we employ the truncated Wigner approximation (TWA) [51, 52]. TWA is a phase-space method that approximates quantum dynamics by an ensemble of classical trajectories, incorporating leading-order quantum fluctuations through statistical sampling of the initial state while remaining computationally tractable beyond exact approaches. Expectation values of observables, such as magnetisation, as well as correlation functions are obtained as statistical averages over this ensemble of trajectories. TWA has been widely and successfully applied to isolated quantum many-body systems [53–55], and more recently, it has also been increasingly used for dissipative spin systems [56–68]. Such framework is particularly suitable here, as the essential physics of dissipative Rydberg arrays is well captured by interacting spin models with coherent driving and local dissipation. This paper is organised as follows. In section 2 we specify our dissipative spin model under consideration, and in section 3 we outline TWA for its dynamics. In section 4 we compare TWA to exact dynamics for small system sizes and also to the fully classical description of [29] in the strongly dissipative limit. Subsequently, we look at the dynamics for large system sizes in both , alongside spatial and temporal correlation functions. We conclude in section 5 with a discussion of the results and an outlook. 2. Model The essential physics of Rydberg systems is captured by paradigmatic quantum spin models, for instance the long-range interacting transverse-field Ising model. We consider a system of atoms at the sites of a regular lattice in and dimensions with uniform lattice spacing and periodic boundary conditions. Each atom is modelled as a two-level system with a ground state and a highly excited Rydberg state , as schematically shown in figure 1 for a two-dimensional lattice. At each lattice site , we can map this onto a quantum spin- particle. In this representation, the states and correspond to the eigenstates of the Pauli operator , and we introduce the usual spin operators acting on site . The local Rydberg excitation number operator is given by which projects onto the excited state. We consider two different initial states. The first is the fully polarised state which is the natural initial condition for the relaxation dynamics since all spins start in the ground state. The second one is the Néel state which represents a staggered configuration on the two sublattices. In the limit of strong interactions and in the absence of dissipation, this state exhibits quantum scarring: it belongs to a low-dimensional manifold of states and the ensuing dynamics displays persistent coherent oscillations [69–71]. Considering this initial state therefore allows to illuminate the role of quantum coherence during relaxation. When coupled to an external environment, the evolution of the system’s density matrix is described by the master equation Here, the first term describes the coherent evolution due to the Hamiltonian , and the subsequent dissipator describes the incoherent processes. Each dissipator is associated with a specific jump process, described by a set of jump operators . The coherent part of the system’s evolution is governed by a Hamiltonian, where an external laser field couples the ground and Rydberg states, and the Rydberg–Rydberg interactions are described by a long-range potential, with In the following we consider van der Waals interactions with , and we regard dephasing as the primary source of dissipation. In practice this arises from environmental noise sources, such as laser phase fluctuations or thermal motion, which randomise the relative phase between the ground and Rydberg states, effectively destroying quantum coherence [15]. At site , the jump operator is where is the dephasing rate. 3. Truncated Wigner approach To access system sizes beyond the reach of exact methods, particularly in two spatial dimensions, we employ the TWA. In this approach, the quantum master equation is mapped onto a stochastic evolution of classical phase-space variables, which can be obtained from the Keldysh path-integral representation [68]. A detailed derivation of this mapping and that of the TWA equations used in this work is provided in appendix with the Poisson bracket . In practice, observables are obtained by evolving an ensemble of classical trajectories according to equation (9), with initial conditions sampled from a phase-space representation of the quantum state. The corresponding classical observable is evaluated along each stochastic trajectory and averaged over the ensemble. This procedure approximates the expectation values computed from the density matrix evolved with the Lindblad equation (4). The fields appearing in (9) are defined as where are complex Gaussian noise fields with correlations in this case for local dephasing at rate . The overline in the above equation denotes the average over noise realisations. The Kronecker delta implies that the noise is local in space such that fluctuations associated with different jump operators are uncorrelated. The Dirac delta corresponds to Markovian white noise in time. More general noise structures, including correlated noise, can be incorporated by replacing these delta functions with an appropriate correlation kernel. In the discrete TWA [55, 72–74], the initial phase-space distribution is replaced by a discrete sampling scheme. For the fully polarised state (2), all spins are aligned along the negative -direction, such that and . For the Néel initial state (3), spins on the two sublattices are polarised in opposite directions, with for and for . In both cases the transverse spin components have vanishing expectation values, while their fluctuations are maximal, . Accordingly, the classical spin variables are sampled from the discrete distribution where for the fully polarised state and () for sites on sublattice () in the Néel state. This sampling exactly reproduces the first and second moments of the corresponding quantum spin operators. The sampled spin configurations are subsequently evolved according to the stochastic classical equations of motion (see appendix where and is given by equation (7). For local dephasing with the Hermitian jump operator , the complex noise fields can be recast as a single real noise field acting on the transverse spin components, as only the combination appears in equation (9), giving The stochastic equations are interpreted in the Stratonovich sense, consistent with their derivation from the Keldysh path-integral formalism [68]. For the numerical integration, we employ a rotation-based integrator that evolves each spin by a finite element of . In this scheme, the deterministic and stochastic contributions define an effective field about which the spin precesses, and the discrete update is implemented as the corresponding finite rotation on the Bloch sphere. This construction preserves the spin length exactly and ensures consistency with Stratonovich calculus. Throughout this work we express all lengths in units of the lattice spacing and time is expressed in units of . The time stepping used in the simulations is . The numerical implementation is described in detail in appendix 4. Relaxation dynamics 4.1. Benchmark and kinetic constraints in small systems We first benchmark the TWA approach against numerically exact simulations obtained via stochastic unravelling of the Lindblad master equation based on the quantum-jump Monte Carlo method [75]. This allows us to test how well TWA captures quantum coherent effects in the relaxation dynamics. To this end, we consider the relaxation of the mean magnetisation, where the expectation value is obtained by averaging over an ensemble of quantum trajectories in the stochastic unravelling. The quantum-jump simulations are performed using the package QuantumToolbox.jl [76]. In TWA, the same observable is evaluated by averaging over ensembles of classical trajectories. In this work, we average over 5000 trajectories for one-dimensional, and 3000 trajectories for two-dimensional results. Beyond mean values, the ensemble directly provides the uncertainty of each observable: the standard deviation across trajectories, , which estimates the physical shot-to-shot fluctuations, while the statistical error of the mean is . For the mean magnetisation this error is of order for the trajectory numbers used here, i.e. smaller than the linewidths in the figures shown below. 4.1.1. Coherent and classical relaxation regimes A particularly interesting regime arises when the dephasing rate is much larger than the Rabi frequency, . In this strongly dissipative limit, coherences decay rapidly and can be adiabatically eliminated. The dynamics is then described by an effective classical master equation [29] This describes stochastic spin flips with configuration-dependent rates given by with the interaction parameter , where is given by equation (7). The resulting stochastic dynamics can be simulated efficiently using kinetic Monte Carlo (kMC) as studied in [29], allowing access to much larger system sizes than in the full quantum evolution. The validity of this approach has also been confirmed experimentally [77–79]. Representative data for a spin chain of length , and a spin lattice of with interaction strength are shown for the fully polarised as well as the Néel initial states in figures 2(a) and (b), respectively. In the limit of strong dephasing (, blue curve), the exact, TWA and kMC relaxation curves agree for all timescales for both initial conditions. For smaller dephasing rates (, yellow and red curves), the dynamics changes qualitatively and coherent oscillations are present at early to intermediate time scales. For the Néel initial state these oscillations originate from the quantum many-body scar manifold, which identifies the time window where coherent quantum dynamics dominates before noise-driven relaxation takes over. In this regime, kMC systematically deviates from the exact quantum dynamics, as equation (19) is no longer applicable. On the other hand, TWA remains in good qualitative and quantitative agreement with the exact results at early and late times, although it underestimates the amplitude of the intermediate oscillations. We emphasise that the coherent oscillations are not a perturbative correction to the classical dynamics, but rather, they reflect a regime where coherent many-body dynamics governs the evolution, inaccessible to classical rate equations by construction. Download figure: Standard image High-resolution image However, at long times, the final stage of relaxation is governed by an exponential decay towards the stationary state, which is well captured by both TWA and kMC for all values of , even though the classical rate equation underlying kMC is formally derived only for . In this late-time regime, the system predominantly occupies configurations with magnetisation close to zero, which are strongly interacting and individual spin flips involve a large interaction-induced energy cost or gain. These transitions are therefore strongly off-resonant and cannot be efficiently driven by coherent processes, but instead rely on dissipative spin flips. Consequently, the dynamics is governed by slow, dissipative configuration changes, which are well captured by the classical rate equation (19), even when coherent processes affect the intermediate-time behaviour. The relaxation dynamics therefore involves a crossover from an early regime dominated by coherent many-body dynamics to a late-time regime where dissipation governs the slow approach to the stationary state. As seen in the next subsection, strong interactions also impose kinetic constraints that restrict the allowed configuration changes. Understanding this crossover and the role of kinetic constraints in the weakly dissipative scenario is a central aspect of the behaviour studied here. 4.1.2. Emergence of kinetic constraints at weak dissipation Figure 3 shows the relaxation of for increasing interaction strength in dimensions for the fully polarised initial state. For weak interactions, , TWA quantitatively reproduces the exact dynamics for all times, while for stronger interactions, , the amplitude of intermediate oscillations is gradually underestimated. For larger interaction strengths, , there is a qualitative disagreement between TWA and the exact dynamics. Here, the relaxation is increasingly influenced by kinetic constraints, which originate from the interplay between coherent laser excitation and strong interactions between excited atoms. In the absence of interactions, the drive induces resonant spin flips between the ground and excited states. However, when an atom is already excited, it shifts the energy levels of nearby atoms via the van der Waals interaction, rendering the excitation of a neighbouring atom off-resonant when [29], a manifestation of the Rydberg blockade. As a result, the local configuration determines whether a transition is energetically accessible, imposing an effective kinetic constraint on the dynamics. Since the van der Waals interaction decays as , nearest-neighbour interactions strongly exceed the coherent drive, , while next-nearest-neighbour contributions can become comparable to . Consequently, spin flips in the immediate vicinity of an excitation are strongly off-resonant. This kinetic constraint manifests in two ways: the emergence of a transient plateau in the magnetisation at intermediate times, and a slow, dissipation-dominated final approach to the stationary state. Download figure: Standard image High-resolution image For the fully polarised initial state in one dimension, this occurs at in the exact numerics. This value corresponds to that of a one-dimensional lattice gas subject to a hard-dimer constraint, where the presence of an excitation on a site suppresses the excitation of its nearest neighbours [80–82]. The magnetisation plateau thus suggests that the Rydberg blockade effectively imposes a nearest-neighbour exclusion constraint on the dynamics: once the system reaches a configuration where the simultaneous excitation of neighbouring sites are blocked, further relaxation is strongly suppressed. While TWA does not capture the extended plateau arising from this at , it does reproduce the key signatures of kinetic constraints at : the onset of a plateau in the magnetisation and the slowed, dissipation-dominated final relaxation, as seen in figures 2(a) and 3. For the Néel initial state, no intermediate plateau develops as seen in figure 2(b). In contrast to the fully polarised case, this configuration already contains excitations on alternating sites, so the relaxation proceeds mainly through their rearrangement rather than through a gradual increase in the excitation density that leads to the intermediate plateau. For weaker dissipation, resonant transitions lead to pronounced oscillations, causing the magnetisation to develop a short-lived minimum at intermediate times before slowly relaxing towards the stationary state. At even stronger interactions, the dynamics is effectively captured by a PXP model with dephasing in which Hilbert-space fragmentation gives rise to a large set of long-lived metastable states fixed by the initial configuration [83]. For strong dissipation, the magnetisation changes only slightly during the evolution (figure 2(b), blue curve). This behaviour is consistent with relaxation dominated by noise-induced off-resonant processes. These observations confirm a clear separation between the quantum and classical regimes at weak dissipation and the role of kinetic constraints. At early times, the dynamics is dominated by coherent, resonant transitions for both initial conditions. For the fully polarised initial state, this regime leads to the emergence of the blockade-induced plateau discussed above. The classical rate equation overestimates its extent (red curve in figure 2(a) for ) because it cannot capture the coherent oscillations present at weak dissipation. By contrast, for the Néel initial state the system does not pass through such an intermediate jammed configuration. Instead, coherent effects, that are not captured by the classical rate equation lead to pronounced oscillations that can reach below the minimum predicted by the classical rate equation, before relaxing to the same stationary state. For the two-dimensional lattice, the relaxation dynamics differs qualitatively from the one-dimensional case. For the polarised initial state, only the onset of a plateau is visible, which becomes more pronounced with both increasing interaction strength and as shown in the following section, with system size. This indicates that, while local blockade effects still suppress nearby excitations, they do not generate a sufficiently sharp separation of time scales to stabilise a long-lived plateau. However, the onset of the plateau observed here is itself a signature of quantum-coherent dynamics, as no such feature was observed in the strongly dissipative regime [29] in . In the weakly dissipative case, coherent dynamics gives rise to this plateau, followed by a slower, dissipation-dominated relaxation. The Néel initial state in suffers from finite-size effects for the small lattice. To investigate this further, we apply TWA in the following to larger systems. 4.2. Kinetic constraints in large systems A key advantage of TWA is its favourable scaling with system size compared to exact quantum simulations, where the Hilbert space grows exponentially. We therefore make use of this scalability to investigate substantially larger systems beyond the reach of numerically exact stochastic unravellings of the Lindblad equation. In one dimension, we consider spin chains of length , while in two dimensions we study square lattices of size . We focus on the regime and , where TWA was found to reproduce the exact quantum dynamics for small system sizes. 4.2.1. Magnetisation relaxation We consider different dephasing strengths in figure 4(a) at fixed interaction for (left) and (right) for the polarised initial state. As seen on the left of figure 4(a) for the chain, the relaxation approaches a plateau at , consistent with the hard-dimer constraint discussed in section 4.1. We also observe the onset of an intermediate plateau in two dimensions (right hand side of figure 4(a)), at magnetisation . By analogy with the one-dimensional case, one might expect this corresponds to that of a hard-square lattice gas [84–86], the two-dimensional equivalent of the hard-dimer constraint, where the presence of an excitation suppresses the excitation of its four nearest neighbours. However, this would correspond to a magnetisation [80], which is substantially higher than the observed value. This suggests that in two dimensions the dynamics in the coherent regime is subject to constraints extending beyond simple nearest-neighbour exclusion, indicating that coherent processes generate longer-ranged constraints that suppress spin flips over a larger surrounding region than in one dimension. Download figure: Standard image High-resolution image The relaxation curves exhibit distinct stages as the dynamics crosses over from coherent to dissipation-dominated behaviour. In this process, the magnetisation first approaches the plateau before slowly relaxing towards the stationary state. To analyse these regimes more quantitatively, we extract two characteristic timescales from the dynamics, shown in figures 4(b) and (c). The plateau relaxation time , defined as the time to reach in and in , scales as in both dimensions. The final relaxation time , extracted by fitting in the regime , obeys in both dimensions, indicating that the asymptotic relaxation is governed by the dissipative timescale. Consequently, decreasing simultaneously accelerates the approach to the plateau via enhanced coherent dynamics, while slowing down the final relaxation towards the stationary state. The scaling holds in the regime , where the relaxation is rate-limited by dissipation. At the largest dephasing strengths , deviates from this behaviour and turns over, signalling the crossover to the strong-dephasing (quantum Zeno [87, 88]) regime: there, dephasing suppresses the coherent spin flips driving relaxation, and adiabatic elimination of the coherences gives an effective incoherent flip rate [29], so that . The two regimes cross near , where the relaxation time is minimal. We now turn to the investigation of the Néel initial condition in large systems, as shown in figure 5 for in (a) and in (b). In one dimension, the dynamics is qualitatively similar to that of the chain discussed in section 4.1. In two dimensions, the time scales are more clearly separated. Coherent oscillations dominate the early-time dynamics, while the minimum in the magnetisation develops only later. The minimum reached in one dimension is noticeably lower than in two dimensions. The opposite trend is observed for the fully polarised initial state, where the emergent plateau occurs at a lower magnetisation in . Starting from a lattice with no excitations, each excitation in two dimensions suppresses the creation of additional excitations on a larger number of neighbouring sites. By contrast, starting from the Néel state, changing the excitation pattern leads to larger interaction shifts in two dimensions, which constrains the dynamics more strongly. Download figure: Standard image High-resolution image We now consider the effect of varying the interaction strength. Figure 6 shows the TWA relaxation curves for the polarised initial state at fixed dephasing , in (left) and (right) while varying the interaction strength. Since kMC accurately describes the late-time relaxation even in the weakly dissipative regime, we characterise this slowdown by analysing the approach to the stationary state as a function of the effective interaction parameter . For van der Waals interactions (), interaction-induced detuning shifts enter the classical rate equation (19) through a squared term in the denominator of the rates (20), leading to a suppression of spin flips that scales with the square of the interaction strength. In the strongly interacting regime , the dominant contribution scales as , which leads to an effective rate that scales as , implying a relaxation time . Rescaling time by collapses the late-time dynamics onto a common curve, as shown in the insets of figure 6. This rescaling works similarly when considering the Néel initial state. Download figure: Standard image High-resolution image 4.2.2. Spatial correlations Spatial correlations provide additional insight into the emergence of collective behaviour during the dynamics. To quantify these, we introduce the Mandel -parameter where denotes the variance of the magnetisation. As a reference value for the upcoming analysis, we consider the variance of a non-interacting lattice system with uncorrelated spins, given by The parameter therefore measures deviations from the behaviour of an uncorrelated random lattice gas: corresponds to uncorrelated spins, while indicates the presence of correlations. Figure 7 shows data for the spin chain and the spin lattice at different for the polarised and the Néel initial states. Download figure: Standard image High-resolution image For the fully polarised initial state in the strongly dissipative regime, remains close to zero for most of the evolution, with a small negative value developing at intermediate times. In this limit, strong dephasing suppresses coherent processes and interactions do not generate significant correlations. As the dephasing strength is reduced, develops a pronounced minimum with negative values during the relaxation process. This reflects the build-up of anticorrelations between excitations as the dynamics becomes dominated by blockade-induced kinetic constraints. Correspondingly, the recovery of towards zero becomes progressively slower for weaker dephasing, reflecting the increasingly collective and constrained character of the relaxation. At late times, approaches zero as the system relaxes towards the stationary state. This is consistent with equations (21)–(23), since the stationary state corresponds to an uncorrelated configuration of spins, where . Showing this more explicitly, we first take the magnetisation per site as . The local spin variable has the same mean value at every site in this case, . Since , we also have . The variance of the magnetisation then evaluates to . This gives . The fluctuations thus follow the random lattice gas result , and therefore . In the stationary state specifically, so that . For the Néel initial state, shows a qualitatively different behaviour. Since excitations initially occupy only one sublattice, nearest-neighbour blockade constraints are already satisfied at , leading to strongly suppressed magnetisation fluctuations. Accordingly, starts close to its minimal value and then increases towards zero as the initially ordered configuration is gradually lost during the relaxation. The relaxation again depends on the dephasing strength. For strong dissipation, the increase of is smooth and monotonic, reflecting the incoherent decay of the initially ordered configuration. In the weakly dissipative regime, however, the relaxation proceeds through multiple stages. The system first rapidly loses the initially ordered configuration. This is followed by a long-lived transient regime where remains close to , indicating that strong anticorrelations persist even after the initial Néel order has been lost. In this regime the evolution is constrained by the blockade-induced kinetic constraints, which restrict the allowed spin flips and thereby temporarily limit the accessible configurations. At later times incoherent processes dominate and the system gradually approaches the uncorrelated stationary state, for which , consistent with the random lattice gas. Overall, we also observe that weaker dissipation initially accelerates the approach to this intermediate flat region but slows the loss of correlations at later times, for both initial conditions. This corroborates the observations in figures 4(b) and (c), where smaller values of lead to an earlier onset of the plateau in the magnetisation relaxation and a slower subsequent approach to the stationary state. We analyse the behaviour of temporal correlation functions in appendix 5. Summary and outlook In this work, we have investigated the relaxation dynamics of a dissipative Rydberg gas in one- and two-dimensional systems, focusing on the regime where coherent interactions and dissipation compete. To study this regime, we employed the TWA, which we first benchmarked against numerically exact simulations for small system sizes and then applied to large lattices beyond the reach of exact methods. Within this framework, we analysed the dynamics starting from two distinct initial states; a fully polarised state and a Néel state belonging to a quantum-scar manifold. For the fully polarised initial state in one dimension, the dynamics evolves through an initial regime governed by resonant spin flips that increase the excitation density until the nearest-neighbour blockade restricts further spin flips, leading to the emergence of an intermediate plateau in the magnetisation. The late-time relaxation towards the stationary state is then governed by slow dissipative configuration changes, for both . Remarkably, a similar plateau is also observed in two dimensions in the weakly dissipative regime, which has not been found in the strongly dissipative limit. This occurs at a magnetisation considerably lower than expected from a simple nearest-neighbour exclusion picture. This suggests the presence of effective constraints extending beyond nearest neighbours, whose role will be explored in more detail in future work. Pronounced spatial anticorrelations, which decay to zero when the stationary state is reached further corroborate the presence of blockade-induced constraints. There are also characteristic time scales from the crossover between coherent and dissipative dynamics. Reducing the dissipation shortens the time required to reach the plateau in the coherent regime, while increasing the final relaxation times. For the Néel initial state coherent oscillations dominate the early-time dynamics, however, by contrast the magnetisation develops a pronounced intermediate minimum rather than a plateau in the weakly dissipative regime. During this regime strong anticorrelations persist even after the initial Néel order has been lost, reflecting the continued influence of blockade-induced constraints. At later times incoherent processes dominate and the system gradually approaches the uncorrelated stationary state. The behaviour of the spatial correlations exhibits time scales similar to those found in the relaxation dynamics of the polarised initial state. In particular, a long-lived transient regime of anticorrelations emerges: for weaker dissipation this regime is reached more quickly, but the correlations also decay more slowly. More broadly, our results demonstrate that the TWA offers a computationally efficient and scalable framework for probing the relaxation dynamics of open many-body systems in regimes that are currently inaccessible to exact methods, in particular large lattices and two spatial dimensions at weak to intermediate dissipation, where coherent driving and dissipation compete. The results are directly applicable to the large-scale Rydberg arrays now routinely realised in experiments, and provide a practical tool both for interpreting their nonequilibrium dynamics and for characterising the role of dissipation in Rydberg-based quantum simulation and computation, where decoherence and loss are intrinsic. While TWA enables access to large system sizes and in higher dimensions, its limitations are apparent in regimes where kinetic constraints dominate. At the same time, owing to the linear scaling of the TWA with system size, an extension of the present study to three-dimensional lattices is computationally feasible without much modification: only the lattice setup needs to be adapted. In three dimensions one would expect the blockade-induced constraints to be even more pronounced, as each excitation suppresses six nearest neighbours, although a quantitative benchmark of the method would be limited to small clusters accessible to exact methods. In the one- and two-dimensional systems studied here, for the polarised initial state, only the onset of the intermediate magnetisation plateau associated with a nearest-neighbour excitation exclusion constraint is captured within TWA. Extending the current approach, for instance via BBGKY hierarchies [89] or cluster extensions of TWA [90] could allow for a more systematic incorporation of quantum correlations beyond the present treatment. Tensor-network approaches could also provide an alternative route to explore this regime beyond phase-space methods. This would be especially relevant for accurately describing the constraint-dominated regime in two dimensions, where the formation of the intermediate plateau could be studied at stronger interactions than in this work. Lastly, we note that while the parameter regime explored here is, in principle, accessible in current Rydberg array experiments, the long-time dynamics considered in this work remains challenging to realise in practice. Additionally, we have only considered dephasing as the source of dissipation. A finite Rydberg lifetime would introduce an additional decay channel at rate , setting an observation window for this dynamics. The constraint-dominated features reported here are most pronounced around and thus remain accessible provided , whereas the full relaxation to the stationary state at is more demanding. At longer times, decay competes with excitations, and creates a correlated steady state, as shown in reference [29] in the strongly dissipative limit. The observation of the slow, constraint-dominated relaxation requires coherence and stability over extended times, whereas dissipation in typical setups leads to heating and loss processes that can obscure these effects. Achieving such regimes would therefore require dissipative quantum platforms that remain stable at long times. Promising candidates in this direction include recently proposed Rydberg-ion systems, where engineered dissipation and longer coherence times may enable access to this dynamics [91]. Acknowledgments We thank M. Cech for fruitful discussions, comments on the manuscript and for helpful guidance with quantum jump Monte Carlo simulations. This work is supported by the ERC grant OPEN-2QS (Grant No. 101164443). We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Research Unit FOR 5413/1, Grant No. 465199066 and through the Research Unit FOR 5522/1, Grant No. 499180199. We also acknowledge funding through JST-DFG 2024: Japanese-German Joint Call for Proposals on “Quantum Technologies” (Japan-JST-DFG-ASPIRE 2024) under DFG Grant No. 55456179. The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the DFG through grant no INST 40/575-1 FUGG (JUSTUS 2 cluster). ChatGPT (OpenAI, GPT-3) was used for minor improvements to figure presentation and language editing. Data availability statement The simulation code used to generate the data can be found at [92]. The data that support the findings of this study are openly available at the following URL/DOI: https://zenodo.org/records/19468166 [93]. Appendix A: Truncated Wigner approach to dissipative spins In this appendix, we summarise the necessary steps for the derivation of the TWA equations used in the main text. First, we briefly outline the general Lindblad-to-TWA procedure and then apply it specifically to the spin model considered in this work. The presentation follows from the Keldysh path-integral formulation of Lindblad dynamics and its semiclassical truncation, as detailed in reference [68, 94]. A.1. Keldysh path-integral formulation The Keldysh path integral is defined on a closed time contour running from to and back to [95, 96]. In the absence of dissipation, for the fields on the forward and backward branches, the generating functional reads where is the functional integration measure over the forward and backward fields. The Keldysh action is given by where denotes the Lagrangian density of the system. Expectation values are obtained from where encodes the phase-space representation of the initial density matrix. For an open system described by a Lindblad master equation (4), the action receives an additional dissipative contribution, with [94] where denote the jump operators evaluated on the forward and backward branches of the contour. Here and throughout, an overbar on a single field or operator, such as denotes complex conjugation. We use this instead of to avoid the cluttering of symbols with the Keldysh superscripts . More generally, equation (A5) follows from a bosonic coherent-state path integral. For an phase space, however, an additional subtlety arises: the jump operators must be combined using the star product rather than ordinary multiplication. In the case of pure dephasing, which is what we consider in this case, this distinction is inconsequential, whereas for decay it leads to corrections. We will discuss this point in detail in future work [97]. A.2. Keldysh rotation Introducing the Keldysh rotation, the fields can be written in terms of “classical” and “quantum” components, given by In this basis, the action admits an expansion in powers of , where the linear term follows by expanding the coherent contribution (A2) in , so that with Here, the index labels internal and site degrees of freedom, for instance, for spins, one has with , where is the site index, and are for the three spin components The quadratic term originates from the dissipative contribution (A5). Writing one finds Substituting into (A5) gives Altogether, we have The kernel is determined by the dissipative term (A5) through the quadratic form . Higher-order terms in represent further quantum corrections. A.3. Semiclassical approximation and emergence of the stochastic equations The TWA consists of truncating the action (A13) at quadratic order in the quantum fields . The weight in the path integral now reads The quadratic term can be represented by a Gaussian identity. Introducing auxiliary complex Gaussian fields , the action becomes linear in , with which implies The functional integral over can thus be performed explicitly, where denotes a functional delta, i.e. . The generating functional therefore reduces to which is the functional representation of a stochastic process. The delta functional enforces the stochastic equations In the absence of dissipation (), equation (A21) reduces to the classical equations These equations can be written in Hamiltonian form as This Poisson structure corresponds to the classical limit of the quantum commutator (Dirac correspondence). For an open system, the kernel originates from the quadratic form in the dissipative action (A5). In a phase-space formulation, the quantum fields play the role of response fields and generate infinitesimal canonical transformations of the classical variables. Consequently, to leading order in , the difference of any functional between the two contour branches is generated by the structure so that and similarly for . Substituting this into yields a quadratic form in with kernel which is positive semidefinite by construction. After the Hubbard–Stratonovich decoupling (A15), the stochastic equations involve a Gaussian noise with covariance . Since factorises in terms of jump-operator brackets as in (A25), it is convenient to parameterise in terms of independent Gaussian noises associated with each jump operator As the Lindblad dissipator generates both deterministic dissipative drift terms and stochastic fluctuations, it is convenient to combine these contributions by introducing auxiliary fields where produces the deterministic dissipative contribution and represents the Gaussian noise associated with the jump process. Thus, the stochastic contribution can be expressed as in the more general equations of motion which reproduces the stochastic equations (9) in the main text. We note that the resulting dynamics can equivalently be interpreted as the semiclassical limit of a quantum Langevin equation [98]. A.4. Application to the long-range interacting, dissipative Rydberg gas For spin- degrees of freedom, operators are replaced by classical vectors with Poisson algebra where is the totally antisymmetric Levi–Civita tensor, with . In particular, the fundamental Poisson brackets reproduce the classical algebra, The classical Hamiltonian is with long-range interaction (7), and the classical dephasing jump operator is The coherent equations of motion follow from . For the driving term, one obtains while for the interaction term, the contributions are with given by equation (16). For the local dephasing, the stochastic contributions are where the real Gaussian noise satisfies equation (17). Appendix B: Numerical implementation B.1. Stochastic trajectories of spins In this appendix, we summarise the numerical implementation used to simulate the dissipative spin dynamics discussed in the main text. The dynamics of each site in our lattice is represented by a classical spin vector of fixed length , whose orientation evolves due to coherent precession and stochastic noise. The time evolution is governed by a Stratonovich stochastic differential equation of the form where the symbol ‘’ denotes the Stratonovich interpretation. The deterministic drift is given by which describes coherent precession in the effective field , where is the interaction field generated by the surrounding spins, while the noise vector is generating stochastic rotations about the -axis. The increments are independent Wiener processes satisfying Both the deterministic and stochastic contributions are constructed such that they are tangent to the Bloch sphere, i.e. for all . This property follows from the fact that both the coherent mean-field evolution and the Lindblad jump operators generate infinitesimal rotations of the local spin vector. As a consequence, the dynamics can be written entirely in terms of an effective angular velocity vector via where is the sum of a deterministic contribution and a stochastic contribution , where denotes Gaussian white noise in equation (17). Equation (B6) implies that each realisation of the dynamics is a trajectory on the sphere of radius : Thus, in exact arithmetic, every trajectory preserves the spin length . B.2. Rotation-based Stratonovich integrator In numerical implementations, conventional SDE solvers (such as Euler–Maruyama or Heun-type schemes) perform linear additive updates, rather than finite rotations on the sphere. As a consequence, the discrete update is not an orthogonal transformation and does not strictly preserve the constraint . Although the per-step deviation is small, it accumulates over the long simulation times considered in the main text and leads to a gradual drift of the spin norm. A simple remedy is to project the updated spin back onto the sphere after each step, so that which restores the constraint at the discrete level. However, this projection slightly modifies the stochastic flow. This motivates the use of structure-preserving integrators that evolve each timestep as a finite rotation on the Bloch sphere, thereby maintaining the constraint intrinsically. Equation (B6) generates a flow on the Lie group and therefore preserves exactly [99]. To maintain this structure at the discrete level, we evolve each timestep as a composition of finite rotations. We employ a symmetric Strang-type splitting scheme for the Stratonovich SDE in the spirit of [100], consisting of: (i) a deterministic half-step rotation, (ii) stochastic rotations generated by the Wiener increments, and (iii) a second deterministic half-step rotation evaluated at the updated configuration. Here, the drift in the second half-step is recomputed from the intermediate configuration obtained after the stochastic rotations. Letting denote the spin at time , the full timestep therefore corresponds to the composition which preserves up to machine precision. Since the stochastic rotations are applied symmetrically between two deterministic half-steps, the resulting composition constitutes a symmetric Strang splitting of the Stratonovich SDE on . As rotations on do not commute in general, the symmetric placement of the stochastic rotation between the two deterministic half-steps is important for maintaining consistency with Stratonovich calculus. For the deterministic part, the drift (using ) determines an instantaneous angular velocity and the spin is rotated by an angle about the axis . The stochastic part is implemented as exact finite rotations about the corresponding axes. For example, dephasing generates rotations about with angle where are Wiener increments. Appendix C: Temporal correlations To probe the temporal correlations underlying the relaxation dynamics, we analyse the connected two-time autocorrelation function, where denotes the time difference after the waiting time , i.e. . This is shown for decreasing dephasing strength for the polarised initial state (top row) as well as the Néel initial state (bottom row) for the spin chain in figure 8 and for the lattice in figure 9. Download figure: Standard image High-resolution image Download figure: Standard image High-resolution image The contour plots reveal a similar behaviour of the relaxation on the waiting time in both and , and for both initial states. For strong dephasing the autocorrelation decays rapidly in the time difference , and the contour lines are nearly vertical, indicating that the dynamics quickly becomes stationary and loses memory of the initial configuration. As the dephasing strength is reduced, correlations persist over increasingly longer time differences and a pronounced dependence on the waiting time emerges. This appears as a bending of the contour lines at intermediate , reflecting that the effective relaxation rate evolves during the approach to the stationary state. However, for a fixed waiting time , the decay in remains exponential throughout the evolution. The nontrivial connected two-time correlations arise because the interacting dynamics retains memory of earlier spin configurations during relaxation. Since the spin-flip rates depend on the local excitation pattern, the configuration at time influences the subsequent evolution, leading to a waiting-time dependent autocorrelation until the stationary state is reached. While the same qualitative behaviour is observed in one and two dimensions, the regime of pronounced waiting-time dependence extends over longer times in , indicating slower relaxation and more persistent temporal correlations.