Physical swap dynamics, shortcuts to relaxation and entropy production in dissipative Rydberg gases

TL;DR

This paper investigates how physical swap processes in dissipative Rydberg gases act as shortcuts to relaxation, influence entropy production, and alter the stationary state, revealing analogies with Monte Carlo methods and glassy dynamics.

Contribution

It introduces the concept of physical swaps as dynamical shortcuts in Rydberg gases and quantifies their impact on relaxation, stationary states, and entropy production.

Findings

Swaps accelerate relaxation and approach to stationarity.

Swaps increase entropy production and irreversibility.

Stationary states are significantly altered by swap dynamics.

Abstract

Dense Rydberg gases are out-of-equilibrium systems where strong density-density interactions give rise to effective kinetic constraints. They cause dynamic arrest associated with highly-constrained many-body configurations, leading to slow relaxation and glassy behavior. Multi-component Rydberg gases feature additional long-range interactions such as excitation-exchange. These are analogous to particle swaps used to artificially accelerate relaxation in simulations of atomistic models of classical glass formers. In Rydberg gases, however, swaps are real physical processes, which provide dynamical shortcuts to relaxation. They permit the accelerated approach to stationarity in experiment and at the same time have an impact on the non-equilibrium stationary state. In particular their interplay with radiative decay processes amplifies irreversibility of the dynamics, an effect which we quantify via the entropy production at stationarity. Our work highlights an intriguing analogy between real dynamical processes in Rydberg gases and artificial dynamics underlying advanced Monte Carlo methods. Moreover, it delivers a quantitative characterization of the dramatic effect swaps have on the structure and dynamics of their stationary state.

Equations

$$H_0=\sum _{k<l}\sum _{s=1}^p\left[V_{kl}^s{n}_s^{(k)}{n}_s^{(l)}+\sum _{s^{\prime }\ne s}{V}_{kl}^{s{s}^{\prime }}{n}_s^{(k)}{n}_{s^{\prime }}^{(l)}\right],$$

$$H_1=\sum _{k=1}^N\sum _{s=1}^p{\Omega }_s{\sigma }_{sx}^{(k)}+\sum _{k<l}\sum _{s<s^{\prime }}{E}_{kl}^{s{s}^{\prime }}\left({\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)}+{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)}\right).$$

$${\mathcal{L}}_0\rho =-i[H_0,\rho ]+\sum _{s=1}^p{\gamma }_s\sum _{k=1}^N\left(n_s^{(k)}\rho n_s^{(k)}-\frac{1}{2}\left\{n_s^{(k)},\rho \right\}\right),$$

$${\mathcal{L}}_1\rho =-i[H_1,\rho ]+\sum _{s=1}^p{\kappa }_s\sum _{k=1}^N\left({\sigma }_{-s}^{(k)}\rho {\sigma }_{+s}^{(k)}-\frac{1}{2}\left\{n_s^{(k)},\rho \right\}\right),$$

$$\underset{l\ne k}{⨂}{e}^{{\mathcal{L}}_{0,d}^{l}t}\left(\begin{array}{ccc}{\rho }_{pp}^{(k)}& e^{-\frac{({\gamma }_{p-1}+{\gamma }_p)t}{2}}{\rho }_{p(p-1)}^{(k)}& \cdots \\ e^{-\frac{({\gamma }_{p-1}+{\gamma }_p)t}{2}}{\rho }_{(p-1)p}^{(k)}& {\rho }_{(p-1)(p-1)}^{(k)}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right).$$

$${\partial }_t\mu =\mathcal{P}{\mathcal{L}}_1\mu +{\int }_0^{\infty }dt\mathcal{P}{\mathcal{L}}_1\mathcal{Q}{e}^{{\mathcal{L}}_{0}t}\mathcal{Q}{\mathcal{L}}_1\mu .$$

$$\mathcal{P}{\mathcal{L}}_1\mu =\sum _{s=1}^p{\kappa }_s\sum _{k=1}^N\left({\sigma }_{-s}^{(k)}\mu {\sigma }_{+s}^{(k)}-n_s^{(k)}\mu \right).$$

$$\mathcal{P}{\mathcal{L}}_1\mathcal{Q}{e}^{{\mathcal{L}}_{0}t}\mathcal{Q}{\mathcal{L}}_1\mu =\underset{D[\mu ,t]}{\underbrace{-\sum _s\sum _k{\Omega }_s^2\mathcal{P}[{\sigma }_{sx}^{(k)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{sx}^{(k)},\mu ]]}}\underset{E[\mu ,t]}{\underbrace{-\sum _{s<s^{\prime }}\sum _{k<l}{\left(E_{kl}^{s{s}^{\prime }}\right)}^2\mathcal{P}[{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)}+{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)}+{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},\mu ]].}}$$

$${\partial }_t\mu$$

$${\int }_0^{\infty }dtD[\mu ,t]=\sum _{s=1}^p\frac{4{\Omega }_s^2}{{\gamma }_s}\sum _k{\Gamma }_s^{(k)}\left[{\sigma }_{sx}^{(k)}\mu {\sigma }_{sx}^{(k)}-{\mathcal{I}}_s^{(k)}\mu \right],$$

$${\Gamma }_s^{(k)}=\frac{1}{1+{\left(\frac{2{\mathcal{V}}_s^k}{{\gamma }_s}\right)}^2}.$$

$${\int }_0^{\infty }dtE[\mu ,t]=\sum _{s<s^{\prime }}\frac{2({\mathcal{E}}_{\beta }^{s{s}^{\prime }}{)}^2}{{\gamma }_s+{\gamma }_s^{\prime }}\sum _{k<l}\frac{1}{|{\mathbf{r}}_k-{\mathbf{r}}_l{|}^{2\beta }}{\Gamma }_{s{s}^{\prime }}^{(kl)}\times$$

$$\times \left({\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)}\mu {\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)}+{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)}\mu {\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)}-{\mathcal{I}}_{s{s}^{\prime }}^{(k,l)}\mu \right)$$

$${{\rm Y}}_{s{s}^{\prime }}^{(kl)}=\frac{2({\mathcal{E}}_{\beta }^{s{s}^{\prime }}{)}^2}{{\gamma }_s+{\gamma }_s^{\prime }}\frac{1}{|{\mathbf{r}}_k-{\mathbf{r}}_l{|}^{2\beta }}{\Gamma }_{s{s}^{\prime }}^{(kl)}$$

$$\frac{1}{{\Gamma }_s^{(k)}}=1+{\left[\frac{2}{{\gamma }_s{a}^{\alpha }}\sum _{m\ne k}\left(\frac{C_{\alpha }^s{n}_s^{(m)}+{\sum }_{s^{\prime }\ne s}{C}_{\alpha }^{s{s}^{\prime }}{n}_{s^{\prime }}^{(m)}}{|{\hat{\mathbf{r}}}_k-{\hat{\mathbf{r}}}_m{|}^{\alpha }}\right)\right]}^2.$$

$${{\rm Y}}_{s{s}^{\prime }}^{(kl)}=\frac{\frac{2({\mathcal{E}}_{\beta }^{s{s}^{\prime }}{)}^2}{{\gamma }_s+{\gamma }_s^{\prime }}}{|{\mathbf{r}}_k-{\mathbf{r}}_l{|}^{2\beta }}\{1+[\frac{a^{-\alpha }}{{\gamma }_s+{\gamma }_s^{\prime }}\sum _{m\ne k,l}\left(\frac{1}{|{\hat{\mathbf{r}}}_k-{\hat{\mathbf{r}}}_m{|}^{\alpha }}-\frac{1}{|{\hat{\mathbf{r}}}_l-{\hat{\mathbf{r}}}_m{|}^{\alpha }}\right)$$

$${{\left(C_{\alpha }^{s^{\prime }}{n}_{s^{\prime }}^{(m)}+\sum _{s^{\prime \prime }\ne s^{\prime }}{C}_{\alpha }^{s^{\prime }{s}^{\prime \prime }}{n}_{s^{\prime \prime }}^{(m)}-C_{\alpha }^s{n}_s^{(m)}-\sum _{s^{\prime \prime }\ne s}{C}_{\alpha }^{s{s}^{\prime \prime }}{n}_{s^{\prime \prime }}^{(m)}\right)]}^2\}}^{-1}.$$

$${\partial }_t\mu$$

$${\Gamma }_{1,2}^{(k)}={\left(1+{\left[\sum _{m\ne k}\frac{R^6{n}_{1,2}^{(m)}+R_c^6{n}_{2,1}^{(m)}}{|{\hat{\mathbf{r}}}_k-{\hat{\mathbf{r}}}_m{|}^6}\right]}^2\right)}^{-1},$$

$${{\rm Y}}_{12}^{(kl)}=\frac{U}{|{\hat{\mathbf{r}}}_k-{\hat{\mathbf{r}}}_l{|}^{12}}{\left(1+{\left\{\frac{1}{4}\sum _{m\ne k,l}\left[\frac{\left(R_c^6-R^6\right)n_1^{(m)}+\left(R^6-R_c^6\right)n_2^{(m)}}{|{\hat{\mathbf{r}}}_k-{\hat{\mathbf{r}}}_m{|}^6}+\frac{\left(R_c^6-R^6\right)n_2^{(m)}+\left(R^6-R_c^6\right)n_1^{(m)}}{|{\hat{\mathbf{r}}}_l-{\hat{\mathbf{r}}}_m{|}^6}\right]\right\}}^2\right)}^{-1}.$$

$$\Pi =\frac{1}{2}\sum _{\mu ,\nu }\left({\Gamma }_{\mu \to \nu }{p}_{\mu }-{\Gamma }_{\nu \to \mu }{p}_{\nu }\right)\log \frac{{\Gamma }_{\mu \to \nu }{p}_{\mu }}{{\Gamma }_{\nu \to \mu }{p}_{\nu }}.$$

$$\Pi =\sum _{\mu }\sum _{k_0}\sum _{e=1,2}\left(\kappa p_{\mu e^{(k_0)}}(U,\kappa )-{\Gamma }_e^{(k_0)}[p_{\mu }(U,\kappa )-p_{\mu e^{(k_0)}}(U,\kappa )]\right)\log \frac{{\Gamma }_e^{(k_0)}+\kappa }{{\Gamma }_e^{(k_0)}}.$$

$${\Pi }_{U\to \infty }=\frac{2N\kappa }{3}\overline{\log (1+\kappa /{\Gamma }_e^{(k_0)})},$$

$$\mathcal{P}{\mathcal{L}}_1\mathcal{Q}{e}^{{\mathcal{L}}_{0}t}\mathcal{Q}{\mathcal{L}}_1\mu =$$

$$D[\mu ,t]$$

$${\sigma }_{1x}^{(k)}{e}^{{\mathcal{L}}_{0}t}\left({\sigma }_{1x}^{(k)}\mu \right)$$

$$\displaystyle-\Omega_{1}^{2}\left(\begin{array}[]{ccc}\ddots&\vdots&\vdots\\ \cdots&2\,e^{-\frac{1}{2}\gamma_{1}t}\cos\left({\mathcal{V}^{k}_{1}t}\right)\left[\rho_{11}^{(k)}-\rho_{00}^{(k)}\right]&0\\ \cdots&0&2\,e^{-\frac{1}{2}\gamma_{1}t}\cos\left(\mathcal{V}^{k}_{1}t\right)\left[\rho_{00}^{(k)}-\rho_{11}^{(k)}\right]\end{array}\right),$$

$$D[\mu ,t]=-\sum _{s=1}^p{\Omega }_s^2\sum _{k}2e^{-\frac{{\gamma }_s}{2}t}\cos \left({\mathcal{V}}_s^{k}t\right)\left[{\mathcal{I}}_s^{(k)}\mu -{\sigma }_{sx}^{(k)}\mu {\sigma }_{sx}^{(k)}\right],$$

$${\int }_0^{\infty }dtD[\mu ,t]=\sum _{s=1}^p\sum _k\frac{4{\Omega }_s^2/{\gamma }_s}{1+(2{\mathcal{V}}_s^k/{\gamma }_s{)}^2}\left[{\sigma }_{sx}^{(k)}\mu {\sigma }_{sx}^{(k)}-{\mathcal{I}}_s^{(k)}\mu \right].$$

$$E[\mu ,t]=-\sum _{s<s^{\prime }}\sum _{k<l}{\left(E_{kl}^{s{s}^{\prime }}\right)}^2\mathcal{P}([{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)},\mu ]]$$

$$+[{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},\mu ]]+[{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{s{s}^{\prime }}^{(k)}{\sigma }_{s^{\prime }s}^{(l)},\mu ]]$$

$$+[{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},e^{{\mathcal{L}}_{0}t}[{\sigma }_{s^{\prime }s}^{(k)}{\sigma }_{s{s}^{\prime }}^{(l)},\mu ]]).$$

Full text

I Introduction

Glasses are out of equilibrium systems typically produced by lowering the temperature of a liquid until the relaxation time exceeds experimentally accessible timescales Binder and Kob (2011); Berthier and Biroli (2011); Biroli and Garrahan (2013). Close to the glass transition, the physics is different from that of a high-temperature fluid, displaying dynamic heterogeneity (i.e. the coexistence of space-time regions with remarkably different timescales) arising from effective kinetic constraints Ritort and Sollich (2003); Chandler and Garrahan (2010). What makes the dynamics interesting, however, renders the simulation of glassy fluids challenging, as the computational equilibration times also increase dramatically. A method to reduce these times in Monte Carlo simulations is to artificially introduce swaps between particles of different species. This approach was proposed in Ref. Grigera and Parisi (2001), and later refined in e.g. Refs. Brumer and Reichman (2004); Fernandez et al. (2007). It enables the numerical equilibration of strongly arrested systems Ninarello et al. (2017), providing insight into lengthscales characterizing glassiness at low temperatures Biroli et al. (2008); Gutiérrez et al. (2015a), the jamming transition at high densities Berthier et al. (2016) or the configurational entropy of deeply supercooled states Berthier et al. (2017).

In the domain of cold atomic gases Rydberg atoms have emerged as a platform Gallagher (2005); Saffman et al. (2010); Löw et al. (2012) for studying collective dynamical behaviors Urvoy et al. (2015); Valado et al. (2016); Gutiérrez et al. (2017); Schempp et al. (2014); Gutiérrez et al. (2016); Letscher et al. (2017) which are closely linked to that of glasses Lesanovsky and Garrahan (2013, 2014); Gutiérrez et al. (2015b). The physics of interacting Rydberg atoms is governed by blockade effects Jaksch et al. (2000); Lukin et al. (2001) analogous to the excluded volume interactions underlying the dynamics of glass-forming liquids and jammed systems of hard objects Binder and Kob (2011); Berthier and Biroli (2011); Sanders et al. (2015). Also, in both Rydberg gases Pérez-Espigares et al. (2018) and glassy soft-matter models Garrahan et al. (2007, 2009); Hedges et al. (2009); Speck et al. (2012) dynamical phase transitions have been shown to be at the heart of the observed collective dynamical effects.

In this paper we investigate the impact of excitation swaps on the relaxation timescales and on the structure and dynamics of the stationary state of Rydberg gases. Swaps are caused by resonant dipole-dipole interactions Saffman et al. (2010); de Léséleuc et al. (2017) which effectuate state exchanges over long distances. This is strongly reminiscent of the above-mentioned advanced Monte Carlo approaches in which particle swaps accelerate the approach of highly arrested systems to stationarity. In Rydberg gases swaps are physical and their strength and range can be controlled. We show that their inclusion indeed shortens the characteristic timescales by several orders of magnitude, which allows one to experimentally introduce shortcuts to relaxation. We also show that swaps drive the system further away from equilibrium conditions, in particular in conjunction with radiative decay processes. The resulting irreversibility of the dynamics is quantified via the entropy production. We outline how this quantity may be experimentally accessed and deduce the existence of an irreversibility bound for large swap rate.

II Multi-component Rydberg gases

II.1 Description of the model

We consider a system of $N$ atoms. The ground state $\left|0\right>$ of each atom is resonantly coupled by a laser field to the Rydberg states $\left|1\right>,\left|2\right>,\ldots,\left|p\right>$ (with energies $E_{0}<E_{1}<\cdots<E_{p}$). Such systems are called multi-component Rydberg gases, as they include different kinds of Rydberg excitations (i.e. atoms excited to different Rydberg levels), see Fig 1 (a) for an illustration of the $p=2$ case. The Hamiltonian $H=H_{0}+H_{1}$ includes the density-density interactions in $H_{0}$ and the exchange interactions and the driving in $H_{1}$. Atoms in the Rydberg states $|s\rangle$ and $|s^{\prime}\rangle$ at positions ${\bf r}_{k}$ and ${\bf r}_{l}$ interact through a power-law potential $V_{kl}^{ss^{\prime}}=C_{\alpha}^{ss^{\prime}}/|{\bf r}_{k}-{\bf r}_{l}|^{\alpha}$ with exponent $\alpha$. For simplicity we denote the intra-level interactions by $V_{kl}^{s}$ instead of $V_{kl}^{ss}$. We thus have

[TABLE]

where $n_{s}^{(k)}=\left|s\right>_{k}\!\left<s\right|$, which is diagonal in the eigenbasis $\left|00\ldots 0\right>$, $\left|00\ldots 1\right>$, etc. The off-diagonal part is

[TABLE]

Here, $\sigma_{sx}^{(k)}=|0\rangle_{k}\langle s|+|s\rangle_{k}\langle 0|$ is the Pauli matrix for a given transition and $\sigma_{ss^{\prime}}^{(k)}=|s\rangle_{k}\langle s^{\prime}|$. The first term in Eq. (2) includes the driving in the rotating wave approximation, where $\Omega_{s}$ is the Rabi frequency of the coupling between $\left|s\right>$ and $\left|0\right>$. The second term contains the exchange interactions, with a power-law dependence on the distance between atoms $E_{kl}^{ss^{\prime}}=\mathcal{E}_{\beta}^{ss^{\prime}}/|{\bf r}_{k}-{\bf r}_{l}|^{\beta}$, with exponent $\beta$.

The full dynamics is governed by a quantum master equation in Lindblad form, which includes as incoherent processes a dephasing term (arising from the laser linewidth, thermal effects, etc. Urvoy et al. (2015); Valado et al. (2016); Gutiérrez et al. (2017)) and a spontaneous decay term for each Rydberg level. This master equation can be written as $\partial_{t}\rho=\mathcal{L}\rho=\mathcal{L}_{0}\rho+\mathcal{L}_{1}\rho$. Here, the diagonal Liouvillian $\mathcal{L}_{0}$ contains not only the interaction Hamiltonian $H_{0}$ but also a dissipative term that governs the dephasing process, and is defined as

[TABLE]

where $\gamma_{s}$ is the dephasing rate of $\left|s\right>$ w.r.t. $\left|0\right>$. The superoperator $\mathcal{L}_{0}$ therefore consists of a Hamiltonian part and a dissipator whose individual terms commute. On the other hand, the Liouvillian $\mathcal{L}_{1}$ includes the evolution due to the off-diagonal part of the Hamiltonian $H_{1}$ and a dissipator that accounts for spontaneous decay in the system. It takes the form

[TABLE]

where $\sigma_{-s}^{(k)}=|0\rangle_{k}\langle s|$, $\sigma_{+s}^{(k)}={\sigma_{-s}^{(k)}}^{{\dagger}}$, and $\kappa_{s}$ is the decay rate of excited level $\left|s\right>$.

II.2 Perturbative analysis based on projection operators

When different dynamical processes evolve on vastly different timescales, one sometimes can adiabatically eliminate fast-evolving degrees of freedom, and focus on the evolution of the remaining (slow) ones. In this section, we derive the effective dynamics in the limit of strong dissipation and strong (density-density) interactions, where the non-diagonal terms of the full Liouvillian $\mathcal{L}$, which are contained in $\mathcal{L}_{1}$ can be treated perturbatively: $\Omega_{s},\kappa_{s},\mathcal{E}_{\beta}^{ss^{\prime}}\ll\gamma_{s},C_{\alpha}^{ss^{\prime}}$. In this strongly-dissipative limit, coherent superpositions of local atomic states dephase exponentially fast, allowing for a description of the dynamics that only includes (slow-evolving) diagonal elements of the density matrix on timescales larger than $\gamma_{s}^{-1}$ (for the atomic level with the smallest dephasing rate $\gamma_{s}$ in the system). An analogous approach has been previously used to explore multi-component Rydberg gases in the absence of exchange processes and decay Gutiérrez et al. (2016), and also to study one-component systems Lesanovsky and Garrahan (2013); Cai and Barthel (2013); Petrosyan et al. (2013); Marcuzzi et al. (2014); Schönleber et al. (2014); Hoening et al. (2014).

We start by working out the effect of dephasing on the dynamics. Using the notation, $\mathcal{L}_{0,d}^{k}=\sum_{s=1}^{p}\left[\gamma_{s}\left(n_{s}^{(k)}\rho\,n_{s}^{(k)}-\frac{1}{2}\left\{n_{s}^{(k)},\rho\right\}\right)\right]$, the action of the Liouvillian $\mathcal{L}_{0}$ leads to the evolution $e^{\mathcal{L}_{0}t}\rho=e^{-iH_{0}t}\,[\bigotimes_{k}e^{\mathcal{L}_{0,d}^{k}t}\rho]\,e^{iH_{0}t}$ where the term in brackets gives

[TABLE]

Here, $\rho_{mn}^{(k)}$ are the $p^{N-1}\times p^{N-1}$ matrices defined by $\rho_{mn}^{(k)}={}_{k}{\left<n\right|}\rho\left|m\right>_{k}$, using as basis states $\left|0\right>_{k}$, $\left|1\right>_{k}$, $\left|2\right>_{k}$, $\ldots,$ $\left|p\right>_{k}$. The off-diagonal entries of the density matrix are seen to decay exponentially, a fact that is not altered by the action of the coherent dynamics $e^{-iH_{0}t}$. Therefore, the evolution under $\mathcal{L}_{0}$ becomes, at times considerably larger than the inverse of the dephasing rates, a projector $\mathcal{P}$ on the diagonal of $\rho$, $\mathcal{P}\rho=\lim_{t\to\infty}e^{\mathcal{L}_{0}t}\rho=\textrm{diag}(\rho)$, as happens in the case of just one Rydberg level Lesanovsky and Garrahan (2013). The removal of all coherences leads to a diagonal density matrix, where each classically accessible configuration (e.g. $\left|0010203\cdots\right>$) is given a certain probability of occurrence.

By using the operator $\mathcal{P}$ (which projects on the slow-evolving diagonal) and its complement $\mathcal{Q}=1-\mathcal{P}$ (which projects on the rest of the Liouville space, where rapidly-decaying quantum coherences reside), we can formulate the effective evolution equation for the diagonal density matrix $\mu=\mathcal{P}\rho$, which describes the slow evolution of the excitation dynamics. To second order in $\mathcal{L}_{1}$, the general Nakajima-Zwanzig expression Breuer and Petruccione (2002) is given by

[TABLE]

See Marcuzzi et al. (2014) for a detailed treatment of the perturbative expansion in the context of Rydberg gases. In the present case, the first term simplifies to

[TABLE]

As $H_{1}$ takes the elements of $\mu$ off the main diagonal, its action vanishes under the projection $\mathcal{P}$. On the other hand, the integrand in Eq. (6) reduces to

[TABLE]

A detailed derivation can be found in Appendix A. For simplicity the first term in Eq. (8), arising from the driving, is denoted as $D[\mu,t]$, and the second, related to exchange processes, as $E[\mu,t]$. Additionally, there is the term that accounts for the decay processes, Eq. (7). Putting it all together, Eq. (6) can be rewritten as

[TABLE]

II.3 Effective dynamics: driving

The laser driving term $\int_{0}^{\infty}dt\,D[\mu,t]$ in Eq. (9), gives rise to the excitations and de-excitations in the system. While this has already been calculated before Gutiérrez et al. (2016), we include the details in Appendix B for completeness. It yields a classical stochastic generator

[TABLE]

where the projection operator $\mathcal{I}_{s}^{(k)}=n_{s}^{(k)}+|0\rangle_{k}\langle 0|$ cancels all the elements in $\mu$ that do not correspond to the ground state or $|s\rangle$ at site $k$. The rates for a transition $|0\rangle\to|s\rangle$ or $|s\rangle\to|0\rangle$ at site $k$ are

[TABLE]

The increase in the interaction energy that is required for the excitation of atom $k$ to level $|s\rangle$ has been denoted as $\mathcal{V}^{k}_{s}=\sum_{m}\left[V^{s}_{km}n_{s}^{(m)}+\sum_{s^{\prime}\neq s}V^{ss^{\prime}}_{km}n_{s^{\prime}}^{(m)}\right]$ for convenience.

II.4 Effective dynamics: exchange

The term $\int_{0}^{\infty}\!dt\,E[\mu,t]$ in the effective dynamics, Eq. (9), which is due to the exchange interactions in $H_{1}$, is analyzed in Appendix C. It can also be written as classical stochastic generator

[TABLE]

where $\mathcal{I}_{ss^{\prime}}^{(k,l)}=n_{s}^{(k)}n_{s^{\prime}}^{(l)}+n_{s^{\prime}}^{(k)}n_{s}^{(l)}$. The rate for an exchange $s\leftrightarrow s^{\prime}$ between sites $k$ and $l$, i.e. the excitation-swap rate, is

[TABLE]

with $\Gamma_{ss^{\prime}}^{(kl)}=\left({1+\left[\frac{1}{\gamma_{s}\!+\!\gamma_{s}^{\prime}}\left(\mathcal{V}^{k}_{s^{\prime}}\!-\!\mathcal{V}^{k}_{s}\!+\!\mathcal{V}^{l}_{s}\!-\!\mathcal{V}^{l}_{s^{\prime}}\right)\right]^{2}}\right)^{-1}\!.$

II.5 Transition rates

We conclude that the equation governing the effective dynamics, Eq. (9), is a classical master equation with three different processes (excitations, swaps and decay) occuring with different rates, see Fig 1 (a). Specifically as follows.

(i) Excitations and de-excitations between $|0\rangle$ and $|s\rangle$, which, according to Eqs. (10, 11), occur with an inverse timescale $\frac{4\Omega_{s}^{2}}{\gamma_{s}}$, and with a configuration-dependent rate

[TABLE]

We have made explicit the power-law dependence on the distance between atoms by giving the positions in a lattice in terms of reduced position vectors $\hat{\bf r}_{k}={\bf r}_{k}/a$, where $a$ is the lattice constant.

(ii) Excitation swaps $|s\rangle_{k}|s^{\prime}\rangle_{l}\leftrightarrow|s^{\prime}\rangle_{k}|s\rangle_{l}$ (12). The power-law suppression $\frac{1}{|{\bf r}_{k}-{\bf r}_{l}|^{2\beta}}$ in Eq. (12) severely limits the possibility of swaps between highly distant atoms. The configuration-dependent part of the swap rates in (13), $\Gamma_{ss^{\prime}}^{(kl)}$, contains in the denominator the difference between the energy associated with the presence of an excitation $|s\rangle$ at site $k$ and an excitation $|s^{\prime}\rangle$ at site $l$ and that resulting from an exchange of the excited levels of $k$ and $l$ (i.e. the energy before and after a swap). If we make explicit the interactions, the rates appear as

[TABLE]

(iii) Decay processes leading from a given excited level $|s\rangle$ to the ground state $|0\rangle$ with a constant rate $\kappa_{s}$, which are accounted for by the last term in Eq. (9).

II.6 Simplifying parameter choices

The transition rates in Eqs. (14, 15) contain many different parameters associated to processes and interactions involving different Rydberg levels. Due to the high dimensionality of the parameter space, we need to make some simplifying parameter choices in order to explore the dynamics that emerges from excitation, swap and decay processes in numerical simulations. The first simplifying assumption that we make is that all transitions are driven with the same Rabi frequency $\Omega_{1}=\Omega_{2}=\cdots=\Omega_{p}\equiv\Omega$ and dephase with the same rate $\gamma_{1}=\gamma_{2}=\cdots=\gamma_{p}\equiv\gamma$. Moreover, the decay rates are considered to be equal too, $\kappa_{1}=\kappa_{2}=\cdots=\kappa_{p}\equiv\kappa$. Following Ref. Gutiérrez et al. (2016), we also define an intra-level interaction parameter $R_{s}=a^{-1}[2C^{s}_{\alpha}/\gamma]^{1/\alpha}$ (for interactions between atoms in the same excited level $s$), and an inter-level interaction parameter $R_{ss^{\prime}}=a^{-1}[2C^{ss^{\prime}}_{\alpha}/\gamma]^{1/\alpha}$ (for interactions between atoms in different levels, $s$ and $s^{\prime}$). These are microscopic (reduced) lengthscales (normalized by the lattice constant) that parametrize the effective range of intra-level and inter-level interactions, respectively, and they determine the dynamical regime of the system Pérez-Espigares et al. (2018). Moreover, we focus on van der Waals interactions, so the exponents in (14, 15) satisfy $\alpha=\beta=6$.

II.7 Full equations of motion for a two-component Rydberg gas

For simplicity, our numerical results focus on the case of two components (i.e. two Rydberg levels), $p=2$, with a single intra-level interaction parameter $R$ ($R_{1}=R_{2}\equiv R$, $C_{6}^{1}=C_{6}^{2}$), and an inter-level interaction parameter which we denote by $R_{c}$ ($R_{12}\equiv R_{c}$) for convenience. A simpler version of this system without excitation swaps or decay was explored in Ref. Gutiérrez et al. (2016).

The effective dynamics is governed by the following classical master equation

[TABLE]

The time variable can be rescaled by $4\Omega^{2}/\gamma$, which is the common factor appearing in all excitation and de-excitation processes (i.e. the first two terms above). Then the factor $(\mathcal{E}_{\beta}^{12})^{2}/\gamma$ in Eq. (15) is rescaled to $(\mathcal{E}_{\beta}^{12})^{2}/4\Omega^{2}\equiv U$, which we take as the overall swap strength.

The resulting Markovian dynamics proceeds along classical (diagonal) states, e.g. $|012100\cdots\rangle$, and includes the following processes [see Fig. 1 (a)]: (i) driven (de-)excitations ($|0\rangle\leftrightarrow|1\rangle$, $|0\rangle\leftrightarrow|2\rangle$) with rates [see (14)]

[TABLE]

(ii) excitation swaps $|12\rangle\leftrightarrow|21\rangle$ with rates

[TABLE]

[which is the form adopted by (15) after introducing the simplifications of this and the previous section] and (iii) spontaneous decay from $|1\rangle$ or $|2\rangle$ to $|0\rangle$, with (constant) rate $\kappa$. Similarly derived effective dynamics involving a single Rydberg level (thus excluding excitation swaps) were recently studied experimentally, yielding an excellent agreement with theoretical predictions Urvoy et al. (2015); Valado et al. (2016); Gutiérrez et al. (2017).

According to Eq. (17), transitions between the ground state and each excited state – type (i) processes – slow down in the vicinity of previously excited atoms. This excitation blockade leads to dynamic bottlenecks Lesanovsky and Garrahan (2013); Valado et al. (2016); Gutiérrez et al. (2016), analogous to kinetic constraints in glassy systems Ritort and Sollich (2003); Chandler and Garrahan (2010). Following the rationale behind the swap Monte Carlo studies discussed in the Introduction, we expect that type (ii) processes may unblock otherwise highly constrained configurations, as illustrated in Fig. 1 (a). In the following section we show compelling numerical results that confirm these expectations.

III Impact of excitation swaps on dynamics

III.1 Accelerating effect of swaps

To gain insight into the role of excitation swaps, we study a chain of $N$ atoms with periodic boundary conditions via continuous-time Monte Carlo (CTMC) Bortz et al. (1975); Newman and Barkema (1999). Atoms $1\leq k\leq N/2$ are initially in state $|1\rangle$ while atoms $N/2<k\leq N$ are initially in state $|2\rangle$. For intra-level interactions appreciably stronger than the inter-level interactions, i.e. if $R$ is considerably larger than $R_{c}$, such an initial configuration is highly arrested, see Eq.  (17). This means that, in the absence of exchange processes ($U=0$), it takes a long time for the system to ‘forget’ its initial configuration, see Fig. 1 (b). When exchange interactions are included (for swap strength $U>0$), the relaxation time is drastically reduced — see a representative trajectory for $U=1$ in Fig. 1 (c). Swaps start being effective at the domain boundaries, due to the proximity of different types of excitations [notice the power-law decay of the swap rates with distance in (18)]. Then they propagate through the chain, eventually erasing all trace of the initial configuration. By contrast, trajectories of the same duration with $U=0$ hardly show any departure from the initial configuration.

To quantify the acceleration caused by swaps, we use the persistence function $P(t)$, defined as the fraction of atoms that have not changed state since $t=0$. It is displayed in Fig. 2 (a) for several values of the intra-level interaction parameter $R>1$ and the swap rate $U$ (the inter-level interaction parameter is $R_{c}=1$, and for the moment decay is excluded, $\kappa=0$). For $R=4$ the initial state is highly arrested, and swaps dramatically shorten the time that the system needs to move away from it. This effect is attenuated when the dynamical arrest is reduced, as in the curve for $R=1.5$, which is only moderately larger than $R_{c}$. By inspecting these and many other persistence curves with different sets of parameters, we conclude that the accelerating role of swaps for a given swap strength $U$ is dramatic when the interaction parameters, $R$ and $R_{c}$, are sufficiently large so the dynamics is glassy Pérez-Espigares et al. (2018), and one is considerably larger than the other. In fact, if the intra-level parameter $R$ and the inter-level parameter $R_{c}$ are very similar, an excitation swap is not more effective than a swap of particles with very similar diameters in a supercooled liquid: the swap would not make the local relaxation significantly easier, as the initial and final configuration would be similarly constrained. In this sense, the role that the interaction parameters $R$ and $R_{c}$ play in the swap mechanism is analogous to that of particle diameters in glassy liquids Berthier et al. (2019). The effect of swaps is enhanced in Rydberg gases with markedly different interaction parameters, as in highly polydisperse atomistic mixtures Berthier et al. (2016) or in coarse-grained glassy models with wide local softness distributions Gutierrez et al. (2019).

In Fig. 2 (a), for either value of the intra-level interaction parameter $R$, the decay of $P(t)$ is approximately logarithmic for large swap strength $U$. This dependence can by understood from the statistics of local persistence times $\tau$, i.e. the time it takes an atom to undergo a state change. The distribution of $\tau$ in logarithmic scale $\pi[\log_{10}(\tau)]$ is highly peaked for $U=0$ but becomes wider and flatter for larger $U$, see Fig. 2 (b). As the persistence is $P(t)=\int_{t}^{\infty}d\tau\,p[\tau]=\int_{t}^{\infty}d\tau\,\tau^{-1}\,\pi[\log_{10}(\tau)]$, for $\log_{10}(\tau)$ uniformly distributed across a $\tau$ interval, we obtain $P(t)\sim-\log_{10}(t)$ for $t$ within that interval [for shorter (longer) times $P(t)=1$ ([math])]. This suggests that the logarithmic decay is connected to the broad distribution of local persistence times caused by excitation swaps.

Further information on the relaxation dynamics can be obtained from the global persistence time $t_{p}$, defined as the time it takes for the persistence function $P(t)$ to reach zero, which corresponds to the longest local persistence time $\tau$ in the system. These are shown in Fig. 2 (c) using the same color code as in panel (a). A non-zero $U$ is seen to cause a departure from the persistence line without swaps, $U=0$, if a threshold value of the intra-level interactions $R$ is exceeded. This threshold value becomes smaller with increasing swap strength $U$, which means that if $R$ is sufficiently large, a very small swap strength $U$ has a strong effect, while a larger $U$ is needed to accelerate a less constrained dynamics. Above the threshold, swaps lead to shorter persistence times and a slower (almost linear) growth with $R$. Reductions of timescales of up to four orders of magnitude are observed, and larger reductions could be in principle achieved (although at a high computational cost) by making the dynamics more constrained or increasing the swap strength or both.

III.2 Interplay of swap and decay processes

The presence of spontaneous decay, which annihilates excitations irrespective of the state of their neighbours, imposes a further timescale which, despite its non-collective nature, is relevant for the effectiveness of swaps. In fact, swaps can only be effective if they take place on timescales that are shorter than the typical decay time (otherwise, constrained configurations automatically relax by the effect of decay). This is illustrated in Fig. 2 (c), where persistence times $t_{p}$ for a decay rate $\kappa=10^{-6}$ are shown in gray, and can be compared to the persistence times in a situation without decay ($\kappa=0$). When the typical decay time $\kappa^{-1}$ is longer than $t_{p}$ for $\kappa=0$, the persistence time is not affected by decay processes. This is why the role of decay in the figure (where $\kappa=10^{-6}$) is negligible for vanishing swap strength $U=0$ when the intra-level interaction parameter satisfies $R\lessapprox 3$, and also for $U=0.01$ when $R\lessapprox 4$ and for $U=1$ throughout the range of $R$ that has been considered. If, as $R$ is increased, the persistence time in the absence of decay starts to exceed $\kappa^{-1}$, then the actual $t_{p}$ becomes independent of $R$ (and essentially equal to $\kappa^{-1}$) when the decay is present [gray symbols in Fig. 2 (c)]. From this we conclude that excitation swaps accelerate the dynamics as long as the typical decay time is longer than the reduced persistence times caused by the swaps. Otherwise the relaxation does not depend on collective effects as it is driven by spontaneous decay.

IV Irreversibility and entropy production

IV.1 Irreversible dynamics in the presence of spontaneous decay

When the decay rate $\kappa$ is zero, the (uncorrelated) stationary state is such that each atom can be in any energy level $|0\rangle,\,|1\rangle,\,|2\rangle$ with the same probability Gutiérrez et al. (2016). As the rate between any two configurations connected by a transition is the same in either direction, see Eqs. (17) and (18), the dynamics trivially satisfies detailed balance, and therefore is reversible Kelly (2011). This is expected to change for non-vanishing decay ($\kappa>0$), as then the excitation of atom $k$ from the ground state $|0\rangle$ to e.g. state $|1\rangle$ occurs with a rate $\Gamma_{1}^{(k)}$ [Eq. (17)], while a larger rate $\Gamma_{1}^{(k)}+\kappa$ is associated with a de-excitation $|1\rangle$ to $|0\rangle$, and the probabilities of the many-body configurations are different.

To investigate the violation of detailed balance in the presence of decay we employ the Kolmogorov criterion. It states that a stationary Markov process is reversible if and only if the product of its transition rates along any cycle in configuration space is the same whether the path is traced in one or the other direction Zia and Schmittmann (2007); Kelly (2011). In Fig. 3 (a) we show a simple cycle in configuration space along which two atoms are successively de-excited and excited back to their initial state. The pictorial notation for the excitation rates depicts the states of the first and second nearest neighbors to either side of the atom undergoing a transition (atoms immediately to the left and right of those shown are assumed to be excited). While the Kolmogorov criterion is seen to be satisfied for $\kappa=0$, in the presence of spontaneous decay, $\kappa>0$, the product of the transition rates depends on the direction taken along the cycle, see Fig. 3 (a). Thus, detailed balance is violated for $\kappa>0$: the dynamics becomes irreversible due to the presence of decay processes.

IV.2 Entropy production rate

The degree of irreversibility of a stochastic dynamics is captured by the entropy production rate Schnakenberg (1976); Zia and Schmittmann (2007),

[TABLE]

Here $\mu$ and $\nu$ are two configurations (e.g. $|0121\cdots\rangle$ and $|1121\cdots\rangle$), $\Gamma_{\mu\to\nu}$ is the rate of the transition $\mu\to\nu$, and $p_{\mu}$ is the probability of configuration $\mu$ in the stationary state. The entropy production rate $\Pi$ is zero for reversible dynamics and positive otherwise. Eq. (19) can be simplified by removing the probabilities $p_{\mu}$ from the logarithm, as they account for the internal entropy production rate, which vanishes in the stationary state Schnakenberg (1976). This yields a somewhat simpler expression: $\Pi=\frac{1}{2}\sum_{\mu,\nu}\left(\Gamma_{\mu\to\nu}\,p_{\mu}-\Gamma_{\nu\to\mu}\,p_{\nu}\right)\log(\Gamma_{\mu\to\nu}/\Gamma_{\nu\to\mu})$. The entropy production rate is accessible experimentally by continuously monitoring the states of each atom in a small cold-atomic lattice system or tweezer array Schauss et al. (2012); Saffman (2016); Labuhn et al. (2016); Jau et al. (2016); Bernien et al. (2017); Kim et al. (2018); Shi (2018), and inferring the probabilities $p_{\mu}$, and also the rates $\Gamma_{\mu\to\nu}$ via the so-called empirical currents, which count the number of transitions in a given time window.

The transition rates $\Gamma_{\mu\to\nu}$ in this context are: (A) $\Gamma_{1,2}^{(k)}$ if $\nu$ has an excited atom that is in the ground state in $\mu$, (B) $\Gamma_{1,2}^{(k)}+\kappa$ if $\nu$ has a ground state atom that is excited in $\mu$, and (C) $\Upsilon_{12}^{(kl)}$ if $\mu$ and $\nu$ are connected by an excitation swap. Case (C) does not contribute directly to the entropy production rate, as $\log(\Upsilon_{12}^{(kl)}/\Upsilon_{21}^{(kl)})=0$ ($\Upsilon_{12}^{(kl)}=\Upsilon_{21}^{(kl)}$). But excitation swaps do play a crucial role in a different way, as they modify the stationary distribution $p_{\mu}$ (see below). Eq. (19) becomes

[TABLE]

The terms in brackets have been rearranged slightly for convenience. Here, $\mu$ is a given configuration and $\mu e^{(k_{0})}$ is a configuration that is reached from $\mu$ by the excitation of site $k_{0}$ to level $|e\rangle$. (Strictly speaking, $\mu$ must have at least one atom in the ground state $|0\rangle$ or the corresponding term will be zero.) This transition occurs at a rate which we denote $\Gamma_{e}^{(k_{0})}$. As the stationary probability distribution depends on the values of $U$ and $\kappa$, we explicitly write this dependence in $p_{\mu}(U,\kappa)$ and $p_{\mu e^{(k_{0})}}(U,\kappa)$. In the sums, $k_{0}$ goes over all atoms in the ground state in $\mu$, and $e=1,2$ denotes excitations to state $|1\rangle$ or $|2\rangle$. For $\kappa=0$, the logarithm vanishes, and so does $\Pi$; otherwise the logarithm is positive, and $p_{\mu}(U,\kappa)\geq p_{\mu e^{(k_{0})}}(U,\kappa)$.

IV.3 The role of swaps in the entropy production

In Fig. 3 (b), we show the entropy production as a function of the decay rate $\kappa$ for different values of the swap strength $U$. To compute the entropy production, we obtain numerically the stationary distribution $p_{\mu}$ by running a CTMC simulation for a sufficiently long time so that the system has visited every configuration at least $10^{4}$ times. Then we simply apply Eq. (20), as the rates for each $\mu$ are defined in Eq. (17). As estimating $p_{\mu}$ for each configuration $\mu$ is computationally demanding, we base our results on a small ring of $N=4$ atoms.

As expected, the entropy production $\Pi$ is zero for $\kappa=0$ (reversible dynamics) and increases for larger $\kappa$. More interestingly, for $U>0$ the existence of excitation swaps leads to a significant increase in $\Pi$, which becomes more pronounced as $U$ gets larger. Somewhat counterintuitively, this is connected to the appearance of “reversible” shortcuts [in the sense that their associated probability currents, $\Gamma_{\mu\to\nu}\,p_{\mu}-\Gamma_{\nu\to\mu}\,p_{\nu}$, vanish Zia and Schmittmann (2007), and they do not contribute directly to the entropy production, as pointed out above]. These shortcuts move the system between configurations with several excitations that are not connected by a single transition when swaps are absent.

While swaps by themselves do not produce entropy, they open alternative routes to highly inaccessible states. For example, a state with many excitations may be practically unreachable from the highly constrained states that have just one excitation less, given the very small rates $\Gamma_{1,2}^{(k)}$, but it may be reached through an appropriate sequence of excitation swaps. Exchange processes thus lead to a stationary distribution where on average the populations $p_{\mu e^{(k_{0})}}(U,\kappa)$ are larger and closer to their associated $p_{\mu}(U,\kappa)$. The basic mechanism is illustrated by a simple model in Appendix D. This means that, as the swap strength $U$ grows, the positive term in brackets in Eq. (20) becomes larger and the absolute value of the negative term becomes smaller. As the stationary probabilities are the only variables in Eq.  (20) that are affected by swaps, the entropy production thus grows with $U$. To give numerical support to this observation, we show in Fig. 3 (c) the average number of excitations in the stationary state $\langle n\rangle=\sum_{\mu}p_{\mu}(U,\kappa)\,n_{\mu}$ as a function of the decay rate $\kappa$ for different values of the swap strength $U$. Here the sum goes over all configurations and $n_{\mu}$ gives the number of sites that are in states $|1\rangle$ or $|2\rangle$ in $\mu$. Spontaneous decay suppresses the probability of configurations with a large number of excitations, but swaps counteract this effect to a significant extent 111The value that the average number of excitations $\left<n\right>$ takes for $\kappa=0$ in Fig. 3 (c) is the value corresponding to a uniform stationary distribution. This is $\sum_{n=0}^{N}n\,{N\choose n}\,2^{n}/3^{N}=2N/3$ in a system of $N$ atoms, as can be easily seen from the binomial identity $x\frac{d}{dx}(1+x)^{N}=\sum_{n=0}^{N}n\,{N\choose n}\,x^{n}$..

To conclude, we derive an upper bound for the entropy production for large swap strength $U$. In that limit, $p_{\mu}\approx p_{\mu e^{(k_{0})}}\approx 1/3^{N}$, i.e. the probability distribution approaches that of a completely mixed (infinite temperature) stationary state, and the term in brackets in Eq.  (20) reaches its maximum, $\kappa/3^{N}$. The entropy production then approaches

[TABLE]

which includes the average value of $\log(1+\kappa/\Gamma_{e}^{(k_{0})})$ taken over all $2N3^{N-1}$ excitation processes. This bound increases with the decay rate $\kappa$ and can be easily computed for small systems [see the red dotted line in Fig. 3 (b)].

V Outlook

Excitation swaps dramatically shorten the relaxation timescales of dissipative Rydberg gases. This is what particle swaps in glassy soft-matter models also do, but there is one crucial difference: while swaps in those systems are just a numerical trick to reduce computational times, in atomic ensembles they are part of the physical dynamics. The accelerating effect of swaps in Rydberg gases is robust and does not require any fine-tuning of parameters, and therefore should be observable in modern cold atomic settings. In fact, the observation of this effect only requires sufficiently strong interactions (which give rise to relaxation timescales much longer than the excitation timescale of isolated atoms), and sufficiently different interaction strengths for different excited levels. (Analogously, particle swaps are effective in liquids with glassy dynamics and if the model includes particles with sufficiently different diameters.) Additionally, there is the somewhat trivial requirement that decay processes must not take place predominantly on timescales that are shorter than the effect of swaps.

We have also shown that, apart from accelerating the dynamics, excitation swaps increase the entropy production rate, effectively moving the system further away from equilibrium conditions. This increased irreversibility is due to a reconfiguration of the stationary state distribution caused by the appearance of shortcuts in configuration space that link highly constrained configurations. The entropy production rate is in principle accessible in Rydberg quantum simulators and may enable a new handle for the characterization of non-equilibrium states in open quantum systems.

Acknowledgements.

We are grateful to M. Marcuzzi and P. Rotondo for insightful discussions. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. 335266 (ESCQUMA) and the EPSRC Grants No. EP/M014266/1, EP/R04421X/1, EP/R04340X/1. RG acknowledges the funding received from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 703683. IL gratefully acknowledges funding through the Royal Society Wolfson Research Merit Award. We are also grateful for access to the University of Nottingham High Performance Computing Facility, and for the computing resources and related technical support provided by CRESCO/ENEAGRID High Performance Computing infrastructure (funded by ENEA, the Italian National Agency for New Technologies, Energy and Sustainable Economic Development and by Italian and European research programmes) and its staff et al. (2014).

Appendix A Derivation of Eq. (8)

The integral in Eq. (6) contains all the (de-)excitation processes as well as the excitation swaps in the effective dynamics. The decay processes appear in the first term $\mathcal{P}\mathcal{L}_{1}\mu$, but not in the integral, as the action of the corresponding dissipator on a matrix previously acted upon with $\mathcal{P}$ remains in the diagonal, and vanishes when projected under $\mathcal{Q}$. Consequently, the integrand in (6) is

[TABLE]

The second and third terms on the right hand side of the equation contain operators such as $\sigma_{ss^{\prime}}^{(k)}\sigma_{s^{\prime}s}^{(l)}$ that move off the diagonal the terms of $\mu$ corresponding to two different sites. Those operators are preceded or followed by operators such as $\sigma_{s^{\prime}x}^{(k)}$ that act on just one site, so the final result vanishes when acted upon by $\mathcal{P}$. A similar reasoning helps to simplify the two remaining terms. The first term simplifies because the action of $\sigma_{sx}^{(k)}=\left|s\right>_{k}\!\left<0\right|+\left|0\right>_{k}\!\left<s\right|$ combined with that of $\sigma_{s^{\prime}x}^{(m)}=\left|s^{\prime}\right>_{m}\!\left<0\right|+\left|0\right>_{m}\!\left<s^{\prime}\right|$ can only produce non-zero diagonal elements if $s^{\prime}=s$ and $m=k$. Analogously, in the fourth term only summands for which $k=m$, $l=n$, $s=s^{\prime\prime}$ and $s^{\prime}=s^{\prime\prime\prime}$ can yield nonzero contributions. Taking these facts into account simplifies the expression in Eq. (22) to that in Eq. (8).

Appendix B Driving term $D[\mu,t]$

The driving term $D[\mu,t]$ in Eq. (9) is

[TABLE]

For concreteness, we focus on the contribution of level $\left|1\right>$ on the right hand side. The first term yields

[TABLE]

The other terms can be similarly worked out. In the following, we use $\mathcal{V}^{k}_{s}=\sum_{m}\left[V^{s}_{km}n_{s}^{(m)}+\sum_{s^{\prime}\neq s}V^{ss^{\prime}}_{km}n_{s^{\prime}}^{(m)}\right]$ as shorthand for the increment in the interaction energy corresponding to the excitation of atom $k$ to level $|s\rangle$. The term corresponding to $s=1$ in Eq. (23) becomes

[TABLE]

and those contributions due to the other levels take an analogous form. We thus obtain

[TABLE]

where the projection operator $\mathcal{I}_{s}^{(k)}=n_{s}^{(k)}+|0\rangle_{k}\langle 0|$. The effective dynamics due to the driving in $\mathcal{H}_{1}$ is given by

[TABLE]

Appendix C Exchange term: $E[\mu,t]$

The excitation swap term $E[\mu,t]$ in Eq. (9) is

[TABLE]

Cancellations occur as diagonal terms are moved off and not brought back on the diagonal, which then vanish under the action of the projector $\mathcal{P}$, resulting in

[TABLE]

Each term of the type $\sigma_{ss^{\prime}}^{(k)}\sigma_{s^{\prime}s}^{(l)}\ e^{\mathcal{L}_{0}t}\left(\sigma_{s^{\prime}s}^{(k)}\sigma_{ss^{\prime}}^{(l)}\mu\right)$ can be unravelled as follows

[TABLE]

Here, $\left[\rho_{s^{\prime}s^{\prime},ss}^{(k,l)}\right]_{ss^{\prime},s^{\prime}s}$ denotes the $p^{N-2}\times p^{N-2}$ matrix defined by $\rho_{s^{\prime}s^{\prime},ss}^{(k,l)}=\left({}_{k}{\left<s^{\prime}\right|}\otimes{}_{l}{\left<s\right|}\right)\rho\left(\left|s^{\prime}\right>_{k}\otimes\left|s\right>_{l}\right)$ located in the position of the (off-diagonal, and therefore rapidly decaying due to the dephasing) matrix element $\rho_{ss^{\prime},s^{\prime}s}^{(k,l)}$, with all the other projected matrix elements being zero. We next work out the second term in Eq. (44):

[TABLE]

The third term in Eq. (44) yields

[TABLE]

and the fourth term yields

[TABLE]

The remaining four terms in (44) can be easily worked out from (45,46,47,48) by swapping the indices of the excited levels $s\leftrightarrow s^{\prime}$. Putting all these into Eq. (44),

[TABLE]

The exchange term in Eq. (9) can thus be written as

[TABLE]

where the projector $\mathcal{I}_{ss^{\prime}}^{(k,l)}=n_{s}^{(k)}n_{s^{\prime}}^{(l)}+n_{s^{\prime}}^{(k)}n_{s}^{(l)}$ and $\Gamma_{ss^{\prime}}^{(kl)}=\left({1+\left[\frac{1}{\gamma_{s}\!+\!\gamma_{s}^{\prime}}\left(\mathcal{V}^{k}_{s^{\prime}}\!-\!\mathcal{V}^{k}_{s}\!+\!\mathcal{V}^{l}_{s}\!-\!\mathcal{V}^{l}_{s^{\prime}}\right)\right]^{2}}\right)^{-1}\!.$

Appendix D Excitation swaps and redistribution of stationary-state probabilities

We consider a simple model that illustrates how excitation swaps lead to a redistribution of probabilities in the stationary state. The configuration space comprises two “ground” states and two “excited” states, denoted as $|g_{1}\rangle$, $|g_{2}\rangle$, $|e_{1}\rangle$ and $|e_{2}\rangle$. From the perspective of the analogy between this abstract model and a Rydberg gas, the “ground” states may actually have excitations, as long as their number is smaller than the number of excitations in the corresponding “excited” states (we will not further consider the configuration of these states). The transitions between states are as follows [see Fig. 4 (a)].

• •

The system can excite from $|g_{1}\rangle$ to $|e_{1}\rangle$ with a rate $\Gamma_{1}$, and from $|g_{2}\rangle$ to $|e_{2}\rangle$ with a rate $\Gamma_{2}$.

• •

The system can de-excite from $|e_{1}\rangle$ to $|g_{1}\rangle$ with a rate $\Gamma_{1}+\kappa$, and from $|e_{2}\rangle$ to $|g_{2}\rangle$ with a rate $\Gamma_{2}+\kappa$. Here $\kappa$ plays the role of a decay rate.

• •

The system can undergo swap processes from $|e_{1}\rangle$ to $|e_{2}\rangle$ with a rate $U_{e}$, and from $|g_{1}\rangle$ to $|g_{2}\rangle$ with a rate $U_{g}$. The swap strength is smaller in the excited states $U_{e}\leq U_{g}$, but both are assumed to be proportional to an overall swap strength $U$.

The probabilities of each of the four states, which we denote as $g_{1}$, $g_{2}$, $e_{1}$ and $e_{2}$, evolve in time as follows

[TABLE]

These equations are not independent, since they satisfy $\dot{g}_{1}+\dot{g}_{2}+\dot{e}_{1}+\dot{e}_{2}=0$ by probability conservation. By setting three of them to zero, we obtain the stationary state values of e.g. $g_{1}$, $g_{2}$ and $e_{1}$ as functions of $e_{2}$. The normalization, $g_{1}+g_{2}+e_{1}+e_{2}=1$, provides the remaining condition, yielding

[TABLE]

For $\kappa=0$ we obtain a fully mixed state $g_{1}^{ss}=g_{2}^{ss}=e_{1}^{ss}=e_{2}^{ss}=\frac{1}{4}$ (and detailed balance is satisfied).

We now assume that $\Gamma_{1}\gg\kappa\gg\Gamma_{2}$. This means that transitions between $|g_{1}\rangle$ and $|e_{1}\rangle$ occur almost as frequently in one direction as in the opposite direction ($\Gamma_{1}+\kappa\approx\Gamma_{1}$), while excitations from $|g_{2}\rangle$ to $|e_{2}\rangle$ occur much less frequently than the opposite (de-excitation) process ($\Gamma_{2}+\kappa\approx\kappa$). This assumption allows us to reproduce (in a simplified form) the dynamic heterogeneity of the full Rydberg dynamics, where both highly mobile ($|g_{1}\rangle$ and $|e_{1}\rangle$) and highly arrested ($|g_{2}\rangle$ and $|e_{2}\rangle$) configurations exist.

We first consider the situation when swap processes are absent, $U=0$. The stationary state probability distribution is then obtained by setting $U_{e}=U_{g}=0$ in Eqs. (51, 52,53,D) [or one can directly simplify Eqs. (54, 55, 56, 57)], which yields

[TABLE]

These are two uncoupled two-level systems with very different stationary distributions: in one case the probabilities of “ground” state and “excited” state are almost identical, in the other the population of the “ground” state is much larger.

We then focus on the opposite situation, i.e. when swap processes occur very frequently. If both $U_{e}$ and $U_{g}$ are very large, the dominant contributions in Eqs. (54, 55, 56, 57) are those of second order in $U$, leading to

[TABLE]

We conclude that an enhancement of exchange processes eventually leads to a vanishing probability imbalance between “ground” and “excited” states, $e_{2}^{ss}\approx g_{2}^{ss}$.

So far we have only considered the two limiting cases of vanishing swap strength or very large swap strength. In Fig. 4 (b), we show $g_{1}^{ss}$, $g_{2}^{ss}$, $e_{1}^{ss}$ and $e_{2}^{ss}$ as functions of $U$ for parameter values satisfying $\Gamma_{1}\gg\kappa\gg\Gamma_{2}$ (similar results are obtained for other parameter choices that fulfill this condition). This illustrates the gradual change of the stationary state probabilities as the overall swap strength $U$ is increased. The curves smoothly join the solution in Eq. (58) to that in Eqs. (59,60) as $U$ is increased. These results provide a simple picture of the basic mechanism underlying the much more complex situation considered in Section IV.C in general, and more specifically in the numerical results shown in Fig. 3 (c).