Non-equilibrium absorbing state phase transitions in discrete-time quantum cellular automaton dynamics on spin lattices

Recent years have witnessed breakthroughs in the realization of quantum simulator platforms based on cold atomic systems [1–5]. One of the most recent generations of these quantum simulators is based on Rydberg atoms and offers freely programmable and addressable spin arrays [2–4, 6]. When excited to (high-lying) Rydberg states atoms interact strongly, thereby offering a versatile platform for the study of quantum matter in an out of equilibrium.

Strong interactions between Rydberg atoms are moreover at the heart of implementations of quantum information processing protocols [7] where they allow the realization of conditional gates [8, 9] that generate entangling operations. Digital quantum simulators [10] employ such gates—similar to the circuit-based approach to quantum computing—and represent a route towards emulating quantum dynamics with exotic interactions. The possibility of digitally simulating open and closed many-body systems with Rydberg lattice systems was theoretically explored in [11] and the capability of this platform for preparing exotic many-body systems and states was highlighted. While their experimental realization has not yet been achieved, first proof-of-principle demonstrations of the feasibility of this idea were demonstrated within a trapped ion quantum simulator [12] and superconducting circuits [13, 14].

Here we introduce a class of spin models with discrete-time quantum dynamics that lends itself rather naturally to the implementation on a Rydberg quantum simulator. The dynamics takes place within a 1 + 1-dimensional lattice in which the directions can be thought of representing time and space, respectively. Propagation between time slices proceeds via the successive application of three-body gates that perform conditional unitary rotations. Despite the fact that the dynamics of the whole system is unitary, entanglement between time slices leads to the convergence of the reduced state on the final time slice. This stationary state may display a non-equilibrium phase transition and features non-classical correlations that become long-ranged in the vicinity of the transition point. We illustrate our idea using an example that is efficiently solvable in the sense that it permits the mapping onto the non-equilibrium process of a classical cellular automaton for site percolation.

Our work highlights the emergence of stationary behavior in closed quantum systems and introduces a new aspect in extending the concept of a cellular automaton into the quantum domain—for a few examples see [15–21]. Moreover, it connects to quantum generalizations of perceptrons in neural networks [22, 23]. Our proposed setting provides a natural test bed for assessing the capabilities of current Rydberg lattice quantum simulators: it possesses non-trivial features, such as a phase transition and long-ranged quantum correlations, but yet can be efficiently solved. It can thus be used for the certification of a Rydberg simulator in a regime (two dimensions, strong interactions, long times) which is usually numerically intractable.

The two-dimensional lattice system we are considering is depicted in figure 1. Each row consists of N sites, with a spin-1/2 degree of freedom per site. The horizontal and vertical directions we consider as space and time, respectively. The dynamics starts from a state where all spins are in the state $\left|\downarrow ,\rangle \right\rangle$ except for the first time slice (first row) which is prepared in the desired initial configuration. The evolution then proceeds by applying a sequence of elementary gates linking the time slice at time t to the time slice at time t + 1.

Zoom In Zoom Out Reset image size

For the case we are mainly interested in this work these elementary gates are unitary operators that act on three spins simultaneously—two consecutive ones on time slice t (control spins) and one on time slice t + 1 (target spin), as shown in figure 1. These gates perform a rotation of the state of the target spin, conditioned on the presence or absence of excited spins (in state $\left|\uparrow ,\rangle \right\rangle$) among the two-control spins. We consider this type of gate here because it can be rather naturally implemented in Rydberg lattice quantum simulators as is discussed further below. Formally, we can write the gate as

$$\begin{eqnarray}&&{G}_{m}^{(t)}={P}_{m,m+1}^{(t-1)}\otimes {U}_{m}^{(t)}+{Q}_{m,m+1}^{(t-1)}\otimes {{\mathbb{I}}}_{m}.\end{eqnarray} \tag{ 1 }$$

Here ${P}_{m,m+1}^{(t-1)}$ and ${Q}_{m,m+1}^{(t-1)}$ are projection operators, which act on the control spins on time slice $t-1$ (with indices m and m + 1) and obey ${P}_{m,m+1}^{(t-1)}+{Q}_{m,m+1}^{(t-1)}={\mathbb{I}}$. To be specific we use for now

$$\begin{eqnarray}&&{P}_{m,m+1}=1-\left(1-{n}_{m}\right)\left(1-{n}_{m+1}\right),\end{eqnarray} \tag{ 2 }$$

where ${n}_{m}=(1+{\sigma }_{z}^{m})/2$, projects onto the excited state of the mth spin on time slice t − 1 and ${\sigma }_{z}^{m}$ is a Pauli matrix. The projector ${P}_{m,m+1}$ returns a non-zero value only if at least one of the control spins is in the excited state. When this is the case the unitary operator U^(t) acts on the target atom on time slice t and performs a spin rotation about the y-axis by an angle α: $U=\exp \left(-{\rm{i}}\tfrac{\alpha }{2}{\sigma }_{y}\right)$. Note, that we dropped the time slice index t in the explicit forms of both the projectors and the unitary in order not to make the notation too contrived.

The rule (1) can be considered as imposing a kinetic constraint in the dynamics, reminiscent of facilitated models of glasses [24]. Constrained dynamics—which has been experimentally shown to take place e.g. in interacting Rydberg gases [25]—can give rise to complex evolution both in classical [26] and closed and open quantum systems [21, 27–32]. In particular, a rule akin to (1), of at least one nearest neighbour in the excited state required to allow for local evolution, is known in classical facilitated models to lead to an effective dynamics of the reaction-diffusion kind [33, 34].

In our model the propagation from time slice t − 1 to t is achieved via the concatenation of gates, ${G}_{N}^{(t)}...{G}_{2}^{(t)}{G}_{1}^{(t)}$, where we assume periodic boundary conditions when applying ${G}_{N}^{(t)}$. Note, that due to the specific choice made in equation (2) the actual order of the gates is not important since the projectors commute. The successive application of the gate ${{ \mathcal G }}^{(t)}$ to subsequent time slices propagates the initial state and creates a pure state (provided that the initial state has been pure) on the entire lattice.

The reduced state ρ_t on time slice t is linked to the reduced state of the previous time slice by a recurrence relation:

$$\begin{eqnarray}&&{\rho }_{t}=\displaystyle \sum _{{i}_{1},\ldots {i}_{N}=0,1}\mathrm{Tr}\left[{X}_{1}^{({i}_{1})}...{X}_{N}^{({i}_{N})}{\rho }_{t-1}\right]{\rho }^{({i}_{1})}\,\otimes \,...\,\otimes \,{\rho }^{({i}_{N})},\end{eqnarray} \tag{ 3 }$$

with ${X}_{m}^{(0)}={P}_{m,m+1}$ and ${X}_{m}^{(1)}={Q}_{m,m+1}$ as well as ${\rho }^{(0)}=U\left|\downarrow ,\rangle \right\rangle \left\langle \downarrow \right|{U}^{\dagger }$ and ${\rho }^{(1)}=\left|\downarrow ,\rangle \right\rangle \left\langle \downarrow \right|$. The state ρ_t is separable as it is formed by a convex superposition of product states of the form ${\rho }^{({i}_{1})}\,\otimes \,...\,\otimes \,{\rho }^{({i}_{N})}$. The weight of each state is given by the expectation value of the product of projection operators taken in the state of the previous time slice, ${\rho }_{t-1}$. In our protocol local quantum operations, such as ${U}_{m}^{(t)}$, are conditioned by a measurement result i_m, that can be communicated 'classically'. Such scheme cannot produce an entangled state on time slice t. Nevertheless, ρ_t can exhibit non-classical correlations as we show later.

For a first understanding of the discrete-time dynamics we conduct a mean field study. To this end we consider the evolution of the local density on site m under the gate (1), which yields

$$\begin{eqnarray*}&&\langle {n}_{m}^{(t)}{\rangle }_{t}=\langle \downarrow | {{G}_{m}^{(t)}}^{\dagger }{n}_{m}^{(t)}{G}_{m}^{(t)}| \downarrow {\rangle }_{t}=x\,{P}_{m,m+1}^{(t-1)},\end{eqnarray*}$$

with $x= \langle \,\downarrow | {{U}_{m}^{(t)}}^{\dagger }{n}_{m}^{(t)}{U}_{m}^{(t)}| \downarrow \,{ \rangle }_{t}={\sin }^{2}\left(\tfrac{\alpha }{2}\right),x\in [0,1]$. We take the expectation value over the t − 1-time slice, make use of the form (2) of the projector ${P}_{m,m+1}$ and perform the mean field approximation (decoupling of pair correlators and assumption of homogeneity). This yields a recurrence relation, connecting the mean field densities ν at time slices t and t − 1:

$$\begin{eqnarray}&&{\nu }^{(t)}=x\,{\nu }^{(t-1)}\left(2-{\nu }^{(t-1)}\right).\end{eqnarray} \tag{ 4 }$$

To make progress we turn the recurrence relation into a differential equation $[{\nu }^{(t)}\to \nu (t)$, ${\nu }^{(t)}-{\nu }^{(t-1)}\to {\partial }_{t}\nu (t)$]. Choosing the initial condition n(0) = 1, we obtain the solution

$$\begin{eqnarray}&&\nu (t)=\displaystyle \frac{2x-1}{x+\exp \left(t\left[1-2x\right]\right)\left(x-1\right)},\end{eqnarray} \tag{ 5 }$$

which has an interesting limiting behavior at long times: for $x\lt {x}_{\mathrm{crit}}=1/2$ we find ${\mathrm{lim}}_{t\to \infty }\nu (t)=0$, while for $x\gt {x}_{\mathrm{crit}}$ the excitation density assumes the non-zero stationary value ${\mathrm{lim}}_{t\to \infty }\nu (t)=(x-{x}_{\mathrm{crit}})/(2x)$. Thus, x_crit defines a critical rotation angle ${\alpha }_{\mathrm{crit}}=\pi /4$ which in the limit $t\to \infty$ separates two qualitatively different states. At x = x_crit we find $\nu (t)=2/(2+t)$ and thus the density displays an algebraic approach to stationarity. This result is reminiscent of mean field calculations of classical reaction-diffusion problems that feature absorbing state phase transitions [35].

Further insight into this phase transition behavior is obtained by exploiting a link to a classical stochastic process: due to the separability of the reduced density matrices ρ_t and the structure of the projectors (2), the probabilities $\mathrm{Tr}\left[{X}_{1}^{({i}_{1})}...{X}_{N}^{({i}_{N})}{\rho }_{t-1}\right]$, which appear in the reduced state (3), can be generated via a classical discrete time dynamics. This process takes place on a two-dimensional lattice, as depicted in figure 1, that contains classical spins (either up or down), initially prepared in the state $\parallel \downarrow \rangle\rangle$. The discrete time evolution proceeds via the classical maps

$$\begin{eqnarray}{W}_{m}^{(t)}={P}_{m,m+1}^{(t-1)}\otimes {\left(\begin{array}{cc}1-x & x\\ x & 1-x\end{array}\right)}_{m}+{Q}_{m,m+1}^{(t-1)}\otimes {{\mathbb{I}}}_{m}\end{eqnarray} \tag{ 6 }$$

which are applied on a probability vector in order to propagate the system between time slices. This dynamics implements an instance of the so-called Domany–Kinzel (DK) cellular automaton [35, 36] and it performs a flip of the target spin (time slice t) with probability x, provided that the projection operator ${P}_{m,m+1}^{(t-1)}$ yields a non-zero value when applied to the control spin on time slice t − 1. Under this dynamics the reduced probability vector $\parallel {\bf{p}}{ \rangle\rangle }_{t}$ of time slice t evolves according to

$$\begin{eqnarray}&&\parallel {\bf{p}}{ \rangle\rangle }_{t}=\displaystyle \sum _{{i}_{1},\ldots {i}_{N}=0,1} \langle\langle +\parallel {X}_{1}^{({i}_{1})}...{X}_{N}^{({i}_{N})}\parallel {\bf{p}}{ \rangle\rangle }_{t-1}\parallel {s}_{{i}_{1}} \rangle\rangle \,\otimes \,...\,\otimes \,\parallel {s}_{{i}_{N}} \rangle\rangle \end{eqnarray} \tag{ 7 }$$

with $\parallel {s}_{0} \rangle\rangle =\parallel \uparrow \rangle\rangle$ and $\parallel {s}_{1} \rangle\rangle =(1-x)\parallel \downarrow \rangle\rangle +x\parallel \uparrow \rangle\rangle$. Note, that instead of taking the trace, expectation values in this classical description are calculated by applying the desired operator to the probability vector and multiplying from the left with a (flat) reference state: for N spins this is $\parallel + \rangle\rangle ={\displaystyle \bigotimes }_{m=1}^{N}\left[\parallel \downarrow { \rangle\rangle }_{m}+\parallel \uparrow { \rangle\rangle }_{m}\right]$.

The structural resemblance between the reduced state (3) and the probability vector (7) is evident. The local quantum states ρ^(k) and classical states $\parallel {s}_{k} \rangle\rangle$ are constructed such that they yield the same expectation values for classical observables, e.g. $\mathrm{Tr}\left(n\,{\rho }^{(m)}\right)= \langle\langle +\parallel n\parallel {s}_{m} \rangle\rangle =x\,{\delta }_{m,1}$. Thus, also the states (3) and (7) yield identical expectation values of classical observables, and in this sense the discrete time quantum dynamics is mapped onto a classical process.

The connection to the DK cellular automaton provides an explanation for the phase transition behavior observed in the mean field calculation: it is known that the cellular automaton dynamics (6) leads to a non-equilibrium stationary state which displays a continuous (absorbing state) phase transition between a so-called inactive phase—in which the expectation value of the average density $\langle n\rangle =\tfrac{1}{N}{\sum }_{m}^{N}\langle {n}_{m}\rangle$ is zero—and an active phase in which $\langle n\rangle \ne 1$. This transition occurs at x ≈ 0.7 and is in the directed percolation universality class. The corresponding numerical data is shown in figures 2(a)–(c).

Zoom In Zoom Out Reset image size

Despite being separable and related to a classical dynamics, the state (3) possesses non-classical correlations, as we show now. Furthermore, by exploiting the mapping to the DK cellular automaton dynamics we find that it is possible to extract quantum correlations from the measurement of classical observables, which are straight-forwardly accessible on Rydberg quantum simulators [3].

As a measure for quantum correlations we employ the LQU put forward in [37] which is a variant of bipartite quantum discord [38–40]. It quantifies how much of the fluctuations of a local measurement is due to the non-commutativity between the state and the measured local observable, which is caused only by the states coherence, not its mixedness. By minimising over the choice of the local observable, only non-local coherences—necessarily corresponding to quantum correlations—are captured. In our model, the non-local coherence in the state (3) can be thought of as being a result of an effectively classical communication of a measurement outcome on time slice t − 1 and conditioned local coherent preparation, via ${U}_{m}^{(t)}$ or ${I}_{m}^{(t)}$, of atoms on time slice t, see (1) and (3).

For the reduced state ρ_ij of two spins the LQU is defined as ${{\ell }}_{{ij}}=1-{\lambda }_{\max }\left\{{W}^{{ij}}\right\}$, where ${\lambda }_{\max }\left\{{W}^{{ij}}\right\}$ is the largest eigenvalue of the matrix W^ij with components ${W}_{\alpha \beta }^{{ij}}=\mathrm{Tr}\left({\rho }_{{ij}}^{1/2}{\sigma }_{\alpha }^{i}{\rho }_{{ij}}^{1/2}{\sigma }_{\beta }^{i}\right)$. The reduced density matrix ρ_ij can be obtained entirely from measuring the local density and density–density correlations between sites i and j. To see this we exploit the special structure of the reduced state (3): each term of the sum contains a product of pure states which allows to relate expectation values of off-diagonal operators to those of diagonal observables, e.g. $\mathrm{Tr}\left({\sigma }_{i}^{\pm }{\rho }_{t}\right)=\sqrt{\tfrac{1-x}{x}}\mathrm{Tr}\left({n}_{i}{\rho }_{t}\right)=\sqrt{\tfrac{1-x}{x}}\left\langle {n}_{i}\right\rangle$. Using this property, and assuming translation invariance ($\left\langle n\right\rangle =\left\langle {n}_{i}\right\rangle =\left\langle {n}_{j}\right\rangle$), one obtains

$$\begin{eqnarray}{\rho }_{{ij}}=\left(\begin{array}{cc}{c}_{{ij}} & x\,{c}_{{ij}}\\ x\,{c}_{{ij}} & 1-{c}_{{ij}}\end{array}\right)\otimes \left(\begin{array}{cc}{c}_{{ij}} & x\,{c}_{{ij}}\\ x\,{c}_{{ij}} & 1-{c}_{{ij}}\end{array}\right)+\left[\left\langle n\right\rangle -{c}_{{ij}}\right]\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 1 & 0 & x\\ 0 & 0 & 1 & x\\ 0 & x & x & -2\end{array}\right),\end{eqnarray} \tag{ 8 }$$

with ${c}_{{ij}}=\sqrt{\left\langle {n}_{i}{n}_{j}\right\rangle }$ being the square root of the density–density correlation function.

In the absence of correlations one has ${c}_{{ij}}=\left\langle n\right\rangle$. Here, the second term in equation (8) vanishes and ρ_ij becomes a product state without quantum correlations. This is the case away from a phase transition where correlations between two sites are decaying rapidly as a function of their distance. Near a phase transition, however, correlations are long-ranged, as is shown in figure 2(c), where we display the connected density–density correlation function ${C}_{{ij}}={c}_{{ij}}^{2}-{\left\langle n\right\rangle }^{2}$. Here also finite and long-ranged quantum correlations, characterized through the LQU, emerge, as can be seen in figure 2(d).

Note, that the entries of the density matrix (8) can be determined experimentally via state tomography, i.e. the measurement of the expectation values of products of the three Pauli matrices and the identity matrix for the atoms on sites i and j. The LQU can thus be inferred directly from experimental measurements, and its scaling, see figure 2(d), requires a number of measurements which scales only linearly in the system size, as it is exactly the case when measuring classical correlations.

The open cellular automaton model discussed here can be implemented on Rydberg quantum simulators [3, 4]. The three-body gates underlying the gate (1) are implemented by employing the blockade interaction [8] which yields conditional unitaries [41] discriminating between the cases in which at least one or none of the source atom is excited, in direct analogy with gate (1). For the experimental investigation of the non-equilibrium dynamics it is moreover not necessary to have a two-dimensional lattice. Two parallel one-dimensional arrays (or concentric rings if periodic boundary conditions are required) are sufficient for the following protocol (see also figure 3(a)): (i) the initial state is prepared on the first chain (ring) and all sites of the second chain (ring) are prepared in the state $\left|\downarrow ,\rangle \right\rangle$. (ii) The discrete-time propagation is performed from the first to the second chain (ring). (iii) The first chain (ring) is reset, so that all sites are in state $\left|\downarrow ,\rangle \right\rangle$. (iv) The process is repeated but the role of the chains (rings) is interchanged.

Zoom In Zoom Out Reset image size

An interesting practical aspect of this protocol might be that it allows to investigate collective phenomena—such as absorbing state phase transitions—that are commonly explored in the context of open quantum systems, see e.g. [42]. The discrete-time propagation presented is fully coherent, it does not involve processes such as radiative decay photon scattering. This removes typically detrimental sources of heating and thus may facilitate the observation of dynamical processes of long times.

Generalizations of the dynamics presented here can be achieved by extending the fundamental gate (1) to more source/target atoms and/or by introducing more conditional spin rotations. One possible extension of the gate to K source atoms and one target atom is given by

$$\begin{eqnarray}&&{F}_{m}^{(t)}=\displaystyle \sum _{k=0}^{K}{{\rm{\Pi }}}_{m}^{(t-1)}(k,K)\otimes {U}_{m}^{(t)}({\alpha }_{k}).\end{eqnarray} \tag{ 9 }$$

Here the operators Π_m(k, K) project on the subspace containing k excitations among the K source atoms whose state conditions the state change of the mth target atom. The latter is rotated by the unitary $U({\alpha }_{k})=\exp \left(-{\rm{i}}\tfrac{{\alpha }_{k}}{2}{\sigma }_{y}\right)$. We anticipate two interesting cases here:

(i) K = 2 source sites and rotation angles are given by α₂ = π, α₁ = α and α₀ = 0: the corresponding non-equilibrium process has the two absorbing states $\left|\downarrow \downarrow ...\downarrow ,\rangle \right\rangle$ and $\left|\uparrow \uparrow ...\uparrow ,\rangle \right\rangle$. At α = π/2 the stationary state switches between these two possibilities and displays a phase transition that is in the directed compact percolation universality class [43].

(ii) α₀ = 0, which ensures the presence of the absorbing state $\left|\downarrow \downarrow ...\downarrow ,\rangle \right\rangle$: here the mean field density follows the recurrence relation

$$\begin{eqnarray}&&{\nu }^{(t)}=\displaystyle \sum _{k=1}^{K}{x}_{k}\left(\displaystyle \genfrac{}{}{0em}{}{K}{k}\right){\left({\nu }^{(t-1)}\right)}^{k}{\left(1-{\nu }^{(t-1)}\right)}^{K-k},\end{eqnarray} \tag{ 10 }$$

where ${x}_{k}={\sin }^{2}\left(\tfrac{{\alpha }_{k}}{2}\right)$. This process features a host of absorbing state phase transitions, coexistence regions and critical lines. Moreover, a suitable choice of the rotation angles ${\alpha }_{k}$ allows to set all terms of order smaller than K to zero which tunes the system to a multi-critical point (similar to tri-critical directed percolation [44]): ${\nu }^{(t)}-{\nu }^{(t-1)}\propto -{\left({\nu }^{(t-1)}\right)}^{K}$. Here the mean field density displays a power-law behavior on approach to stationarity: $\nu (t)\sim {t}^{1/(1-K)}$. In figure 3(b) we illustrate the case K = 4.

We studied examples of quantum non-equilibrium processes that can be mapped onto classical cellular automata and therefore efficiently solved. Beyond being of conceptual interesting these findings can be applied for scrutinizing current quantum simulation platforms under challenging, yet numerically tractable, conditions. An interesting subject for future investigations is the realization of non-equilibrium processes with absorbing (dark) states [45] that feature entanglement and/or phase coherence between different sites. Those can be achieved by employing projectors in the fundamental gate (9) that project for example on two-site entangled states, in conjunction with unitary operations acting on two and more target sites. Future studies may also investigate the role of imperfections in the fundamental gates and how they influence the statics and dynamics of the absorbing state phase transition discussed here.

We thank M Marcuzzi and M Müller for useful discussions. The research leading to these results has received funding from the European Research Council under the European Unions Seventh Framework Programme (FP/2007-2013)/ERC [grant agreement number 335266 (ESCQUMA)], and was supported by the Engineering and Physical Sciences Council [grant numbers EP/M014266/1 and EP/R04340X/1] as well as the Leverhulme Trust [grant number RPG-2018-181]. IL gratefully acknowledges funding through the Royal Society Wolfson Research Merit Award.